ARTS 2.5.11 (git: 725533f0)
psd.cc
Go to the documentation of this file.
1
9#include "psd.h"
10
11/*===========================================================================
12 === External declarations
13 ===========================================================================*/
14#include <algorithm>
15#include <cmath>
16#include <ctime>
17#include <limits>
18#include <random>
19#include <stdexcept>
20
21#include "arts.h"
22#include "arts_constants.h"
23#include "check_input.h"
24#include "cloudbox.h"
25#include "lin_alg.h"
26#include "logic.h"
27#include "math_funcs.h"
28#include "mc_antenna.h"
29#include "messages.h"
30#include "physics_funcs.h"
31#include "ppath.h"
32#include "rng.h"
33#include "sorting.h"
34
35inline constexpr Numeric PI=Constant::pi;
38
39void psd_cloudice_MH97(Vector& psd,
40 const Vector& diameter,
41 const Numeric& iwc,
42 const Numeric& t,
43 const bool noisy) {
44 Index nD = diameter.nelem();
45 psd.resize(nD);
46 psd = 0.;
47
48 ARTS_ASSERT(t > 0.);
49
50 // skip calculation if IWC is 0.0
51 if (iwc == 0.0) {
52 return;
53 }
54 ARTS_ASSERT(iwc > 0.);
55
56 // convert m to microns
57 Vector d_um(nD);
58 for (Index iD = 0; iD < nD; iD++) d_um[iD] = 1e6 * diameter[iD];
59 //convert T from Kelvin to Celsius
60 Numeric Tc = t - 273.15;
61
62 //[kg/m3] -> [g/m3] as used by parameterisation
63 Numeric ciwc = iwc * 1e3;
64 Numeric cdensity = DENSITY_OF_ICE * 1e3;
65
66 Numeric sig_a = 0., sig_b1 = 0.;
67 Numeric sig_b2 = 0., sig_m = 0.;
68 Numeric sig_aamu = 0., sig_bamu = 0.;
69 Numeric sig_abmu = 0., sig_bbmu = 0.;
70 Numeric sig_aasigma = 0., sig_basigma = 0;
71 Numeric sig_absigma = 0., sig_bbsigma = 0.;
72
73 if (noisy) {
74 RandomNumberGenerator<> rng; //Random Number generator
75
76 sig_a = 0.068, sig_b1 = 0.054;
77 sig_b2 = 5.5e-3, sig_m = 0.0029;
78 sig_aamu = 0.02, sig_bamu = 0.0005;
79 sig_abmu = 0.023, sig_bbmu = 0.5e-3;
80 sig_aasigma = 0.02, sig_basigma = 0.5e-3;
81 sig_absigma = 0.023, sig_bbsigma = 4.7e-4;
82
83 sig_a = rng.get<std::normal_distribution>(0.0, sig_a)();
84 sig_b1 = rng.get<std::normal_distribution>(0.0, sig_b1)();
85 sig_b2 = rng.get<std::normal_distribution>(0.0, sig_b2)();
86 sig_m = rng.get<std::normal_distribution>(0.0, sig_m)();
87 sig_aamu = rng.get<std::normal_distribution>(0.0, sig_aamu)();
88 sig_bamu = rng.get<std::normal_distribution>(0.0, sig_bamu)();
89 sig_abmu = rng.get<std::normal_distribution>(0.0, sig_abmu)();
90 sig_bbmu = rng.get<std::normal_distribution>(0.0, sig_bbmu)();
91 sig_aasigma = rng.get<std::normal_distribution>(0.0, sig_aasigma)();
92 sig_basigma = rng.get<std::normal_distribution>(0.0, sig_basigma)();
93 sig_absigma = rng.get<std::normal_distribution>(0.0, sig_absigma)();
94 sig_bbsigma = rng.get<std::normal_distribution>(0.0, sig_bbsigma)();;
95 }
96
97 // Split IWC in small and large size modes
98
99 // Calculate IWC in each mode
100 Numeric a = 0.252 + sig_a; //g/m^3
101 Numeric b1 = 0.837 + sig_b1;
102 Numeric IWCs100 = min(ciwc, a * pow(ciwc, b1));
103 Numeric IWCl100 = ciwc - IWCs100;
104
105 // Gamma distribution component (small mode)
106
107 Numeric b2 = -4.99e-3 + sig_b2; //micron^-1
108 Numeric m = 0.0494 + sig_m; //micron^-1
109 Numeric alphas100 = b2 - m * log10(IWCs100); //micron^-1
110
111 // alpha, and hence dNdD1, becomes NaN if IWC>0.
112 // this should be ensured to not happen before.
113 //
114 // alpha, and hence dNdD1, becomes negative for large IWC.
115 // towards this limit, particles anyway get larger than 100um, i.e., outside
116 // the size region the small-particle mode gamma distrib is intended for.
117 // hence, leave dNdD1 at 0 for those cases.
118 Vector dNdD1(nD, 0.);
119 if (alphas100 > 0.) {
120 Numeric Ns100 = 6 * IWCs100 * pow(alphas100, 5.) /
121 (PI * cdensity * tgamma(5.)); //micron^-5
122 for (Index iD = 0; iD < nD; iD++)
123 dNdD1[iD] = 1e18 * Ns100 * d_um[iD] *
124 exp(-alphas100 * d_um[iD]); //micron^-4 -> m^-3 micron^-1
125 }
126
127 // Log normal distribution component (large mode)
128
129 // for small IWC, IWCtotal==IWC<100 & IWC>100=0.
130 // this will give dNdD2=NaN. avoid that by explicitly setting to 0
131 Vector dNdD2(nD, 0.);
132 if (IWCl100 > 0.) {
133 //FIXME: Do we need to ensure mul100>0 and sigmal100>0?
134
135 Numeric aamu = 5.20 + sig_aamu;
136 Numeric bamu = 0.0013 + sig_bamu;
137 Numeric abmu = 0.026 + sig_abmu;
138 Numeric bbmu = -1.2e-3 + sig_bbmu;
139 Numeric amu = aamu + bamu * Tc;
140 Numeric bmu = abmu + bbmu * Tc;
141 Numeric mul100 = amu + bmu * log10(IWCl100);
142
143 Numeric aasigma = 0.47 + sig_aasigma;
144 Numeric basigma = 2.1e-3 + sig_basigma;
145 Numeric absigma = 0.018 + sig_absigma;
146 Numeric bbsigma = -2.1e-4 + sig_bbsigma;
147 Numeric asigma = aasigma + basigma * Tc;
148 Numeric bsigma = absigma + bbsigma * Tc;
149 Numeric sigmal100 = asigma + bsigma * log10(IWCl100);
150
151 if ((mul100 > 0.) & (sigmal100 > 0.)) {
152 Numeric a1 = 6 * IWCl100; //g/m^3
153 Numeric a2_fac = pow(PI, 3. / 2.) * cdensity * sqrt(2) *
154 exp(3 * mul100 + 9. / 2. * pow(sigmal100, 2)) *
155 sigmal100 * pow(1., 3);
156 //a2 = a2_fac * d_um; //g/m^3/micron^4
157 for (Index iD = 0; iD < nD; iD++)
158 dNdD2[iD] = 1e18 * a1 / (a2_fac * d_um[iD]) *
159 exp(-0.5 * pow((log(d_um[iD]) - mul100) / sigmal100, 2));
160 //micron^-4 -> m^-3 micron^-1
161 }
162 }
163
164 for (Index iD = 0; iD < nD; iD++) {
165 // FIXME: Do we still need this check here? Non-NaN of each mode should
166 // now be ensure by the checks/if-loops for each of the modes separately.
167 // I, JM, think (and hope).
168 //if ( !std::isnan(dNdD1[iD]) && !std::isnan(dNdD2[iD]) )
169 psd[iD] = (dNdD1[iD] + dNdD2[iD]) * 1e6; // m^-3 m^-1
170 }
171}
172
173void psd_mgd_mass_and_something(Matrix& psd_data,
174 Tensor3& dpsd_data_dx,
175 const String& something,
176 const Vector& psd_size_grid,
177 const Vector& pnd_agenda_input_t,
178 const Matrix& pnd_agenda_input,
179 const ArrayOfString& pnd_agenda_input_names,
180 const ArrayOfString& dpnd_data_dx_names,
181 const Numeric& scat_species_a,
182 const Numeric& scat_species_b,
183 const Numeric& n0,
184 const Numeric& mu,
185 const Numeric& la,
186 const Numeric& ga,
187 const Numeric& t_min,
188 const Numeric& t_max,
189 const Index& picky,
190 const Verbosity&) {
191 // Standard checks
193
194 // Additional (basic) checks
195 if (nin < 1 || nin > 4)
196 throw runtime_error(
197 "The number of columns in *pnd_agenda_input* must "
198 "be 2, 3 or 4.");
199 if (scat_species_a <= 0)
200 throw runtime_error("*scat_species_a* should be > 0.");
201 if (scat_species_b <= 0 || scat_species_b >= 5)
202 throw runtime_error("*scat_species_b* should be > 0 and < 5.");
203
204 // Check and determine dependent and fixed parameters
205 const Index n0_depend = (Index)n0 == -999;
206 const Index mu_depend = (Index)mu == -999;
207 const Index la_depend = (Index)la == -999;
208 const Index ga_depend = (Index)ga == -999;
209 //
210 if (n0_depend + mu_depend + la_depend + ga_depend != 2)
211 throw runtime_error(
212 "Two (but only two) of n0, mu, la and ga must be NaN, "
213 "to flag that these parameters are the ones dependent of "
214 "mass content and mean particle size.");
215 if (mu_depend || ga_depend)
216 throw runtime_error(
217 "Sorry, mu and la are not yet allowed to be a "
218 "dependent parameter.");
219 //
220 const Index n0_fixed = (Index) !(n0_depend || std::isnan(n0));
221 const Index mu_fixed = (Index) !(mu_depend || std::isnan(mu));
222 const Index la_fixed = (Index) !(la_depend || std::isnan(la));
223 const Index ga_fixed = (Index) !(ga_depend || std::isnan(ga));
224 //
225 if (nin + n0_fixed + mu_fixed + la_fixed + ga_fixed != 4)
226 throw runtime_error(
227 "This PSD has four free parameters. This means that "
228 "the number\nof columns in *pnd_agenda_input* and the "
229 "number of numerics\n(i.e. not -999 or NaN) and among "
230 "the GIN arguments n0, mu, la and\nga must add up to "
231 "four. And this was found not to be the case.");
232
233 // Create vectors to hold the four MGD and the "extra" parameters
234 Vector mgd_pars(4), ext_pars(2);
235 ArrayOfIndex mgd_i_pai = {-1, -1, -1, -1}; // Position in pnd_agenda_input
236 const ArrayOfIndex ext_i_pai = {0, 1}; // Position in pnd_agenda_input
237 {
238 Index nhit = 2; // As mass and Dm always occupy first two positions
239 if (n0_fixed) {
240 mgd_pars[0] = n0;
241 } else if (!n0_depend) {
242 mgd_i_pai[0] = nhit++;
243 }
244 if (mu_fixed) {
245 mgd_pars[1] = mu;
246 } else if (!mu_depend) {
247 mgd_i_pai[1] = nhit++;
248 }
249 if (la_fixed) {
250 mgd_pars[2] = la;
251 } else if (!la_depend) {
252 mgd_i_pai[2] = nhit++;
253 }
254 if (ga_fixed) {
255 mgd_pars[3] = ga;
256 } else if (!ga_depend) {
257 mgd_i_pai[3] = nhit++;
258 }
259 }
260
261 // Determine what derivatives to do and their positions
262 ArrayOfIndex mgd_do_jac = {
263 0,
264 0,
265 0,
266 0,
267 };
268 ArrayOfIndex ext_do_jac = {0, 0};
269 ArrayOfIndex mgd_i_jac = {
270 -1, -1, -1, -1}; // Position among jacobian quantities
271 ArrayOfIndex ext_i_jac = {-1, -1}; // Position among jacobian quantities
272 //
273 for (Index i = 0; i < ndx; i++) {
274 if (dx2in[i] == 0) // That is, mass is a derivative
275 {
276 ext_do_jac[0] = 1;
277 ext_i_jac[0] = i;
278 } else if (dx2in[i] == 1) // That is, "something" is a derivative
279 {
280 ext_do_jac[1] = 1;
281 ext_i_jac[1] = i;
282 } else // Otherwise, either n0, mu, la or ga
283 {
284 for (Index j = 0; j < 4; j++) {
285 if (dx2in[i] == mgd_i_pai[j]) {
286 mgd_do_jac[j] = 1;
287 mgd_i_jac[j] = i;
288 break;
289 }
290 }
291 }
292 }
293
294 // Loop input data and calculate PSDs
295 for (Index ip = 0; ip < np; ip++) {
296
297 // Extract mass
298 ext_pars[0] = pnd_agenda_input(ip, ext_i_pai[0]);
299
300 // No calc needed if mass==0 and no jacobians requested.
301 if ((ext_pars[0] == 0.) && (!ndx)) {
302 continue;
303 } // If here, we are ready with this point!
304
305 // Extract "something"
306 ext_pars[1] = pnd_agenda_input(ip, ext_i_pai[1]);
307 if (ext_pars[1] <= 0) {
308 ostringstream os;
309 os << "Negative or zero " << something << " found in a position where "
310 << "this is not allowed.";
311 if (ndx)
312 os << "\nNote that for retrievals, " << something
313 << " must be set to a value > 0 even where\nmass is zero.";
314 throw std::runtime_error(os.str());
315 }
316
317 // Extract core MGD parameters
318 for (Index i = 0; i < 4; i++) {
319 if (mgd_i_pai[i] >= 0) {
320 mgd_pars[i] = pnd_agenda_input(ip, mgd_i_pai[i]);
321 }
322 }
323 Numeric t = pnd_agenda_input_t[ip];
324
325 // Outside of [t_min,tmax]?
326 if (t < t_min || t > t_max) {
327 if (picky) {
328 ostringstream os;
329 os << "Method called with a temperature of " << t << " K.\n"
330 << "This is outside the specified allowed range: [ max(0.," << t_min
331 << "), " << t_max << " ]";
332 throw runtime_error(os.str());
333 } else {
334 continue;
335 } // If here, we are ready with this point!
336 }
337
338 // Derive the dependent parameters (see ATD)
339 //
340 Numeric mu1 = 0, mub1 = 0, eterm = 0, gterm = 0;
341 Numeric scfac1 = 0, scfac2 = 0, gab = 0;
342 //
343 // *** Mean size ***
344 if (something == "mean size") {
345 if (n0_depend && la_depend) {
346 mub1 = mgd_pars[1] + scat_species_b + 1;
347 if (mub1 <= 0)
348 throw runtime_error("Bad MGD parameter detected: mu + b + 1 <= 0");
349 eterm = mub1 / mgd_pars[3];
350 // Start by deriving la
351 gterm = tgamma(eterm);
352 scfac2 = pow(tgamma(eterm+1/ga)/gterm, mgd_pars[3]);
353 mgd_pars[2] = scfac2 * pow(ext_pars[1], -mgd_pars[3]);
354 // We can now derive n0
355 scfac1 =
356 (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
357 mgd_pars[0] = scfac1 * ext_pars[0];
358 } else {
359 ARTS_ASSERT(0);
360 }
361 }
362
363 // *** Median size ***
364 else if (something == "median size") {
365 if (n0_depend && la_depend) {
366 mub1 = mgd_pars[1] + scat_species_b + 1;
367 if (mub1 <= 0)
368 throw runtime_error("Bad MGD parameter detected: mu + b + 1 <= 0");
369 eterm = mub1 / mgd_pars[3];
370 // Start by deriving la
371 scfac2 = (mgd_pars[1] + 1 + scat_species_b - 0.327 * mgd_pars[3]) /
372 mgd_pars[3];
373 mgd_pars[2] = scfac2 * pow(ext_pars[1], -mgd_pars[3]);
374 // We can now derive n0
375 gterm = tgamma(eterm);
376 scfac1 =
377 (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
378 mgd_pars[0] = scfac1 * ext_pars[0];
379 } else {
380 ARTS_ASSERT(0);
381 }
382 }
383
384 // *** Mean particle size ***
385 else if (something == "mean particle mass") {
386 if (n0_depend && la_depend) {
387 mu1 = mgd_pars[1] + 1;
388 if (mu1 <= 0)
389 throw runtime_error("Bad MGD parameter detected: mu + 1 <= 0");
390 eterm = (mgd_pars[1] + scat_species_b + 1) / mgd_pars[3];
391 gterm = tgamma(eterm);
392 // Start by deriving la
393 gab = mgd_pars[3] / scat_species_b;
394 scfac2 = pow(scat_species_a * gterm / tgamma(mu1 / mgd_pars[3]), gab);
395 mgd_pars[2] = scfac2 * pow(ext_pars[1], -gab);
396 // We can now derive n0
397 scfac1 =
398 (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
399 mgd_pars[0] = scfac1 * ext_pars[0];
400 } else {
401 ARTS_ASSERT(0);
402 }
403 }
404
405 else if (something == "Ntot") {
406 if (n0_depend && la_depend) {
407 mu1 = mgd_pars[1] + 1;
408 if (mu1 <= 0)
409 throw runtime_error("Bad MGD parameter detected: mu + 1 <= 0");
410 eterm = (mgd_pars[1] + scat_species_b + 1) / mgd_pars[3];
411 gterm = tgamma(eterm);
412 // Start by deriving la
413 gab = mgd_pars[3] / scat_species_b;
414 scfac2 = pow(scat_species_a * gterm / tgamma(mu1 / mgd_pars[3]), gab);
415 mgd_pars[2] = scfac2 * pow(ext_pars[1] / ext_pars[0], gab);
416 // We can now derive n0
417 scfac1 =
418 (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
419 mgd_pars[0] = scfac1 * ext_pars[0];
420 } else {
421 ARTS_ASSERT(0);
422 }
423 }
424
425 // String something not recognised
426 else {
427 ARTS_ASSERT(0);
428 }
429
430 // Now when all four MGD parameters are set, check that la and ga are OK
431 if (mgd_pars[2] <= 0)
432 throw runtime_error("Bad MGD parameter detected: la <= 0");
433 if (mgd_pars[3] <= 0)
434 throw runtime_error("Bad MGD parameter detected: ga <= 0");
435
436 // Calculate PSD and derivatives
437 Matrix jac_data(4, nsi);
438 mgd_with_derivatives(psd_data(ip, joker),
439 jac_data,
440 psd_size_grid,
441 mgd_pars[0],
442 mgd_pars[1],
443 mgd_pars[2],
444 mgd_pars[3],
445 (bool)mgd_do_jac[0] || n0_depend,
446 (bool)mgd_do_jac[1] || mu_depend,
447 (bool)mgd_do_jac[2] || la_depend,
448 (bool)mgd_do_jac[3] || ga_depend);
449
450 // Derivatives for mass and something
451 if (ext_do_jac[0] | ext_do_jac[1]) {
452 // *** Mean size ***
453 if (something == "mean size") {
454 if (n0_depend && la_depend) {
455 // Derivative with respect to mass
456 if (ext_do_jac[0]) {
457 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
458 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
459 }
460 // Derivative with respect to mean size
461 if (ext_do_jac[1]) {
462 // 1. Term associated with n0
463 // Calculated as dpsd/dn0 * dn0/dla * dla/dXm
464 dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
465 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
466 ext_pars[0] * mgd_pars[3] * eterm *
467 pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
468 // 2. Term associated with la
469 // Calculated as dpsd/dla * dla/dXm
470 dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
471 // Apply dla/dXm to sum
472 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
473 -mgd_pars[3] * scfac2 * pow(ext_pars[1], -(mgd_pars[3] + 1));
474 }
475 } else {
476 ARTS_ASSERT(0);
477 }
478 }
479
480 // *** Median size ***
481 else if (something == "median size") {
482 if (n0_depend && la_depend) {
483 // Derivative with respect to mass
484 if (ext_do_jac[0]) {
485 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
486 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
487 }
488 // Derivative with respect to median size
489 if (ext_do_jac[1]) {
490 // 1. Term associated with n0
491 // Calculated as dpsd/dn0 * dn0/dla * dla/dXm
492 dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
493 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
494 ext_pars[0] * mgd_pars[3] * eterm *
495 pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
496 // 2. Term associated with la
497 // Calculated as dpsd/dla * dla/dXm
498 dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
499 // Apply dla/dXm to sum
500 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
501 -mgd_pars[3] * scfac2 * pow(ext_pars[1], -(mgd_pars[3] + 1));
502 }
503 } else {
504 ARTS_ASSERT(0);
505 }
506 }
507
508 // *** Mean particle size ***
509 else if (something == "mean particle mass") {
510 if (n0_depend && la_depend) {
511 // Derivative with respect to mass
512 if (ext_do_jac[0]) {
513 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
514 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
515 }
516 // Derivative with respect to mean particle size
517 if (ext_do_jac[1]) {
518 // 1. Term associated with n0
519 // Calculated as dpsd/dn0 * dn0/dla * dla/dMm
520 dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
521 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
522 ext_pars[0] * mgd_pars[3] * eterm *
523 pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
524 // 2. Term associated with la
525 // Calculated as dpsd/dla * dla/dMm
526 dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
527 // Apply dla/dMm to sum
528 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
529 scfac2 * (-mgd_pars[3] / scat_species_b) *
530 pow(ext_pars[1], -(gab + 1));
531 } else {
532 ARTS_ASSERT(0);
533 }
534 }
535 }
536
537 else if (something == "Ntot") {
538 if (n0_depend && la_depend) {
539 // Term part of both derivatives
540 const Numeric dn0dla = ext_pars[0] * mgd_pars[3] * eterm *
541 pow(mgd_pars[2], eterm - 1) /
542 (scat_species_a * gterm);
543 // Derivative with respect to mass
544 if (ext_do_jac[0]) {
545 // Repeated term
546 const Numeric dladw = scfac2 * pow(ext_pars[1], gab) *
547 (-mgd_pars[3] / scat_species_b) *
548 pow(ext_pars[0], -(gab + 1));
549 // 1. Term associated with n0
550 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
551 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1 + dn0dla * dladw;
552 // 2. Term associated with la
553 Vector term2{jac_data(2, joker)};
554 term2 *= dladw;
555 // Sum up
556 dpsd_data_dx(ext_i_jac[0], ip, joker) += term2;
557 }
558 // Derivative with respect to Ntot
559 if (ext_do_jac[1]) {
560 // 1. Term associated with n0
561 // Calculated as dpsd/dn0 * dn0/dla * dla/dNtot
562 dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
563 dpsd_data_dx(ext_i_jac[1], ip, joker) *= dn0dla;
564 // 2. Term associated with la
565 // Calculated as dpsd/dla * dla/dNtot
566 dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
567 // Apply dla/dNtot to sum
568 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
569 scfac2 * pow(ext_pars[0], -gab) *
570 (mgd_pars[3] / scat_species_b) * pow(ext_pars[1], gab - 1);
571 }
572 } else {
573 ARTS_ASSERT(0);
574 }
575 }
576
577 // String something not recognised
578 else {
579 ARTS_ASSERT(0);
580 }
581 }
582
583 // Derivatives for non-dependent native parameters
584 for (Index i = 0; i < 4; i++) {
585 if (mgd_do_jac[i]) {
586 dpsd_data_dx(mgd_i_jac[i], ip, joker) = jac_data(i, joker);
587 }
588 }
589 }
590}
591
592void psd_mono_common(Matrix& psd_data,
593 Tensor3& dpsd_data_dx,
594 const String& type,
595 const Vector& pnd_agenda_input_t,
596 const Matrix& pnd_agenda_input,
597 const ArrayOfString& pnd_agenda_input_names,
598 const ArrayOfString& dpnd_data_dx_names,
599 const ArrayOfArrayOfScatteringMetaData& scat_meta,
600 const Index& species_index,
601 const Numeric& t_min,
602 const Numeric& t_max,
603 const Index& picky,
604 const Verbosity&) {
605 // Standard checcks
606 const Vector psd_size_grid(1, 0); // As this WSV is not input for thse WSM
608
609 // Extra checks for this PSD
610 const Index nss = scat_meta.nelem();
611 if (nss == 0) throw runtime_error("*scat_meta* is empty!");
612 if (nss < species_index + 1) {
613 ostringstream os;
614 os << "Selected scattering species index is " << species_index
615 << " but this "
616 << "is not allowed since *scat_meta* has only " << nss << " elements.";
617 throw runtime_error(os.str());
618 }
619 if (scat_meta[species_index].nelem() != 1) {
620 ostringstream os;
621 os << "This method only works with scattering species consisting of a\n"
622 << "single element, but your data do not match this demand.\n"
623 << "Selected scattering species index is " << species_index << ".\n"
624 << "This species has " << scat_meta[species_index].nelem()
625 << " elements.";
626 throw runtime_error(os.str());
627 }
628 //
629 if (pnd_agenda_input.ncols() != 1)
630 throw runtime_error("*pnd_agenda_input* must have one column.");
631 if (nsi != 1)
632 throw runtime_error(
633 "This method demands that length of *psd_size_grid* is 1.");
634
635 // Extract particle mass
636 Numeric pmass = 0;
637 if (type == "mass") {
638 pmass = scat_meta[species_index][0].mass;
639 }
640
641 for (Index ip = 0; ip < np; ip++) {
642 // Extract the input variables
643 Numeric x = pnd_agenda_input(ip, 0);
644 Numeric t = pnd_agenda_input_t[ip];
645
646 // No calc needed if n==0 and no jacobians requested.
647 if ((x == 0.) && (!ndx)) {
648 continue;
649 } // If here, we are ready with this point!
650
651 // Outside of [t_min,tmax]?
652 if (t < t_min || t > t_max) {
653 if (picky) {
654 ostringstream os;
655 os << "Method called with a temperature of " << t << " K.\n"
656 << "This is outside the specified allowed range: [ max(0.," << t_min
657 << "), " << t_max << " ]";
658 throw runtime_error(os.str());
659 } else {
660 continue;
661 } // If here, we are ready with this point!
662 }
663
664 // Set PSD
665 //
666 if (type == "ntot") {
667 psd_data(ip, 0) = x;
668 //
669 if (ndx) {
670 dpsd_data_dx(0, ip, 0) = 1;
671 }
672 } else if (type == "mass") {
673 psd_data(ip, 0) = x / pmass;
674 //
675 if (ndx) {
676 dpsd_data_dx(0, ip, 0) = 1 / pmass;
677 }
678 } else {
679 ARTS_ASSERT(0);
680 }
681 }
682}
683
684void psd_rain_W16(Vector& psd, const Vector& diameter, const Numeric& rwc) {
685 Index nD = diameter.nelem();
686 psd.resize(nD);
687 psd = 0.;
688
689 // skip calculation if RWC is 0.0
690 if (rwc == 0.0) {
691 return;
692 }
693 ARTS_ASSERT(rwc > 0.);
694
695 // a and b relates N0 to lambda N0 = a*lambda^b
696 Numeric a = 0.000141;
697 Numeric b = 1.49;
698 Numeric c1 = DENSITY_OF_WATER * PI / 6;
699 Numeric base = c1 / rwc * a * tgamma(4);
700 Numeric exponent = 1. / (4 - b);
701 Numeric lambda = pow(base, exponent);
702 Numeric N0 = a * pow(lambda, b);
703
704 //psd_general_MGD( psd, N0, 0, lambda, 1 );
705 N0 *= 1e8;
706 lambda *= 100;
707 for (Index iD = 0; iD < nD; iD++) {
708 psd[iD] = N0 * exp(-lambda * diameter[iD]);
709 }
710}
711
712void psd_mgd_smm_common(Matrix& psd_data,
713 Tensor3& dpsd_data_dx,
714 const String& psd_name,
715 const Vector& psd_size_grid,
716 const Vector& pnd_agenda_input_t,
717 const Matrix& pnd_agenda_input,
718 const ArrayOfString& pnd_agenda_input_names,
719 const ArrayOfString& dpnd_data_dx_names,
720 const Numeric& scat_species_a,
721 const Numeric& scat_species_b,
722 const Numeric& n_alpha_in,
723 const Numeric& n_b_in,
724 const Numeric& mu_in,
725 const Numeric& gamma_in,
726 const Numeric& t_min,
727 const Numeric& t_max,
728 const Index& picky,
729 const Verbosity&) {
730 // Standard checks
732
733 // All PSDs are defined in terms of hydrometeor mass content
734 if (pnd_agenda_input.ncols() != 1)
735 throw runtime_error("*pnd_agenda_input* must have one column.");
736
737 // Extra checks for rain PSDs which should be consistent with
738 // spherical liquid drops
739 if (psd_name == "Abel12" || psd_name == "Wang16"){
740 if (scat_species_b < 2.9 || scat_species_b > 3.1) {
741 ostringstream os;
742 os << "This PSD treats rain, using Dveq as size grid.\n"
743 << "This means that *scat_species_b* should be close to 3,\n"
744 << "but it is outside of the tolerated range of [2.9,3.1].\n"
745 << "Your value of *scat_species_b* is: " << scat_species_b;
746 throw runtime_error(os.str());
747 }
748 if (scat_species_a < 500 || scat_species_a > 540) {
749 ostringstream os;
750 os << "This PSD treats rain, using Dveq as size grid.\n"
751 << "This means that *scat_species_a* should be close to 520,\n"
752 << "but it is outside of the tolerated range of [500,540].\n"
753 << "Your value of *scat_species_a* is: " << scat_species_a;
754 throw runtime_error(os.str());
755 }
756 }
757 // Extra checks for graupel/hail PSDs which assume constant effective density
758 //
759 if (psd_name == "Field19"){
760 if (scat_species_b < 2.8 || scat_species_b > 3.2) {
761 ostringstream os;
762 os << "This PSD treats graupel, assuming a constant effective density.\n"
763 << "This means that *scat_species_b* should be close to 3,\n"
764 << "but it is outside of the tolerated range of [2.8,3.2].\n"
765 << "Your value of *scat_species_b* is: " << scat_species_b;
766 throw runtime_error(os.str());
767 }
768 }
769
770 for (Index ip = 0; ip < np; ip++) {
771 // Extract the input variables
772 Numeric water_content = pnd_agenda_input(ip, 0);
773 Numeric t = pnd_agenda_input_t[ip];
774
775 // No calc needed if water_content==0 and no jacobians requested.
776 if ((water_content == 0.) && (!ndx)) {
777 continue;
778 } // If here, we are ready with this point!
779
780 // Outside of [t_min,tmax]?
781 if (t < t_min || t > t_max) {
782 if (picky) {
783 ostringstream os;
784 os << "Method called with a temperature of " << t << " K.\n"
785 << "This is outside the specified allowed range: [ max(0.," << t_min
786 << "), " << t_max << " ]";
787 throw runtime_error(os.str());
788 } else {
789 continue;
790 } // If here, we are ready with this point!
791 }
792
793 // Negative wc?
794 Numeric psd_weight = 1.0;
795 if (water_content < 0) {
796 psd_weight = -1.0;
797 water_content *= -1.0;
798 }
799
800 // PSD settings for different parametrizations
801 Numeric gamma = 0.0;
802 Numeric n_alpha = 0.0;
803 Numeric n_b = 0.0;
804 Numeric mu = 0.0;
805 if (psd_name == "Abel12"){
806 n_alpha = 0.22;
807 n_b = 2.2;
808 mu = 0.0;
809 gamma = 1.0;
810 }
811 else if (psd_name == "Wang16"){
812 // Wang 16 parameters converted to SI units
813 n_alpha = 14.764;
814 n_b = 1.49;
815 mu = 0.0;
816 gamma = 1.0;
817 }
818 else if (psd_name == "Field19"){
819 n_alpha = 7.9e9;
820 n_b = -2.58;
821 mu = 0.0;
822 gamma = 1.0;
823 }
824 else if (psd_name == "generic"){
825 n_alpha = n_alpha_in;
826 n_b = n_b_in;
827 mu = mu_in;
828 gamma = gamma_in;
829 }
830 else {
831 ARTS_ASSERT(0);
832 }
833
834 // Calculate PSD
835 // Calculate lambda for modified gamma distribution from mass density
836 Numeric k = (scat_species_b + mu + 1 - gamma)/gamma;
837 Numeric expo = 1.0 / (n_b - k - 1);
838 Numeric denom = scat_species_a * n_alpha * tgamma(k + 1);
839 Numeric lam = pow(water_content*gamma/denom, expo);
840 Numeric n_0 = n_alpha * pow(lam, n_b);
841 Vector psd_1p(nsi);
842 Matrix jac_data(4, nsi);
843
844 psd_1p = 0.0;
845 if (water_content != 0) {
846 mgd_with_derivatives(psd_1p, jac_data, psd_size_grid, n_0, mu, lam, gamma,
847 true, // n_0 jacobian
848 false,// mu jacobian
849 true, // lambda jacobian
850 false); // gamma jacobian
851 } else {
852 ARTS_ASSERT(0);
853 }
854 //
855 for (Index i = 0; i < nsi; i++) {
856 psd_data(ip, i) = psd_weight * psd_1p[i];
857 }
858
859 // Calculate derivative with respect to water content
860 if (ndx) {
861 const Numeric dlam_dwc = pow(gamma/denom, expo) * expo * pow(water_content, expo-1);
862 const Numeric dn_0_dwc = n_alpha * n_b * pow(lam, n_b-1) * dlam_dwc;
863 for (Index i = 0; i < nsi; i++) {
864 dpsd_data_dx(0, ip, i) = psd_weight * (jac_data(0,i)*dn_0_dwc +
865 jac_data(2,i)*dlam_dwc);
866 }
867 }
868 }
869}
870
871void psd_snow_F07(Vector& psd,
872 const Vector& diameter,
873 const Numeric& swc,
874 const Numeric& t,
875 const Numeric alpha,
876 const Numeric beta,
877 const String& regime) {
878 Index nD = diameter.nelem();
879 psd.resize(nD);
880 psd = 0.;
881
882 ARTS_ASSERT(t > 0.);
883
884 // skip calculation if SWC is 0.0
885 if (swc == 0.0) {
886 return;
887 }
888 ARTS_ASSERT(swc > 0.);
889
890 Numeric An, Bn, Cn;
891 Numeric M2, Mn, M2Mn;
892 Numeric base, pp;
893 Numeric x, phi23;
894 //Numeric dN;
895
896 Vector q(5);
897
898 ARTS_ASSERT((regime == "TR") || (regime == "ML"));
899
900 //factors of phi23
901 if (regime == "TR")
902 q = {152., -12.4, 3.28, -0.78, -1.94};
903 else // if( regime=="ML" )
904 q = {141., -16.8, 102., 2.07, -4.82};
905
906 //Factors of factors of the moment estimation parametrization
907 Vector Aq{13.6, -7.76, 0.479};
908 Vector Bq{-0.0361, 0.0151, 0.00149};
909 Vector Cq{0.807, 0.00581, 0.0457};
910
911 //convert T from Kelvin to Celsius
912 Numeric Tc = t - 273.15;
913
914 // estimate second moment
915 M2 = swc / alpha;
916 if (beta != 2) {
917 // calculate factors of the moment estimation parametrization
918 An = exp(Aq[0] + Aq[1] * beta + Aq[2] * pow(beta, 2));
919 Bn = Bq[0] + Bq[1] * beta + Bq[2] * pow(beta, 2);
920 Cn = Cq[0] + Cq[1] * beta + Cq[2] * pow(beta, 2);
921
922 base = M2 * exp(-Bn * Tc) / An;
923 pp = 1. / (Cn);
924
925 M2 = pow(base, pp);
926 }
927
928 // order of the moment parametrization
929 const Numeric n = 3;
930
931 // calculate factors of the moment estimation parametrization
932 An = exp(Aq[0] + Aq[1] * n + Aq[2] * pow(n, 2));
933 Bn = Bq[0] + Bq[1] * n + Bq[2] * pow(n, 2);
934 Cn = Cq[0] + Cq[1] * n + Cq[2] * pow(n, 2);
935
936 // moment parametrization
937 Mn = An * exp(Bn * Tc) * pow(M2, Cn);
938
939 M2Mn = pow(M2, 4.) / pow(Mn, 3.);
940
941 for (Index iD = 0; iD < nD; iD++) {
942 // define x
943 x = diameter[iD] * M2 / Mn;
944
945 // characteristic function
946 phi23 = q[0] * exp(q[1] * x) + q[2] * pow(x, q[3]) * exp(q[4] * x);
947
948 // distribution function
949 //dN = phi23*M2Mn;
950
951 //if ( !std::isnan(psd[dN]) )
952 // psd[iD] = dN;
953
954 // set psd directly. Non-NaN should be (and is, hopefully) ensured by
955 // checks above (if we encounter further NaN, that should be handled above.
956 // which intermediate quantities make problems? at which parametrisation
957 // values, ie WC, T, alpha, beta?).
958 psd[iD] = phi23 * M2Mn;
959 }
960}
961
962void psd_SB06(Vector& psd,
963 Matrix& dpsd,
964 const Vector& mass,
965 const Numeric& N_tot,
966 const Numeric& WC,
967 const String& hydrometeor_type) {
968 Numeric N0;
969 Numeric Lambda;
970 Numeric arg1;
971 Numeric arg2;
972 Numeric brk;
973 Numeric mu;
974 Numeric gamma;
975 Numeric xmin;
976 Numeric xmax;
977 Numeric M0min;
978 Numeric M0max;
979 Numeric M0;
980 Numeric M1;
981 Numeric c1;
982 Numeric c2;
983 Numeric L1;
984 Numeric mMu;
985 Numeric mGamma;
986 Numeric brkMu1;
987
988 // Get the coefficients for the right hydrometeor
989 if (hydrometeor_type == "cloud_ice") //Cloud ice water
990 {
991 mu = 0.;
992 gamma = 1. / 3.;
993 xmin = 1e-12;
994 xmax = 1e-5;
995 } else if (hydrometeor_type == "rain") //Rain
996 {
997 mu = 0.;
998 gamma = 1. / 3.;
999 xmin = 2.6e-10;
1000 xmax = 3e-6;
1001 } else if (hydrometeor_type == "snow") //Snow
1002 {
1003 mu = 0.;
1004 gamma = 1. / 2.;
1005 xmin = 1e-10;
1006 xmax = 2e-5;
1007 } else if (hydrometeor_type == "graupel") //Graupel
1008 {
1009 mu = 1.;
1010 gamma = 1. / 3.;
1011 xmin = 1e-9;
1012 xmax = 5e-4;
1013 } else if (hydrometeor_type == "hail") //Hail
1014 {
1015 mu = 1.;
1016 gamma = 1. / 3.;
1017 xmin = 2.6e-10;
1018 xmax = 5e-4;
1019 } else if (hydrometeor_type == "cloud_water") //Cloud liquid water
1020 {
1021 mu = 1;
1022 gamma = 1;
1023 xmin = 4.2e-15;
1024 xmax = 2.6e-10;
1025 } else {
1026 ostringstream os;
1027 os << "You use a wrong tag! ";
1028 throw runtime_error(os.str());
1029 }
1030
1031 M0 = N_tot;
1032 M1 = WC;
1033
1034 Index nD = mass.nelem();
1035 psd.resize(nD);
1036 psd = 0.;
1037
1038 dpsd.resize(nD, 2);
1039 dpsd = 0.;
1040
1041 if (M1 > 0.0) {
1042 // lower and upper limit check is taken from the ICON code of the two moment
1043 //scheme
1044
1045 M0max = M1 / xmax;
1046 M0min = M1 / xmin;
1047
1048 //check lower limit of the scheme
1049 if (M0 > M0min) {
1050 M0 = M0min;
1051 }
1052
1053 //check upper limit of the scheme
1054 if (M0 < M0max) {
1055 M0 = M0max;
1056 }
1057
1058 //arguments for Gamma function
1059 arg2 = (mu + 2) / gamma;
1060 arg1 = (mu + 1) / gamma;
1061
1062 // results of gamma function
1063 c1 = tgamma(arg1);
1064 c2 = tgamma(arg2);
1065
1066 // variable to shorten the formula
1067 brk = M0 / M1 * c2 / c1;
1068 brkMu1 = pow(brk, (mu + 1));
1069
1070 //Lambda (parameter for modified gamma distribution)
1071 Lambda = pow(brk, gamma);
1072
1073 L1 = pow(Lambda, arg1);
1074
1075 //N0
1076 N0 = M0 * gamma / tgamma(arg1) * L1;
1077
1078 // Calculate distribution function
1079 for (Index iD = 0; iD < nD; iD++) {
1080 //Distribution function
1081 psd[iD] = mod_gamma_dist(mass[iD], N0, Lambda, mu, gamma);
1082
1083 if (std::isnan(psd[iD])) psd[iD] = 0.0;
1084 if (std::isinf(psd[iD])) psd[iD] = 0.0;
1085
1086 //Calculate derivatives analytically
1087 mMu = pow(mass[iD], mu);
1088 mGamma = pow(mass[iD], gamma);
1089
1090 // dpsd/dM1
1091 dpsd(iD, 0) = gamma / c1 * M0 / M1 * mMu * exp(-Lambda * mGamma) *
1092 brkMu1 * (-1 - mu + gamma * mGamma * Lambda);
1093
1094 // dpsd/dM0
1095 dpsd(iD, 1) = -gamma / c1 * mMu * exp(-Lambda * mGamma) * brkMu1 *
1096 (-2 - mu - gamma * mGamma * Lambda);
1097 }
1098 } else {
1099 return;
1100 }
1101}
1102
1103void psd_MY05(Vector& psd,
1104 Matrix& dpsd,
1105 const Vector& diameter_max,
1106 const Numeric N_tot,
1107 const Numeric WC,
1108 const String psd_type) {
1109 Numeric N0;
1110 Numeric Lambda;
1111 Numeric arg1;
1112 Numeric arg2;
1113 Numeric temp;
1114 Numeric mu;
1115 Numeric gamma;
1116 Numeric alpha;
1117 Numeric beta;
1118 Numeric M0;
1119 Numeric M1;
1120 Numeric c1;
1121 Numeric c2;
1122 Numeric Lmg;
1123 Numeric DMu;
1124 Numeric DGamma;
1125
1126 // Get the coefficients for the right hydrometeor
1127 if (psd_type == "cloud_ice") //Cloud ice water
1128 {
1129 mu = 0.;
1130 gamma = 1.;
1131 alpha = 440.; //[kg]
1132 beta = 3;
1133 } else if (psd_type == "rain") //Rain
1134 {
1135 mu = 0.;
1136 gamma = 1;
1137 alpha = 523.5988; //[kg]
1138 beta = 3;
1139 } else if (psd_type == "snow") //Snow
1140 {
1141 mu = 0.;
1142 gamma = 1;
1143 alpha = 52.35988; //[kg]
1144 beta = 3;
1145 } else if (psd_type == "graupel") //Graupel
1146 {
1147 mu = 0.;
1148 gamma = 1;
1149 alpha = 209.4395; //[kg]
1150 beta = 3;
1151 } else if (psd_type == "hail") //Hail
1152 {
1153 mu = 0.;
1154 gamma = 1;
1155 alpha = 471.2389; //[kg]
1156 beta = 3;
1157 } else if (psd_type == "cloud_water") //Cloud liquid water
1158 {
1159 mu = 1;
1160 gamma = 1;
1161 alpha = 523.5988; //[kg]
1162 beta = 3;
1163 } else {
1164 ostringstream os;
1165 os << "You use a wrong tag! ";
1166 throw runtime_error(os.str());
1167 }
1168
1169 M0 = N_tot;
1170 M1 = WC;
1171
1172 Index nD = diameter_max.nelem();
1173 psd.resize(nD);
1174 psd = 0.;
1175
1176 dpsd.resize(nD, 2);
1177 dpsd = 0.;
1178
1179 if (M1 > 0.0 && M0 > 0) {
1180 //arguments for Gamma function
1181 arg2 = (mu + beta + 1) / gamma;
1182 arg1 = (mu + 1) / gamma;
1183
1184 // results of gamma function
1185 c1 = tgamma(arg1);
1186 c2 = tgamma(arg2);
1187
1188 //base of lambda
1189 temp = alpha * M0 / M1 * c2 / c1;
1190
1191 //Lambda (parameter for modified gamma distribution)
1192 Lambda = pow(temp, gamma / beta);
1193
1194 Lmg = pow(Lambda, arg1);
1195
1196 //N0
1197 N0 = M0 * gamma / c1 * Lmg;
1198
1199 //Distribution function
1200
1201 // Calculate distribution function
1202 for (Index iD = 0; iD < nD; iD++) {
1203 psd[iD] = mod_gamma_dist(diameter_max[iD], N0, Lambda, mu, gamma);
1204
1205 if (std::isnan(psd[iD])) psd[iD] = 0.0;
1206 if (std::isinf(psd[iD])) psd[iD] = 0.0;
1207
1208 //Calculate derivatives analytically
1209 DMu = pow(diameter_max[iD], mu);
1210 DGamma = pow(diameter_max[iD], gamma);
1211
1212 // dpsd/dM1
1213 dpsd(iD, 0) = (DMu * exp(-DGamma * Lambda) * gamma * M0 * Lmg *
1214 (-1 - mu + DGamma * gamma * Lambda) / (M1 * beta * c1));
1215
1216 // dpsd/dM0
1217 dpsd(iD, 1) = (DMu * exp(-DGamma * Lambda) * gamma * Lmg *
1218 (1 + beta + mu - DGamma * gamma * Lambda) / (beta * c1));
1219 }
1220
1221 } else {
1222 return;
1223 }
1224}
1225
1226Numeric dm_from_iwc_n0(Numeric iwc, Numeric n0, Numeric rho) {
1227 if (iwc == 0.0) {
1228 return 1e-9;
1229 } else {
1230 return pow(256.0 * iwc / PI / rho / n0, 0.25);
1231 }
1232}
1233
1234Numeric n0_from_iwc_dm(Numeric iwc, Numeric dm, Numeric rho) {
1235 if (dm > 1e-9) {
1236 return 256.0 * iwc / PI / rho / pow(dm, 4.0);
1237 } else {
1238 return 0.0;
1239 }
1240}
1241
1242Numeric n0_from_t(Numeric t) { return exp(-0.076586 * (t - 273.15) + 17.948); }
base min(const Array< base > &x)
Min function.
Definition array.h:144
The global header file for ARTS.
Constants of physical expressions as constexpr.
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
A C++ standards dependent random number generator class.
Definition rng.h:22
auto get(Ts &&...x) const
Returns a random number generator of some random distribution.
Definition rng.h:103
Internal cloudbox functions.
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
Numeric mod_gamma_dist(Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma)
Generalized Modified Gamma Distribution.
void mgd_with_derivatives(VectorView psd, MatrixView jac_data, const Vector &x, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga, const bool &do_n0_jac, const bool &do_mu_jac, const bool &do_la_jac, const bool &do_ga_jac)
Workspace functions for the solution of cloud-box radiative transfer by Monte Carlo methods....
Declarations having to do with the four output streams.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric density_of_ice_at_0c
Global constant, Density of water ice at 0C [kg/m3] source: http://en.wikipedia.org/wiki/Ice.
constexpr Numeric denity_of_water_at_4c
Global constant, Density of liquid water +4C [kg/m3] source: http://en.wikipedia.org/wiki/Water.
This file contains declerations of functions of physical character.
Propagation path structure and functions.
void psd_rain_W16(Vector &psd, const Vector &diameter, const Numeric &rwc)
The Wang16 rain DSD DEPRECATED BY NEW MGD_SMM_COMMON Only included for compatibility with "old" pnd_f...
Definition psd.cc:684
void psd_snow_F07(Vector &psd, const Vector &diameter, const Numeric &swc, const Numeric &t, const Numeric alpha, const Numeric beta, const String &regime)
The F07 snow PSD.
Definition psd.cc:871
void psd_MY05(Vector &psd, Matrix &dpsd, const Vector &diameter_max, const Numeric N_tot, const Numeric WC, const String psd_type)
Definition psd.cc:1103
void psd_mgd_mass_and_something(Matrix &psd_data, Tensor3 &dpsd_data_dx, const String &something, const Vector &psd_size_grid, const Vector &pnd_agenda_input_t, const Matrix &pnd_agenda_input, const ArrayOfString &pnd_agenda_input_names, const ArrayOfString &dpnd_data_dx_names, const Numeric &scat_species_a, const Numeric &scat_species_b, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga, const Numeric &t_min, const Numeric &t_max, const Index &picky, const Verbosity &)
Code common to MGD PSD involving the integrated mass.
Definition psd.cc:173
Numeric n0_from_t(Numeric t)
Sets N0star based on temperature.
Definition psd.cc:1242
void psd_cloudice_MH97(Vector &psd, const Vector &diameter, const Numeric &iwc, const Numeric &t, const bool noisy)
The MH97 cloud ice PSD.
Definition psd.cc:39
void psd_mono_common(Matrix &psd_data, Tensor3 &dpsd_data_dx, const String &type, const Vector &pnd_agenda_input_t, const Matrix &pnd_agenda_input, const ArrayOfString &pnd_agenda_input_names, const ArrayOfString &dpnd_data_dx_names, const ArrayOfArrayOfScatteringMetaData &scat_meta, const Index &species_index, const Numeric &t_min, const Numeric &t_max, const Index &picky, const Verbosity &)
Code common to PSDs of mono type.
Definition psd.cc:592
constexpr Numeric DENSITY_OF_ICE
Definition psd.cc:36
Numeric n0_from_iwc_dm(Numeric iwc, Numeric dm, Numeric rho)
Derives N0star from IWC and Dm.
Definition psd.cc:1234
Numeric dm_from_iwc_n0(Numeric iwc, Numeric n0, Numeric rho)
Derives Dm from IWC and N0star.
Definition psd.cc:1226
constexpr Numeric PI
Definition psd.cc:35
void psd_SB06(Vector &psd, Matrix &dpsd, const Vector &mass, const Numeric &N_tot, const Numeric &WC, const String &hydrometeor_type)
Definition psd.cc:962
constexpr Numeric DENSITY_OF_WATER
Definition psd.cc:37
void psd_mgd_smm_common(Matrix &psd_data, Tensor3 &dpsd_data_dx, const String &psd_name, const Vector &psd_size_grid, const Vector &pnd_agenda_input_t, const Matrix &pnd_agenda_input, const ArrayOfString &pnd_agenda_input_names, const ArrayOfString &dpnd_data_dx_names, const Numeric &scat_species_a, const Numeric &scat_species_b, const Numeric &n_alpha_in, const Numeric &n_b_in, const Numeric &mu_in, const Numeric &gamma_in, const Numeric &t_min, const Numeric &t_max, const Index &picky, const Verbosity &)
Code common to a number of modified gamma PSDs used with single-moment mass schemes.
Definition psd.cc:712
Internal functions associated with size distributions.
#define START_OF_PSD_METHODS()
Definition psd.h:24
Contains sorting routines.
#define a
#define b