-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy patha_ps.m
87 lines (77 loc) · 3.64 KB
/
a_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
function[a]=a_ps(p,s)
% speed of sound as a function of pressure and entropy (old tgas7.f)
p0 = 101330;
s0 = 6779.2;
r0 = 1.292;
a0 = 331.3613;
gascon = 287.06;
x=(log(p/p0))/(log(10));
snon=(log(s/gascon))/(log(10));
z = x - snon;
if(snon < 1.23)
delts = s-s0;
asqlog = log(1.4*p0/r0)+(log(p/p0)+delts/gascon)/3.5;
a = exp(asqlog/2.0);
return
elseif(snon <= 1.4)
y = (-.138377-8.84138*snon) + ((2.6105-3.16535*snon)*z) +...
((11.0866+.988389*z-3.25761*snon)*snon*snon) + ...
((-.100224+.0662193*snon+.00082061*z)*z*z);
elseif(snon < 1.595)
y = (131.057-288.847*snon) + ((-5.04887+7.73862*snon)*z) +...
((210.147-2.88963*z-50.396*snon)*snon*snon) + ...
((.0548031-.0439459*snon-.000210202*z)*z*z) + ...
( (-133.465+284.739*snon) + ((7.57389-10.7749*snon)*z) ...
+ ((-202.362+3.8313*z+47.9075*snon)*snon*snon) + ...
((-.153453+.108531*snon+.00097931*z)*z*z))/(1.0+...
(exp(-15.0*(x+54.179-86.947*snon+33.583*snon*snon))));
elseif(snon < 1.693)
zm = -9.842*snon+14.19;
if(z <= zm)
y = (-61.3548+78.0742*snon) + ((2.08524-1.21609*snon)*z) ...
+ (-24.3686*snon*snon) + ((.0877563-.0546311*snon)*z*z) ...
+ ( (20.7952-27.1591*snon) + ((-.743673+.408312*snon)*z) ...
+ (8.68124*snon*snon) + ((-.093592+.0532328*snon)*z*z))/...
(1.0+(exp(-2.0*(z+38.785*snon-57.157))));
else
y = (3.37056-4.87016*snon) + ((-.385754+.287192*snon)*z) + ...
(2.02041*snon*snon) + ((-.00463144+.00830832*snon)*z*z);
end
elseif(snon <1.80)
zm = -1.917*snon+0.092;
if( z <= zm)
y = (-80.4927+76.3739*snon) + ((-9.381+5.72104*snon)*z) + ...
( -16.3435*snon*snon) + ((-.748578+.450043*snon)*z*z) + ...
((83.4054-85.8837*snon) + ((4.84197-3.11188*snon)*z) + ...
(21.1196*snon*snon) + ((.233945-.159099*snon)*z*z))/(1.0+...
(exp(-2.0*(z+7.874*snon-7.569))));
else
y = (-4.73308+4.69363*snon) + ((.0943798+.00354953*snon)*z) +...
(-.798293*snon*snon) + ((.0202561-.00873036*snon)*z*z);
end
elseif(snon < 1.9)
y = (-660.574+738.042*snon) + ((-87.7589+97.894*snon)*z) + ...
((-206.156-27.3753*z)*snon*snon) + (-.0214028*z*z) + ...
((665.014-743.416*snon) + ((88.7679-99.0508*snon)*z) + ...
((208.117+27.7374*z)*snon*snon) + (.0281148*z*z))/(1.0+...
(exp(-2.0*(z+35.275*snon-58.624))));
elseif(snon < 2.0)
y = (-5.93554-7.79929*snon) + ((-7.23618+3.31162*snon)*z) + ...
(5.06381*snon*snon) + ((-.5377735+.246865*snon)*z*z) + ...
((39.526-29.0994*snon) + ((6.24136-2.70007*snon)*z) + ...
(5.42786*snon*snon) + ((.595235-.267771*snon)*z*z))/(1.0+...
(exp(-2.0*(z+13.6751*snon-20.1676))));
elseif(snon < 2.1)
y = (70.2453-73.4732*snon) + (-.586844*z) + ((19.5548+.181521*z)*...
snon*snon) + (.00736786*z*z) + ((-16.961+18.0169*snon) + ...
( -.365779*z) + ((-4.75063+.107583*z)*snon*snon) + ...
(.00201803*z*z))/(1.0+(exp(-5.0*(z+28.284*snon-53.185))));
else
y = (-112.793+135.296*snon) + (2.74295*z) + ((-40.8171-1.29257*z)*...
snon*snon) + (-.225652*z*z) + ((257.057-283.199*snon) + ...
( -5.84656*z) + ((79.1574+2.17312*z)*snon*snon) + ...
(.335658*z*z))/(1.0+(exp(-2.0*(z+15.048*snon-26.307))));
end
a =(10^y)*a0;
return;
end