source: issm/trunk-jpl/src/m/materials/cuffey.py@ 26358

Last change on this file since 26358 was 26358, checked in by jdquinn, 4 years ago

CHG: Completed MATLAB -> Python updates for SE; archive updates now that GMSH can be used on macOS and Linux; various minor bug fixes; formatting; cleanup
File size: 3.5 KB
Line 
1import numpy as np
2
3
4def cuffey(temperature):
5 """CUFFEY - calculates ice rigidity as a function of temperature
6
7 rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law
8 sigma = B * e(1/3) (Cuffey and Paterson, p75).
9
10 temperature is in Kelvin degrees
11
12 Usage:
13 rigidity = cuffey(temperature)
14 """
15
16 if np.any(temperature < 0.0):
17 raise RuntimeError('input temperature should be in Kelvin (positive)')
18
19 if np.ndim(temperature) == 2:
20 #T = temperature.reshape(-1, ) - 273.15
21 T = temperature.flatten() - 273.15
22 elif isinstance(temperature, float) or isinstance(temperature, int):
23 T = np.array([temperature]) - 273.15
24 else:
25 T = temperature - 273.15
26
27 rigidity = np.zeros_like(T)
28 pos = np.nonzero(T <= -45)
29 if len(pos):
30 rigidity[pos] = pow(10, 8) * (-0.000396645116301 * (T[pos] + 50)**3 + 0.013345579471334 * (T[pos] + 50)**2 - 0.356868703259105 * (T[pos] + 50) + 7.272363035371383)
31 pos = np.nonzero(np.logical_and(-45 <= T, T < -40))
32 if len(pos):
33 rigidity[pos] = pow(10, 8) * (-0.000396645116301 * (T[pos] + 45)**3 + 0.007395902726819 * (T[pos] + 45)**2 - 0.253161292268336 * (T[pos] + 45) + 5.772078366321591)
34 pos = np.nonzero(np.logical_and(-40 <= T, T < -35))
35 if len(pos):
36 rigidity[pos] = pow(10, 8) * (0.000408322072669 * (T[pos] + 40)**3 + 0.001446225982305 * (T[pos] + 40)**2 - 0.208950648722716 * (T[pos] + 40) + 4.641588833612773)
37 pos = np.nonzero(np.logical_and(-35 <= T, T < -30))
38 if len(pos):
39 rigidity[pos] = pow(10, 8) * (-0.000423888728124 * (T[pos] + 35)**3 + 0.007571057072334 * (T[pos] + 35)**2 - 0.163864233449525 * (T[pos] + 35) + 3.684031498640382)
40 pos = np.nonzero(np.logical_and(-30 <= T, T < -25))
41 if len(pos):
42 rigidity[pos] = pow(10, 8) * (0.000147154327025 * (T[pos] + 30)**3 + 0.001212726150476 * (T[pos] + 30)**2 - 0.119945317335478 * (T[pos] + 30) + 3.001000667185614)
43 pos = np.nonzero(np.logical_and(-25 <= T, T < -20))
44 if len(pos):
45 rigidity[pos] = pow(10, 8) * (-0.000193435838672 * (T[pos] + 25)**3 + 0.003420041055847 * (T[pos] + 25)**2 - 0.096781481303861 * (T[pos] + 25) + 2.449986525148220)
46 pos = np.nonzero(np.logical_and(-20 <= T, T < -15))
47 if len(pos):
48 rigidity[pos] = pow(10, 8) * (0.000219771255067 * (T[pos] + 20)**3 + 0.000518503475772 * (T[pos] + 20)**2 - 0.077088758645767 * (T[pos] + 20) + 2.027400665191131)
49 pos = np.nonzero(np.logical_and(-15 <= T, T < -10))
50 if len(pos):
51 rigidity[pos] = pow(10, 8) * (-0.000653438900191 * (T[pos] + 15)**3 + 0.003815072301777 * (T[pos] + 15)**2 - 0.055420879758021 * (T[pos] + 15) + 1.682390865739973)
52 pos = np.nonzero(np.logical_and(-10 <= T, T < -5))
53 if len(pos):
54 rigidity[pos] = pow(10, 8) * (0.000692439419762 * (T[pos] + 10)**3 - 0.005986511201093 * (T[pos] + 10)**2 - 0.066278074254598 * (T[pos] + 10) + 1.418983411970382)
55 pos = np.nonzero(np.logical_and(-5 <= T, T < -2))
56 if len(pos):
57 rigidity[pos] = pow(10, 8) * (-0.000132282004110 * (T[pos] + 5)**3 + 0.004400080095332 * (T[pos] + 5)**2 - 0.074210229783403 * (T[pos] + 5) + 1.024485188140279)
58 pos = np.nonzero(-2 <= T)
59 if len(pos):
60 rigidity[pos] = pow(10, 8) * (-0.000132282004110 * (T[pos] + 2)**3 + 0.003209542058346 * (T[pos] + 2)**2 - 0.051381363322371 * (T[pos] + 2) + 0.837883605537096)
61
62 # Now make sure that rigidity is positive
63 pos = np.nonzero(rigidity < 0)
64 rigidity[pos] = pow(1, 6)
65
66 return rigidity
Note: See TracBrowser for help on using the repository browser.