function [Nbeta]=kh_N_beta_sz(si,sf,ti,tf,pi,pf)
% function [Nbeta]=kh_N_beta_sz(si,sf,ti,tf,pi,pf)
%
% Caluclates N'*beta*S|z (movement of isopycnal*alpha*Salinity change on 
% pressure surfaces
% All terms need to be entered on neutral density levels.
%
% si=salinity initial (neutral density levels*profiles)
% sf=salinity final etc.
% ti= initial potnetial temperature (needs to be potential temperature)
%
% N=difference in height between two surfaces where + value =upward
% movement of isopycnal
%
% Designed for use in calculating terms from Bindoff and McDougall (1994).
% Equation from B&M (1994) p1143 (adapted equation 14)
%
% Written by Kieran Helm (17-5-2007)

% % Approach 1: Difference approach (EQN 11 - bindoff and mcd)
% [sn]=kh_beta_sn(si,sf,ti,tf,pi,pf);
% [sz]=kh_beta_sz(si,sf,ti,tf,pi,pf);
% Nbeta=sn-nz;

%% Approach 2: Full approach (Given in EQN 14)
N=pi-pf; % + = shoaling of isopycnal (to agree with BM1994 p1143) - added 25-9-07
[sm tm pm]=kh_mid_points(si,sf,ti,tf,pi,pf); % Get midpoints
be=sw_beta(sm,tm,pm,'ptmp'); % Calculate alpha

    %% Calculate sfz-siz
    warning off % for divide by zero
    [m n]=size(si);
    gi=NaN(m,n); gf=NaN(m,n); % gradients of temp/depth slope

    % gradient (ds/dz) intial and final profiles
    % note: p*-1 is to reverse xaxis (as looking from depth to surface)
    for k=1:n
         gi(:,k)=gradient(si(:,k),pi(:,k)*-1); 
         gf(:,k)=gradient(sf(:,k),pf(:,k)*-1);
    end

%% If gradient too steep then replace with NaN
gi(find(gi>0.5 | gi<-0.5))=NaN;
gf(find(gf>0.5 | gf<-0.5))=NaN;
    
    
%% combine all elements of equation 14
Nbeta=0.5*N.*be.*[gf+gi];