6 static std::complex<double> t1(0.368017112855051e-1, 0.510878114959572e-1);
7 static std::complex<double> r1(0.447050716285388e2, 0.656876847463481e2);
8 static std::complex<double> t2(0.337315741065416, 0.335449415919309);
9 static std::complex<double> r20(-0.725974574329220e2, -0.781008427112870e2);
10 static std::complex<double> r21(-0.557107698030123e-4, 0.464578634580806e-4);
11 static std::complex<double> r22(0.234801409215913e-10, -0.285651142904972e-10);
16 static double T_t = 273.16,
21 static double g00 = -0.632020233449497e6;
22 static double g01 = 0.655022213658955;
23 static double g02 = -0.189369929326131e-7;
24 static double g03 = 0.339746123271053e-14;
25 static double g04 = -0.556464869058991e-21;
26 static double s0 = -0.332733756492168e4;
38 double a[] = {0, -0.212144006e2, 0.273203819e2, -0.610598130e1};
39 double b[] = {0, 0.333333333e-2, 0.120666667e1, 0.170333333e1};
40 double summer = 0, theta;
42 for (
int i = 1; i <= 3; i++) {
43 summer += a[i] * pow(theta, b[i]);
45 return p_t * exp(1 / theta * summer);
51 double g_Ice(
double T,
double p) {
53 std::complex<double> r2, term1, term2;
54 double g0, theta, pi, pi_0;
58 g0 = g00 * pow(pi - pi_0, 0.0) + g01 * pow(pi - pi_0, 1.0) + g02 * pow(pi - pi_0, 2.0) + g03 * pow(pi - pi_0, 3.0) + g04 * pow(pi - pi_0, 4.0);
59 r2 = r20 * pow(pi - pi_0, 0.0) + r21 * pow(pi - pi_0, 1.0) + r22 * pow(pi - pi_0, 2.0);
61 term1 = r1 * ((t1 - theta) * log(t1 - theta) + (t1 + theta) * log(t1 + theta) - 2.0 * t1 * log(t1) - theta * theta / t1);
62 term2 = r2 * ((t2 - theta) * log(t2 - theta) + (t2 + theta) * log(t2 + theta) - 2.0 * t2 * log(t2) - theta * theta / t2);
63 return g0 - s0 * T_t * theta + T_t * real(term1 + term2);
71 std::complex<double> r2_p;
72 double g0_p, theta, pi, pi_0;
76 g0_p = g01 * 1.0 / p_t * pow(pi - pi_0, 1 - 1.0) + g02 * 2.0 / p_t * pow(pi - pi_0, 2 - 1.0) + g03 * 3.0 / p_t * pow(pi - pi_0, 3 - 1.0)
77 + g04 * 4.0 / p_t * pow(pi - pi_0, 4 - 1.0);
78 r2_p = r21 * 1.0 / p_t * pow(pi - pi_0, 1 - 1.0) + r22 * 2.0 / p_t * pow(pi - pi_0, 2 - 1.0);
79 return g0_p + T_t * real(r2_p * ((t2 - theta) * log(t2 - theta) + (t2 + theta) * log(t2 + theta) - 2.0 * t2 * log(t2) - theta * theta / t2));
87 std::complex<double> r2_pp;
88 double g0_pp, theta, pi, pi_0;
92 g0_pp = g02 * 2.0 * (2.0 - 1.0) / p_t / p_t * pow(pi - pi_0, 2.0 - 2.0) + g03 * 3.0 * (3.0 - 1.0) / p_t / p_t * pow(pi - pi_0, 3.0 - 2.0)
93 + g04 * 4.0 * (4.0 - 1.0) / p_t / p_t * pow(pi - pi_0, 4 - 2.0);
94 r2_pp = r22 * 2.0 / p_t / p_t;
95 return g0_pp + T_t * real(r2_pp * ((t2 - theta) * log(t2 - theta) + (t2 + theta) * log(t2 + theta) - 2.0 * t2 * log(t2) - theta * theta / t2));
103 std::complex<double> r2, term1, term2;
104 double theta, pi, pi_0;
108 r2 = r20 * pow(pi - pi_0, 0.0) + r21 * pow(pi - pi_0, 1.0) + r22 * pow(pi - pi_0, 2.0);
110 term1 = r1 * (-log(t1 - theta) + log(t1 + theta) - 2.0 * theta / t1);
111 term2 = r2 * (-log(t2 - theta) + log(t2 + theta) - 2.0 * theta / t2);
112 return -s0 + real(term1 + term2);