% INSITU2POT  Convert depth, salinity and insitu temperature to potential 
%       temperature. Assume reference pressure of 0 (which also helps 
%       reduce memory requirements.)
%
% Most of code cut straight from the SEAWATER library, compiled by Phil
% Morgan & Lindsay Pender.             - JRD 9/2/06
%
% INPUT:  t    [ncast ndep]
%         s    [ncast ndep]
%         lat  [ncast 1]
%         dep  [ncast ndep] or [1 ndep] if appropriate
%
% *** DO NOT READ INPUT DETAILS ABOVE     (and get what you deserve) ***
%
% USAGE: pt = insitu2pot(t,s,lat,dep);

function pt = insitu2pot(T,S,lat,dep)

[ncst,ndep] = size(S);

X = latcor(lat);
C1 = 5.92E-3 + (X(:,ones(1,ndep)).^2 * 5.25E-3);

if min(dims(dep))==1 & ncst>1
   P = ((1-C1)-sqrt(((1-C1).^2)-(8.84E-6*dep(ones(ncst,1),:))))/4.42E-6;
else
   P = ((1-C1)-sqrt(((1-C1).^2)-(8.84E-6*dep)))/4.42E-6;
end
clear C1 X dep

% theta1
del_th = -P.*adtg_sw(S,T*1.00024,P);
th     = T * 1.00024 + 0.5*del_th;
q      = del_th;

% theta2
del_th = -P.*adtg_sw(S,th,0.5*P);
th     = th + (1 - 1/sqrt(2))*(del_th - q);
q      = (2-sqrt(2))*del_th + (-2+3/sqrt(2))*q;

% theta3
del_th = -P.*adtg_sw(S,th,0.5*P);
th     = th + (1 + 1/sqrt(2))*(del_th - q);
q      = (2 + sqrt(2))*del_th + (-2-3/sqrt(2))*q;

% theta4
del_th = -P.*adtg_sw(S,th,0);
pt     = (th + (del_th - 2*q)/6)/1.00024;

return
%-----------------------------------------------------------------------
% Just the SEAWATER sw_adtg reformulated to handle bigger arrays without 
% swamping memory

function ADTG = adtg_sw(S,T,P)

a0 =  3.5803E-5;
a1 = +8.5258E-6;
a2 = -6.836E-8;
a3 =  6.6228E-10;

b0 = +1.8932E-6;
b1 = -4.2393E-8;

c0 = +1.8741E-8;

c1 = -6.7795E-10;
c2 = +8.733E-12;
c3 = -5.4481E-14;

d0 = -1.1351E-10;
d1 =  2.7759E-12;

e0 = -4.6206E-13;
e1 = +1.8676E-14;
e2 = -2.1687E-16;

[m,n] = size(T);
if m*n > 10E6

   ADTG = zeros(m,n);
   for jj = 1:100:size(T,1)
      kk = jj:(min([jj+99 m])); 
      ADTG(kk,:) = (a2 + a3.*T(kk,:)).*T(kk,:);
      ADTG(kk,:) = a0 + (a1 + ADTG(kk,:)).*T(kk,:);
      tmp = (b0 + b1.*T(kk,:)).*(S(kk,:)-35);
      ADTG(kk,:) = ADTG(kk,:) + tmp;
      if P ~= 0
	 tmp = (c2 + c3.*T(kk,:)).*T(kk,:);
	 tmp = c0 + (c1 + tmp).*T(kk,:);
	 tmp = tmp + (d0 + d1.*T(kk,:)).*(S(kk,:)-35);
	 ADTG(kk,:) = ADTG(kk,:) + tmp.*P(kk,:);
	 tmp = e0 + (e1 + e2.*T(kk,:)).*T(kk,:);
	 ADTG(kk,:) = ADTG(kk,:) + tmp.*P(kk,:).*P(kk,:);
      end
   end
else   
   ADTG = a0 + (a1 + (a2 + a3.*T).*T).*T ...
	  + (b0 + b1.*T).*(S-35);
   if P ~= 0 
      ADTG = ADTG + ( (c0 + (c1 + (c2 + c3.*T).*T).*T) ...
		      + (d0 + d1.*T).*(S-35) ).*P ...
	     + (  e0 + (e1 + e2.*T).*T ).*P.*P;
   end
end

%-------------------------------------------------------------------------
