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