« back to all changes in this revision
Viewing changes to gsw/gibbs/basic_thermodynamic_t.py
- browse files at revision 5
- compare with another revision
- download diff
- download tarball
- view history from revision 5
- Committer: Package Import Robot
- Author(s): Alastair McKinstry
- Date: 2013-07-18 07:26:49 UTC
- mfrom: (1.1.1)
- Revision ID: package-import@ubuntu.com-20130718072649-s2yct1q8we4z1w37
Tags: 3.0.2-1
* New upstream release.
* Move to DH 9.
* Move to DH 9.
- files added:
- .pc/python3.patch
- .pc/python3.patch/gsw
- .pc/python3.patch/gsw/gibbs
- .pc/python3.patch/gsw/test
- .pc/python3.patch/gsw/utilities
- files removed:
- .pc/remove-setuptools.patch
- gsw.egg-info
- files modified:
- .pc/applied-patches
added
removed
gsw/gibbs/basic_thermodynamic_t.py
4
4
5
5
import numpy as np
6
6
7
from library import gibbs
8
from gsw.utilities import match_args_return, strip_mask
9
from constants import Kelvin, db2Pascal, P0, SSO, cp0, R, sfac
10
from conversions import pt_from_CT, pt_from_t, pt0_from_t, molality_from_SA
11
from absolute_salinity_sstar_ct import CT_from_t
7
from .library import gibbs
8
from .absolute_salinity_sstar_ct import CT_from_t
9
from ..utilities import match_args_return, strip_mask
10
from .conversions import pt_from_CT, pt_from_t, pt0_from_t
11
from .constants import Kelvin, db2Pascal, P0, SSO, cp0, R, sfac, M_S
12
12
13
__all__ = [
13
14
'rho_t_exact',
14
15
'pot_rho_t_exact',
42
43
'osmotic_pressure_t_exact'
43
44
]
44
45
46
n0, n1, n2 = 0, 1, 2
47
45
48
46
49
@match_args_return
47
50
def Helmholtz_energy_t_exact(SA, t, p):
96
99
2011-03-29. Trevor McDougall
97
100
"""
98
101
99
n0, n1 = 0, 1
100
g000 = gibbs(n0, n0, n0, SA, t, p)
101
g001 = gibbs(n0, n0, n1, SA, t, p)
102
return (g000 - (db2Pascal * p + P0) * g001)
102
return (gibbs(n0, n0, n0, SA, t, p) -
103
(db2Pascal * p + P0) * gibbs(n0, n0, n1, SA, t, p))
103
104
104
105
105
106
@match_args_return
150
151
2011-03-29. Paul Barker, David Jackett and Trevor McDougal
151
152
"""
152
153
153
n0, n1 = 0, 1
154
rho = 1. / gibbs(n0, n0, n1, SA, t, p)
155
156
return rho
154
return 1. / gibbs(n0, n0, n1, SA, t, p)
157
155
158
156
159
157
@match_args_return
222
220
g08 = x2 * (-3310.49154044839 +
223
221
x * (199.459603073901 +
224
222
x * (-54.7919133532887 +
225
x * 36.0284195611086 -
226
y * 22.6683558512829) +
223
x * 36.0284195611086 -
224
y * 22.6683558512829) +
227
225
y * (-175.292041186547 +
228
226
y * (383.058066002476 +
229
227
y * (-460.319931801257 +
230
y * 234.565187611355)))) +
228
y * 234.565187611355)))) +
231
229
y * (729.116529735046 +
232
230
y * (-860.764303783977 +
233
231
y * (694.244814133268 +
239
237
240
238
return 100000000. / (g03 + g08) - 1000.0
241
239
240
242
241
@match_args_return
243
242
def enthalpy_t_exact(SA, t, p):
244
243
r"""Calculates the specific enthalpy of seawater.
292
291
2011-03-29. David Jackett, Trevor McDougall and Paul Barker.
293
292
"""
294
293
295
n0, n1 = 0, 1
296
g000 = gibbs(n0, n0, n0, SA, t, p)
297
g010 = gibbs(n0, n1, n0, SA, t, p)
298
return g000 - (t + Kelvin) * g010
294
return (gibbs(n0, n0, n0, SA, t, p) -
295
(t + Kelvin) * gibbs(n0, n1, n0, SA, t, p))
299
296
300
297
301
298
@match_args_return
345
342
2011-03-23. David Jackett and Paul Barker.
346
343
"""
347
344
348
n0, n1 = 0, 1
349
345
return gibbs(n0, n0, n1, SA, t, p)
350
346
351
347
405
401
2011-03-29. David Jackett, Trevor McDougall and Paul Barker.
406
402
"""
407
403
408
n0, n1 = 0, 1
409
404
return -gibbs(n0, n1, n0, SA, t, p)
410
405
411
406
454
449
2011-03-29. David Jackett, Trevor McDougall and Paul Barker
455
450
"""
456
451
457
n0, n2 = 0, 2
458
459
452
return -(t + Kelvin) * gibbs(n0, n2, n0, SA, t, p)
460
453
461
454
520
513
2011-03-29. David Jackett, Paul Barker and Trevor McDougall.
521
514
"""
522
515
523
n0, n1, n2 = 0, 1, 2
524
g020 = gibbs(n0, n2, n0, SA, t, p)
525
g011 = gibbs(n0, n1, n1, SA, t, p)
526
g001 = gibbs(n0, n0, n1, SA, t, p)
527
g002 = gibbs(n0, n0, n2, SA, t, p)
528
return g001 * np.sqrt(g020 / (g011 ** 2 - g020 * g002))
516
return (gibbs(n0, n0, n1, SA, t, p) * np.sqrt(gibbs(n0, n2, n0, SA, t, p) /
517
(gibbs(n0, n1, n1, SA, t, p) ** 2 - gibbs(n0, n2, n0, SA, t, p) *
518
gibbs(n0, n0, n2, SA, t, p))))
529
519
530
520
531
521
@match_args_return
579
569
2011-03-23. Trevor McDougall and Paul Barker
580
570
"""
581
571
582
n0, n1 = 0, 1
583
584
572
pt_zero = pt_from_CT(SSO, 0)
585
586
573
t_zero = pt_from_t(SSO, pt_zero, 0, p)
587
588
return (gibbs(n0, n0, n1, SA, t, p) - gibbs(n0, n0, n1, SSO, t_zero, p))
574
return (gibbs(n0, n0, n1, SA, t, p) -
575
gibbs(n0, n0, n1, SSO, t_zero, p))
589
576
590
577
591
578
@match_args_return
635
622
2011-03-29. Trevor McDougall and Paul Barker
636
623
"""
637
624
638
n0, n1 = 0, 1
639
625
return gibbs(n1, n0, n0, SA, t, p)
640
626
641
627
696
682
2011-03-29. Trevor McDougall
697
683
"""
698
684
699
n0, n1 = 0, 1
700
g000 = gibbs(n0, n0, n0, SA, t, p)
701
g010 = gibbs(n0, n1, n0, SA, t, p)
702
g001 = gibbs(n0, n0, n1, SA, t, p)
703
return g000 - (Kelvin + t) * g010 - (db2Pascal * p + P0) * g001
685
return (gibbs(n0, n0, n0, SA, t, p) -
686
(Kelvin + t) * gibbs(n0, n1, n0, SA, t, p) -
687
(db2Pascal * p + P0) * gibbs(n0, n0, n1, SA, t, p))
704
688
705
689
706
690
@match_args_return
757
741
2011-03-29. David Jackett, Trevor McDougall and Paul Barker
758
742
"""
759
743
760
n0, n1, n2 = 0, 1, 2
761
g002 = gibbs(n0, n0, n2, SA, t, p)
762
g001 = gibbs(n0, n0, n1, SA, t, p)
763
return -g002 / g001
744
return -gibbs(n0, n0, n2, SA, t, p) / gibbs(n0, n0, n1, SA, t, p)
764
745
765
746
766
747
@match_args_return
816
797
2011-03-29. David Jackett, Trevor McDougall and Paul Barker
817
798
"""
818
799
819
n0, n1 = 0, 1
820
821
g011 = gibbs(n0, n1, n1, SA, t, p)
822
g001 = gibbs(n0, n0, n1, SA, t, p)
823
return g011 / g001
800
return gibbs(n0, n1, n1, SA, t, p) / gibbs(n0, n0, n1, SA, t, p)
824
801
825
802
826
803
@match_args_return
870
847
2011-03-29. Trevor McDougall
871
848
"""
872
849
873
n0, n1, n2 = 0, 1, 2
874
875
g020 = gibbs(n0, n2, n0, SA, t, p)
876
g011 = gibbs(n0, n1, n1, SA, t, p)
877
g002 = gibbs(n0, n0, n2, SA, t, p)
878
879
return -(Kelvin + t) * (g020 - g011 * g011 / g002)
850
return (-(Kelvin + t) * (gibbs(n0, n2, n0, SA, t, p) -
851
gibbs(n0, n1, n1, SA, t, p) ** 2 / gibbs(n0, n0, n2, SA, t, p)))
880
852
881
853
882
854
@match_args_return
944
916
2011-03-23. David Jackett, Trevor McDougall and Paul Barker
945
917
"""
946
918
947
n0, n1, n2 = 0, 1, 2
948
949
g020 = gibbs(n0, n2, n0, SA, t, p)
950
g011 = gibbs(n0, n1, n1, SA, t, p)
951
g002 = gibbs(n0, n0, n2, SA, t, p)
952
g001 = gibbs(n0, n0, n1, SA, t, p)
953
954
return (g011 ** 2 - g020 * g002) / (g001 * g020)
919
return ((gibbs(n0, n1, n1, SA, t, p) ** 2 - gibbs(n0, n2, n0, SA, t, p) *
920
gibbs(n0, n0, n2, SA, t, p)) / (gibbs(n0, n0, n1, SA, t, p) *
921
gibbs(n0, n2, n0, SA, t, p)))
955
922
956
923
957
924
@match_args_return
1017
984
2011-03-28. Trevor McDougall and Paul Barker.
1018
985
"""
1019
986
1020
n0, n1 = 0, 1
1021
987
v_lab = np.ones_like(rho) / rho
1022
988
v_0 = gibbs(n0, n0, n1, 0, t, p)
1023
989
v_120 = gibbs(n0, n0, n1, 120, t, p)
1024
SA = 120 * (v_lab - v_0) / (v_120 - v_0) # Initial estimate of SA.
1025
Ior = (SA < 0) | (SA > 120)
1026
v_SA = (v_120 - v_0) / 120 # Initial estimate of v_SA, SA derivative of v
990
991
# Initial estimate of SA.
992
SA = 120 * (v_lab - v_0) / (v_120 - v_0)
993
Ior = np.logical_or(SA < 0, SA > 120)
994
995
# Initial estimate of v_SA, SA derivative of v
996
v_SA = (v_120 - v_0) / 120
1027
997
1028
998
for k in range(0, 2):
1029
999
SA_old = SA
1030
1000
delta_v = gibbs(n0, n0, n1, SA_old, t, p) - v_lab
1031
# Half way through the modified N-R method.
1001
# Half way the mod. N-R method (McDougall and Wotherspoon, 2012)
1032
1002
SA = SA_old - delta_v / v_SA
1033
1003
SA_mean = 0.5 * (SA + SA_old)
1034
1004
v_SA = gibbs(n1, n0, n1, SA_mean, t, p)
1133
1103
1134
1104
I_salty = alpha_freezing > alpha_limit
1135
1105
1136
t_diff = 40*ones(size(I_salty)) - t_freezing(I_salty);
1106
t_diff = 40. * np.ones_like(I_salty) - t_freezing(I_salty)
1137
1107
1138
1108
top = (rho_40[I_salty] - rho_freezing[I_salty] +
1139
1109
rho_freezing[I_salty] * alpha_freezing[I_salty] * t_diff)
1140
1110
1141
a = top / (t_diff ** 2)
1111
a = top / (t_diff ** 2)
1142
1112
b = -rho_freezing[I_salty] * alpha_freezing[I_salty]
1143
1113
c = rho_freezing[I_salty] - rho[I_salty]
1144
1114
sqrt_disc = np.sqrt(b ** 2 - 4 * a * c)
1201
1171
I_quad = ~np.isnan(t_a)
1202
1172
t[I_quad] = t_a[I_quad]
1203
1173
1204
I_quad = ~isnan(t_b)
1174
I_quad = ~np.isnan(t_b)
1205
1175
t_multiple[I_quad] = t_b[I_quad]
1206
1176
1207
1177
# After three iterations of this modified Newton-Raphson iteration,
1209
1179
1210
1180
return t, t_multiple
1211
1181
1182
1212
1183
@match_args_return
1213
1184
def pot_rho_t_exact(SA, t, p, p_ref=0):
1214
1185
r"""Calculates potential density of seawater.
1314
1285
2011-03-29. Trevor McDougall and Paul Barker
1315
1286
"""
1316
1287
1317
n0, n1, n2 = 0, 1, 2
1318
1319
1288
pt0 = pt0_from_t(SA, t, p)
1320
1321
1289
factor = -cp0 / ((Kelvin + pt0) * gibbs(n0, n2, n0, SA, t, p))
1322
1323
g011 = gibbs(n0, n1, n1, SA, t, p)
1324
g001 = gibbs(n0, n0, n1, SA, t, p)
1325
1326
return factor * (g011 / g001)
1290
return factor * (gibbs(n0, n1, n1, SA, t, p) / gibbs(n0, n0, n1, SA, t, p))
1327
1291
1328
1292
1329
1293
@match_args_return
1374
1338
2011-03-29. David Jackett, Trevor McDougall and Paul Barker
1375
1339
"""
1376
1340
1377
n0, n1, n2 = 0, 1, 2
1378
1379
1341
pt0 = pt0_from_t(SA, t, p)
1380
1381
1342
factor = gibbs(n0, n2, n0, SA, pt0, 0) / gibbs(n0, n2, n0, SA, t, p)
1382
g011 = gibbs(n0, n1, n1, SA, t, p)
1383
g001 = gibbs(n0, n0, n1, SA, t, p)
1384
1385
return factor * (g011 / g001)
1343
return factor * (gibbs(n0, n1, n1, SA, t, p) / gibbs(n0, n0, n1, SA, t, p))
1386
1344
1387
1345
1388
1346
@match_args_return
1435
1393
1436
1394
# TODO: Original GSW-V3 re-implements gibbs, check what to do here!
1437
1395
1438
n0, n1, n2 = 0, 1, 2
1439
1440
1396
pt0 = pt0_from_t(SA, t, p)
1441
1397
1442
g001 = gibbs(n0, n0, n1, SA, t, p)
1443
g110 = gibbs(n1, n1, n0, SA, t, p)
1444
g100 = gibbs(n1, n0, n0, SA, pt0, 0)
1445
1446
factora = g110 - g100 / (Kelvin + pt0)
1447
g020 = gibbs(n0, n2, n0, SA, t, p)
1448
factor = factora / (g001 * g020)
1449
1450
g011 = gibbs(n0, n1, n1, SA, t, p)
1451
g101 = gibbs(n1, n0, n1, SA, t, p)
1452
1453
return g011 * factor - g101 / g001
1398
factora = (gibbs(n1, n1, n0, SA, t, p) - gibbs(n1, n0, n0, SA, pt0, 0) /
1399
(Kelvin + pt0))
1400
factor = (factora / (gibbs(n0, n0, n1, SA, t, p) *
1401
gibbs(n0, n2, n0, SA, t, p)))
1402
1403
return (gibbs(n0, n1, n1, SA, t, p) * factor -
1404
gibbs(n1, n0, n1, SA, t, p) / gibbs(n0, n0, n1, SA, t, p))
1454
1405
1455
1406
1456
1407
@match_args_return
1500
1451
Modifications:
1501
1452
2011-04-10. Trevor McDougall and Paul Barker
1502
1453
"""
1503
# TODO: Original GSW-V3 re-implements gibbs, check what to do here!
1504
1505
n0, n1, n2 = 0, 1, 2
1454
# NOTE: The original Matlab toolbox re-implement some code here. Why?
1506
1455
1507
1456
pt0 = pt0_from_t(SA, t, p)
1508
1457
1509
g001 = gibbs(n0, n0, n1, SA, t, p)
1510
g110 = gibbs(n1, n1, n0, SA, t, p)
1511
factora = g110 - gibbs(n1, n1, n0, SA, pt0, 0)
1512
1513
g020 = gibbs(n0, n2, n0, SA, t, p)
1514
factor = factora / (g001 * g020)
1515
1516
g011 = gibbs(n0, n1, n1, SA, t, p)
1517
g101 = gibbs(n1, n0, n1, SA, t, p)
1518
1519
return g011 * factor - g101 / g001
1458
factora = gibbs(n1, n1, n0, SA, t, p) - gibbs(n1, n1, n0, SA, pt0, 0)
1459
1460
factor = (factora / (gibbs(n0, n0, n1, SA, t, p) *
1461
gibbs(n0, n2, n0, SA, t, p)))
1462
1463
return (gibbs(n0, n1, n1, SA, t, p) * factor -
1464
gibbs(n1, n0, n1, SA, t, p) / gibbs(n0, n0, n1, SA, t, p))
1520
1465
1521
1466
1522
1467
@match_args_return
1567
1512
2011-03-29. David Jackett, Trevor McDougall and Paul Barker
1568
1513
"""
1569
1514
1570
n0, n1 = 0, 1
1571
1572
g101 = gibbs(n1, n0, n1, SA, t, p)
1573
g001 = gibbs(n0, n0, n1, SA, t, p)
1574
return -g101 / g001
1515
return -gibbs(n1, n0, n1, SA, t, p) / gibbs(n0, n0, n1, SA, t, p)
1575
1516
1576
1517
1577
1518
@match_args_return
1798
1739
2011-03-29. Trevor McDougall and Paul Barker
1799
1740
"""
1800
1741
1801
n0, n1, n2 = 0, 1, 2
1802
g011 = gibbs(n0, n1, n1, SA, t, p)
1803
g020 = gibbs(n0, n2, n0, SA, t, p)
1804
return - g011 / g020
1742
return -gibbs(n0, n1, n1, SA, t, p) / gibbs(n0, n2, n0, SA, t, p)
1805
1743
1806
1744
1807
1745
@match_args_return
1848
1786
UNESCO (English), 196 pp.
1849
1787
1850
1788
Modifications:
1851
2011-04-01. Trevor McDougall and Paul Barker
1789
2011-04-01. Trevor McDougall and Paul Barker.
1790
2012-11-15. Trevor McDougall and Paul Barker.
1852
1791
"""
1853
1792
1854
n0 = 0
1855
1856
# Molality of seawater in mol/kg.
1857
molal = molality_from_SA(SA)
1858
part = molal * R * (Kelvin + t)
1859
1860
SAzero = 0
1861
g000 = gibbs(n0, n0, n0, SAzero, t, p)
1862
return (g000 - chem_potential_water_t_exact(SA, t, p)) / part
1793
SA = np.maximum(SA, 0)
1794
k = M_S / R
1795
part = k * (1000 - SA) / (Kelvin + t)
1796
1797
x2 = sfac * SA
1798
x = np.sqrt(x2)
1799
y = t * 0.025
1800
# Note that the input pressure (p) is sea pressure in units of dbar.
1801
z = p / db2Pascal
1802
1803
oc = (7.231916621570606e1, 1.059039593127674e1, -3.025914794694813e1,
1804
5.040733670521486e1, -4.074543321119333e1, 1.864215613820487e1,
1805
-3.022566485046178, -6.138647522851840, 1.353207379758663e1,
1806
-7.316560781114737, 1.829232499785750, -5.358042980767074e-1,
1807
-1.705887283375562, -1.246962174707332e-1, 1.228376913546017,
1808
1.089364009088042e-2, -4.264828939262248e-1, 6.213127679460041e-2,
1809
2.481543497315280, -1.363368964861909, -5.640491627443773e-1,
1810
1.344724779893754, -2.180866793244492, 4.765753255963401,
1811
-5.726993916772165, 2.918303792060746, -6.506082399183509e-1,
1812
-1.015695507663942e-1, 1.035024326471108, -6.742173543702397e-1,
1813
8.465642650849419e-1, -7.508472135244717e-1, -3.668086444057845e-1,
1814
3.189939162107803e-1, -4.245629194309487e-2)
1815
1816
tl = (oc[0] + oc[1] * y + x * (oc[2] + x * (oc[3] + x * (oc[4] + x *
1817
(oc[5] + oc[6] * x))) + y * (oc[7] + x * (oc[8] + x *
1818
(oc[9] + oc[10] * x)) + y * (oc[11] + oc[12] * x + y * (oc[13] +
1819
oc[14] * x + y * (oc[15] + x * (oc[16] + oc[17] * y))))) + z *
1820
(oc[18] + x * (oc[19] + oc[20] * y + oc[21] * x) + y * (oc[22] + y *
1821
(oc[23] + y * (oc[24] + oc[25] * y))) + z * (oc[26] + oc[27] * x + y *
1822
(oc[28] + oc[29] * y) + z * (oc[30] + oc[31] * x + y * (oc[32] +
1823
oc[33] * y) + oc[34] * z)))))
1824
1825
return tl * part
1863
1826
1864
1827
1865
1828
@match_args_return
1959
1922
2011-04-03. Trevor McDougall and Paul Barker
1960
1923
"""
1961
1924
1962
n0, n1 = 0, 1
1963
1925
# The temperature increment for calculating the gibbs_PTT derivative.
1964
1926
dt = 0.001
1965
1927
t = 3.978 - 0.22072 * SA # The initial guess of t_maxden.
1968
1930
for Number_of_iterations in range(0, 3):
1969
1931
t_old = t
1970
1932
gibbs_PT = gibbs(n0, n1, n1, SA, t_old, p)
1971
# This is half way through the modified method
1933
# Half way through the mod. method (McDougall and Wotherspoon, 2012)
1972
1934
t = t_old - gibbs_PT / gibbs_PTT
1973
1935
t_mean = 0.5 * (t + t_old)
1974
1936
gibbs_PTT = (gibbs(n0, n1, n1, SA, t_mean + dt, p) -
2021
1983
Modifications:
2022
1984
2011-05-26. Trevor McDougall and Paul Barker
2023
1985
"""
2024
1986
SA = np.maximum(SA, 0)
2025
1987
gibbs_pure_water = gibbs(0, 0, 0, 0, t, pw)
2026
1988
2027
1989
# Initial guess of p, in dbar.
2028
1990
p = pw + 235.4684
2029
1991
2030
1992
# Initial guess of df/dp.
2031
df_dp = -db2Pascal * (gibbs(0, 0, 1, SA, t, p) -
2032
SA * gibbs(1, 0, 1, SA, t, p))
1993
df_dp = -db2Pascal * (gibbs(n0, n0, n1, SA, t, p) -
1994
SA * gibbs(n1, n0, n1, SA, t, p))
2033
1995
2034
1996
for Number_of_iterations in range(0, 2):
2035
1997
p_old = p
Loggerhead is a web-based interface for Breezy
Version: 2.0.1