%% Compute the degree of orthogonality based on gamman
% Define the gradients arrays
grad_gamman = zeros(npi,npj,npk);
grad_salt = zeros(npi,npj,npk);
grad_spice = zeros(npi,npj,npk);
gradi_spicity = zeros(npi,npj,npk);

% Define the gradients 
[grad_gamman(:,:,:,1),grad_gamman(:,:,:,2),grad_gamman(:,:,:,3)] = ...
    grad_density(gamman,longitude,latitude,zlev,ra);
[grad_spice(:,:,:,1),grad_spice(:,:,:,2),grad_spice(:,:,:,3)] = ...
    grad_density(tau_jmd,longitude,latitude,zlev,ra);
[grad_spicity(:,:,:,1),grad_spicity(:,:,:,2),grad_spicity(:,:,:,3)] = ...
    grad_density(tau_huang,longitude,latitude,zlev,ra);
[grad_salt(:,:,:,1),grad_salt(:,:,:,2),grad_salt(:,:,:,3)] = ...
    grad_density(salt,longitude,latitude,zlev,ra);

% Define the angles
sin_gamma_salt = sin_angle(grad_gamman,grad_salt);
sin_gamma_spice = sin_angle(grad_gamman,grad_spice);
sin_gamma_spicity = sin_angle(grad_gamman,grad_spicity);

%% Compute the degree of orthogonality based on sigma1
% Define the gradients arrays
grad_sigma1 = zeros(npi,npj,npk);
grad_salt = zeros(npi,npj,npk);
grad_spice = zeros(npi,npj,npk);
gradi_spicity = zeros(npi,npj,npk);

% Define the gradients 
[grad_sigma1(:,:,:,1),grad_sigma1(:,:,:,2),grad_sigma1(:,:,:,3)] = ...
    grad_density(gamman,longitude,latitude,zlev,ra);
[grad_spice(:,:,:,1),grad_spice(:,:,:,2),grad_spice(:,:,:,3)] = ...
    grad_density(tau_jmd,longitude,latitude,zlev,ra);
[grad_spicity(:,:,:,1),grad_spicity(:,:,:,2),grad_spicity(:,:,:,3)] = ...
    grad_density(tau_huang,longitude,latitude,zlev,ra);
[grad_salt(:,:,:,1),grad_salt(:,:,:,2),grad_salt(:,:,:,3)] = ...
    grad_density(salt,longitude,latitude,zlev,ra);

% Define the angles
sin_sigma1_salt = sin_angle(grad_sigma1,grad_salt);
sin_sigma1_spice = sin_angle(grad_sigma1,grad_spice);
sin_sigma1_spicity = sin_angle(grad_sigma1,grad_spicity);
