001 function rho = wdensity(S,T,P)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036 error(nargchk(0,3,nargin))
037 if nargin<1|isempty(S), S = 35; end
038 if nargin<2|isempty(T), T = 4; end
039 if nargin<3|isempty(P), P = 0; end
040 [ec, S,T,P] = comnsize(S,T,P);
041 if ec>0,
042 error('S, T and P must be scalar or of common size')
043 end
044
045
046
047 rho0 = 999.842594;
048 rho = rho0(ones(size(S)));
049
050 k = find(T~=0);
051 if any(k),
052 Tk = T(k);
053 a = [6.536332e-9, -1.120083e-6, 1.001685e-4, -9.095290e-3, 6.793952e-2];
054 rho(k) = rho(k)+(a(5)+(a(4)+(a(3)+(a(2)+a(1)*Tk).*Tk).*Tk).*Tk).*Tk;
055
056 end
057
058
059
060 k1 = find(S~=0);
061 if any(k1),
062 c3 = -5.72466e-3;
063 b5 = 8.24493e-1;
064 tmp = b5(ones(size(k1)));
065 tmp1 = c3(ones(size(k1)));
066 Tk = T(k1);
067 k12 = find(Tk~=0);
068 if any(k12),
069 Tk = Tk(k12);
070 b = [ 5.3875e-9, -8.2467e-7, 7.6438e-5, -4.0899e-3, 0];
071 c = [ -1.6546e-6, 1.0227e-4, 0];
072 tmp(k12) = tmp(k12)+(b(4) + (b(3) + (b(2) + b(1)*Tk).*Tk).*Tk).*Tk;
073 tmp1(k12) = tmp1(k12)+(c(2) + c(1)*Tk).*Tk;
074 clear Tk
075 end
076 d = 4.8314e-4;
077 Sk = S(k1);
078 rho(k1) = rho(k1) + (tmp + tmp1.*sqrt(Sk) + d*Sk).*Sk;
079 end
080
081 k2 = find(P~=0);
082 if any(k2),
083 Pk = P(k2);
084 K = seck(S(k2),T(k2),Pk);
085 Pk = Pk/10;
086 rho(k2) = rho(k2)./(1-Pk./K);
087 end
088 return
089
090
091 function K = seck(S,T,P)
092
093
094
095
096
097
098
099
100
101
102
103
104
105
106 error(nargchk(3,3,nargin))
107
108
109
110
111
112 K0 = 19652.21;
113 K = K0(ones(size(S)));
114
115 k1 = find(T~=0);
116 if any(k1),
117 Tk = T(k1);
118 e = [-5.155288E-5, 1.360477E-2, -2.327105, 148.4206,0];
119 K(k1) = K(k1) + (e(4) + (e(3) + (e(2) + e(1)*Tk).*Tk).*Tk).*Tk;
120 end
121
122
123
124
125
126
127 SR = sqrt(S);
128 k1 = find(S~=0);
129 if any(k1),
130 Tk = T(k1);
131 f = [-6.1670E-5, 1.09987E-2,-0.603459,54.6746];
132 g = [-5.3009E-4, 1.6483E-2,7.944E-2];
133
134 K(k1) = K(k1) + ( f(4) + (f(3) + (f(2) + f(1)*Tk).*Tk).*Tk ...
135 + (g(3) + (g(2) + g(1)*Tk).*Tk).*SR(k1)).*S(k1);
136 end
137
138
139 k2 = find(P~=0);
140 if any(k2),
141 A0 = 3.239908; B0 = 8.50935E-5;
142 A = A0(ones(size(k2)));
143 B = B0(ones(size(k2)));
144 Tk = T(k2);
145
146 k12 = find(Tk~=0);
147 if any(k12)
148 Tk1 = Tk(k12);
149 k12 = k(k12);
150 h = [-5.77905E-7, 1.16092E-4, 1.43713E-3, 0];
151 A(k12) = A(k12) + (h(3) + (h(2) + h(1).*Tk1).*Tk1).*Tk1;
152
153 b = [ 5.2787E-8, -6.12293E-6,0];
154 B(k12) = B(k12) + (b(2) + b(1)*Tk1).*Tk1;
155 end
156 J = 1.91075E-4;
157 I =[ -1.6078E-6, -1.0981E-5 2.2838E-3];
158 A = A + (I(3) + (I(2) + I(1)*Tk).*Tk + J*SR(k2)).*S(k2);
159
160 m = [ 9.1697E-10, 2.0816E-8, -9.9348E-7];
161 B = B + (m(3) + (m(2) + m(1)*Tk).*Tk).*S(k2);
162
163 Pk = P(k2)/10;
164 K(k2) = K(k2) + (A + B.*Pk).*Pk;
165 end
166 return
167
168
169
170
171
172
173
174
175
176