Canoe
Comprehensive Atmosphere N' Ocean Engine
attenuation_nh3.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 // nh3 data headers
17 #include "mwr_absorber_nh3.hpp"
18 
19 // cli, 190801
20 inline double SQUARE(double x) { return x * x; }
21 inline double atm_from_bar(double P) { return 0.9869232667160128 * P; }
22 static double const pi = 3.14159265358979;
23 // end
24 
25 double attenuation_NH3_Bellotti_switch(double freq, double P, double P_idl,
26  double T, double XH2, double XHe,
27  double XNH3, double XH2O)
28 // sum the contributions of all lines, first the Inversion transitions,
29 // then the Rotational transitions
30 
31 {
32 #if defined(FAKE_CUTOFF_TO_AVOID_SINGULARITY_IN_BELLOTTI_ABS)
33  // Bellotti's parameterization suffers from two problems at high pressure:
34  // (1) for P~7000 bar, absorption becomes negative!!!
35  // (2) slope of opacity wrt pressure changes sign around P~1000 bar.
36  // We patch this by readjusting opacity for T>T_cutoff. (We choose a
37  // temperature cutoff instead of a pressure cutoff, since the opacity is more
38  // strongly dependent on temperature.) Update: this anomalous behavior has
39  // been corrected by zeroing zhe
40  constexpr double T_Bellotti_cutoff = 1250;
41  if (T > T_Bellotti_cutoff) {
42  T = T_Bellotti_cutoff;
43  }
44 #endif
45 
46  // TEMPORARY KLUDGE to compare to Hanley, fig19!!!
47  // T=295;
48  // XNH3=0.004366812227074;
49  bool use_mm_parameters(freq >= 30 ? true : false);
50  // miscellaneous constants:
51  double const ppscale = 2.99792458e18; // converts to cm^2*cm-1 from nm^2*MHz
52  double const gcon = 25.923;
53  double const cmghz = 29.9792458; // converts to GHz from cm-1
54  double const hckt0 = 0.0047959232; // hc/kT0, where T0=300K
55  double const cdens =
56  7.242971565e21; // 1/kB * (1e5 Pa/bar) * (1e-6 m^-3/cm^-3)
57 
58  // partial pressures in bars
59  double const PNH3 = P_idl * XNH3;
60  double const PH2 = P_idl * XH2;
61  double const PHe = P_idl * XHe;
62  double const PH2O = P_idl * XH2O;
63 
64  double const dens = cdens * (P_idl * XNH3) / T;
65  // double const dens = cdens*PNH3/T;
66  double const beta300 = 300.0 / T;
67  double const beta295 = 295.0 / T;
68  double const beta300_25 = pow(beta300, 2.5);
69  double gh2, zh2, ghe, zhe, gnh3, znh3, gh2o, zh2o, dd, Dinv;
70 
71  if (freq > 30) {
72  gh2 = 1.7465 * PH2 * pow(beta300, 0.8202);
73  ghe = 0.9779 * PHe * pow(beta300, 1.0);
74  gnh3 = 0.7298 * PNH3 * pow(beta295, 1.0);
75  gh2o = 4.8901 * PH2O * pow(beta300, 0.5);
76 
77  zh2 = 1.2163 * PH2 * pow(beta300, 0.8873);
78  // 161014: Bellotti says zhe should be zero to prevent anomolous behavior at
79  // high temp
80  // zhe = 0.0291 * PHe * pow(beta300,0.8994);
81  zhe = 0.0;
82  znh3 = 0.5152 * PNH3 * pow(beta295, 0.6667);
83  zh2o = 2.731 * PH2O * pow(beta300, 1.0);
84 
85  dd = -0.0627; // pressure shift factor from linewidth
86  Dinv = 0.9862; // unitless emperical scale factor
87  } else {
88  // coefficients from Bellotti (private correspondence 160407)
89 
90  gh2 = 1.69369709 * PH2 * pow(beta300, 0.80847492);
91  ghe = 0.69968615 * PHe * pow(beta300, 1.);
92  gnh3 = 0.75232086 * PNH3 * pow(beta295, 1.0);
93  // gh2o = 4.8901 * PH2O * pow(beta300,0.5);
94  gh2o = 5.2333 * PH2O * pow(beta300, 0.6224);
95 
96  zh2 = 1.32634109 * PH2 * pow(beta300, 0.81994711);
97  // 161014: Bellotti says zhe should be zero to prevent anomolous behavior at
98  // high temp
99  // zhe = 0.16068662 * PHe * pow(beta300,-0.72691007);
100  zhe = 0.0;
101  znh3 = 0.61622276 * PNH3 * pow(beta295, 1.38318269);
102  // zh2o = 2.731 * PH2O * pow(beta300,1.0);
103  zh2o = 5.2333 * PH2O * pow(beta300, 2.1248);
104 
105  dd = -0.01386851; // pressure shift factor from linewidth
106  Dinv = 0.96190579; // unitless emperical scale factor
107  }
108 
109  // Inversion:
110  double alpha_inv(0.0);
111  int const Ninv = sizeof(NH3_inv_JPL5) / sizeof(NH3_inv_JPL5[0]);
112  for (int i = 0; i < Ninv; i++) {
113  LineParameters_NH3_inversion_JPL5 const& x = NH3_inv_JPL5[i];
114  double f = x.f0; // convert to GHz
115  double fx = freq / f;
116  double sjkt = x.I0 * beta300_25 * exp(hckt0 * (x.E0) * (1.0 - beta300));
117 
118  // pressure-broadening widths:
119  double gam = gh2 + ghe + gnh3 * x.gamma_self;
120  double zet = zh2 + zhe + znh3 * x.gamma_self;
121 
122  if (PH2O != 0.0) { // Broadening line width's due to water Vapor:
123  gam += gh2o;
124  zet += zh2o;
125  }
126 
127  // Ben-Reuven line profile (eqn 2.6):
128  double num =
129  (gam - zet) * SQUARE(freq) +
130  (gam + zet) * (SQUARE(f + dd * gam) + SQUARE(gam) - SQUARE(zet));
131  double den = SQUARE(SQUARE(freq) - SQUARE(f + dd * gam) - SQUARE(gam) +
132  SQUARE(zet)) +
133  SQUARE(2.0 * freq * gam);
134  double fline = cmghz * 2.0 * fx * num / (pi * den);
135 
136  alpha_inv += Dinv * dens * sjkt * fx * fline;
137  }
138 
139  double alpha_final;
140 
141  if (freq > 30) {
142  // Rotational:
143  double alpha_rot(0.0);
144  double const gH2_rot = 0.2984 * PH2 * pow(beta300, 0.8730);
145  double const gHe_rot = 0.75 * PHe * pow(beta300, 0.6667);
146  double const gNH3_rot = 3.1789 * PNH3 * pow(beta300, 1.0);
147  double const Drot = 2.4268;
148 
149  int const Nrot = sizeof(NH3_rot_JPL5) / sizeof(NH3_rot_JPL5[0]);
150  for (int i = 0; i < Nrot; i++) {
151  LineParameters_NH3_rotational_JPL5 const& x = NH3_rot_JPL5[i];
152  double f = x.f0; // convert to GHz
153  double fx = freq / f;
154  // pressure-broadening width:
155  double gam = gH2_rot * (x.gamma_H2) + gHe_rot * (x.gamma_He) +
156  gNH3_rot * (x.gamma_self);
157 
158  if (PH2O != 0.0) {
159  gam += gh2o;
160  }
161  // Gross line profile (from Joiner at al.):
162  double fline = 4.0 / pi * cmghz * fx * gam /
163  (SQUARE(f - fx * freq) + SQUARE(2.0 * fx * gam));
164  // the rest is as for Inversion case:
165  alpha_rot += Drot * dens * (x.I0) * fx * beta300_25 *
166  exp(hckt0 * (x.E0) * (1.0 - beta300)) * fline;
167  }
168 
169  // Rotovibrational:
170  double alpha_rv(0.0);
171  double const gH2_nu2 = 1.4 * PH2 * pow(beta300, 0.73);
172  double const gHe_nu2 = 0.68 * PHe * pow(beta300, 0.5716);
173  double const gNH3_nu2 = 9.5 * PNH3 * pow(beta300, 1.0);
174  double const Drv = 1.1206;
175 
176  int const Nrv = sizeof(NH3_rv_JPL5) / sizeof(NH3_rv_JPL5[0]);
177  for (int i = 0; i < Nrv; i++) {
178  LineParameters_NH3_rotovibrational_JPL5 const& x = NH3_rv_JPL5[i];
179  double f = x.f0; // convert to GHz
180  double fx = freq / f;
181  // pressure-broadening width:
182  double gam = gH2_nu2 + gHe_nu2 + gNH3_nu2;
183 
184  if (PH2O != 0.0) {
185  gam += gh2o;
186  }
187  // Gross line profile (from Joiner at al.):
188  double fline = 4.0 / pi * cmghz * fx * gam /
189  (SQUARE(f - fx * freq) + SQUARE(2.0 * fx * gam));
190  // the rest is as for Inversion case:
191  alpha_rv += Drv * dens * (x.I0) * fx * beta300_25 *
192  exp(hckt0 * (x.E0) * (1.0 - beta300)) * fline;
193  }
194  alpha_final = alpha_inv + alpha_rot + alpha_rv; // in cm^-1
195 
196  } else {
197  alpha_final = alpha_inv;
198  }
199 
200  return alpha_final;
201 }
202 
203 double attenuation_NH3_Bellotti(double freq, double P, double P_idl, double T,
204  double XH2, double XHe, double XNH3,
205  double XH2O)
206 // sum the contributions of all lines, first the Inversion transitions,
207 // then the Rotational transitions
208 
209 {
210 #if defined(FAKE_CUTOFF_TO_AVOID_SINGULARITY_IN_BELLOTTI_ABS)
211  // Bellotti's parameterization suffers from two problems at high pressure:
212  // (1) for P~7000 bar, absorption becomes negative!!!
213  // (2) slope of opacity wrt pressure changes sign around P~1000 bar.
214  // We patch this by readjusting opacity for T>T_cutoff. (We choose a
215  // temperature cutoff instead of a pressure cutoff, since the opacity is more
216  // strongly dependent on temperature.) Update: this anomalous behavior has
217  // been corrected by zeroing zhe
218  constexpr double T_Bellotti_cutoff = 1250;
219  if (T > T_Bellotti_cutoff) {
220  T = T_Bellotti_cutoff;
221  }
222 #endif
223 
224  // TEMPORARY KLUDGE to compare to Hanley, fig19!!!
225  // T=295;
226  // XNH3=0.004366812227074;
227 
228  // miscellaneous constants:
229  double const ppscale = 2.99792458e18; // converts to cm^2*cm-1 from nm^2*MHz
230  double const gcon = 25.923;
231  double const cmghz = 29.9792458; // converts to GHz from cm-1
232  double const hckt0 = 0.0047959232; // hc/kT0, where T0=300K
233  double const cdens =
234  7.242971565e21; // 1/kB * (1e5 Pa/bar) * (1e-6 m^-3/cm^-3)
235 
236  // partial pressures in bars
237  double const PNH3 = P_idl * XNH3;
238  double const PH2 = P_idl * XH2;
239  double const PHe = P_idl * XHe;
240  double const PH2O = P_idl * XH2O;
241 
242  double const dens = cdens * (P_idl * XNH3) / T;
243  // double const dens = cdens*PNH3/T;
244  double const beta300 = 300.0 / T;
245  double const beta295 = 295.0 / T;
246  double const beta300_25 = pow(beta300, 2.5);
247  double gh2, zh2, ghe, zhe, gnh3, znh3, gh2o, zh2o, dd, Dinv;
248 
249  // coefficients from Bellotti (private correspondence 160407)
250  gh2 = 1.6948 * PH2 * pow(beta300, 1.0);
251  ghe = 0.6385 * PHe * pow(beta300, 0.7995);
252  gnh3 = 0.7 * PNH3 * pow(beta295, 1.0);
253  // gh2o = 4.8901 * PH2O * pow(beta300,0.5);
254  gh2o = 5.2333 * PH2O * pow(beta300, 0.6224);
255 
256  zh2 = 1.2718 * PH2 * pow(beta300, 0.9988);
257  zhe = 0.1316 * PHe * pow(beta300, -0.9129);
258  znh3 = 0.4879 * PNH3 * pow(beta295, 0.7885);
259  // zh2o = 2.731 * PH2O * pow(beta300,1.0);
260  zh2o = 5.2333 * PH2O * pow(beta300, 2.1248);
261 
262  dd = 0.; // pressure shift factor from linewidth
263  Dinv = 1.0323; // unitless emperical scale factor
264 
265  // Inversion:
266  double alpha_inv(0.0);
267  int const Ninv = sizeof(NH3_inv_JPL5) / sizeof(NH3_inv_JPL5[0]);
268  for (int i = 0; i < Ninv; i++) {
269  LineParameters_NH3_inversion_JPL5 const& x = NH3_inv_JPL5[i];
270  double f = x.f0; // convert to GHz
271  double fx = freq / f;
272  double sjkt = x.I0 * beta300_25 * exp(hckt0 * (x.E0) * (1.0 - beta300));
273 
274  // pressure-broadening widths:
275  double gam = gh2 + ghe + gnh3 * x.gamma_self;
276  double zet = zh2 + zhe + znh3 * x.gamma_self;
277 
278  if (PH2O != 0.0) { // Broadening line width's due to water Vapor:
279  gam += gh2o;
280  zet += zh2o;
281  }
282 
283  // Ben-Reuven line profile (eqn 2.6):
284  double num =
285  (gam - zet) * SQUARE(freq) +
286  (gam + zet) * (SQUARE(f + dd * gam) + SQUARE(gam) - SQUARE(zet));
287  double den = SQUARE(SQUARE(freq) - SQUARE(f + dd * gam) - SQUARE(gam) +
288  SQUARE(zet)) +
289  SQUARE(2.0 * freq * gam);
290  double fline = cmghz * 2.0 * fx * num / (pi * den);
291 
292  alpha_inv += Dinv * dens * sjkt * fx * fline;
293  }
294  // Rotational:
295  double alpha_rot(0.0);
296 
297  double const gH2_rot = 1.7761 * PH2 * pow(beta300, 0.5);
298  double const gHe_rot = 0.6175 * PHe * pow(beta300, 0.5663);
299  double const gNH3_rot = 3.1518 * PNH3 * pow(beta300, 1.0);
300  double const Drot = 2.7252;
301 
302  int const Nrot = sizeof(NH3_rot_JPL5) / sizeof(NH3_rot_JPL5[0]);
303  for (int i = 0; i < Nrot; i++) {
304  LineParameters_NH3_rotational_JPL5 const& x = NH3_rot_JPL5[i];
305  double f = x.f0; // convert to GHz
306  double fx = freq / f;
307  // pressure-broadening width:
308  double gam = gH2_rot * (x.gamma_H2) + gHe_rot * (x.gamma_He) +
309  gNH3_rot * (x.gamma_self);
310  if (PH2O != 0.0) {
311  gam += gh2o;
312  }
313  // Gross line profile (from Joiner at al.):
314  double fline = 4.0 / pi * cmghz * fx * gam /
315  (SQUARE(f - fx * freq) + SQUARE(2.0 * fx * gam));
316  // the rest is as for Inversion case:
317  alpha_rot += Drot * dens * (x.I0) * fx * beta300_25 *
318  exp(hckt0 * (x.E0) * (1.0 - beta300)) * fline;
319  }
320 
321  // Rotovibrational:
322  double alpha_rv(0.0);
323 
324  double const gH2_nu2 = 0.5982 * PH2 * pow(beta300, 0.5);
325  double const gHe_nu2 = 0.6175 * PHe * pow(beta300, 0.5505);
326  double const gNH3_nu2 = 5.0894 * PNH3 * pow(beta300, 0.9996);
327  double const Drv = 0.7286;
328 
329  int const Nrv = sizeof(NH3_rv_JPL5) / sizeof(NH3_rv_JPL5[0]);
330  for (int i = 0; i < Nrv; i++) {
331  LineParameters_NH3_rotovibrational_JPL5 const& x = NH3_rv_JPL5[i];
332  double f = x.f0; // convert to GHz
333  double fx = freq / f;
334  // pressure-broadening width:
335  double gam = gH2_nu2 + gHe_nu2 + gNH3_nu2;
336  if (PH2O != 0.0) {
337  gam += gh2o;
338  }
339  // Gross line profile (from Joiner at al.):
340  double fline = 4.0 / pi * cmghz * fx * gam /
341  (SQUARE(f - fx * freq) + SQUARE(2.0 * fx * gam));
342  // the rest is as for Inversion case:
343  alpha_rv += Drv * dens * (x.I0) * fx * beta300_25 *
344  exp(hckt0 * (x.E0) * (1.0 - beta300)) * fline;
345  }
346 
347  return alpha_inv + alpha_rot + alpha_rv; // in cm^-1
348 }
349 
350 double attenuation_NH3_Devaraj(double freq, double P, double P_idl, double T,
351  double XH2, double XHe, double XNH3, double XH2O,
352  int version)
353 // sum the contributions of all lines, first the Inversion transitions,
354 // then the Rotational transitions
355 
356 {
357  // TEMPORARY KLUDGE to compare to Hanley, fig19!!!
358  // T=295;
359  // XNH3=0.004366812227074;
360 
361  // versioning
362  bool use_v2011(false);
363  bool use_high_pressure_parameters(P_idl >= 15 ? true : false);
364  bool use_low_pressure_parameters(!use_high_pressure_parameters);
365 
366  switch (version) {
367  case 0: // use defaults; do not override
368  break;
369  case -1: // use version 2011
370  use_v2011 = true;
371  use_high_pressure_parameters = false;
372  use_low_pressure_parameters = false;
373  break;
374  case 1: // use high pressure parameters
375  use_high_pressure_parameters = true;
376  use_low_pressure_parameters = false;
377  break;
378  case 2: // use low pressure parameters
379  use_high_pressure_parameters = false;
380  use_low_pressure_parameters = true;
381  break;
382  }
383 
384  // miscellaneous constants:
385  double const ppscale = 2.99792458e18; // converts to cm^2*cm-1 from nm^2*MHz
386  double const gcon = 25.923;
387  double const cmghz = 29.9792458; // converts to GHz from cm-1
388  double const hckt0 = 0.0047959232; // hc/kT0, where T0=300K
389  double const cdens =
390  7.242971565e21; // 1/kB * (1e5 Pa/bar) * (1e-6 m^-3/cm^-3)
391 
392  // partial pressures in bars
393  double const PNH3 = P_idl * XNH3;
394  double const PH2 = P_idl * XH2;
395  double const PHe = P_idl * XHe;
396  double const PH2O = P_idl * XH2O;
397 
398  double const dens = cdens * (P_idl * XNH3) / T;
399  // double const dens = cdens*PNH3/T;
400  double const beta300 = 300.0 / T;
401  double const beta295 = 295.0 / T;
402  double const beta300_25 = pow(beta300, 2.5);
403 
404  // Inversion:
405  double alpha_inv(0.0);
406  double gh2, zh2, ghe, zhe, gnh3, znh3, gh2o, zh2o, dd, Dinv;
407  if (use_v2011) {
408  gh2 = 1.7947 * PH2 * pow(beta300, 0.8357);
409  ghe = 0.75 * PHe * pow(beta300, 0.6667);
410  gnh3 = 0.7719 * PNH3 * pow(beta295, 1.0);
411  gh2o = 4.8901 * PH2O * pow(beta300, 0.5);
412 
413  zh2 = 1.2031 * PH2 * pow(beta300, 0.8610);
414  zhe = 0.3 * PHe * pow(beta300, 0.6667);
415  znh3 = 0.5620 * PNH3 * pow(beta295, 0.6206);
416  zh2o = 2.731 * PH2O * pow(beta300, 1.0);
417 
418  dd = -0.0404; // pressure shift factor from linewidth
419  Dinv = 0.9903; // unitless emperical scale factor
420  } else {
421  if (use_low_pressure_parameters) {
422  gh2 = 1.7465 * PH2 * pow(beta300, 0.8202);
423  ghe = 0.9779 * PHe * pow(beta300, 1.0);
424  gnh3 = 0.7298 * PNH3 * pow(beta295, 1.0);
425  gh2o = 4.8901 * PH2O * pow(beta300, 0.5);
426 
427  zh2 = 1.2163 * PH2 * pow(beta300, 0.8873);
428  zhe = 0.0291 * PHe * pow(beta300, 0.8994);
429  znh3 = 0.5152 * PNH3 * pow(beta295, 0.6667);
430  zh2o = 2.731 * PH2O * pow(beta300, 1.0);
431 
432  dd = -0.0627; // pressure shift factor from linewidth
433  Dinv = 0.9862; // unitless emperical scale factor
434  } else if (use_high_pressure_parameters) {
435  gh2 = 1.6361 * PH2 * pow(beta300, 0.8);
436  ghe = 0.4555 * PHe * pow(beta300, 0.5);
437  gnh3 = 0.7298 * PNH3 * pow(beta295, 1.0);
438  gh2o = 4.8901 * PH2O * pow(beta300, 0.5);
439 
440  zh2 = 1.1313 * PH2 * pow(beta300, 0.6234);
441  zhe = 0.1 * PHe * pow(beta300, 0.5);
442  znh3 = 0.5152 * PNH3 * pow(beta295, 0.6667);
443  zh2o = 2.731 * PH2O * pow(beta300, 1.0);
444 
445  dd = 0.2; // pressure shift factor from linewidth
446  Dinv = 1.3746; // unitless emperical scale factor
447  }
448  }
449 
450  int const Ninv = sizeof(NH3_inv_JPL5) / sizeof(NH3_inv_JPL5[0]);
451  for (int i = 0; i < Ninv; i++) {
452  LineParameters_NH3_inversion_JPL5 const& x = NH3_inv_JPL5[i];
453  double f = x.f0; // convert to GHz
454  double fx = freq / f;
455  double sjkt = x.I0 * beta300_25 * exp(hckt0 * (x.E0) * (1.0 - beta300));
456 
457  // pressure-broadening widths:
458  double gam = gh2 + ghe + gnh3 * x.gamma_self;
459  double zet = zh2 + zhe + znh3 * x.gamma_self;
460  if (PH2O != 0.0) { // Broadening line width's due to water Vapor:
461  gam += gh2o;
462  zet += zh2o;
463  }
464 
465  // Ben-Reuven line profile (eqn 2.6):
466  double num =
467  (gam - zet) * SQUARE(freq) +
468  (gam + zet) * (SQUARE(f + dd * gam) + SQUARE(gam) - SQUARE(zet));
469  double den = SQUARE(SQUARE(freq) - SQUARE(f + dd * gam) - SQUARE(gam) +
470  SQUARE(zet)) +
471  SQUARE(2.0 * freq * gam);
472  double fline = cmghz * 2.0 * fx * num / (pi * den);
473 
474  alpha_inv += Dinv * dens * sjkt * fx * fline;
475  }
476 
477  // Rotational:
478  double alpha_rot(0.0);
479  double const gH2_rot = 0.2984 * PH2 * pow(beta300, 0.8730);
480  double const gHe_rot = 0.75 * PHe * pow(beta300, 0.6667);
481  double const gNH3_rot = 3.1789 * PNH3 * pow(beta300, 1.0);
482  double const Drot = 2.4268;
483 
484  int const Nrot = sizeof(NH3_rot_JPL5) / sizeof(NH3_rot_JPL5[0]);
485  for (int i = 0; i < Nrot; i++) {
486  LineParameters_NH3_rotational_JPL5 const& x = NH3_rot_JPL5[i];
487  double f = x.f0; // convert to GHz
488  double fx = freq / f;
489  // pressure-broadening width:
490  double gam = gH2_rot * (x.gamma_H2) + gHe_rot * (x.gamma_He) +
491  gNH3_rot * (x.gamma_self);
492  if (PH2O != 0.0) {
493  gam += gh2o;
494  }
495  // Gross line profile (from Joiner at al.):
496  double fline = 4.0 / pi * cmghz * fx * gam /
497  (SQUARE(f - fx * freq) + SQUARE(2.0 * fx * gam));
498  // the rest is as for Inversion case:
499  alpha_rot += Drot * dens * (x.I0) * fx * beta300_25 *
500  exp(hckt0 * (x.E0) * (1.0 - beta300)) * fline;
501  }
502 
503  // Rotovibrational:
504  double alpha_rv(0.0);
505  double const gH2_nu2 = 1.4 * PH2 * pow(beta300, 0.73);
506  double const gHe_nu2 = 0.68 * PHe * pow(beta300, 0.5716);
507  double const gNH3_nu2 = 9.5 * PNH3 * pow(beta300, 1.0);
508  double const Drv = 1.1206;
509 
510  int const Nrv = sizeof(NH3_rv_JPL5) / sizeof(NH3_rv_JPL5[0]);
511  for (int i = 0; i < Nrv; i++) {
512  LineParameters_NH3_rotovibrational_JPL5 const& x = NH3_rv_JPL5[i];
513  double f = x.f0; // convert to GHz
514  double fx = freq / f;
515  // pressure-broadening width:
516  double gam = gH2_nu2 + gHe_nu2 + gNH3_nu2;
517  if (PH2O != 0.0) {
518  gam += gh2o;
519  }
520  // Gross line profile (from Joiner at al.):
521  double fline = 4.0 / pi * cmghz * fx * gam /
522  (SQUARE(f - fx * freq) + SQUARE(2.0 * fx * gam));
523  // the rest is as for Inversion case:
524  alpha_rv += Drv * dens * (x.I0) * fx * beta300_25 *
525  exp(hckt0 * (x.E0) * (1.0 - beta300)) * fline;
526  }
527 
528  return alpha_inv + alpha_rot + alpha_rv; // in cm^-1
529 }
530 
531 double attenuation_NH3_Hanley(double freq, double P, double P_idl, double T,
532  double XH2, double XHe, double XNH3, double XH2O,
533  double power) {
534  // KLUDGE to compare to Hanley, fig19!!!
535  // P=P_idl;
536  // T=295;
537  // XNH3=0.004366812227074;
538 
539  static double hif0[NINV], hie0[NINV], hiint[NINV], hrf0[NROT], hre0[NROT],
540  hrint[NROT];
541 
542  int hij[NINV], hik[NINV], hrj[NROT], hrk[NROT];
543  static bool first = true;
544 
545  // for line broadening (Table 4.7):
546  double const gh2 = 1.64;
547  double const cgh2 = 0.7756;
548  double const zh2 = 1.262;
549  double const czh2 = 0.7964;
550 
551  double const ghe = 0.75;
552  double const cghe = 0.6667;
553  double const zhe = 0.3;
554  double const czhe = 0.6667;
555 
556  double const gnh3 = 0.852;
557  double const cgnh3 = 1.0;
558  double const znh3 = 0.5296;
559  double const cznh3 = 1.554;
560 
561  // miscellaneous constants:
562  double const ppscale = 3.00e18; // converts to cm^2*cm-1 from nm^2*MHz
563  double const gcon = 25.923;
564  double const cmghz = 29.979; // converts to GHz from cm-1
565 
566  double const dd = -0.0498; // pressure shift factor from linewidth
567  double const ddd = 0.9301; // unitless emperical scale factor
568  double const hckt0 = 0.004796; // hc/kT0, where T0=300K
569  double const cdens = 7.244e21; // 1/kB * (1e5 Pa/bar) * (1e-6 m^-3/cm^-3)
570  // double const cdens = 7.242971565e21; // 1/kB * (1e5 Pa/bar) * (1e-6
571  // m^-3/cm^-3)
572 
573  // on first call, load the arrays from the tables:
574  if (first || 1) {
575  for (int i = 0; i < NINV; i++) {
576  hif0[i] = 1.0e-3 * hireals[i][0]; // convert from MHz to GHz
577  // intensities from Poynter Pickett catalog:
578  hiint[i] = pow(10.0, hireals[i][2]) / ppscale;
579  hie0[i] = hireals[i][3];
580  hij[i] = hiints[i][3];
581  hik[i] = hiints[i][4];
582  }
583  for (int i = 0; i < NROT; i++) {
584  hrf0[i] = 1.0e-3 * hrreals[i][0]; // convert from MHz to GHz
585  // intensities from Poynter Pickett catalog:
586  hrint[i] = pow(10.0, hrreals[i][2]) / ppscale;
587  hre0[i] = hrreals[i][3];
588  hrj[i] = hrints[i][3];
589  hrk[i] = hrints[i][4];
590  }
591  first = false;
592  }
593 
594  // partial pressures in bars
595  double const PNH3 = P_idl * XNH3;
596  double const PH2 = P_idl * XH2;
597  double const PHe = P_idl * XHe;
598  double const PH2O = P_idl * XH2O;
599 
600  // double const dens = cdens*PNH3/T;
601  double const dens = cdens * (P_idl * XNH3) / T;
602  double const beta300 = 300.0 / T;
603  double const beta295 = 295.0 / T;
604 
605  double const beta300_25 = pow(beta300, 2.5);
606  double const beta300_cgh2 = pow(beta300, cgh2);
607  double const beta300_cghe = pow(beta300, cghe);
608  double const beta300_czh2 = pow(beta300, czh2);
609  double const beta300_czhe = pow(beta300, czhe);
610  double const beta295_cgnh3 = pow(beta295, cgnh3);
611  double const beta295_cznh3 = pow(beta295, cznh3);
612 
613  // water broadening parameters
614  double const gh2o = 4.8901;
615  double const cgh2o = 0.5;
616  double const zh2o = 2.731;
617  double const czh2o = 1;
618  double const beta300_cgh2o = pow(beta300, cgh2o);
619  double const beta300_czh2o = pow(beta300, czh2o);
620 
621  // sum the contributions of all lines, first innversion then rotational
622  double dtau(0.0);
623  double dtau_inv(0.0);
624  double dtau_rot(0.0);
625 
626  // Inversion:
627  for (int i = 0; i < NINV; ++i) {
628  double fx = freq / hif0[i];
629  double sjkt =
630  hiint[i] * beta300_25 * exp(hckt0 * hie0[i] * (1.0 - beta300));
631  // pressure-broadening widths:
632  double gam = gh2 * PH2 * beta300_cgh2 + ghe * PHe * beta300_cghe;
633  double zet = zh2 * PH2 * beta300_czh2 + zhe * PHe * beta300_czhe;
634  if (PH2O != 0.0) {
635  gam += gh2o * PH2O * beta300_cgh2o;
636  zet += zh2o * PH2O * beta300_czh2o;
637  }
638 
639  // self-broadening line width:
640  double gam0 = gcon * double(hik[i]) / sqrt(double(hij[i] * (hij[i] + 1)));
641  gam += gnh3 * gam0 * PNH3 * beta295_cgnh3;
642  zet += znh3 * gam0 * PNH3 * beta295_cznh3;
643 
644  // Ben-Reuven line profile (eqn 2.6):
645  double num =
646  (gam - zet) * SQUARE(freq) +
647  (gam + zet) * (SQUARE(hif0[i] + dd * gam) + SQUARE(gam) - SQUARE(zet));
648  double den = SQUARE(SQUARE(freq) - SQUARE(hif0[i] + dd * gam) -
649  SQUARE(gam) + SQUARE(zet)) +
650  SQUARE(2.0 * freq * gam);
651  double fline = cmghz * 2.0 * fx * num / (pi * den);
652  double alpha = ddd * dens * sjkt * fx * fline;
653  dtau_inv += alpha;
654  dtau += alpha;
655  }
656 
657  double const beta296 = 296 / T;
658  double const gH2_rot = PH2 * pow(beta296, cgh2);
659  double const gHe_rot = PHe * pow(beta296, cghe);
660  double const gNH3_rot = PNH3 * pow(beta296, cgnh3);
661 
662  // Rotational:
663  for (int i = 0; i < NROT; i++) {
664  LineParameters_NH3_rotational_JPL5 const& x = NH3_rot_Hanley[i];
665  double fx = freq / hrf0[i];
666 
667  double gam = atm_from_bar(gH2_rot * (x.gamma_H2) + gHe_rot * (x.gamma_He) +
668  gNH3_rot * (x.gamma_self));
669  if (PH2O != 0.0) {
670  gam += gh2o * PH2O * beta300_cgh2o;
671  }
672 
673  // Gross line profile (from Joiner at al.):
674  double fline = 4.0 / pi * cmghz * fx * gam /
675  (SQUARE(hrf0[i] - fx * freq) + SQUARE(2.0 * fx * gam));
676  // the rest is as for Inversion case:
677  double dtaua = ddd * dens * hrint[i] * fx * beta300_25 *
678  exp(hckt0 * hre0[i] * (1.0 - beta300)) * fline;
679  dtau_rot += dtaua;
680  dtau += dtaua;
681  }
682 
683 #ifdef ABSORPTION_KLUDGE
684  constexpr double P_top = 150;
685  constexpr double P_bot = 4000;
686  constexpr double scale_bot = 2.0;
687  if (P > P_bot)
688  dtau *= scale_bot;
689  else if (P > P_top)
690  dtau *= ((P - P_top) * scale_bot + (P_bot - P)) / (P_bot - P_top);
691 #endif
692 
693  if (T > 750.) dtau *= pow(750. / T, power);
694 
695  return dtau; // in cm^-1
696 }
697 
698 // helper function for attenuation_NH3_radtran() -- array ordering
699 inline int Index_NH3(int I, int J) {
700  int INDX = J;
701 
702  if (I > 1) {
703  for (int L = 1; L <= I - 1; L++) {
704  INDX += L;
705  }
706  }
707  return INDX;
708 }
709 
710 // helper function for attenuation_NH3_radtran() -- computes Ben-Reuven
711 // lineshape
712 inline double shape(double NU0, double GAMMA, double NU, double DELTA,
713  double ZETA) {
714  double const NUSQ = SQUARE(NU);
715  double const NU0SQ = SQUARE(NU0);
716  double const GSQ = SQUARE(GAMMA);
717  double const ZSQ = SQUARE(ZETA);
718  double const T = (NU0 + DELTA) * (NU0 + DELTA);
719  return 2.0 * GAMMA * (NUSQ / NU0SQ) *
720  (((GAMMA - ZETA) * NUSQ + (GAMMA + ZETA) * (T + GSQ - ZSQ)) /
721  ((NUSQ - T - GSQ + ZSQ) * (NUSQ - T - GSQ + ZSQ) + 4.0 * NUSQ * GSQ));
722 }
723 
724 double attenuation_NH3_radtran(double freq, double P, double T, double XH2,
725  double XHe, double XNH3) {
726  /*
727  SUBROUTINE NH3ABS CALCULATES THE ABSORPTION DUE TO AMMONIA
728  INCLUDING THE PRESSURE BROADENING DUE TO HYDROGEN.
729 
730  PNH3 ..... PARTIAL PRESSURE OF NH3 IN ATM
731  PH2 ...... PARTIAL PRESSURE OF H2 IN ATM
732  T ........ TEMPERATURE IN DEGREES KELVIN
733  FQ ....... FREQUENCY IN GHZ AT WHICH ABSORPTION
734  IS TO BE CALCULATED
735  DTAU ..... ABSORPTION IN 1/CM AT A PARTICULAR FREQ DUE TO ALL
736  NH3 LINES.
737  */
738 
739  double DTAU = 0.0, DTAUA = 0.0;
740 
741  double const PNH3 = atm_from_bar(P * XNH3);
742  double const PH2 = atm_from_bar(P * XH2);
743  double const PHe = atm_from_bar(P * XHe);
744 
745  double const TESTPR = 1.0e-4;
746  double const EXPON = 4.0e0;
747  double const GAMFAC = 1.0e0;
748  double const BOLTZ = 2.084e4;
749  double const XMU = 1.476;
750  double const DROT = -1.09107e5;
751  double const B = 2.98107e5;
752 
753  // COLLSW IS A FUDGE FACTOR WHICH INCREASES THE SELF-BROADENING
754  // OF NH3 INVERSION LINES ... ... FACNH3 = COLLSW*FACNH3
755  double const collsw = 1.0;
756  double const ZFAC = 1.0;
757  double const DFAC = 1.0;
758 
759  // COMMON/ATMCA/WT,TAU,DTAU,DTB,freq,ANH3SV,AH2OSV,acldsv
760 
761  if (T <= 0.0) throw "Non-positive temperature";
762 
763  double const PCROSS = pow(((PNH3 + PH2 / 8.0) / TESTPR), EXPON);
764  double const SWITCH = ((PCROSS <= 9.0) ? 1.0 - exp(-PCROSS) : 1.0);
765  double const FACT = pow(T, 3.5);
766  double const TBOLTZ = T * BOLTZ;
767  double const FACH2 = 2.318 * PH2 * pow((300.0 / T), (2. / 3.));
768  double const FACHE = (35.403 * PHe) / pow(T, (2. / 3.));
769  double const FNH3 = collsw * 224.2 * PNH3 / T;
770  double const DELTA = -SWITCH * 0.45 * PNH3 * DFAC;
771 
772  // FACMP IS A FUDGE FACTOR DETERMINED BY COMPARISON OF THIS ROUTINE
773  // WITH THE RESULTS OF MORRIS AND PARSONS.
774  double const FACMP =
775  1.00745 + 3.08301e-2 * (PH2 / T) + 5.52505e-2 * pow((PH2 / T), 2.0);
776 
777  for (int J = 1; J <= 16; J++) {
778  double ENJ = B * J * (J + 1);
779  double FJMUSQ = (2 * J + 1) * XMU * XMU / (J * (J + 1));
780 
781  for (int K = 1; K <= J; K++) {
782  int I = Index_NH3(J, K) - 1;
783  if (FI[I] <= 1.0 || GAMIN[I] <= 1.0) continue;
784 
785  double EN = -(ENJ + DROT * K * K) / TBOLTZ;
786  double FXMUSQ = FJMUSQ * K * K;
787  double FACNH3 = FNH3 * GAMIN[I];
788  double GAMMA = GAMFAC * (FACNH3 + FACH2 + FACHE);
789  double ZETA =
790  SWITCH * (0.655 * FACNH3 + 0.83 * FACH2 + 0.3797 * FACHE) * ZFAC;
791  double GK = (K % 3 == 0 ? 3.0 : 1.5);
792 
793  double TEXP = 1.132e3 * GK * FXMUSQ * FI[I] * FI[I] * PNH3 * exp(EN) /
794  (GAMMA * FACT);
795  double R = shape(FI[I], GAMMA, freq, DELTA, ZETA);
796 
797  // MAKE FURTHER ADJUSTMENT BECAUSE ONLY 1/2 OF MOLECULES ABSORB
798  DTAUA = 0.5 * FACMP * TEXP * R;
799  DTAU += DTAUA;
800  }
801  }
802 
803  // disable this check, as it excludes these transitions for our freqs:
804  // if (freq <= FI[6]) return DTAU;
805 
806  // INCLUDE THE ABSORPTION DUE TO ROTATIONAL TRANSITIONS
807  double ROTFAC = 1.0;
808  double FACNH3 = 20.0 * FNH3;
809  double GAMMA = GAMFAC * (FACNH3 + FACH2 + FACHE);
810  double ZETA =
811  SWITCH * (0.655 * FACNH3 + 0.83 * FACH2 + 0.3797 * FACHE) * ZFAC;
812 
813  for (int L = 0; L < 49; L++) {
814  int const J = JROT[L];
815  int const K = KROT[L];
816  double const F = FROT[L];
817 
818  double R = shape(F, GAMMA, freq, DELTA, ZETA);
819  double FXMUSQ = XMU * XMU * ((J + 1) * (J + 1) - K * K) / (J + 1);
820  double EN = -(B * J * (J + 1) + DROT * K * K) / TBOLTZ;
821 
822  double GK = (K % 3 == 0 ? 3.0 : 1.5);
823  if (K == 0) GK *= 0.5;
824 
825  double TEXP =
826  1.132e3 * GK * FXMUSQ * F * F * PNH3 * exp(EN) / (GAMMA * FACT);
827  DTAUA = 0.5 * ROTFAC * FACMP * TEXP * R;
828  DTAU += DTAUA;
829  }
830 
831  return DTAU;
832 }
atm_from_bar
double atm_from_bar(double P)
Definition: attenuation_nh3.cpp:21
pi
static double const pi
Definition: attenuation_nh3.cpp:22
shape
double shape(double NU0, double GAMMA, double NU, double DELTA, double ZETA)
Definition: attenuation_nh3.cpp:712
attenuation_NH3_Bellotti
double attenuation_NH3_Bellotti(double freq, double P, double P_idl, double T, double XH2, double XHe, double XNH3, double XH2O)
Definition: attenuation_nh3.cpp:203
attenuation_NH3_radtran
double attenuation_NH3_radtran(double freq, double P, double T, double XH2, double XHe, double XNH3)
Definition: attenuation_nh3.cpp:724
attenuation_NH3_Hanley
double attenuation_NH3_Hanley(double freq, double P, double P_idl, double T, double XH2, double XHe, double XNH3, double XH2O, double power)
Definition: attenuation_nh3.cpp:531
SQUARE
double SQUARE(double x)
Definition: attenuation_nh3.cpp:20
attenuation_NH3_Bellotti_switch
double attenuation_NH3_Bellotti_switch(double freq, double P, double P_idl, double T, double XH2, double XHe, double XNH3, double XH2O)
Definition: attenuation_nh3.cpp:25
Index_NH3
int Index_NH3(int I, int J)
Definition: attenuation_nh3.cpp:699
attenuation_NH3_Devaraj
double attenuation_NH3_Devaraj(double freq, double P, double P_idl, double T, double XH2, double XHe, double XNH3, double XH2O, int version)
Definition: attenuation_nh3.cpp:350
LineParameters_NH3_inversion_JPL5
Definition: mwr_absorber_nh3.hpp:729
LineParameters_NH3_inversion_JPL5::f0
double f0
Definition: mwr_absorber_nh3.hpp:731
LineParameters_NH3_inversion_JPL5::I0
double I0
Definition: mwr_absorber_nh3.hpp:731
LineParameters_NH3_inversion_JPL5::gamma_self
double gamma_self
Definition: mwr_absorber_nh3.hpp:731
LineParameters_NH3_inversion_JPL5::E0
double E0
Definition: mwr_absorber_nh3.hpp:731
LineParameters_NH3_rotational_JPL5
Definition: mwr_absorber_nh3.hpp:734
LineParameters_NH3_rotational_JPL5::gamma_He
double gamma_He
Definition: mwr_absorber_nh3.hpp:736
LineParameters_NH3_rotational_JPL5::gamma_self
double gamma_self
Definition: mwr_absorber_nh3.hpp:736
LineParameters_NH3_rotational_JPL5::gamma_H2
double gamma_H2
Definition: mwr_absorber_nh3.hpp:736
LineParameters_NH3_rotational_JPL5::E0
double E0
Definition: mwr_absorber_nh3.hpp:736
LineParameters_NH3_rotational_JPL5::I0
double I0
Definition: mwr_absorber_nh3.hpp:736
LineParameters_NH3_rotational_JPL5::f0
double f0
Definition: mwr_absorber_nh3.hpp:736
LineParameters_NH3_rotovibrational_JPL5
Definition: mwr_absorber_nh3.hpp:739
LineParameters_NH3_rotovibrational_JPL5::E0
double E0
Definition: mwr_absorber_nh3.hpp:741
LineParameters_NH3_rotovibrational_JPL5::I0
double I0
Definition: mwr_absorber_nh3.hpp:741
LineParameters_NH3_rotovibrational_JPL5::f0
double f0
Definition: mwr_absorber_nh3.hpp:741
mwr_absorber_nh3.hpp
JROT
static int const JROT[49]
Definition: mwr_absorber_nh3.hpp:670
hrints
static double const hrints[NROT][9]
Definition: mwr_absorber_nh3.hpp:473
hrreals
static double const hrreals[NROT][4]
Definition: mwr_absorber_nh3.hpp:153
NH3_inv_JPL5
static LineParameters_NH3_inversion_JPL5 const NH3_inv_JPL5[]
Definition: mwr_absorber_nh3.hpp:1006
FI
static double const FI[136]
Definition: mwr_absorber_nh3.hpp:691
NH3_rv_JPL5
static LineParameters_NH3_rotovibrational_JPL5 const NH3_rv_JPL5[]
Definition: mwr_absorber_nh3.hpp:4442
NINV
int const NINV
Definition: mwr_absorber_nh3.hpp:16
hireals
static double const hireals[NINV][4]
Definition: mwr_absorber_nh3.hpp:23
NROT
int const NROT
Definition: mwr_absorber_nh3.hpp:17
NH3_rot_Hanley
static LineParameters_NH3_rotational_JPL5 const NH3_rot_Hanley[]
Definition: mwr_absorber_nh3.hpp:746
FROT
static double const FROT[49]
Definition: mwr_absorber_nh3.hpp:679
NH3_rot_JPL5
static LineParameters_NH3_rotational_JPL5 const NH3_rot_JPL5[]
Definition: mwr_absorber_nh3.hpp:1838
KROT
static int const KROT[49]
Definition: mwr_absorber_nh3.hpp:674
GAMIN
static double const GAMIN[136]
Definition: mwr_absorber_nh3.hpp:715
hiints
static double const hiints[NINV][9]
Definition: mwr_absorber_nh3.hpp:327
K
Real K
Definition: straka.cpp:66