function sigma = sigma_p(t,s,p)

% SIGMA = SIGMA_P(T,S,P) computes sigma (in situ density minus 1000) 
% knowing p (db),t (deg C), and s (psu).
% Coefficients from Millero & Poisson, DSR 28 (6a),625 (1981).
%
% checkvalues: s (psu)  t (deg C)  p (dB)  sigma +1000 (kg m^-3)
%                0         5         0         999.96675
%                0         5       10000      1044.12802
%                0         25        0         997.04796
%                0         25      10000      1037.90204
%                35        5         0        1027.67547
%                35        5       10000      1069.48914
%                35        25        0        1023.34306
%                35        25      10000      1062.53817
%
% Converted to a matlab function from a fortran subroutine
% June 30, 1993 by TDM.
% 2018-04-19, RGL ammended description by indication density - 1000
% Function checked with their check values.

a00=-0.157406; a10=6.793952e-02; a20=-9.095290e-03; a30=1.001685e-04;
a40=-1.120083e-06; a50=6.536336e-09; a01=0.824493; a11=-4.0899e-03;
a21=7.6438e-05; a31=-8.2467e-07; a41=5.3875e-09; a02=-5.72466e-03;
a12=1.02270e-04; a22=-1.65460e-06; a03=4.8314e-04;

e00=19652.21; e10=148.4206; e20=-2.327105; e30=1.360477e-02;
e40=-5.155288e-5; e01=54.6746; e11=-.603459; e21=1.09987e-2;
e31=-6.167e-5; e02=7.944e-02; e12=1.6483e-02; e22=-5.3009e-04;

h00=3.239908; h10=1.43713e-03; h20=1.16092e-04;
h30=-5.77905e-07; h01=2.2838e-03; h11=-1.0981e-5; h21=-1.6078e-06;
h02=1.91075e-04;

m00=8.50935e-05; m10=-6.12293e-6; m20=5.2787e-08;
m01=-9.9348e-07; m11=2.0816e-08; m21=9.1697e-10;

p       = p/10.;

sigma0  = a00 + a10*t + a20*t.^2 + a30*t.^3 + a40*t.^4 + a50*t.^5 ...
        + s.*(a01 + a11*t + a21*t.^2 + a31*t.^3 + a41*t.^4) ...
        + sqrt(s.^3).*(a02 + a12*t + a22*t.^2) ...
        + s.^2*a03;

k       = e00 + e10*t + e20*t.^2 + e30*t.^3 + e40*t.^4 ... 
        + s.*(e01 + e11*t + e21*t.^2 + e31*t.^3) ...
        + sqrt(s.^3).*(e02 + e12*t + e22*t.^2) ...
        + (h00 + h10*t + h20*t.^2 + h30*t.^3 ...
        + s.*(h01 + h11*t + h21*t.^2) ...
        + sqrt(s.^3)*h02).*p ...
        + (m00 + m10*t + m20*t.^2 ...
        + s.*(m01 + m11*t + m21*t.^2)).*p.*p;

sigma   = (sigma0 + 1000*p./k)./(1. - p./k);
