-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgsw_specvol_from_pot_enthalpy_ice_poly_work.m
115 lines (97 loc) · 3.85 KB
/
gsw_specvol_from_pot_enthalpy_ice_poly_work.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
pot_h_ices = pot_h_ice.*1e-5;
ps = p.*1e-4;
v00 = 1.090843882189585e-3;
v0_fit =[-2.993503247981132e-05
-2.187697844100046e-05
-6.597339965467078e-06
-1.092035946035247e-06
-9.164174131511919e-08
-2.874852704079762e-09];
v_fit = [ -1.558856898604863e-05
-1.323880019968672e-06
-2.685780994854041e-07
-6.212785565052150e-08
-4.797551830947803e-09
2.953039321948178e-07
-2.403547886994993e-08
-1.461069600352674e-08];
v_poly = v00 + v0_fit(1).*pot_h_ices + v0_fit(2).*pot_h_ices.^2 ...
+ v0_fit(3).*pot_h_ices.^3 + v0_fit(4).*pot_h_ices.^4 ...
+ v0_fit(5).*pot_h_ices.^5 + v0_fit(6).*pot_h_ices.^6 ...
+ v_fit(1).*ps ...
+ v_fit(2).*pot_h_ices.*ps + v_fit(3).*pot_h_ices.^2.*ps ...
+ v_fit(4).*pot_h_ices.^3.*ps + v_fit(5).*pot_h_ices.^4.*ps ...
+ v_fit(6).*ps.^2 ...
+ v_fit(7).*pot_h_ices.*ps.^2 ...
+ v_fit(8).*ps.^3;
v_h_Ih_poly = v01 + 2.*v02.*pot_h_ices ...
+ 3.*v03.*pot_h_ices.^2 + 4.*v04.*pot_h_ices.^3 ...
+ 5.*v05.*pot_h_ices.^4 + 6.*v06.*pot_h_ices.^5 ...
+ v08.*ps + 2.*v09.*pot_h_ices.*ps ...
+ 3.*v10.*pot_h_ices.^2.*ps + 4.*v11.*pot_h_ices.^3.*ps ...
+ v13.*ps.^2;
v00 = 1.090843882189585e-3;
v01 = -2.993503247981132e-05*(1e-5);
v02 = -2.187697844100046e-05*(1e-5).^2;
v03 = -6.597339965467078e-06*(1e-5).^3;
v04 = -1.092035946035247e-06*(1e-5).^4;
v05 = -9.164174131511919e-08*(1e-5).^5;
v06 = -2.874852704079762e-09*(1e-5).^6;
v07 = -1.558856898604863e-05*(1e-4);
v08 = -1.323880019968672e-06*(1e-4)*(1e-5);
v09 = -2.685780994854041e-07*(1e-4)*(1e-5).^2;
v10 = -6.212785565052150e-08*(1e-4)*(1e-5).^3;
v11 = -4.797551830947803e-09*(1e-4)*(1e-5).^4;
v12 = 2.953039321948178e-07*(1e-4).^2;
v13 = -2.403547886994993e-08*(1e-4).^2*(1e-5);
v14 = -1.461069600352674e-08*(1e-4).^3;
v_poly2 = v00 + pot_h_ice.*(v01 + pot_h_ice.*(v02 + pot_h_ice.*(v03 ...
+ pot_h_ice.*(v04 + pot_h_ice.*(v05 + v06.*pot_h_ice))))) + p.*(v07 ...
+ pot_h_ice.*(v08 + pot_h_ice.*(v09 + pot_h_ice.*(v10 + v11.*pot_h_ice))) ...
+ p.*(v12 + v13.*pot_h_ice + v14.*p));
plot(v_poly(:) -v_poly2(:))
h00 = -3.333548730778702e5;
h0_fit =[ -1.158533659785955e+11
3.088110638980851e+11
-1.574635299523583e+11
-2.315779173662935e+11
2.790387951606607e+11
-8.296646858006085e+10];
h_fit =[2.315408333461980e+09
-4.021450585146860e+09
-5.705978545141118e+08
4.048693790458912e+09
-1.769028853661338e+09
-2.307016488735996e+06
2.135637957175861e+06
8.593828063620491e+03];
h_poly = h00 ...
+ h0_fit(1).*(v_ices) + h0_fit(2).*(v_ices.^2) ...
+ h0_fit(3).*(v_ices.^3) + h0_fit(4).*(v_ices.^4) ...
+ h0_fit(5).*(v_ices.^5) + h0_fit(6).*(v_ices.^6) ...
+ h_fit(1).*(ps) ...
+ h_fit(2).*(v_ices.*ps) + h_fit(3).*(v_ices.^2.*ps) ...
+ h_fit(4).*(v_ices.^3.*ps) + h_fit(5).*(v_ices.^4.*ps) ...
+ h_fit(6).*(ps.^2) ...
+ h_fit(7).*(v_ices.*ps.^2) ...
+ h_fit(8).*(ps.^3);
h00 = -3.333548730778702e5;
h01 = -1.158533659785955e+11*(1e3);
h02 = 3.088110638980851e+11*(1e3).^2;
h03 = -1.574635299523583e+11*(1e3).^3;
h04 = -2.315779173662935e+11*(1e3).^4;
h05 = 2.790387951606607e+11*(1e3).^5;
h06 = -8.296646858006085e+10*(1e3).^6;
h07 = 2.315408333461980e+09*(1e-4);
h08 = -4.021450585146860e+09*(1e-4)*(1e3);
h09 = -5.705978545141118e+08*(1e-4)*(1e3).^2;
h10 = 4.048693790458912e+09*(1e-4)*(1e3).^3;
h11 = -1.769028853661338e+09*(1e-4)*(1e3).^4;
h12 = -2.307016488735996e+06*(1e-4).^2;
h13 = 2.135637957175861e+06*(1e-4).^2*(1e3);
h14 = 8.593828063620491e+03*(1e-4).^3;
h_poly2 = h00 + specvol_ice.*(h01 + specvol_ice.*(h02 + specvol_ice.*(h03 ...
+ specvol_ice.*(h04 + specvol_ice.*(h05 + h06.*specvol_ice))))) + p.*(h07 ...
+ specvol_ice.*(h08 + specvol_ice.*(h09 + specvol_ice.*(h10 + h11.*specvol_ice))) ...
+ p.*(h12 + h13.*specvol_ice + h14.*p));
plot(h_poly(:) - h_poly2(:))