-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathr_ps.m
127 lines (110 loc) · 4.38 KB
/
r_ps.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
function[r]=r_ps(p,s)
% Determines density as a function or pressure and entropy (old tgas5.f)
p0 = 101330.0;
r0 = 1.292;
s0 = 6779.2;
gascon = 287.6;
x = log(p/p0)/log(10);
snon = log(s/gascon)/log(10);
z = x - snon;
if (snon <= 1.23)
delts = s - s0;
rratio = (log(p/p0))/(log(10.))/1.4-delts/(3.5*gascon);
r = r0*exp(rratio);
return;
end
if (snon <= 1.42)
y = -17.2119+54.9354*snon;
y = y+(-1.99776+3.17884*snon)*z;
y = y+(-46.9831-0.866580*z+12.1069*snon)*snon*snon;
y = y+(0.158567-0.103055*snon-0.00152322*z)*z*z;
r = (10^y)*r0;
return;
end
if (snon <= 1.592 )
y = 280.393-595.834*snon+(-17.5934+24.9706*snon)*z+...
(421.767-8.85461*z-99.4515*snon)*snon*snon+...
(0.36189-0.255458*snon-0.00296892*z)*z*z;
y = y/(1.0+exp(-15.0*(x+54.179-86.947*snon+33.583*snon*snon)));
y = y+(-278.074+611.791*snon+(13.7528-19.2394*snon)*z+...
(-442.909+7.10425*z+105.869*snon)*snon*snon+...
(-0.197269+0.149708*snon+0.00119153*z)*z*z);
r = (10^y)*r0;
return;
end
if (snon<=1.70)
zm = 7.5269*snon-14.9366;
if (z>=zm)
y = -50.1859+63.1564*snon+(2.39925-1.47883*snon)*z...
-19.8852*snon*snon+(0.466791-0.306926*snon)*z*z;
y = y/(1.0+exp(-2.0*(z+11.16*snon-18.203)));
y = y + (39.7149-48.0988*snon+(-1.17821+1.27522*snon)...
*z+14.5552*snon*snon+(-0.250279+0.158399*snon)*z*z);
else
y = -68.5292+82.3834*snon+(-1.24942+0.831615*snon)*...
z-24.524*snon*snon+(-0.113019+0.0746719*snon)*z*z;
y = y/(1.0+exp(-2.0*(z+33.976*snon-49.659)));
y = y+(110.732-133.968*snon+(0.37583+0.277887*snon)*z+...
40.3018*snon*snon+(0.118506-0.0698812*snon)*z*z);
end
r = (10^y)*r0;
return;
end
if (snon<=1.80)
zm = 7.5269*snon-14.9366;
if(z<zm)
y = -72.8463+94.2375*snon+(4.04979-2.45467*snon)*...
z-30.2865*snon*snon+(0.12872-0.0863721*snon)*z*z;
y = y/(1.0+exp(-2.0*(z+13.4576*snon-18.36)));
y =y+(62.2639-78.884*snon+(-2.19123+1.84533*snon)*...
z+24.8413*snon*snon+(-0.0876105+0.0555047*snon)*z*z);
else
y = 60.5677-70.8307*snon+(-2.95622+1.70504*snon)*...
z+20.7549*snon*snon+(-0.270962+0.180057*snon)*z*z;
y = y/(1.0+exp(-2.0*(z+16.822*snon-29.139)));
y = y+(-29.9578+38.9998*snon+(2.90256-1.25088*snon)*...
z-12.6641*snon*snon+(0.0613867-0.0517059*snon)*z*z);
end
r=(10^y)*r0;
return;
end
if(snon<=1.90)
y=-639.514+697.458*snon+(-100.154+106.701*snon)*...
z+(-190.745-28.6323*z)*snon*snon+(-0.481471+...
0.200348*snon-0.00643371*z)*z*z;
y=y/(1.0+exp(-2.0*(z+35.275*snon-58.624)));
y=y+(623.124-677.571*snon+(91.2811-95.548*snon)*...
z+(184.603+25.4274*z)*snon*snon+(0.697635-0.327916*...
snon+0.00490838*z)*z*z);
r=(10^y)*r0;
return;
end
if(snon<=2.00)
y=-113.217+124.304*snon+(4.41505-2.4853*snon)*...
z-34.2370*snon*snon+(-0.0784297+0.0121459*snon)*z*z;
y=y/(1.0+exp(-2.00*(z+20.884*snon-36.54)));
y=y+(14.0088-15.8855*snon+(0.171245+0.432352*snon)*...
z+4.42836*snon*snon+(-0.00954417+0.00892335*snon)*z*z);
r=(10^y)*r0;
return;
end
if(snon<=2.10)
y=(65.2920-62.8154*snon)+(2.10906-1.11759*snon)*z+15.0026*...
snon*snon+(0.271716-0.14431*snon)*z*z;
y=y/(1.0+exp(-5.00*(z+28.284*snon-53.185)));
y=y+(-36.2767+38.634*snon+(1.48507-0.341824*snon)*...
z-10.4718*snon*snon+(-0.0325437+0.0122032*snon)*z*z);
r=(10^y)*r0;
return;
end
if(snon>2.10)
y = -127.867+85.8453*snon+(-22.8808+11.3034*snon)*z-11.234*...
snon*snon+(-1.56005+0.779397*snon)*z*z;
y=y/(1.0+exp(-2.00*(z+15.048*snon-26.307)));
y=y+(186.938-145.261*snon+(22.3883-10.4974*snon)*z...
+26.0284*snon*snon+(1.85949-0.898218*snon)*z*z);
r = (10^y)*r0;
return;
end
return;
end