function [betasw,beta90sw,bsw]= betasw_hs_ZH2009(lambda,Tc,theta,SA,p,delta)
% Xiaodong Zhang and Lianbo Hu (2009). Scattering by pure seawater at high salinity, Optics Express,
% 17(15), 12685-12691.
%
% lambda [nm]: wavelength, scalar or vector  
% Tc [deg C]: temperature in Celsius, scalar                   
% theta [deg]: scattering angle, scalar or vector      
% SA [g kg^-1]: absolute salinity. SA = S*1.0047, where S is PSS. The reason using absolute salinity
%               is that PSS is not defined for S > 42                      
% p [dbar]: sea (gauge) pressure. Default = 0, which is under normal condition at sea surface;                    
% delta: depolarization ratio. Default = 0.039.
%
% betasw: volume scattering at angles defined by theta. Its size is [x y], where x is the number of 
%         angles (x = length(theta)) and y is the number of wavelengths in lambda (y = length(lambda))
% beta90sw: volume scattering at 90 degree. Its size is [1 y]
% bw: total scattering coefficient. Its size is [1 y]
%     for backscattering coefficients, divide total scattering by 2
%
% Xiaodong Zhang & Lianbo Hu, June 10, 2009

% values of the constants
Kbz = 1.3806503e-23 ;    %  Boltzmann constant
Tk = Tc+273.15 ;         %  Absolute tempearture

betasw = zeros(length(theta),length(lambda));
if ~isscalar(Tc) || ~isscalar (SA)
    error('Both Tc and S need to be scalar variable');
end

if nargin == 4
    p = 0;
    delta = 0.039;
elseif nargin == 5
    delta = 0.039;
end

rad = theta(:)*pi/180; % angle in radian as a colum variable

% nsw: absolute refractive index of seawater
S = SA/1.0047;
nsw = RInw(lambda,Tc,S);

% density derivative of refractive index from PMH model
DFRI = PMH(nsw);  %% PMH model

% partial derivative of g w.r.t p
gp = mygsw(0,1,SA,Tc,p);

% partial derivative of g w.r.t p and p
gpp = mygsw(0,2,SA,Tc,p);

% partial derivative of g w.r.t. p and S
gps = mygsw(1,1,SA,Tc,p);

% partial derivative of g w.r.t. S and S
gss = mygsw(2,0,SA,Tc,p);

% total volume scattering at 90 degree
if SA == 0
    beta90sw =  pi*pi/2*((lambda*1e-9).^(-4))*Kbz*Tk.*DFRI.^2* ...
           (-gpp/gp)*(6+6*delta)/(6-7*delta);
else
    beta90sw = pi*pi/2*((lambda*1e-9).^(-4))*Kbz*Tk.*DFRI.^2* ...
           (-gpp/gp+gps^2/gp/gss)*(6+6*delta)/(6-7*delta);
end

% volume scattering angle distribution
for i=1:length(lambda)
    betasw(:,i)=beta90sw(i)*(1+((cos(rad)).^2).*(1-delta)/(1+delta));
end
% total volume scattering 
bsw=8*pi/3*beta90sw*(2+delta)/(1+delta);



function nsw = RInw(lambda,Tc,S)
% this function is used to calculate the refractive index of seawater 
% from Quan and Fry (1994, Applied Optics)
% lambda: wavelength (nm)
%     Tc: temperature in degree
%      S: salinity (pss)
%    nsw: refractive index of seawater

% refractive index of air is from Ciddor (1996,Applied Optics)
n_air = 1.0+(5792105.0./(238.0185-1./(lambda/1e3).^2)+167917.0./......
        (57.362-1./(lambda/1e3).^2))/1e8;

n0 = 1.31405; n1 = 1.779e-4 ; n2 = -1.05e-6 ; n3 = 1.6e-8 ; n4 = -2.02e-6 ;
n5 = 15.868; n6 = 0.01155;  n7 = -0.00423;  n8 = -4382 ; n9 = 1.1455e6;

nsw = n0+(n1+n2*Tc+n3*Tc^2)*S+n4*Tc^2+(n5+n6*S+n7*Tc)./......
                lambda+n8./lambda.^2+n9./lambda.^3; 
nsw = nsw.*n_air;


function n_density_derivative=PMH(nw)
% this function is used to calculate the density fluctuation of the
% refractive index (DFRI) based on Proutiere et al (1992,J. Phys. Chem.)
% nw: refractive index 

nw2 = nw.^2;
n_density_derivative=(nw2-1).*(1+2/3*(nw2+2).*(nw/3-1/3./nw).^2);


function [n0,gama0,ci,di,ti,ni,alfa,beta,gama,cigma,ai,Ai,bi,Bi,Ci,Di]=......
         IAPWS95_Const
% ----------------------------------constant---------------------------
% Constants coefficients from IAPWS95
% n0 and n1 are updated according to Feistel and Marion 2007 Table1
 n0 = [-8.320446483749693022547976804526 6.683210527593225852567901929071......
       3.00632 0.012436 0.97315 1.27950 0.96956 0.24873];
 gama0 = [0 0 0 1.28728967 3.53734222 7.74073708 9.24437796 27.5075105];
 
 ci = zeros(1,54);
 di = zeros(1,54);
 ti = zeros(1,54); 
 ni = zeros(1,56);

 ci(1:7) = 0; 
 ci(8:22) = 1; 
 ci(23:42) = 2; 
 ci(43:46) = 3;
 ci(47)= 4;
 ci(48:51) = 6;
 ci(52:54) = 0;
 di=[1 1 1 2 2 3 4 1 1 1 2 2 3 4 4 5 7 9 10 11 13 15 1 2 2 2 3 4 4 4 ......
     5 6 6 7 9 9 9 9 9 10 10 12 3 4 4 5 14 3 6 6 6  3 3 3];
 ti = [-0.5 0.875 1 0.5 0.75 .375 1 4 6 12 1 5 4 2 13 9 3 4 11 4 13 ......
       1 7 1 9 10 10 3 7 10 10 6 10 10 1 2 3 4 8 6 9 8 16 22 23 23 ......
       10 50 44 46 50 0 1 4];
 ni = [......
       0.12533547935523d-1  0.78957634722828d1  -0.87803203303561d1  ......
       0.31802509345418    -0.26145533859358    -0.78199751687981d-2 ......
       0.88089493102134d-2 -0.66856572307965     0.20433810950965    ......
      -0.66212605039687d-4 -0.19232721156002    -0.25709043003438    ......
       0.16074868486251    -0.40092828925807d-1  0.39343422603254d-6 ......
      -0.75941377088144d-5  0.56250979351888d-3 -0.15608652257135d-4 ......
       0.11537996422951d-8  0.36582165144204d-6 -0.13251180074668d-11......
      -0.62639586912454d-9 -0.10793600908932     0.17611491008752d-1 ......
       0.22132295167546    -0.40247669763528     0.58083399985759    ......
       0.49969146990806d-2 -0.31358700712549d-1 -0.74315929710341    ......
       0.47807329915480     0.20527940895948d-1 -0.13636435110343    ......
       0.14180634400617d-1  0.83326504880713d-2 -0.29052336009585d-1 ......
       0.38615085574206d-1 -0.20393486513704d-1 -0.16554050063734d-2 ......
       0.19955571979541d-2  0.15870308324157d-3 -0.16388568342530d-4 ......
       0.43613615723811d-1  0.34994005463765d-1 -0.76788197844621d-1 ......
       0.22446277332006d-1 -0.62689710414685d-4 -0.55711118565645d-9 ......
      -0.19905718354408     0.31777497330738    -0.11841182425981    ......
      -0.31306260323435d2   0.31546140237781d2  -0.25213154341695d4 ......
      -0.14874640856724     0.31806110878444                        .......
      ];
  alfa = [20 20 20 ];
  beta = [150 150 250 0.3 0.3 ];
  gama = [1.21 1.21 1.25];
  cigma = [1 1 1];
  ai = [3.5 3.5];
  bi = [0.85 0.95];
  Bi = [0.2 0.2];
  Ci = [28 32];
  Di = [700 800];
  Ai = [0.32 0.32]; 
  %------------------------------------------------------------------------
  
function [x,errvec]=newton(fun, dfun, x0, TOL, DAMPED, nmax, numeric) 
 
 % this fucntion is used toimplementation of the well known Newton method.
 % 
 %      fun: the function to investigate, given as inline function
 %     dfun: the derivative, given as inline function
 %       x0: starting guess
 %      TOL: tolerance for declaration of convergence
 %   DAMPED: 0 for the standard undamped method
 %           1 for the damped method to achieve global convergence
 %     nmax: maximum number of iteration steps
 %  numeric: use numerical differentiation
 %        x: conbergence points
  
 small = 1e-8;
 sigma = 0.0001; % parameter for damped method
 rho = 0.5; % parameter for damped method
 
 xvec = x0;
 x = x0;
 err = abs(feval(fun, x0)); 
 errvec = err;
 errvec2 = [err^2];
 
 nit = 0; 
 while((nit < nmax) && (err > TOL))
   nit = nit + 1; 
   if (numeric)
     dx = max(sqrt(eps)*abs(x), 1e-8); % step width for central difference
     dfx = (feval(fun, x + dx) - feval(fun, x - dx))/(2*dx);
   else
     dfx = feval(dfun,x); % get derivative
   end
 
   if (abs(dfx) < small) % zero derivative -> divergence of the method
     disp('Zero derivative!');
     break;
   else
     d = - feval(fun,x)/dfx; % get search direction
     alpha = 1.0; % damping parameter set to 1
     xold = x;
     x = x + alpha*d; % new iterate
     if (DAMPED),
       f2 = feval(fun,xold).^2;
       f2n = feval(fun,x).^2;
       while (f2n > f2 + sigma*alpha*feval(fun,xold)*d) 
         alpha = rho*alpha; % damping parameter is adjusted
         x = xold + alpha*d; % adjust iterate
         f2n = feval(fun,x).^2;
       end
     end
     xvec = [xvec; x];
     err = abs(feval(fun, x));
     errvec = [errvec; err];
     errvec2 = [errvec2; err^2];
   end
 end
 
function y=fn(x)
% x   : density

% declare global variables
global pG TG

% Call coefficients and parameters used in Helmholtz free energy
 [n0,gama0,ci,di,ti,ni,alfa,beta,gama,cigma,ai,Ai,bi,Bi,Ci,Di]=IAPWS95_Const;

% Reference Constants 
  Tc = 647.096;     % K
  rouC = 322;       % Kg m-3
  R = 0.46151805d3; % J Kg-1 K-1
  
 % ------------------------------------------------------------------------
  tao = Tc/TG;
  dp = x/rouC;
  theta = (1-tao)+Ai.*((dp-1)^2).^(0.5./beta(4:5));
  thea = exp(-Ci.*(dp-1)^2-Di.*(tao-1)^2);
  delta = theta.^2+Bi.*((dp-1).^2).^ai;
  dthea_dp = -2*Ci.*(dp-1).*thea;
  ddelta_dp = (dp-1)*(Ai.*theta*2./beta(4:5).*(dp-1)^2.^......
              (0.5./beta(4:5)-1)+2*Bi.*ai.*((dp-1)^2).^(0.5./beta(4:5)-1));
  ddelta_bi_dp = bi.*delta.^(bi-1).*ddelta_dp;
  
  i=1:7; j=8:51; k=52:54; m=55:56; 
  phi_r = sum(ni(i).*di(i).*dp.^(di(i)-1).*tao.^ti(i))+......
         sum(ni(j).*exp(-dp.^ci(j)).*(dp.^(di(j)-1).*.......
               tao.^ti(j)).*(di(j)-ci(j).*dp.^(ci(j))))+......
         sum(ni(k).*dp.^di(k).*tao.^ti(k).*......
               exp(-alfa.*(dp-cigma).^2-beta(1:3).*(tao-gama).^2).*......
               (di(k)/dp-2*alfa.*(dp-cigma)))+......
         sum(ni(m).*(delta.^bi.*(thea+dp*dthea_dp)+ddelta_bi_dp*dp.*thea));
          
  y = pG-x*R*TG*(1+dp*phi_r);
  
function y=dfn(x)
% x: density

% declare global variables
 global TG

 % Call coefficients and parameters used in Helmholtz free energy
 [n0,gama0,ci,di,ti,ni,alfa,beta,gama,cigma,ai,Ai,bi,Bi,Ci,Di]=IAPWS95_Const;
 
 % Reference Constants 
  Tc = 647.096;     % K
  rouC = 322;       % Kg m-3
  R = 0.46151805d3; % J Kg-1 K-1
  
 % ------------------------------------------------------------------------
  tao = Tc/TG;
  dp = x/rouC;
  theta = (1-tao)+Ai.*((dp-1)^2).^(0.5./beta(4:5));
  thea = exp(-Ci.*(dp-1)^2-Di.*(tao-1)^2);
  delta = theta.^2+Bi.*((dp-1).^2).^ai;
  dthea_dp = -2*Ci.*(dp-1).*thea;
  d2thea_dp = 2*(2*Ci.*(dp-1)^2-1).*Ci.*thea;
  
  
  ddelta_dp = (dp-1)*(Ai.*theta*2./beta(4:5).*(dp-1)^2.^(0.5./......
               beta(4:5)-1)+2*Bi.*ai.*((dp-1)^2).^(0.5./beta(4:5)-1));
  d2delta_dp = 1/(dp-1).*ddelta_dp+(dp-1)^2.*(4*Bi.*ai.*(ai-1).*......
              ((dp-1)^2).^(ai-2)+2*Ai.^2.*(1./beta(4:5)).^2.*......
              (((dp-1)^2)).^(0.5./beta(4:5-1)).^2+Ai.*theta.*4./......
              beta(4:5).*(0.5./beta(4:5)-1).*((dp-1)^2).^(0.5./beta(4:5)-2)) ;
 
 ddelta_bi_dp = bi.*delta.^(bi-1).*ddelta_dp;
 d2delta_bi_dp = bi.*(delta.^(bi-1).*d2delta_dp+(bi-1).*delta.^(bi-2).*......
                (ddelta_dp).^2);
 
  i=1:7; j=8:51; k=52:54; m=55:56;
  
  
  phi_r1 = sum(ni(i).*di(i).*dp.^(di(i)-1).*tao.^ti(i))+......
           sum(ni(j).*exp(-dp.^ci(j)).*(dp.^(di(j)-1).*.......
               tao.^ti(j)).*(di(j)-ci(j).*dp.^(ci(j))))+......
           sum(ni(k).*dp.^di(k).*tao.^ti(k).*......
               exp(-alfa.*(dp-cigma).^2-beta(1:3).*(tao-gama).^2).*......
               (di(k)/dp-2*alfa.*(dp-cigma)))+......
           sum(ni(m).*(delta.^bi.*(thea+dp*dthea_dp)+ddelta_bi_dp*dp.*thea)); 
       
   phi_r2 = sum(ni(i).*di(i).*(di(i)-1).*dp.^(di(i)-2).*tao.^ti(i))+......
            sum(ni(j).*exp(-dp.^ci(j)).*(dp.^(di(j)-2).*tao.^ti(j).*......
               ((di(j)-ci(j).*dp.^ci(j)).*(di(j)-1-ci(j).*dp.^ci(j))-......
                ci(j).^2.*dp.^ci(j))))+......
            sum(ni(k).*tao.^ti(k).*exp(-alfa.*(dp-cigma).^2-beta(1:3).*......
                (tao-gama).^2).*(-2*alfa.*dp.^di(k)+4*alfa.^2.*dp.^di(k).*......
                (dp-cigma).^2-4*di(k).*alfa.*dp.^(di(k)-1).*(dp-cigma)+......
                 di(k).*(di(k)-1).*dp.^(di(k)-2)))+......
            sum(ni(m).*(delta.^bi.*(2*dthea_dp+dp.*d2thea_dp)+2*......
                ddelta_bi_dp.*(thea+dp*dthea_dp)+d2delta_bi_dp.*dp.*thea));
            
  y = -R*TG*(1+dp+2*dp*phi_r1+dp^2*phi_r2);  

function [density,errvec] = IAPWS95_dens(p,T,x)
% this function is used to calculate the density from pressure value using 
% Newton Iteration method

% declare global variables
global pG TG
pG = p;
TG = T;
x0 = x;   % starting point, in Kg.m-3
TOL = 1e-6;  % tolerance for declaration of convergence
DAMPED = 0;  % the standard undamped method
nmax = 100;  % maximum interation times
numeric = 0; % not use numerical differentiation
fun = @fn;
dfun = @dfn;

[x,errvec] = newton(fun, dfun, x0, TOL, DAMPED, nmax, numeric) ;
 
density = x;
        

function res = Gibbs08(ns,np,SA,t,p)

% this function is used to calculate the gibss function based on R.Feistel
% 2003 and R.Feistel,2008,Deep sea research
%
% ns                  : order of s derivative
% np                  : order of p derivative
% SA                  : Absolute Salinity                  [g/kg]
% t                   : temperature                        [deg C]
% p                   : sea (gauge) pressure               [dbar]
% 
% result              : specific Gibbs energy or its derivative

% Gibbs function coefficients
  gcoef = zeros(8,8,7);

  %% coefficients for seawater i>1. Table 17, Feistel 2008
  gcoef(2,1:2,1) = [5812.81456626732 851.226734946706];
  gcoef(2:8,1,1) = [5812.81456626732 1416.27648484197 -2432.14662381794 ......
                    2025.80115603697 -1091.66841042967 374.601237877840 ......
                    -48.5891069025409];
  gcoef(2:7,2,1) = [851.226734946706 168.072408311545 -493.407510141682 ......
                    543.835333000098 -196.028306689776 36.7571622995805];
  gcoef(3:5,3,1) = [880.031352997204 -43.0664675978042 -68.5572509204491]; 
  gcoef(3:5,4,1) = [-225.267649263401 -10.0227370861875 49.3667694856254];
  gcoef(3:5,5,1) = [91.4260447751259 0.875600661808945 -17.1397577419788];
  gcoef(3,6,1) = -21.6603240875311; 
  gcoef(5,6,1) = 2.49697009569508;
  gcoef(3,7,1) = 2.13016970847183; 
  gcoef(3:6,1,2) = [-3310.49154044839 199.459603073901 -54.7919133532887 .....
                    36.0284195611086];
  gcoef(3:5,2,2) = [729.116529735046 -175.292041186547 -22.6683558512829];
  gcoef(3:4,3,2) = [-860.764303783977 383.058066002476];
  gcoef(3:4,4,2) = [694.244814133268 -460.319931801257];
  gcoef(3:4,5,2) = [-297.728741987187 234.565187611355];
  gcoef(3:5,1,3) = [384.794152978599 -52.2940909281335 -4.08193978912261 ];
  gcoef(3:4,2,3) = [-343.956902961561 83.1923927801819];
  gcoef(3:4,3,3) = [337.409530269367 -54.1917262517112 ];
  gcoef(3,4:5,3) = [-204.889641964903 74.726141138756];
  gcoef(3:5,1,4) = [-96.5324320107458 68.0444942726459 -30.1755111971161];
  gcoef(3:4,2,4) = [124.687671116248 -29.483064349429];
  gcoef(3:4,3,4) = [-178.314556207638 25.6398487389914];
  gcoef(3,4:5,4) = [113.561697840594 -36.4872919001588];
  gcoef(3:4,1,5) = [15.8408172766824 -3.41251932441282 ];
  gcoef(3,2:4,5) = [-31.656964386073 44.2040358308 -11.1282734326413 ];
  gcoef(3,1:3,6) = [-2.62480156590992 7.04658803315449 -7.92001547211682];
  
  
  %% Gibbs function
  
  sfac = 0.0248826675584615;                % sfac = 1/(40*(35.16504/35)); 
  x2 = sfac.*SA;
  x = sqrt(x2);
  y = t.*0.025d0; 
  z = p.*1d-4;
  
  if ns==0 && np==0   
     gsw = x^2*log(x)*(gcoef(2,1,1)+y*gcoef(2,2,1));
     for i=2:7
         for j=0:6
           for k=0:5
               gsw = gsw+gcoef(i+1,j+1,k+1)*x.^i.*y.^j.*z.^k;
           end
         end
     end
     res = gsw;
  end
     
  if ns==0 && np==1
     gsw_p = zeros(1,length(t));
     for i=2:7
         for j=0:6
             for k=1:5
                 gsw_p = gsw_p+k*gcoef(i+1,j+1,k+1)*x.^i.*y.^j.*z.^(k-1);
             end
         end
     end
     res = 1d-8*gsw_p;     
  end
    
   % frist oder of salinity and pressure. gsp = dg/dsdp; unit is m3.g-1
   if ns==1 && np==1
       gsw_sp = zeros(1,length(t));
       for i=2:7
           for j=0:6
               for k=1:5
                   gsw_sp = gsw_sp+i*k*gcoef(i+1,j+1,k+1)*......
                            x.^(i-2).*y.^j.*z.^(k-1);
               end
           end
       end
       res = 0.5*sfac*1d-8*gsw_sp;
   end
   
  
  % second order of pressure derivatives. gpp = dg^2/dp^2.; Unit is
  % m3.Kg-1.Pa-1
  if ns==0 && np==2
     gsw_pp = zeros(1,length(t));
     for i=2:7
         for j=0:6
             for k=2:5
                 gsw_pp = gsw_pp+k*(k-1)*gcoef(i+1,j+1,k+1)*......
                          x.^i.*y.^j.*z.^(k-2);
             end
         end
     end
     res = 1d-16*gsw_pp;
  end

  % second order of salinity derivatives. gss = dg^2/d^2s;
  if ns==2 && np==0
     gss = 2*(gcoef(2,1,1)+gcoef(2,2,1)*y)/x^2;
     for i=4:7
       for j=0:6
         for k=0:5
             gss = gss+i*(i-2)*gcoef(i+1,j+1,k+1)*x.^(i-4).*y.^j.*z.^k;
         end
       end
     end
     res = sfac^2/4*gss;
  end
  


function res = mygsw(ns,np,SA,t,p)

% this routine is used to calculate the thermoaldynamic properties of 
% seawater based on IAPWS95 and IAPWS08. The pure water part is from
% IAPWS95,while the saline part is from IAPWS08
%
% ns                  : order of s derivative
% np                  : order of p derivative
% SA                  : Absolute Salinity                  [g/kg]
% t                   : temperature                        [deg C]
% p                   : sea (gauge) pressure               [dbar]
% 
% res                 : specific Gibbs energy or its derivative

% Reference Constant
  rouC = 322;       % Kg m-3
  R = 0.46151805d3; % J Kg-1 K-1
% Convert to SI unit
  Tk = t+273.15;   % in Kalvin
  P = p*1d4+101325;   % in Pa

 % Density of pure water
  x0 = 1000;
  rou = IAPWS95_dens(P,Tk,x0);
  dp = rou/rouC;
  
  [phi0 phir]= IAPWS95(0,0,Tk,rou);
  [phi0_d,phir_d] = IAPWS95(1,0,Tk,rou);
  [phi0_dd,phir_dd] = IAPWS95(2,0,Tk,rou);
   
  if ns<0 || np<0
      error('------- order of derivatives is negative.------');
  end
  
  if ns>0
      gw = 0;
  elseif  np==0
      fh = R*Tk*(phi0+phir);
      fhd = R*Tk*(1+dp*phir_d);
      gw = fh+fhd;
  elseif np==1
      gw = 1/rou;
  elseif np==2
      fhdd = R*Tk/rouC^2*(phi0_dd+phir_dd);
      fhd = P/rou^2;
      gw = -1/rou^3/(2*fhd+rou*fhdd);
  else
      error('The wrong input arguments!');
  end
  gsw = 0;
 if SA>0
    gsw = Gibbs08(ns,np,SA,t,p);
 end
 res = gw+gsw;
 
function [phi_0,phi_r]=IAPWS95(np,nt,T,rou)
% this routin is to implement the IAPWS95 calculating the first and second
% derivatives of Helmholtz free energy 
% 
% np   : order of density ratio derivatives 
% nt   : order of temperature ratio derivatives
% T    : temperature in Kelvin, Scalar
% rou  : density in Kg.m-3,  Scalar
% 
% Author: Lianbo Hu , May 12,2009
%         lianbo.hu@gmail.com

% Call coefficients and parameters used in Helmholtz free energy
[n0,gama0,ci,di,ti,ni,alfa,beta,gama,cigma,ai,Ai,bi,Bi,Ci,Di]=IAPWS95_Const;
% Reference Constants
Tc = 647.096;   % K
rouC = 322;     % Kg m-3
%   R = 0.46151805; % KJ Kg-1 K-1
  
% ------------------------------------------------------------------------
tao = Tc/T;
dp = rou/rouC;
theta = (1-tao)+Ai.*((dp-1)^2).^(0.5./beta(4:5));
thea = exp(-Ci.*(dp-1)^2-Di.*(tao-1)^2);
delta = theta.^2+Bi.*((dp-1).^2).^ai;


dthea_dp = -2*Ci.*(dp-1).*thea;
d2thea_dp = 2*(2*Ci.*(dp-1)^2-1).*Ci.*thea;
dthea_dt = -2*Di.*(tao-1).*thea;
d2thea_dt = 2*(2*Di.*(tao-1).^2-1).*Di.*thea;
d2thea_dpt = 4*Ci.*Di.*(dp-1).*(tao-1).*thea;

ddelta_dp = (dp-1)*(Ai.*theta*2./beta(4:5).*((dp-1)^2).^(0.5./beta(4:5)-1)+......
    2*Bi.*ai.*((dp-1)^2).^(ai-1));
d2delta_dp = 1/(dp-1).*ddelta_dp+(dp-1)^2.*(4*Bi.*ai.*(ai-1).*......
    ((dp-1)^2).^(ai-2)+2*Ai.^2.*(1./beta(4:5)).^2.*......
    (((dp-1)^2).^(0.5./beta(4:5)-1)).^2+Ai.*theta.*4./beta(4:5).*......
    (0.5./beta(4:5)-1).*((dp-1)^2).^(0.5./beta(4:5)-2)) ;

ddelta_bi_dp = bi.*delta.^(bi-1).*ddelta_dp;
d2delta_bi_dp = bi.*(delta.^(bi-1).*d2delta_dp+(bi-1).*delta.^(bi-2).*......
    (ddelta_dp).^2);
ddelta_bi_dt = -2*theta.*bi.*delta.^(bi-1);
d2delta_bi_dt = 2*bi.*delta.^(bi-1)+4*theta.^2.*bi.*(bi-1).*delta.^(bi-2);
d2delta_bi_dpt = -Ai.*bi.*2./beta(4:5).*delta.^(bi-1).*(dp-1).*((dp-1)^2).^......
    (0.5./beta(4:5)-1)-2*theta.*bi.*(bi-1).*delta.^(bi-2).*ddelta_dp;


i=1:7; j=8:51; k=52:54; m=55:56; n=4:8;

% Helmholtz free energy
if np==0 && nt==0
    phi_0 = log(dp)+n0(1)+n0(2)*tao+n0(3)*log(tao)+......
        sum(n0(n).*log(1-exp(-gama0(n)*tao)));
    phi_r = sum(ni(i).*dp.^di(i).*tao.^ti(i))+......
        sum(ni(j).*dp.^di(j).*tao.^ti(j).*exp(-dp.^ci(j)))+......
        sum(ni(k).*dp.^di(k).*tao.^ti(k).*......
        exp(-alfa.*(dp-cigma).^2-beta(1:3).*(tao-gama).^2))+......
        sum(ni(m).*delta.^bi.*dp.*thea);
    
    
end

% first order derivatives of dp
if np==1 && nt==0
    phi_0 = 1/dp;
    
    phi_r = sum(ni(i).*di(i).*dp.^(di(i)-1).*tao.^ti(i))+......
        sum(ni(j).*exp(-dp.^ci(j)).*(dp.^(di(j)-1).*.......
        tao.^ti(j).*(di(j)-ci(j).*dp.^ci(j))))+......
        sum(ni(k).*dp.^di(k).*tao.^ti(k).*......
        exp(-alfa.*(dp-cigma).^2-beta(1:3).*(tao-gama).^2).*......
        (di(k)/dp-2*alfa.*(dp-cigma)))+......
        sum(ni(m).*(delta.^bi.*(thea+dp*dthea_dp)+ddelta_bi_dp*dp.*thea));
    
end

if np==2 && nt==0
    phi_0 = -1/dp^2;
    
    phi_r = sum(ni(i).*di(i).*(di(i)-1).*dp.^(di(i)-2).*tao.^ti(i))+......
        sum(ni(j).*exp(-dp.^ci(j)).*(dp.^(di(j)-2).*tao.^ti(j).*......
        ((di(j)-ci(j).*dp.^ci(j)).*(di(j)-1-ci(j).*dp.^ci(j))-......
        ci(j).^2.*dp.^ci(j))))+......
        sum(ni(k).*tao.^ti(k).*exp(-alfa.*(dp-cigma).^2-beta(1:3).*......
        (tao-gama).^2).*(-2*alfa.*dp.^di(k)+4*alfa.^2.*dp.^di(k).*......
        (dp-cigma).^2-4*di(k).*alfa.*dp.^(di(k)-1).*(dp-cigma)+......
        di(k).*(di(k)-1).*dp.^(di(k)-2)))+......
        sum(ni(m).*(delta.^bi.*(2*dthea_dp+dp.*d2thea_dp)+2*ddelta_bi_dp.*......
        (thea+dp*dthea_dp)+d2delta_bi_dp.*dp.*thea));
    
end

if np==0 && nt==1
    phi_0 = n0(2)+n0(3)/tao+sum(n0(n).*gama0(n).*((1-exp(-gama0(n)*tao)).^(-1)-1));
    phi_r = sum(ni(i).*ti(i).*dp.^di(i).*tao.^(ti(i)-1))+......
        sum(ni(j).*ti(j).*dp.^di(j).*tao.^(ti(j)-1).*exp(-dp.^ci(j)))+......
        sum(ni(k).*dp.^di(k).*tao.^ti(k).*exp(-alfa.*(dp-cigma).^2-......
        beta(1:3).*(tao-gama).^2).*(ti(k)/tao-2*beta(1:3).*......
        (tao-gama)))+......
        sum(ni(m).*dp.*(ddelta_bi_dt.*thea+delta.^bi.*dthea_dt));
    
end

if np==0 && nt==2
    phi_0 = -n0(3)/tao^2-sum(n0(n).*(gama0(n)).^2.*exp(-gama0(n)*tao).*......
        (1-exp(-gama0(n)*tao)).^(-2));
    phi_r = sum(ni(i).*ti(i).*(ti(i)-1).*dp.^di(i).*tao.^(ti(i)-2))+......
        sum(ni(j).*ti(j).*(ti(j)-1).*dp.^di(j).*tao.^(ti(j)-2).*......
        exp(-dp.^ci(j)))+......
        sum(ni(k).*dp.^di(k).*tao.^ti(k).*exp(-alfa.*(dp-cigma).^2-......
        beta(1:3).*(tao-gama).^2).*((ti(k)/tao-2*beta(1:3).*......
        (tao-gama)).^2-ti(k)./tao^2-2*beta(1:3)))+......
        sum(ni(m).*dp.*(d2delta_bi_dt.*thea+2*ddelta_bi_dt.*dthea_dt+......
        delta.^bi.*d2thea_dt));
    
    
end

if np==1 && nt==1
    phi_0 = 0;
    phi_r = sum(ni(i).*di(i).*ti(i).*dp.^(di(i)-1).*tao.^(ti(i)-1))+......
        sum(ni(j).*ti(j).*dp.^(di(j)-1).*tao.^(ti(j)-1).*......
        (di(j)-ci(j).*dp.^ci(j)).*exp(-dp.^ci(j)))+......
        sum(ni(k).*dp.^di(k).*tao.^ti(k).*exp(-alfa.*(dp-cigma).^2-......
        beta(1:3).*(tao-gama).^2).*(di(k)/dp-2*alfa.*(dp-cigma)).*......
        (ti(k)/tao-2*beta(1:3).*(tao-gama)))+......
        sum(ni(m).*(delta.^bi.*(dthea_dt+dp*d2thea_dpt)+dp*......
        ddelta_bi_dp.*dthea_dt+ddelta_bi_dt.*(thea+dp*dthea_dp)+......
        d2delta_bi_dpt.*dp.*thea));
end


return;
  