ARTS  2.4.0(git:4fb77825)
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 
54 extern const Numeric PI;
55 extern const Numeric DENSITY_OF_ICE;
56 extern 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  assert(t > 0.);
68 
69  // skip calculation if IWC is 0.0
70  if (iwc == 0.0) {
71  return;
72  }
73  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 
197  const String& something,
198  const Vector& psd_size_grid,
199  const Vector& pnd_agenda_input_t,
200  const Matrix& pnd_agenda_input,
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  // Extract mass
319  ext_pars[0] = pnd_agenda_input(ip, ext_i_pai[0]);
320  ext_pars[1] = pnd_agenda_input(ip, ext_i_pai[1]);
321  if (ext_pars[1] <= 0) {
322  ostringstream os;
323  os << "Negative " << something << "found.\nThis is not allowed.";
324  throw std::runtime_error(os.str());
325  }
326  // Extract core MGD parameters
327  for (Index i = 0; i < 4; i++) {
328  if (mgd_i_pai[i] >= 0) {
329  mgd_pars[i] = pnd_agenda_input(ip, mgd_i_pai[i]);
330  }
331  }
332  Numeric t = pnd_agenda_input_t[ip];
333 
334  // No calc needed if mass==0 and no jacobians requested.
335  if ((ext_pars[0] == 0.) && (!ndx)) {
336  continue;
337  } // If here, we are ready with this point!
338 
339  // Outside of [t_min,tmax]?
340  if (t < t_min || t > t_max) {
341  if (picky) {
342  ostringstream os;
343  os << "Method called with a temperature of " << t << " K.\n"
344  << "This is outside the specified allowed range: [ max(0.," << t_min
345  << "), " << t_max << " ]";
346  throw runtime_error(os.str());
347  } else {
348  continue;
349  } // If here, we are ready with this point!
350  }
351 
352  // Derive the dependent parameters (see ATD)
353  //
354  Numeric mu1 = 0, mub1 = 0, eterm = 0, gterm = 0;
355  Numeric scfac1 = 0, scfac2 = 0, gab = 0;
356  //
357  // *** Mean size ***
358  if (something == "mean size") {
359  if (n0_depend && la_depend) {
360  mub1 = mgd_pars[1] + scat_species_b + 1;
361  if (mub1 <= 0)
362  throw runtime_error("Bad MGD parameter detected: mu + b + 1 <= 0");
363  eterm = mub1 / mgd_pars[3];
364  // Start by deriving la
365  scfac2 = pow(eterm, mgd_pars[3]);
366  mgd_pars[2] = scfac2 * pow(ext_pars[1], -mgd_pars[3]);
367  // We can now derive n0
368  gterm = tgamma(eterm);
369  scfac1 =
370  (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
371  mgd_pars[0] = scfac1 * ext_pars[0];
372  } else {
373  assert(0);
374  }
375  }
376 
377  // *** Median size ***
378  else if (something == "median size") {
379  if (n0_depend && la_depend) {
380  mub1 = mgd_pars[1] + scat_species_b + 1;
381  if (mub1 <= 0)
382  throw runtime_error("Bad MGD parameter detected: mu + b + 1 <= 0");
383  eterm = mub1 / mgd_pars[3];
384  // Start by deriving la
385  scfac2 = (mgd_pars[1] + 1 + scat_species_b - 0.327 * mgd_pars[3]) /
386  mgd_pars[3];
387  mgd_pars[2] = scfac2 * pow(ext_pars[1], -mgd_pars[3]);
388  // We can now derive n0
389  gterm = tgamma(eterm);
390  scfac1 =
391  (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
392  mgd_pars[0] = scfac1 * ext_pars[0];
393  } else {
394  assert(0);
395  }
396  }
397 
398  // *** Mean particle size ***
399  else if (something == "mean particle mass") {
400  if (n0_depend && la_depend) {
401  mu1 = mgd_pars[1] + 1;
402  if (mu1 <= 0)
403  throw runtime_error("Bad MGD parameter detected: mu + 1 <= 0");
404  eterm = (mgd_pars[1] + scat_species_b + 1) / mgd_pars[3];
405  gterm = tgamma(eterm);
406  // Start by deriving la
407  gab = mgd_pars[3] / scat_species_b;
408  scfac2 = pow(scat_species_a * gterm / tgamma(mu1 / mgd_pars[3]), gab);
409  mgd_pars[2] = scfac2 * pow(ext_pars[1], -gab);
410  // We can now derive n0
411  scfac1 =
412  (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
413  mgd_pars[0] = scfac1 * ext_pars[0];
414  } else {
415  assert(0);
416  }
417  }
418 
419  else if (something == "Ntot") {
420  if (n0_depend && la_depend) {
421  mu1 = mgd_pars[1] + 1;
422  if (mu1 <= 0)
423  throw runtime_error("Bad MGD parameter detected: mu + 1 <= 0");
424  eterm = (mgd_pars[1] + scat_species_b + 1) / mgd_pars[3];
425  gterm = tgamma(eterm);
426  // Start by deriving la
427  gab = mgd_pars[3] / scat_species_b;
428  scfac2 = pow(scat_species_a * gterm / tgamma(mu1 / mgd_pars[3]), gab);
429  mgd_pars[2] = scfac2 * pow(ext_pars[1] / ext_pars[0], gab);
430  // We can now derive n0
431  scfac1 =
432  (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
433  mgd_pars[0] = scfac1 * ext_pars[0];
434  } else {
435  assert(0);
436  }
437  }
438 
439  // String something not recognised
440  else {
441  assert(0);
442  }
443 
444  // Now when all four MGD parameters are set, check that la and ga are OK
445  if (mgd_pars[2] <= 0)
446  throw runtime_error("Bad MGD parameter detected: la <= 0");
447  if (mgd_pars[3] <= 0)
448  throw runtime_error("Bad MGD parameter detected: ga <= 0");
449 
450  // Calculate PSD and derivatives
451  Matrix jac_data(4, nsi);
453  jac_data,
455  mgd_pars[0],
456  mgd_pars[1],
457  mgd_pars[2],
458  mgd_pars[3],
459  (bool)mgd_do_jac[0] || n0_depend,
460  (bool)mgd_do_jac[1] || mu_depend,
461  (bool)mgd_do_jac[2] || la_depend,
462  (bool)mgd_do_jac[3] || ga_depend);
463 
464  // Derivatives for mass and something
465  if (ext_do_jac[0] | ext_do_jac[1]) {
466  // *** Mean size ***
467  if (something == "mean size") {
468  if (n0_depend && la_depend) {
469  // Derivative with respect to mass
470  if (ext_do_jac[0]) {
471  dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
472  dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
473  }
474  // Derivative with respect to mean size
475  if (ext_do_jac[1]) {
476  // 1. Term associated with n0
477  // Calculated as dpsd/dn0 * dn0/dla * dla/dXm
478  dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
479  dpsd_data_dx(ext_i_jac[1], ip, joker) *=
480  ext_pars[0] * mgd_pars[3] * eterm *
481  pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
482  // 2. Term associated with la
483  // Calculated as dpsd/dla * dla/dXm
484  dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
485  // Apply dla/dXm to sum
486  dpsd_data_dx(ext_i_jac[1], ip, joker) *=
487  -mgd_pars[3] * scfac2 * pow(ext_pars[1], -(mgd_pars[3] + 1));
488  }
489  } else {
490  assert(0);
491  }
492  }
493 
494  // *** Median size ***
495  else if (something == "median size") {
496  if (n0_depend && la_depend) {
497  // Derivative with respect to mass
498  if (ext_do_jac[0]) {
499  dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
500  dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
501  }
502  // Derivative with respect to median size
503  if (ext_do_jac[1]) {
504  // 1. Term associated with n0
505  // Calculated as dpsd/dn0 * dn0/dla * dla/dXm
506  dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
507  dpsd_data_dx(ext_i_jac[1], ip, joker) *=
508  ext_pars[0] * mgd_pars[3] * eterm *
509  pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
510  // 2. Term associated with la
511  // Calculated as dpsd/dla * dla/dXm
512  dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
513  // Apply dla/dXm to sum
514  dpsd_data_dx(ext_i_jac[1], ip, joker) *=
515  -mgd_pars[3] * scfac2 * pow(ext_pars[1], -(mgd_pars[3] + 1));
516  }
517  } else {
518  assert(0);
519  }
520  }
521 
522  // *** Mean particle size ***
523  else if (something == "mean particle mass") {
524  if (n0_depend && la_depend) {
525  // Derivative with respect to mass
526  if (ext_do_jac[0]) {
527  dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
528  dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
529  }
530  // Derivative with respect to mean particle size
531  if (ext_do_jac[1]) {
532  // 1. Term associated with n0
533  // Calculated as dpsd/dn0 * dn0/dla * dla/dMm
534  dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
535  dpsd_data_dx(ext_i_jac[1], ip, joker) *=
536  ext_pars[0] * mgd_pars[3] * eterm *
537  pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
538  // 2. Term associated with la
539  // Calculated as dpsd/dla * dla/dMm
540  dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
541  // Apply dla/dMm to sum
542  dpsd_data_dx(ext_i_jac[1], ip, joker) *=
543  scfac2 * (-mgd_pars[3] / scat_species_b) *
544  pow(ext_pars[1], -(gab + 1));
545  } else {
546  assert(0);
547  }
548  }
549  }
550 
551  else if (something == "Ntot") {
552  if (n0_depend && la_depend) {
553  // Term part of both derivatives
554  const Numeric dn0dla = ext_pars[0] * mgd_pars[3] * eterm *
555  pow(mgd_pars[2], eterm - 1) /
556  (scat_species_a * gterm);
557  // Derivative with respect to mass
558  if (ext_do_jac[0]) {
559  // Repeated term
560  const Numeric dladw = scfac2 * pow(ext_pars[1], gab) *
561  (-mgd_pars[3] / scat_species_b) *
562  pow(ext_pars[0], -(gab + 1));
563  // 1. Term associated with n0
564  dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
565  dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1 + dn0dla * dladw;
566  // 2. Term associated with la
567  Vector term2 = jac_data(2, joker);
568  term2 *= dladw;
569  // Sum up
570  dpsd_data_dx(ext_i_jac[0], ip, joker) += term2;
571  }
572  // Derivative with respect to Ntot
573  if (ext_do_jac[1]) {
574  // 1. Term associated with n0
575  // Calculated as dpsd/dn0 * dn0/dla * dla/dNtot
576  dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
577  dpsd_data_dx(ext_i_jac[1], ip, joker) *= dn0dla;
578  // 2. Term associated with la
579  // Calculated as dpsd/dla * dla/dNtot
580  dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
581  // Apply dla/dNtot to sum
582  dpsd_data_dx(ext_i_jac[1], ip, joker) *=
583  scfac2 * pow(ext_pars[0], -gab) *
584  (mgd_pars[3] / scat_species_b) * pow(ext_pars[1], gab - 1);
585  }
586  } else {
587  assert(0);
588  }
589  }
590 
591  // String something not recognised
592  else {
593  assert(0);
594  }
595  }
596 
597  // Derivatives for non-dependent native parameters
598  for (Index i = 0; i < 4; i++) {
599  if (mgd_do_jac[i]) {
600  dpsd_data_dx(mgd_i_jac[i], ip, joker) = jac_data(i, joker);
601  }
602  }
603  }
604 }
605 
608  const String& type,
609  const Vector& pnd_agenda_input_t,
610  const Matrix& pnd_agenda_input,
614  const Index& species_index,
615  const Numeric& t_min,
616  const Numeric& t_max,
617  const Index& picky,
618  const Verbosity&) {
619  // Standard checcks
620  const Vector psd_size_grid(1, 0); // As this WSV is not input for thse WSM
622 
623  // Extra checks for this PSD
624  const Index nss = scat_meta.nelem();
625  if (nss == 0) throw runtime_error("*scat_meta* is empty!");
626  if (nss < species_index + 1) {
627  ostringstream os;
628  os << "Selected scattering species index is " << species_index
629  << " but this "
630  << "is not allowed since *scat_meta* has only " << nss << " elements.";
631  throw runtime_error(os.str());
632  }
633  if (scat_meta[species_index].nelem() != 1) {
634  ostringstream os;
635  os << "This method only works with scattering species consisting of a\n"
636  << "single element, but your data do not match this demand.\n"
637  << "Selected scattering species index is " << species_index << ".\n"
638  << "This species has " << scat_meta[species_index].nelem()
639  << " elements.";
640  throw runtime_error(os.str());
641  }
642  //
643  if (pnd_agenda_input.ncols() != 1)
644  throw runtime_error("*pnd_agenda_input* must have one column.");
645  if (nsi != 1)
646  throw runtime_error(
647  "This method demands that length of "
648  "*psd_size_grid* is 1.");
649 
650  // Extract particle mass
651  Numeric pmass = 0;
652  if (type == "mass") {
653  pmass = scat_meta[species_index][0].mass;
654  }
655 
656  for (Index ip = 0; ip < np; ip++) {
657  // Extract the input variables
658  Numeric x = pnd_agenda_input(ip, 0);
659  Numeric t = pnd_agenda_input_t[ip];
660 
661  // No calc needed if n==0 and no jacobians requested.
662  if ((x == 0.) && (!ndx)) {
663  continue;
664  } // If here, we are ready with this point!
665 
666  // Outside of [t_min,tmax]?
667  if (t < t_min || t > t_max) {
668  if (picky) {
669  ostringstream os;
670  os << "Method called with a temperature of " << t << " K.\n"
671  << "This is outside the specified allowed range: [ max(0.," << t_min
672  << "), " << t_max << " ]";
673  throw runtime_error(os.str());
674  } else {
675  continue;
676  } // If here, we are ready with this point!
677  }
678 
679  // Set PSD
680  //
681  if (type == "ntot") {
682  psd_data(ip, 0) = x;
683  //
684  if (ndx) {
685  dpsd_data_dx(0, ip, 0) = 1;
686  }
687  } else if (type == "mass") {
688  psd_data(ip, 0) = x / pmass;
689  //
690  if (ndx) {
691  dpsd_data_dx(0, ip, 0) = 1 / pmass;
692  }
693  } else {
694  assert(0);
695  }
696  }
697 }
698 
699 void psd_rain_W16(Vector& psd, const Vector& diameter, const Numeric& rwc) {
700  Index nD = diameter.nelem();
701  psd.resize(nD);
702  psd = 0.;
703 
704  // skip calculation if RWC is 0.0
705  if (rwc == 0.0) {
706  return;
707  }
708  assert(rwc > 0.);
709 
710  // a and b relates N0 to lambda N0 = a*lambda^b
711  Numeric a = 0.000141;
712  Numeric b = 1.49;
713  Numeric c1 = DENSITY_OF_WATER * PI / 6;
714  Numeric base = c1 / rwc * a * tgamma(4);
715  Numeric exponent = 1. / (4 - b);
716  Numeric lambda = pow(base, exponent);
717  Numeric N0 = a * pow(lambda, b);
718 
719  //psd_general_MGD( psd, N0, 0, lambda, 1 );
720  N0 *= 1e8;
721  lambda *= 100;
722  for (Index iD = 0; iD < nD; iD++) {
723  psd[iD] = N0 * exp(-lambda * diameter[iD]);
724  }
725 }
726 
729  const String& psd_name,
730  const Vector& psd_size_grid,
731  const Vector& pnd_agenda_input_t,
732  const Matrix& pnd_agenda_input,
735  const Numeric& scat_species_a,
736  const Numeric& scat_species_b,
737  const Numeric& n_alpha_in,
738  const Numeric& n_b_in,
739  const Numeric& mu_in,
740  const Numeric& gamma_in,
741  const Numeric& t_min,
742  const Numeric& t_max,
743  const Index& picky,
744  const Verbosity&) {
745  // Standard checks
747 
748  // All PSDs are defined in terms of hydrometeor mass content
749  if (pnd_agenda_input.ncols() != 1)
750  throw runtime_error("*pnd_agenda_input* must have one column.");
751 
752  // Extra checks for rain PSDs which should be consistent with
753  // spherical liquid drops
754  if (psd_name == "Abel12" || psd_name == "Wang16"){
755  if (scat_species_b < 2.9 || scat_species_b > 3.1) {
756  ostringstream os;
757  os << "This PSD treats rain, using Dveq as size grid.\n"
758  << "This means that *scat_species_b* should be close to 3,\n"
759  << "but it is outside of the tolerated range of [2.9,3.1].\n"
760  << "Your value of *scat_species_b* is: " << scat_species_b;
761  throw runtime_error(os.str());
762  }
763  if (scat_species_a < 500 || scat_species_a > 540) {
764  ostringstream os;
765  os << "This PSD treats rain, using Dveq as size grid.\n"
766  << "This means that *scat_species_a* should be close to 520,\n"
767  << "but it is outside of the tolerated range of [500,540].\n"
768  << "Your value of *scat_species_a* is: " << scat_species_a;
769  throw runtime_error(os.str());
770  }
771  }
772  // Extra checks for graupel/hail PSDs which assume constant effective density
773  //
774  if (psd_name == "Field19"){
775  if (scat_species_b < 2.8 || scat_species_b > 3.2) {
776  ostringstream os;
777  os << "This PSD treats graupel, assuming a constant effective density.\n"
778  << "This means that *scat_species_b* should be close to 3,\n"
779  << "but it is outside of the tolerated range of [2.8,3.2].\n"
780  << "Your value of *scat_species_b* is: " << scat_species_b;
781  throw runtime_error(os.str());
782  }
783  }
784 
785  for (Index ip = 0; ip < np; ip++) {
786  // Extract the input variables
787  Numeric water_content = pnd_agenda_input(ip, 0);
788  Numeric t = pnd_agenda_input_t[ip];
789 
790  // No calc needed if water_content==0 and no jacobians requested.
791  if ((water_content == 0.) && (!ndx)) {
792  continue;
793  } // If here, we are ready with this point!
794 
795  // Outside of [t_min,tmax]?
796  if (t < t_min || t > t_max) {
797  if (picky) {
798  ostringstream os;
799  os << "Method called with a temperature of " << t << " K.\n"
800  << "This is outside the specified allowed range: [ max(0.," << t_min
801  << "), " << t_max << " ]";
802  throw runtime_error(os.str());
803  } else {
804  continue;
805  } // If here, we are ready with this point!
806  }
807 
808  // Negative wc?
809  Numeric psd_weight = 1.0;
810  if (water_content < 0) {
811  psd_weight = -1.0;
812  water_content *= -1.0;
813  }
814 
815  // PSD settings for different parametrizations
816  Numeric gamma = 0.0;
817  Numeric n_alpha = 0.0;
818  Numeric n_b = 0.0;
819  Numeric mu = 0.0;
820  if (psd_name == "Abel12"){
821  n_alpha = 0.22;
822  n_b = 2.2;
823  mu = 0.0;
824  gamma = 1.0;
825  }
826  else if (psd_name == "Wang16"){
827  // Wang 16 parameters converted to SI units
828  n_alpha = 14.764;
829  n_b = 1.49;
830  mu = 0.0;
831  gamma = 1.0;
832  }
833  else if (psd_name == "Field19"){
834  n_alpha = 7.9e9;
835  n_b = -2.58;
836  mu = 0.0;
837  gamma = 1.0;
838  }
839  else if (psd_name == "generic"){
840  n_alpha = n_alpha_in;
841  n_b = n_b_in;
842  mu = mu_in;
843  gamma = gamma_in;
844  }
845  else {
846  assert(0);
847  }
848 
849  // Calculate PSD
850  // Calculate lambda for modified gamma distribution from mass density
851  Numeric k = (scat_species_b + mu + 1 - gamma)/gamma;
852  Numeric expo = 1.0 / (n_b - k - 1);
853  Numeric denom = scat_species_a * n_alpha * tgamma(k + 1);
854  Numeric lam = pow(water_content*gamma/denom, expo);
855  Numeric n_0 = n_alpha * pow(lam, n_b);
856  Vector psd_1p(nsi);
857  Matrix jac_data(4, nsi);
858 
859  psd_1p = 0.0;
860  if (water_content != 0) {
861  mgd_with_derivatives(psd_1p, jac_data, psd_size_grid, n_0, mu, lam, gamma,
862  true, // n_0 jacobian
863  false,// mu jacobian
864  true, // lambda jacobian
865  false); // gamma jacobian
866  } else {
867  assert(0);
868  }
869  //
870  for (Index i = 0; i < nsi; i++) {
871  psd_data(ip, i) = psd_weight * psd_1p[i];
872  }
873 
874  // Calculate derivative with respect to water content
875  if (ndx) {
876  const Numeric dlam_dwc = pow(gamma/denom, expo) * expo * pow(water_content, expo-1);
877  const Numeric dn_0_dwc = n_alpha * n_b * pow(lam, n_b-1) * dlam_dwc;
878  for (Index i = 0; i < nsi; i++) {
879  dpsd_data_dx(0, ip, i) = psd_weight * (jac_data(0,i)*dn_0_dwc +
880  jac_data(2,i)*dlam_dwc);
881  }
882  }
883  }
884 }
885 
887  const Vector& diameter,
888  const Numeric& swc,
889  const Numeric& t,
890  const Numeric alpha,
891  const Numeric beta,
892  const String& regime) {
893  Index nD = diameter.nelem();
894  psd.resize(nD);
895  psd = 0.;
896 
897  assert(t > 0.);
898 
899  // skip calculation if SWC is 0.0
900  if (swc == 0.0) {
901  return;
902  }
903  assert(swc > 0.);
904 
905  Numeric An, Bn, Cn;
906  Numeric M2, Mn, M2Mn;
907  Numeric base, pp;
908  Numeric x, phi23;
909  //Numeric dN;
910 
911  Vector q(5);
912 
913  assert((regime == "TR") || (regime == "ML"));
914 
915  //factors of phi23
916  if (regime == "TR")
917  q = {152., -12.4, 3.28, -0.78, -1.94};
918  else // if( regime=="ML" )
919  q = {141., -16.8, 102., 2.07, -4.82};
920 
921  //Factors of factors of the moment estimation parametrization
922  Vector Aq{13.6, -7.76, 0.479};
923  Vector Bq{-0.0361, 0.0151, 0.00149};
924  Vector Cq{0.807, 0.00581, 0.0457};
925 
926  //convert T from Kelvin to Celsius
927  Numeric Tc = t - 273.15;
928 
929  // estimate second moment
930  M2 = swc / alpha;
931  if (beta != 2) {
932  // calculate factors of the moment estimation parametrization
933  An = exp(Aq[0] + Aq[1] * beta + Aq[2] * pow(beta, 2));
934  Bn = Bq[0] + Bq[1] * beta + Bq[2] * pow(beta, 2);
935  Cn = Cq[0] + Cq[1] * beta + Cq[2] * pow(beta, 2);
936 
937  base = M2 * exp(-Bn * Tc) / An;
938  pp = 1. / (Cn);
939 
940  M2 = pow(base, pp);
941  }
942 
943  // order of the moment parametrization
944  const Numeric n = 3;
945 
946  // calculate factors of the moment estimation parametrization
947  An = exp(Aq[0] + Aq[1] * n + Aq[2] * pow(n, 2));
948  Bn = Bq[0] + Bq[1] * n + Bq[2] * pow(n, 2);
949  Cn = Cq[0] + Cq[1] * n + Cq[2] * pow(n, 2);
950 
951  // moment parametrization
952  Mn = An * exp(Bn * Tc) * pow(M2, Cn);
953 
954  M2Mn = pow(M2, 4.) / pow(Mn, 3.);
955 
956  for (Index iD = 0; iD < nD; iD++) {
957  // define x
958  x = diameter[iD] * M2 / Mn;
959 
960  // characteristic function
961  phi23 = q[0] * exp(q[1] * x) + q[2] * pow(x, q[3]) * exp(q[4] * x);
962 
963  // distribution function
964  //dN = phi23*M2Mn;
965 
966  //if ( !std::isnan(psd[dN]) )
967  // psd[iD] = dN;
968 
969  // set psd directly. Non-NaN should be (and is, hopefully) ensured by
970  // checks above (if we encounter further NaN, that should be handled above.
971  // which intermediate quantities make problems? at which parametrisation
972  // values, ie WC, T, alpha, beta?).
973  psd[iD] = phi23 * M2Mn;
974  }
975 }
976 
977 void psd_SB06(Vector& psd,
978  Matrix& dpsd,
979  const Vector& mass,
980  const Numeric& N_tot,
981  const Numeric& WC,
982  const String& hydrometeor_type) {
983  Numeric N0;
984  Numeric Lambda;
985  Numeric arg1;
986  Numeric arg2;
987  Numeric brk;
988  Numeric mu;
989  Numeric gamma;
990  Numeric xmin;
991  Numeric xmax;
992  Numeric M0min;
993  Numeric M0max;
994  Numeric M0;
995  Numeric M1;
996  Numeric c1;
997  Numeric c2;
998  Numeric L1;
999  Numeric mMu;
1000  Numeric mGamma;
1001  Numeric brkMu1;
1002 
1003  // Get the coefficients for the right hydrometeor
1004  if (hydrometeor_type == "cloud_ice") //Cloud ice water
1005  {
1006  mu = 0.;
1007  gamma = 1. / 3.;
1008  xmin = 1e-12;
1009  xmax = 1e-5;
1010  } else if (hydrometeor_type == "rain") //Rain
1011  {
1012  mu = 0.;
1013  gamma = 1. / 3.;
1014  xmin = 2.6e-10;
1015  xmax = 3e-6;
1016  } else if (hydrometeor_type == "snow") //Snow
1017  {
1018  mu = 0.;
1019  gamma = 1. / 2.;
1020  xmin = 1e-10;
1021  xmax = 2e-5;
1022  } else if (hydrometeor_type == "graupel") //Graupel
1023  {
1024  mu = 1.;
1025  gamma = 1. / 3.;
1026  xmin = 1e-9;
1027  xmax = 5e-4;
1028  } else if (hydrometeor_type == "hail") //Hail
1029  {
1030  mu = 1.;
1031  gamma = 1. / 3.;
1032  xmin = 2.6e-10;
1033  xmax = 5e-4;
1034  } else if (hydrometeor_type == "cloud_water") //Cloud liquid water
1035  {
1036  mu = 1;
1037  gamma = 1;
1038  xmin = 4.2e-15;
1039  xmax = 2.6e-10;
1040  } else {
1041  ostringstream os;
1042  os << "You use a wrong tag! ";
1043  throw runtime_error(os.str());
1044  }
1045 
1046  M0 = N_tot;
1047  M1 = WC;
1048 
1049  Index nD = mass.nelem();
1050  psd.resize(nD);
1051  psd = 0.;
1052 
1053  dpsd.resize(nD, 2);
1054  dpsd = 0.;
1055 
1056  if (M1 > 0.0) {
1057  // lower and upper limit check is taken from the ICON code of the two moment
1058  //scheme
1059 
1060  M0max = M1 / xmax;
1061  M0min = M1 / xmin;
1062 
1063  //check lower limit of the scheme
1064  if (M0 > M0min) {
1065  M0 = M0min;
1066  }
1067 
1068  //check upper limit of the scheme
1069  if (M0 < M0max) {
1070  M0 = M0max;
1071  }
1072 
1073  //arguments for Gamma function
1074  arg2 = (mu + 2) / gamma;
1075  arg1 = (mu + 1) / gamma;
1076 
1077  // results of gamma function
1078  c1 = tgamma(arg1);
1079  c2 = tgamma(arg2);
1080 
1081  // variable to shorten the formula
1082  brk = M0 / M1 * c2 / c1;
1083  brkMu1 = pow(brk, (mu + 1));
1084 
1085  //Lambda (parameter for modified gamma distribution)
1086  Lambda = pow(brk, gamma);
1087 
1088  L1 = pow(Lambda, arg1);
1089 
1090  //N0
1091  N0 = M0 * gamma / tgamma(arg1) * L1;
1092 
1093  // Calculate distribution function
1094  for (Index iD = 0; iD < nD; iD++) {
1095  //Distribution function
1096  psd[iD] = mod_gamma_dist(mass[iD], N0, Lambda, mu, gamma);
1097 
1098  if (std::isnan(psd[iD])) psd[iD] = 0.0;
1099  if (std::isinf(psd[iD])) psd[iD] = 0.0;
1100 
1101  //Calculate derivatives analytically
1102  mMu = pow(mass[iD], mu);
1103  mGamma = pow(mass[iD], gamma);
1104 
1105  // dpsd/dM1
1106  dpsd(iD, 0) = gamma / c1 * M0 / M1 * mMu * exp(-Lambda * mGamma) *
1107  brkMu1 * (-1 - mu + gamma * mGamma * Lambda);
1108 
1109  // dpsd/dM0
1110  dpsd(iD, 1) = -gamma / c1 * mMu * exp(-Lambda * mGamma) * brkMu1 *
1111  (-2 - mu - gamma * mGamma * Lambda);
1112  }
1113  } else {
1114  return;
1115  }
1116 }
1117 
1118 void psd_MY05(Vector& psd,
1119  Matrix& dpsd,
1120  const Vector& diameter_max,
1121  const Numeric N_tot,
1122  const Numeric WC,
1123  const String psd_type) {
1124  Numeric N0;
1125  Numeric Lambda;
1126  Numeric arg1;
1127  Numeric arg2;
1128  Numeric temp;
1129  Numeric mu;
1130  Numeric gamma;
1131  Numeric alpha;
1132  Numeric beta;
1133  Numeric M0;
1134  Numeric M1;
1135  Numeric c1;
1136  Numeric c2;
1137  Numeric Lmg;
1138  Numeric DMu;
1139  Numeric DGamma;
1140 
1141  // Get the coefficients for the right hydrometeor
1142  if (psd_type == "cloud_ice") //Cloud ice water
1143  {
1144  mu = 0.;
1145  gamma = 1.;
1146  alpha = 440.; //[kg]
1147  beta = 3;
1148  } else if (psd_type == "rain") //Rain
1149  {
1150  mu = 0.;
1151  gamma = 1;
1152  alpha = 523.5988; //[kg]
1153  beta = 3;
1154  } else if (psd_type == "snow") //Snow
1155  {
1156  mu = 0.;
1157  gamma = 1;
1158  alpha = 52.35988; //[kg]
1159  beta = 3;
1160  } else if (psd_type == "graupel") //Graupel
1161  {
1162  mu = 0.;
1163  gamma = 1;
1164  alpha = 209.4395; //[kg]
1165  beta = 3;
1166  } else if (psd_type == "hail") //Hail
1167  {
1168  mu = 0.;
1169  gamma = 1;
1170  alpha = 471.2389; //[kg]
1171  beta = 3;
1172  } else if (psd_type == "cloud_water") //Cloud liquid water
1173  {
1174  mu = 1;
1175  gamma = 1;
1176  alpha = 523.5988; //[kg]
1177  beta = 3;
1178  } else {
1179  ostringstream os;
1180  os << "You use a wrong tag! ";
1181  throw runtime_error(os.str());
1182  }
1183 
1184  M0 = N_tot;
1185  M1 = WC;
1186 
1187  Index nD = diameter_max.nelem();
1188  psd.resize(nD);
1189  psd = 0.;
1190 
1191  dpsd.resize(nD, 2);
1192  dpsd = 0.;
1193 
1194  if (M1 > 0.0 && M0 > 0) {
1195  //arguments for Gamma function
1196  arg2 = (mu + beta + 1) / gamma;
1197  arg1 = (mu + 1) / gamma;
1198 
1199  // results of gamma function
1200  c1 = tgamma(arg1);
1201  c2 = tgamma(arg2);
1202 
1203  //base of lambda
1204  temp = alpha * M0 / M1 * c2 / c1;
1205 
1206  //Lambda (parameter for modified gamma distribution)
1207  Lambda = pow(temp, gamma / beta);
1208 
1209  Lmg = pow(Lambda, arg1);
1210 
1211  //N0
1212  N0 = M0 * gamma / c1 * Lmg;
1213 
1214  //Distribution function
1215 
1216  // Calculate distribution function
1217  for (Index iD = 0; iD < nD; iD++) {
1218  psd[iD] = mod_gamma_dist(diameter_max[iD], N0, Lambda, mu, gamma);
1219 
1220  if (std::isnan(psd[iD])) psd[iD] = 0.0;
1221  if (std::isinf(psd[iD])) psd[iD] = 0.0;
1222 
1223  //Calculate derivatives analytically
1224  DMu = pow(diameter_max[iD], mu);
1225  DGamma = pow(diameter_max[iD], gamma);
1226 
1227  // dpsd/dM1
1228  dpsd(iD, 0) = (DMu * exp(-DGamma * Lambda) * gamma * M0 * Lmg *
1229  (-1 - mu + DGamma * gamma * Lambda) / (M1 * beta * c1));
1230 
1231  // dpsd/dM0
1232  dpsd(iD, 1) = (DMu * exp(-DGamma * Lambda) * gamma * Lmg *
1233  (1 + beta + mu - DGamma * gamma * Lambda) / (beta * c1));
1234  }
1235 
1236  } else {
1237  return;
1238  }
1239 }
1240 
1242  if (iwc == 0.0) {
1243  return 1e-9;
1244  } else {
1245  return pow(256.0 * iwc / PI / rho / n0, 0.25);
1246  }
1247 }
1248 
1250  if (dm > 1e-9) {
1251  return 256.0 * iwc / PI / rho / pow(dm, 4.0);
1252  } else {
1253  return 0.0;
1254  }
1255 }
1256 
1257 Numeric n0_from_t(Numeric t) { return exp(-0.076586 * (t - 273.15) + 17.948); }
Matrix
The Matrix class.
Definition: matpackI.h:1193
ran_gaussian
Numeric ran_gaussian(Rng &rng, const Numeric sigma)
ran_gaussian.
Definition: mc_antenna.cc:38
ARTS::Var::scat_meta
ArrayOfArrayOfScatteringMetaData scat_meta(Workspace &ws) noexcept
Definition: autoarts.h:5983
ARTS::Var::scat_species_a
Numeric scat_species_a(Workspace &ws) noexcept
Definition: autoarts.h:6103
psd_MY05
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:1118
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
joker
const Joker joker
temp
#define temp
Definition: legacy_continua.cc:20951
psd_SB06
void psd_SB06(Vector &psd, Matrix &dpsd, const Vector &mass, const Numeric &N_tot, const Numeric &WC, const String &hydrometeor_type)
Definition: psd.cc:977
ARTS::Var::psd_size_grid
Vector psd_size_grid(Workspace &ws) noexcept
Definition: autoarts.h:5460
Rng::seed
void seed(unsigned long int n, const Verbosity &verbosity)
Seeds the Rng with the integer argument.
Definition: rng.cc:54
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
b2
#define b2
Definition: complex.h:59
ARTS::Var::dpnd_data_dx_names
ArrayOfString dpnd_data_dx_names(Workspace &ws) noexcept
Definition: autoarts.h:3235
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
rng.h
member functions of the Rng class and gsl_rng code
Array< String >
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
ARTS::Var::dpsd_data_dx
Tensor3 dpsd_data_dx(Workspace &ws) noexcept
Definition: autoarts.h:3281
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
min
#define min(a, b)
Definition: legacy_continua.cc:20628
psd_rain_W16
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:699
messages.h
Declarations having to do with the four output streams.
my_basic_string< char >
mod_gamma_dist
Numeric mod_gamma_dist(Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma)
Generalized Modified Gamma Distribution.
Definition: math_funcs.cc:653
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
PI
const Numeric PI
n0_from_iwc_dm
Numeric n0_from_iwc_dm(Numeric iwc, Numeric dm, Numeric rho)
Derives N0star from IWC and Dm.
Definition: psd.cc:1249
physics_funcs.h
This file contains declerations of functions of physical character.
mgd_with_derivatives
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:530
START_OF_PSD_METHODS
#define START_OF_PSD_METHODS()
Definition: psd.h:45
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
Verbosity
Definition: messages.h:49
dm_from_iwc_n0
Numeric dm_from_iwc_n0(Numeric iwc, Numeric n0, Numeric rho)
Derives Dm from IWC and N0star.
Definition: psd.cc:1241
math_funcs.h
ARTS::Var::time
Time time(Workspace &ws) noexcept
Definition: autoarts.h:7043
mc_antenna.h
Workspace functions for the solution of cloud-box radiative transfer by Monte Carlo methods....
Matrix::resize
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
a1
#define a1
Definition: complex.h:56
cloudbox.h
Internal cloudbox functions.
n0_from_t
Numeric n0_from_t(Numeric t)
Sets N0star based on temperature.
Definition: psd.cc:1257
psd.h
Internal functions associated with size distributions.
ARTS::Var::scat_species_b
Numeric scat_species_b(Workspace &ws) noexcept
Definition: autoarts.h:6112
psd_snow_F07
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:886
ppath.h
Propagation path structure and functions.
q
#define q
Definition: legacy_continua.cc:21712
DENSITY_OF_WATER
const Numeric DENSITY_OF_WATER
logic.h
Header file for logic.cc.
Rng
Definition: rng.h:554
DENSITY_OF_ICE
const Numeric DENSITY_OF_ICE
ARTS::Var::pnd_agenda_input_t
Vector pnd_agenda_input_t(Workspace &ws) noexcept
Definition: autoarts.h:5038
ARTS::Var::pnd_agenda_input_names
ArrayOfString pnd_agenda_input_names(Workspace &ws) noexcept
Definition: autoarts.h:5023
psd_mgd_mass_and_something
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
psd_mgd_smm_common
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:727
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
ARTS::Var::mc_seed
Index mc_seed(Workspace &ws) noexcept
Definition: autoarts.h:4369
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
check_input.h
psd_cloudice_MH97
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
Vector
The Vector class.
Definition: matpackI.h:860
beta
#define beta
Definition: legacy_continua.cc:21435
sorting.h
Contains sorting routines.
psd_mono_common
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:606
ARTS::Var::pnd_agenda_input
Matrix pnd_agenda_input(Workspace &ws) noexcept
Definition: autoarts.h:5011
ARTS::Group::Verbosity
Verbosity Verbosity
Definition: autoarts.h:114
b1
#define b1
Definition: complex.h:57
ARTS::Var::psd_data
Matrix psd_data(Workspace &ws) noexcept
Definition: autoarts.h:5446
arts.h
The global header file for ARTS.