计算海水密度

根据海水压力，温度，盐度计算海水密度

注：此计算根据联合国1983年文件得出，所以其中温标采用的是68年温标
P（s,t,0）是一个标准大气压下的数据，海面压强设为0

//根据盐度S，温度T，压力P算出海水密度
        public static double getDensity(double S,double T,double P) {

             
            //p(s,t,p)=p(s,t,0)*[1-p/k(s,t,p)]

            //************************************8首先算p(s,t,0) 设为Pst0**********************************************
            //先算Pw
            double a0 = 999.842594;//***是什么意思
            double a1 = 6.793952 * Math.Pow(10,-2);
            double a2 = -9.095290 * Math.Pow(10,-3);
            double a3 = 1.001685 * Math.Pow(10, -4);
            double a4 = -1.120083 * Math.Pow(10,-6);
            double a5 = 6.536332 * Math.Pow(10,-9);
            double Pw = a0 + a1 * T + a2 * T * T + a3 * Math.Pow(T, 3) + a4 * Math.Pow(T, 4) + a5 * Math.Pow(T, 5);

            //再算Pst0
            double b0 = 8.24493 * Math.Pow(10, -1);
            double b1 = -4.0899 * Math.Pow(10, -3);
            double b2 = 7.6438 * Math.Pow(10, -5);
            double b3 = -8.2467 * Math.Pow(10, -7);
            double b4 = 5.3875 * Math.Pow(10, -9);
            double c0 = -5.72466 * Math.Pow(10, -3);
            double c1 = 1.0227 * Math.Pow(10, -4);
            double c2 = -1.6546 * Math.Pow(10, -6);
            double d0 = 4.8314 * Math.Pow(10, -4);
            double Pst0 = Pw + (b0 + b1 * T + b2 * Math.Pow(T, 2) + b3 * Math.Pow(T, 3) + b4 * Math.Pow(T, 4)) * S + (c0 + c1 * T + c2 * Math.Pow(T, 2)) * S*Math.Sqrt(S) + d0 * S * S;


            //算K(s,t,p) K(s,t,p)=K(s,t,0) +A*P +B*P*P    设为Kstp
            //Aw
            double h0 = 3.239908;//带有******什么意思
            double h1 = 1.43713 * Math.Pow(10, -3);
            double h2 = 1.16092 * Math.Pow(10, -4);
            double h3 = -5.77905 * Math.Pow(10, -7);
            double Aw = h0 + h1 * T + h2 * T * T + h3 * T * T * T;
            //Bw
            double k0 = 8.50935*Math.Pow(10,-5);//带有****不知道什么意思
            double k1 = -6.12293 * Math.Pow(10, -6);
            double k2 = 5.2787 * Math.Pow(10,-8);
            double Bw = k0 + k1 * T + k2 * T * T;
            //Kw
            double e0 = 19652.21;//带有**不知道什么意思
            double e1 = 148.4206;
            double e2 = -2.327105;
            double e3 = 1.360477 * Math.Pow(10, -2);
            double e4 = -5.155288 * Math.Pow(10,-5);
            double Kw = e0 + e1 * T + e2 * T * T + e3 * Math.Pow(T, 3) + e4 * Math.Pow(T, 4);
            //B
            double m0 = -9.9348 * Math.Pow(10, -7);
            double m1 = 2.0816 * Math.Pow(10,-8);
            double m2 = 9.1697 * Math.Pow(10,-10);
            double B = Bw + (m0 + m1 * T + m2 * T * T) * S;
            //A
            double i0 = 2.2838 * Math.Pow(10, -3);
            double i1 = -1.0981 * Math.Pow(10, -5);
            double i2 = -1.6078 * Math.Pow(10,-6);
            double j0 = 1.91075 * Math.Pow(10,-4);
            double A = Aw + (i0 + i1 * T + i2 * T * T) * S + j0* S * Math.Sqrt(S);
            //K(s,t,0)
            double f0 = 54.6746;
            double f1 = -0.603459;
            double f2 = 1.09987 * Math.Pow(10, -2);
            double f3 = -6.1670 * Math.Pow(10, -5);
            double g0 = 7.944 * Math.Pow(10,-2);
            double g1 = 1.6483 * Math.Pow(10,-2);
            double g2 = -5.3009 * Math.Pow(10,-4);
            double Kst0 = Kw + (f0 + f1 * T + f2 * T * T + f3 * Math.Pow(T, 3)) * S + (g0 + g1 * T + g2 * T * T) * S * Math.Sqrt(S);
            //Kstp
            double Kstp = Kst0 + A * P + B * P * P;

            //Pstp //p(s,t,p)=p(s,t,0)/[1-p/k(s,t,p)]
            double Pstp = Pst0 /(1-(P/Kstp));

            return Pstp;
        }

这是文件中的验证数据，我把p/10代进去得出的数据才于文件数据一致。
例如：getDensity（35，25，1000）=1062.53817
而getDensity（35，25，10000）得出的数据不对，反正除10带入就对了