Canoe
Comprehensive Atmosphere N' Ocean Engine
attenuation_h2o.cpp
Go to the documentation of this file.
1 /**************************************************************************
2 
3  ===== JUNO codes =====
4 
5  Copyright (C) 2019 California Institute of Technology (Caltech)
6  All rights reserved.
7 
8  Written by Fabiano Oyafuso (fabiano@jpl.nasa.gov)
9  Adapted by Cheng Li (cli@gps.caltech.edu) to SNAP structure
10 
11 **************************************************************************/
12 
13 // C/C++ headers
14 #include <cmath>
15 
16 // cli, 190801
17 inline double SQUARE(double x) { return x * x; }
18 inline double CUBE(double x) { return x * x * x; }
19 inline double atm_from_bar(double P) { return 0.9869232667160128 * P; }
20 inline double bar_from_atm(double P) { return 1.013250 * P; }
21 static double const pi = 3.14159265358979;
22 static double const kBoltz_mks = 1.3806504e-23; // J/K
23 static double invcm_from_dB_per_km = 2.3059e-6; // dB/km to cm^-1
24 // end
25 
26 double attenuation_H2O_Karpowicz(double freq, double P_idl, double T,
27  double XH2, double XHe, double XH2O,
28  double scale) {
29  double const C_H2Ob = 4.36510480961e-7; // (Mbar*GHz)^-2 (km)^-1
30  double const C_H2Oa =
31  3.1e-7; // (Mbar*GHz)^-2 (km)^-1 (latest update: Nov 2012)
32  double C_H2O = C_H2Oa * (1. - scale) + C_H2Ob * scale;
33 
34  double const Cprime_H2Ob = 2.10003048186e-26;
35  double const Cprime_H2Oa = 0; // (latest update: Nov 2012)
36  double Cprime_H2O = Cprime_H2Oa * (1. - scale) + Cprime_H2Ob * scale;
37 
38  double const C_H2 = 5.07722009423e-11; // (Mbar*GHz)^-2 (km)^-1
39  double const C_He = 1.03562010226e-10; // (Mbar*GHz)^-2 (km)^-1
40  double const x_ctmb = 13.3619799812;
41  double const x_ctma = 12; // (latest update: Nov 2012)
42  double x_ctm = x_ctma * (1. - scale) + x_ctmb * scale;
43 
44  double const n_ctm = 6.76418487001;
45  double const xprime_ctm = 0.0435525417274;
46 
47  // GHz
48  static double const line[] = {22.2351, 183.3101, 321.2256, 325.1529,
49  380.1974, 439.1508, 443.0183, 448.0011,
50  470.8890, 474.6891, 488.4911, 556.9360,
51  620.7008, 752.0332, 916.1712};
52 
53  // Np-Hz/cm^2 (corigendum says this should be cm^2-Hz/molecule)
54  static double const intensity[] = {
55  0.1314e-13, 0.2279e-11, 0.8058e-13, 0.2701e-11, 0.2444e-10,
56  0.2185e-11, 0.4637e-12, 0.2568e-10, 0.8392e-12, 0.3272e-11,
57  0.6676e-12, 0.1535e-8, 0.1711e-10, 0.1014e-8, 0.4238e-10};
58 
59  // GHz/bar
60  static double const linewidth_H2O[] = {13.49, 14.66, 10.57, 13.81, 14.54,
61  9.715, 7.88, 12.75, 9.83, 10.95,
62  13.13, 14.05, 11.836, 12.53, 12.75};
63 
64  // GHz/bar
65  static double const linewidth_H2[] = {2.395, 2.400, 2.395, 2.395, 2.390,
66  2.395, 2.395, 2.395, 2.395, 2.395,
67  2.395, 2.395, 2.395, 2.395, 2.395};
68 
69  // GHz/bar
70  static double const linewidth_He[] = {0.67, 0.71, 0.67, 0.67, 0.63,
71  0.67, 0.67, 0.67, 0.67, 0.67,
72  0.67, 0.67, 0.67, 0.67, 0.67};
73 
74  static double const exp_H2O[] = {0.61, 0.85, 0.54, 0.74, 0.89,
75  0.62, 0.50, 0.67, 0.65, 0.64,
76  0.72, 1.0, 0.68, 0.84, 0.78};
77 
78  static double const temp_coeff[] = {2.144, 0.668, 6.179, 1.541, 1.048,
79  3.595, 5.048, 1.405, 3.597, 2.379,
80  2.852, 0.159, 2.391, 0.396, 1.441};
81 
82  static double const exp_H2[] = {0.900, 0.950, 0.900, 0.900, 0.850,
83  0.900, 0.900, 0.900, 0.900, 0.900,
84  0.900, 0.900, 0.900, 0.900, 0.900};
85 
86  static double const exp_He[] = {0.515, 0.490, 0.515, 0.490, 0.540,
87  0.515, 0.515, 0.515, 0.515, 0.515,
88  0.515, 0.515, 0.515, 0.515, 0.515};
89 
90  double const PH2O_idl = P_idl * XH2O;
91  double const PH2_idl = P_idl * XH2;
92  double const PHe_idl = P_idl * XHe;
93 
94  double const TR = 300 / T;
95  double const nH2O = 0.1 * PH2O_idl / (kBoltz_mks * T); // molecules/cm^3
96 
97  // for now, assume only one isotope
98  double abs_lines = 0.0;
99  for (auto i = 0; i < sizeof(line) / sizeof(double); ++i) {
100  // should have units of GHz
101  double gamma = PH2O_idl * pow(TR, exp_H2O[i]) * linewidth_H2O[i] +
102  PH2_idl * pow(TR, exp_H2[i]) * linewidth_H2[i] +
103  PHe_idl * pow(TR, exp_He[i]) * linewidth_He[i];
104 
105  // 1/GHz
106  double sum_nu_line = fabs(freq + line[i]);
107  double diff_nu_line = fabs(freq - line[i]);
108  double cutoff_term = 1.0 / (SQUARE(750.0) + SQUARE(gamma));
109 
110  double plus_term =
111  (sum_nu_line < 750
112  ? 1.0 / (SQUARE(sum_nu_line) + SQUARE(gamma)) - cutoff_term
113  : 0.0);
114 
115  double minus_term =
116  (diff_nu_line < 750
117  ? 1.0 / (SQUARE(diff_nu_line) + SQUARE(gamma)) - cutoff_term
118  : 0.0);
119 
120  // no cutoff
121  // plus_term = 1.0/(SQUARE(sum_nu_line) + SQUARE(gamma));
122  // minus_term = 1.0/(SQUARE(diff_nu_line) + SQUARE(gamma));
123 
124  double lineshape =
125  SQUARE(freq / line[i]) * gamma / pi * (plus_term + minus_term);
126 
127  abs_lines += intensity[i] * exp((1 - TR) * temp_coeff[i]) * lineshape;
128  }
129  abs_lines *= pow(TR, 2.5) * nH2O;
130  abs_lines *= 1e-9; // Hz/GHz --> unitless
131 
132  double mbar_of_bar = 1.0e3;
133  double abs_ctm = 0.0;
134  abs_ctm += C_H2O * SQUARE(mbar_of_bar * PH2O_idl) * pow(TR, x_ctm);
135  abs_ctm +=
136  Cprime_H2O * pow(mbar_of_bar * PH2O_idl, n_ctm) * pow(TR, xprime_ctm);
137  abs_ctm += (C_H2 * PH2_idl + C_He * PHe_idl) * PH2O_idl * CUBE(TR) *
138  SQUARE(mbar_of_bar);
139  abs_ctm *= SQUARE(freq);
140  abs_ctm *= 1e-5; // inv(km) -> inv(cm)
141 
142  double abs_tot = abs_ctm + abs_lines; // 1/cm
143 
144 #ifdef ABSORPTION_KLUDGE
145  constexpr double P_top = 150;
146  constexpr double P_bot = 4000;
147  constexpr double scale_bot = 1.0;
148  if (P_idl > P_bot)
149  abs_tot *= scale_bot;
150  else if (P_idl > P_top)
151  abs_tot *=
152  ((P_idl - P_top) * scale_bot + (P_bot - P_idl)) / (P_bot - P_top);
153 #endif
154 
155  return abs_tot;
156 }
157 
158 double attenuation_H2O_Goodman(double freq, double P, double T, double XH2,
159  double XHe, double XH2O) {
160  double const PH2O = atm_from_bar(P * XH2O);
161  double const PH2 = atm_from_bar(P * XH2);
162  double const PHe = atm_from_bar(P * XHe);
163 
164  double TR = 273.0 / T;
165  double F2 = freq * freq;
166 
167  double DNU1 = 9.88e-2 * pow(TR, 0.6667) * (0.81 * PH2 + 0.35 * PHe);
168  double DNU2 = DNU1 * DNU1;
169  double DNUP = freq / 29.97 + 0.74;
170  double DNUM = freq / 29.97 - 0.74;
171  double DTAUA =
172  PH2O * F2 * pow(TR, 4.3333) *
173  (9.07e-9 * (DNU1 / (DNUM * DNUM + DNU2) + DNU1 / (DNUP * DNUP + DNU2)) +
174  1.45e-7 * DNU1);
175 
176  return DTAUA;
177 }
178 
179 double attenuation_H2O_Waters(double freq, double P, double T, double XH2,
180  double XHe, double XH2O) {
181  double const P_atm = atm_from_bar(P);
182  double const PH2O = atm_from_bar(P) * XH2O;
183 
184  double F2 = freq * freq;
185 
186  double TR = 273.0 / T;
187  double DTAUA = 0.5596 * F2 * P_atm * PH2O * pow(TR, 3.126);
188  // higher-order terms:
189  double X1 = 1.0 + 3.897 * PH2O / P_atm;
190  double DNU1 = 2.96 * X1 * pow((300.0 / T), 0.626) * P_atm;
191  double X2 = pow((494.4019 - F2), 2) + 4 * F2 * DNU1 * DNU1;
192 
193  DTAUA *= X1 * (2.77e-8 + 7.18 * exp(-644.0 / T) / (T * X2));
194 
195  return DTAUA;
196 }
197 
198 double attenuation_H2O_deBoer(double freq, double P, double T, double XH2,
199  double XHe, double XH2O) {
200  static double const fi[10] = {22.23515, 183.31012, 323.0, 325.1538, 380.1968,
201  390.0, 436.0, 438.0, 442.0, 448.0};
202 
203  static double const ei[10] = {644.0, 196.0, 1850.0, 454.0, 306.0,
204  2199.0, 1507.0, 1070.0, 1507.0, 412.0};
205 
206  static double const ai[10] = {1.0, 41.9, 334.4, 115.7, 651.8,
207  127.0, 191.4, 697.6, 590.2, 973.1};
208 
209  static double const c1[10][3] = {{2.395, 0.67, 10.67}, {2.40, 0.71, 11.64},
210  {2.395, 0.67, 9.59}, {2.395, 0.67, 11.99},
211  {2.39, 0.63, 12.42}, {2.395, 0.67, 9.16},
212  {2.395, 0.67, 6.32}, {2.395, 0.67, 8.34},
213  {2.395, 0.67, 6.52}, {2.395, 0.67, 11.57}};
214 
215  static double const c2[10][3] = {{0.9, 0.515, 0.626}, {0.95, 0.49, 0.649},
216  {0.9, 0.515, 0.42}, {0.9, 0.515, 0.619},
217  {0.85, 0.54, 0.63}, {0.9, 0.515, 0.33},
218  {0.9, 0.515, 0.29}, {0.9, 0.515, 0.36},
219  {0.9, 0.515, 0.332}, {0.9, 0.515, 0.51}};
220 
221  double const PH2O = atm_from_bar(P * XH2O);
222  double const PH2 = atm_from_bar(P * XH2);
223  double const PHe = atm_from_bar(P * XHe);
224 
225  double TR = 300.0 / T;
226  double F2 = freq * freq;
227 
228  // convert P(atm) to P(bars) and store in convenient array:
229  double px[3];
230  px[0] = bar_from_atm(PH2);
231  px[1] = bar_from_atm(PHe);
232  px[2] = bar_from_atm(PH2O);
233 
234  double sum = 0.0;
235  for (int i = 0; i < 10; i++) {
236  double gami = 0.0;
237  for (int j = 0; j < 3; j++) {
238  gami += c1[i][j] * px[j] * pow(TR, c2[i][j]);
239  }
240  sum += ai[i] * exp(-ei[i] / T) * gami /
241  (pow((fi[i] * fi[i] - F2), 2) + 4 * F2 * gami * gami);
242  }
243 
244  double DTAUA = 1444.5 * px[2] * pow(TR, 3.5) * F2 * sum;
245  DTAUA += 0.00339 * px[2] * pow(TR, 3.1) * (0.81 * px[0] + 0.35 * px[1]) * F2;
246  // convert dB/km to 1/cm:
247  DTAUA *= invcm_from_dB_per_km;
248 
249  return DTAUA;
250 }
atm_from_bar
double atm_from_bar(double P)
Definition: attenuation_h2o.cpp:19
pi
static double const pi
Definition: attenuation_h2o.cpp:21
CUBE
double CUBE(double x)
Definition: attenuation_h2o.cpp:18
attenuation_H2O_Karpowicz
double attenuation_H2O_Karpowicz(double freq, double P_idl, double T, double XH2, double XHe, double XH2O, double scale)
Definition: attenuation_h2o.cpp:26
attenuation_H2O_deBoer
double attenuation_H2O_deBoer(double freq, double P, double T, double XH2, double XHe, double XH2O)
Definition: attenuation_h2o.cpp:198
bar_from_atm
double bar_from_atm(double P)
Definition: attenuation_h2o.cpp:20
invcm_from_dB_per_km
static double invcm_from_dB_per_km
Definition: attenuation_h2o.cpp:23
attenuation_H2O_Goodman
double attenuation_H2O_Goodman(double freq, double P, double T, double XH2, double XHe, double XH2O)
Definition: attenuation_h2o.cpp:158
kBoltz_mks
static double const kBoltz_mks
Definition: attenuation_h2o.cpp:22
SQUARE
double SQUARE(double x)
Definition: attenuation_h2o.cpp:17
attenuation_H2O_Waters
double attenuation_H2O_Waters(double freq, double P, double T, double XH2, double XHe, double XH2O)
Definition: attenuation_h2o.cpp:179
fit_ammonia.e
e
Definition: fit_ammonia.py:8
make_plots.j
j
Definition: make_plots.py:26