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