function r0 = density(t,s,p) 
%function r0 = density(t,s,p) 
% computes density in sigma units given
%       t = degrees celsius
%       s = psu
%       p = decibars
% from UNESCO 1980 eqn. of state
                                              
      t2 = t.*t;                                                                  
      t3 = t2.*t;                                                                  
      t4 = t3.*t;                                                                
      t5 = t4.*t;
                                                                  
      s15 = s.^(1.5);                                                               
      s2 = s.*s;                                                                 
       p = p/10;   % now in bars
      p2 = p.*p;                                                                  
%----------------------------------------------------------------------         
% EACH COLUMN CONTAINS COEFFICIENTS OF 1, T, T**2, T**3, T**4, T**5;            
% EACH ROW CONTAINS COEFFICIENTS OF:                                            
%                    1             S            S**3/2        S**2              
%-----------------------------------------------------------------------        
      sigt1=    -0.157406    +8.24493e-1*s  -5.72466e-3*s15 +4.8314e-4*s2 ...    
          + t.*( 6.793952e-2  -4.0899e-3*s  +1.0227e-4*s15);             
      sigt2= t2.*(-9.095290e-3  +7.6438e-5*s  -1.6546e-6*s15)  ...       
            +t3.*( 1.001685e-4  -8.2467e-7*s);       
      sigt3= t4.*(-1.120083e-6  +5.3875e-9*s) +t5.*( 6.536332e-9); 
sigt = sigt1 + sigt2 + sigt3;                                                  
                                                                             
      cko=       1.965221e+4  +5.46746e+1*s  +7.944e-2 *s15    ...                 
          + t.*( 1.484206e+2  -6.03459e-1*s  +1.6483e-2*s15);                 
      ck1 =t2.*(-2.327105e0   +1.09987e-2*s  -5.3009e-4*s15)   ...                 
          +t3.*( 1.360477e-2  -6.1670e-5 *s) +t4*(-5.155288e-5); 
cko = cko + ck1;                                                          
                                                                             
      ckp=  p.*(3.239908      +2.2838E-3*s  +1.91075E-4*s15   ...                 
          + t.*(1.43713E-3    -1.0981E-5*s)                   ...                 
          +t2.*(1.16092E-4    -1.6078E-6*s)                   ...                 
          +t3.*(-5.77905E-7));                                                 
                                                                              
      ckp2=p2.*(8.50935e-5    -9.9348e-7*s      ...                              
          + t.*(-6.12293e-6  +2.0816e-8 *s)     ...                               
          +t2.*(5.2787e-8    +9.1697e-10*s));  
ck = cko + ckp + ckp2;                                 
%                                                                               
% NOW EVALUATE RO(P)= RO(0)/(1-P/CK),             
% AND RO= 1000.+SIGT IN KG/M**3.                                                
%                                                                               
      r0=(1000.+sigt)./(1.-p./ck);                                    
      r0=r0-1000;                                                               
