40 const Vector& diameter,
44 Index nD = diameter.nelem();
58 for (Index iD = 0; iD < nD; iD++) d_um[iD] = 1e6 * diameter[iD];
60 Numeric Tc = t - 273.15;
63 Numeric ciwc = iwc * 1e3;
66 Numeric sig_a = 0., sig_b1 = 0.;
67 Numeric sig_b2 = 0., sig_m = 0.;
68 Numeric sig_aamu = 0., sig_bamu = 0.;
69 Numeric sig_abmu = 0., sig_bbmu = 0.;
70 Numeric sig_aasigma = 0., sig_basigma = 0;
71 Numeric sig_absigma = 0., sig_bbsigma = 0.;
76 sig_a = 0.068, sig_b1 = 0.054;
77 sig_b2 = 5.5e-3, sig_m = 0.0029;
78 sig_aamu = 0.02, sig_bamu = 0.0005;
79 sig_abmu = 0.023, sig_bbmu = 0.5e-3;
80 sig_aasigma = 0.02, sig_basigma = 0.5e-3;
81 sig_absigma = 0.023, sig_bbsigma = 4.7e-4;
83 sig_a = rng.
get<std::normal_distribution>(0.0, sig_a)();
84 sig_b1 = rng.
get<std::normal_distribution>(0.0, sig_b1)();
85 sig_b2 = rng.
get<std::normal_distribution>(0.0, sig_b2)();
86 sig_m = rng.
get<std::normal_distribution>(0.0, sig_m)();
87 sig_aamu = rng.
get<std::normal_distribution>(0.0, sig_aamu)();
88 sig_bamu = rng.
get<std::normal_distribution>(0.0, sig_bamu)();
89 sig_abmu = rng.
get<std::normal_distribution>(0.0, sig_abmu)();
90 sig_bbmu = rng.
get<std::normal_distribution>(0.0, sig_bbmu)();
91 sig_aasigma = rng.
get<std::normal_distribution>(0.0, sig_aasigma)();
92 sig_basigma = rng.
get<std::normal_distribution>(0.0, sig_basigma)();
93 sig_absigma = rng.
get<std::normal_distribution>(0.0, sig_absigma)();
94 sig_bbsigma = rng.
get<std::normal_distribution>(0.0, sig_bbsigma)();;
100 Numeric
a = 0.252 + sig_a;
101 Numeric b1 = 0.837 + sig_b1;
102 Numeric IWCs100 =
min(ciwc,
a * pow(ciwc, b1));
103 Numeric IWCl100 = ciwc - IWCs100;
107 Numeric b2 = -4.99e-3 + sig_b2;
108 Numeric m = 0.0494 + sig_m;
109 Numeric alphas100 = b2 - m * log10(IWCs100);
118 Vector dNdD1(nD, 0.);
119 if (alphas100 > 0.) {
120 Numeric Ns100 = 6 * IWCs100 * pow(alphas100, 5.) /
121 (
PI * cdensity * tgamma(5.));
122 for (Index iD = 0; iD < nD; iD++)
123 dNdD1[iD] = 1e18 * Ns100 * d_um[iD] *
124 exp(-alphas100 * d_um[iD]);
131 Vector dNdD2(nD, 0.);
135 Numeric aamu = 5.20 + sig_aamu;
136 Numeric bamu = 0.0013 + sig_bamu;
137 Numeric abmu = 0.026 + sig_abmu;
138 Numeric bbmu = -1.2e-3 + sig_bbmu;
139 Numeric amu = aamu + bamu * Tc;
140 Numeric bmu = abmu + bbmu * Tc;
141 Numeric mul100 = amu + bmu * log10(IWCl100);
143 Numeric aasigma = 0.47 + sig_aasigma;
144 Numeric basigma = 2.1e-3 + sig_basigma;
145 Numeric absigma = 0.018 + sig_absigma;
146 Numeric bbsigma = -2.1e-4 + sig_bbsigma;
147 Numeric asigma = aasigma + basigma * Tc;
148 Numeric bsigma = absigma + bbsigma * Tc;
149 Numeric sigmal100 = asigma + bsigma * log10(IWCl100);
151 if ((mul100 > 0.) & (sigmal100 > 0.)) {
152 Numeric a1 = 6 * IWCl100;
153 Numeric a2_fac = pow(
PI, 3. / 2.) * cdensity * sqrt(2) *
154 exp(3 * mul100 + 9. / 2. * pow(sigmal100, 2)) *
155 sigmal100 * pow(1., 3);
157 for (Index iD = 0; iD < nD; iD++)
158 dNdD2[iD] = 1e18 * a1 / (a2_fac * d_um[iD]) *
159 exp(-0.5 * pow((log(d_um[iD]) - mul100) / sigmal100, 2));
164 for (Index iD = 0; iD < nD; iD++) {
169 psd[iD] = (dNdD1[iD] + dNdD2[iD]) * 1e6;
174 Tensor3& dpsd_data_dx,
176 const Vector& psd_size_grid,
177 const Vector& pnd_agenda_input_t,
178 const Matrix& pnd_agenda_input,
181 const Numeric& scat_species_a,
182 const Numeric& scat_species_b,
187 const Numeric& t_min,
188 const Numeric& t_max,
195 if (nin < 1 || nin > 4)
197 "The number of columns in *pnd_agenda_input* must "
199 if (scat_species_a <= 0)
200 throw runtime_error(
"*scat_species_a* should be > 0.");
201 if (scat_species_b <= 0 || scat_species_b >= 5)
202 throw runtime_error(
"*scat_species_b* should be > 0 and < 5.");
205 const Index n0_depend = (Index)n0 == -999;
206 const Index mu_depend = (Index)mu == -999;
207 const Index la_depend = (Index)la == -999;
208 const Index ga_depend = (Index)ga == -999;
210 if (n0_depend + mu_depend + la_depend + ga_depend != 2)
212 "Two (but only two) of n0, mu, la and ga must be NaN, "
213 "to flag that these parameters are the ones dependent of "
214 "mass content and mean particle size.");
215 if (mu_depend || ga_depend)
217 "Sorry, mu and la are not yet allowed to be a "
218 "dependent parameter.");
220 const Index n0_fixed = (Index) !(n0_depend || std::isnan(n0));
221 const Index mu_fixed = (Index) !(mu_depend || std::isnan(mu));
222 const Index la_fixed = (Index) !(la_depend || std::isnan(la));
223 const Index ga_fixed = (Index) !(ga_depend || std::isnan(ga));
225 if (nin + n0_fixed + mu_fixed + la_fixed + ga_fixed != 4)
227 "This PSD has four free parameters. This means that "
228 "the number\nof columns in *pnd_agenda_input* and the "
229 "number of numerics\n(i.e. not -999 or NaN) and among "
230 "the GIN arguments n0, mu, la and\nga must add up to "
231 "four. And this was found not to be the case.");
234 Vector mgd_pars(4), ext_pars(2);
241 }
else if (!n0_depend) {
242 mgd_i_pai[0] = nhit++;
246 }
else if (!mu_depend) {
247 mgd_i_pai[1] = nhit++;
251 }
else if (!la_depend) {
252 mgd_i_pai[2] = nhit++;
256 }
else if (!ga_depend) {
257 mgd_i_pai[3] = nhit++;
273 for (Index i = 0; i < ndx; i++) {
278 }
else if (dx2in[i] == 1)
284 for (Index j = 0; j < 4; j++) {
285 if (dx2in[i] == mgd_i_pai[j]) {
295 for (Index ip = 0; ip < np; ip++) {
298 ext_pars[0] = pnd_agenda_input(ip, ext_i_pai[0]);
301 if ((ext_pars[0] == 0.) && (!ndx)) {
306 ext_pars[1] = pnd_agenda_input(ip, ext_i_pai[1]);
307 if (ext_pars[1] <= 0) {
309 os <<
"Negative or zero " << something <<
" found in a position where "
310 <<
"this is not allowed.";
312 os <<
"\nNote that for retrievals, " << something
313 <<
" must be set to a value > 0 even where\nmass is zero.";
314 throw std::runtime_error(os.str());
318 for (Index i = 0; i < 4; i++) {
319 if (mgd_i_pai[i] >= 0) {
320 mgd_pars[i] = pnd_agenda_input(ip, mgd_i_pai[i]);
323 Numeric t = pnd_agenda_input_t[ip];
326 if (t < t_min || t > t_max) {
329 os <<
"Method called with a temperature of " << t <<
" K.\n"
330 <<
"This is outside the specified allowed range: [ max(0.," << t_min
331 <<
"), " << t_max <<
" ]";
332 throw runtime_error(os.str());
340 Numeric mu1 = 0, mub1 = 0, eterm = 0, gterm = 0;
341 Numeric scfac1 = 0, scfac2 = 0, gab = 0;
344 if (something ==
"mean size") {
345 if (n0_depend && la_depend) {
346 mub1 = mgd_pars[1] + scat_species_b + 1;
348 throw runtime_error(
"Bad MGD parameter detected: mu + b + 1 <= 0");
349 eterm = mub1 / mgd_pars[3];
351 gterm = tgamma(eterm);
352 scfac2 = pow(tgamma(eterm+1/ga)/gterm, mgd_pars[3]);
353 mgd_pars[2] = scfac2 * pow(ext_pars[1], -mgd_pars[3]);
356 (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
357 mgd_pars[0] = scfac1 * ext_pars[0];
364 else if (something ==
"median size") {
365 if (n0_depend && la_depend) {
366 mub1 = mgd_pars[1] + scat_species_b + 1;
368 throw runtime_error(
"Bad MGD parameter detected: mu + b + 1 <= 0");
369 eterm = mub1 / mgd_pars[3];
371 scfac2 = (mgd_pars[1] + 1 + scat_species_b - 0.327 * mgd_pars[3]) /
373 mgd_pars[2] = scfac2 * pow(ext_pars[1], -mgd_pars[3]);
375 gterm = tgamma(eterm);
377 (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
378 mgd_pars[0] = scfac1 * ext_pars[0];
385 else if (something ==
"mean particle mass") {
386 if (n0_depend && la_depend) {
387 mu1 = mgd_pars[1] + 1;
389 throw runtime_error(
"Bad MGD parameter detected: mu + 1 <= 0");
390 eterm = (mgd_pars[1] + scat_species_b + 1) / mgd_pars[3];
391 gterm = tgamma(eterm);
393 gab = mgd_pars[3] / scat_species_b;
394 scfac2 = pow(scat_species_a * gterm / tgamma(mu1 / mgd_pars[3]), gab);
395 mgd_pars[2] = scfac2 * pow(ext_pars[1], -gab);
398 (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
399 mgd_pars[0] = scfac1 * ext_pars[0];
405 else if (something ==
"Ntot") {
406 if (n0_depend && la_depend) {
407 mu1 = mgd_pars[1] + 1;
409 throw runtime_error(
"Bad MGD parameter detected: mu + 1 <= 0");
410 eterm = (mgd_pars[1] + scat_species_b + 1) / mgd_pars[3];
411 gterm = tgamma(eterm);
413 gab = mgd_pars[3] / scat_species_b;
414 scfac2 = pow(scat_species_a * gterm / tgamma(mu1 / mgd_pars[3]), gab);
415 mgd_pars[2] = scfac2 * pow(ext_pars[1] / ext_pars[0], gab);
418 (mgd_pars[3] * pow(mgd_pars[2], eterm)) / (scat_species_a * gterm);
419 mgd_pars[0] = scfac1 * ext_pars[0];
431 if (mgd_pars[2] <= 0)
432 throw runtime_error(
"Bad MGD parameter detected: la <= 0");
433 if (mgd_pars[3] <= 0)
434 throw runtime_error(
"Bad MGD parameter detected: ga <= 0");
437 Matrix jac_data(4, nsi);
445 (
bool)mgd_do_jac[0] || n0_depend,
446 (
bool)mgd_do_jac[1] || mu_depend,
447 (
bool)mgd_do_jac[2] || la_depend,
448 (
bool)mgd_do_jac[3] || ga_depend);
451 if (ext_do_jac[0] | ext_do_jac[1]) {
453 if (something ==
"mean size") {
454 if (n0_depend && la_depend) {
457 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
458 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
464 dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
465 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
466 ext_pars[0] * mgd_pars[3] * eterm *
467 pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
470 dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
472 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
473 -mgd_pars[3] * scfac2 * pow(ext_pars[1], -(mgd_pars[3] + 1));
481 else if (something ==
"median size") {
482 if (n0_depend && la_depend) {
485 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
486 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
492 dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
493 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
494 ext_pars[0] * mgd_pars[3] * eterm *
495 pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
498 dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
500 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
501 -mgd_pars[3] * scfac2 * pow(ext_pars[1], -(mgd_pars[3] + 1));
509 else if (something ==
"mean particle mass") {
510 if (n0_depend && la_depend) {
513 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
514 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1;
520 dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
521 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
522 ext_pars[0] * mgd_pars[3] * eterm *
523 pow(mgd_pars[2], eterm - 1) / (scat_species_a * gterm);
526 dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
528 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
529 scfac2 * (-mgd_pars[3] / scat_species_b) *
530 pow(ext_pars[1], -(gab + 1));
537 else if (something ==
"Ntot") {
538 if (n0_depend && la_depend) {
540 const Numeric dn0dla = ext_pars[0] * mgd_pars[3] * eterm *
541 pow(mgd_pars[2], eterm - 1) /
542 (scat_species_a * gterm);
546 const Numeric dladw = scfac2 * pow(ext_pars[1], gab) *
547 (-mgd_pars[3] / scat_species_b) *
548 pow(ext_pars[0], -(gab + 1));
550 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
551 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac1 + dn0dla * dladw;
553 Vector term2{jac_data(2, joker)};
556 dpsd_data_dx(ext_i_jac[0], ip, joker) += term2;
562 dpsd_data_dx(ext_i_jac[1], ip, joker) = jac_data(0, joker);
563 dpsd_data_dx(ext_i_jac[1], ip, joker) *= dn0dla;
566 dpsd_data_dx(ext_i_jac[1], ip, joker) += jac_data(2, joker);
568 dpsd_data_dx(ext_i_jac[1], ip, joker) *=
569 scfac2 * pow(ext_pars[0], -gab) *
570 (mgd_pars[3] / scat_species_b) * pow(ext_pars[1], gab - 1);
584 for (Index i = 0; i < 4; i++) {
586 dpsd_data_dx(mgd_i_jac[i], ip, joker) = jac_data(i, joker);
593 Tensor3& dpsd_data_dx,
595 const Vector& pnd_agenda_input_t,
596 const Matrix& pnd_agenda_input,
600 const Index& species_index,
601 const Numeric& t_min,
602 const Numeric& t_max,
606 const Vector psd_size_grid(1, 0);
610 const Index nss = scat_meta.
nelem();
611 if (nss == 0)
throw runtime_error(
"*scat_meta* is empty!");
612 if (nss < species_index + 1) {
614 os <<
"Selected scattering species index is " << species_index
616 <<
"is not allowed since *scat_meta* has only " << nss <<
" elements.";
617 throw runtime_error(os.str());
619 if (scat_meta[species_index].nelem() != 1) {
621 os <<
"This method only works with scattering species consisting of a\n"
622 <<
"single element, but your data do not match this demand.\n"
623 <<
"Selected scattering species index is " << species_index <<
".\n"
624 <<
"This species has " << scat_meta[species_index].
nelem()
626 throw runtime_error(os.str());
629 if (pnd_agenda_input.ncols() != 1)
630 throw runtime_error(
"*pnd_agenda_input* must have one column.");
633 "This method demands that length of *psd_size_grid* is 1.");
637 if (type ==
"mass") {
638 pmass = scat_meta[species_index][0].mass;
641 for (Index ip = 0; ip < np; ip++) {
643 Numeric x = pnd_agenda_input(ip, 0);
644 Numeric t = pnd_agenda_input_t[ip];
647 if ((x == 0.) && (!ndx)) {
652 if (t < t_min || t > t_max) {
655 os <<
"Method called with a temperature of " << t <<
" K.\n"
656 <<
"This is outside the specified allowed range: [ max(0.," << t_min
657 <<
"), " << t_max <<
" ]";
658 throw runtime_error(os.str());
666 if (type ==
"ntot") {
670 dpsd_data_dx(0, ip, 0) = 1;
672 }
else if (type ==
"mass") {
673 psd_data(ip, 0) = x / pmass;
676 dpsd_data_dx(0, ip, 0) = 1 / pmass;
684void psd_rain_W16(Vector& psd,
const Vector& diameter,
const Numeric& rwc) {
685 Index nD = diameter.nelem();
696 Numeric
a = 0.000141;
699 Numeric base = c1 / rwc *
a * tgamma(4);
700 Numeric exponent = 1. / (4 -
b);
701 Numeric lambda = pow(base, exponent);
702 Numeric N0 =
a * pow(lambda,
b);
707 for (Index iD = 0; iD < nD; iD++) {
708 psd[iD] = N0 * exp(-lambda * diameter[iD]);
713 Tensor3& dpsd_data_dx,
715 const Vector& psd_size_grid,
716 const Vector& pnd_agenda_input_t,
717 const Matrix& pnd_agenda_input,
720 const Numeric& scat_species_a,
721 const Numeric& scat_species_b,
722 const Numeric& n_alpha_in,
723 const Numeric& n_b_in,
724 const Numeric& mu_in,
725 const Numeric& gamma_in,
726 const Numeric& t_min,
727 const Numeric& t_max,
734 if (pnd_agenda_input.ncols() != 1)
735 throw runtime_error(
"*pnd_agenda_input* must have one column.");
739 if (psd_name ==
"Abel12" || psd_name ==
"Wang16"){
740 if (scat_species_b < 2.9 || scat_species_b > 3.1) {
742 os <<
"This PSD treats rain, using Dveq as size grid.\n"
743 <<
"This means that *scat_species_b* should be close to 3,\n"
744 <<
"but it is outside of the tolerated range of [2.9,3.1].\n"
745 <<
"Your value of *scat_species_b* is: " << scat_species_b;
746 throw runtime_error(os.str());
748 if (scat_species_a < 500 || scat_species_a > 540) {
750 os <<
"This PSD treats rain, using Dveq as size grid.\n"
751 <<
"This means that *scat_species_a* should be close to 520,\n"
752 <<
"but it is outside of the tolerated range of [500,540].\n"
753 <<
"Your value of *scat_species_a* is: " << scat_species_a;
754 throw runtime_error(os.str());
759 if (psd_name ==
"Field19"){
760 if (scat_species_b < 2.8 || scat_species_b > 3.2) {
762 os <<
"This PSD treats graupel, assuming a constant effective density.\n"
763 <<
"This means that *scat_species_b* should be close to 3,\n"
764 <<
"but it is outside of the tolerated range of [2.8,3.2].\n"
765 <<
"Your value of *scat_species_b* is: " << scat_species_b;
766 throw runtime_error(os.str());
770 for (Index ip = 0; ip < np; ip++) {
772 Numeric water_content = pnd_agenda_input(ip, 0);
773 Numeric t = pnd_agenda_input_t[ip];
776 if ((water_content == 0.) && (!ndx)) {
781 if (t < t_min || t > t_max) {
784 os <<
"Method called with a temperature of " << t <<
" K.\n"
785 <<
"This is outside the specified allowed range: [ max(0.," << t_min
786 <<
"), " << t_max <<
" ]";
787 throw runtime_error(os.str());
794 Numeric psd_weight = 1.0;
795 if (water_content < 0) {
797 water_content *= -1.0;
802 Numeric n_alpha = 0.0;
805 if (psd_name ==
"Abel12"){
811 else if (psd_name ==
"Wang16"){
818 else if (psd_name ==
"Field19"){
824 else if (psd_name ==
"generic"){
825 n_alpha = n_alpha_in;
836 Numeric k = (scat_species_b + mu + 1 - gamma)/gamma;
837 Numeric expo = 1.0 / (n_b - k - 1);
838 Numeric denom = scat_species_a * n_alpha * tgamma(k + 1);
839 Numeric lam = pow(water_content*gamma/denom, expo);
840 Numeric n_0 = n_alpha * pow(lam, n_b);
842 Matrix jac_data(4, nsi);
845 if (water_content != 0) {
855 for (Index i = 0; i < nsi; i++) {
856 psd_data(ip, i) = psd_weight * psd_1p[i];
861 const Numeric dlam_dwc = pow(gamma/denom, expo) * expo * pow(water_content, expo-1);
862 const Numeric dn_0_dwc = n_alpha * n_b * pow(lam, n_b-1) * dlam_dwc;
863 for (Index i = 0; i < nsi; i++) {
864 dpsd_data_dx(0, ip, i) = psd_weight * (jac_data(0,i)*dn_0_dwc +
865 jac_data(2,i)*dlam_dwc);
872 const Vector& diameter,
878 Index nD = diameter.nelem();
891 Numeric M2, Mn, M2Mn;
902 q = {152., -12.4, 3.28, -0.78, -1.94};
904 q = {141., -16.8, 102., 2.07, -4.82};
907 Vector Aq{13.6, -7.76, 0.479};
908 Vector Bq{-0.0361, 0.0151, 0.00149};
909 Vector Cq{0.807, 0.00581, 0.0457};
912 Numeric Tc = t - 273.15;
918 An = exp(Aq[0] + Aq[1] * beta + Aq[2] * pow(beta, 2));
919 Bn = Bq[0] + Bq[1] * beta + Bq[2] * pow(beta, 2);
920 Cn = Cq[0] + Cq[1] * beta + Cq[2] * pow(beta, 2);
922 base = M2 * exp(-Bn * Tc) / An;
932 An = exp(Aq[0] + Aq[1] * n + Aq[2] * pow(n, 2));
933 Bn = Bq[0] + Bq[1] * n + Bq[2] * pow(n, 2);
934 Cn = Cq[0] + Cq[1] * n + Cq[2] * pow(n, 2);
937 Mn = An * exp(Bn * Tc) * pow(M2, Cn);
939 M2Mn = pow(M2, 4.) / pow(Mn, 3.);
941 for (Index iD = 0; iD < nD; iD++) {
943 x = diameter[iD] * M2 / Mn;
946 phi23 = q[0] * exp(q[1] * x) + q[2] * pow(x, q[3]) * exp(q[4] * x);
958 psd[iD] = phi23 * M2Mn;
965 const Numeric& N_tot,
967 const String& hydrometeor_type) {
989 if (hydrometeor_type ==
"cloud_ice")
995 }
else if (hydrometeor_type ==
"rain")
1001 }
else if (hydrometeor_type ==
"snow")
1007 }
else if (hydrometeor_type ==
"graupel")
1013 }
else if (hydrometeor_type ==
"hail")
1019 }
else if (hydrometeor_type ==
"cloud_water")
1027 os <<
"You use a wrong tag! ";
1028 throw runtime_error(os.str());
1034 Index nD = mass.nelem();
1059 arg2 = (mu + 2) / gamma;
1060 arg1 = (mu + 1) / gamma;
1067 brk = M0 / M1 * c2 / c1;
1068 brkMu1 = pow(brk, (mu + 1));
1071 Lambda = pow(brk, gamma);
1073 L1 = pow(Lambda, arg1);
1076 N0 = M0 * gamma / tgamma(arg1) * L1;
1079 for (Index iD = 0; iD < nD; iD++) {
1083 if (std::isnan(psd[iD])) psd[iD] = 0.0;
1084 if (std::isinf(psd[iD])) psd[iD] = 0.0;
1087 mMu = pow(mass[iD], mu);
1088 mGamma = pow(mass[iD], gamma);
1091 dpsd(iD, 0) = gamma / c1 * M0 / M1 * mMu * exp(-Lambda * mGamma) *
1092 brkMu1 * (-1 - mu + gamma * mGamma * Lambda);
1095 dpsd(iD, 1) = -gamma / c1 * mMu * exp(-Lambda * mGamma) * brkMu1 *
1096 (-2 - mu - gamma * mGamma * Lambda);
1105 const Vector& diameter_max,
1106 const Numeric N_tot,
1127 if (psd_type ==
"cloud_ice")
1133 }
else if (psd_type ==
"rain")
1139 }
else if (psd_type ==
"snow")
1145 }
else if (psd_type ==
"graupel")
1151 }
else if (psd_type ==
"hail")
1157 }
else if (psd_type ==
"cloud_water")
1165 os <<
"You use a wrong tag! ";
1166 throw runtime_error(os.str());
1172 Index nD = diameter_max.nelem();
1179 if (M1 > 0.0 && M0 > 0) {
1181 arg2 = (mu + beta + 1) / gamma;
1182 arg1 = (mu + 1) / gamma;
1189 temp = alpha * M0 / M1 * c2 / c1;
1192 Lambda = pow(temp, gamma / beta);
1194 Lmg = pow(Lambda, arg1);
1197 N0 = M0 * gamma / c1 * Lmg;
1202 for (Index iD = 0; iD < nD; iD++) {
1203 psd[iD] =
mod_gamma_dist(diameter_max[iD], N0, Lambda, mu, gamma);
1205 if (std::isnan(psd[iD])) psd[iD] = 0.0;
1206 if (std::isinf(psd[iD])) psd[iD] = 0.0;
1209 DMu = pow(diameter_max[iD], mu);
1210 DGamma = pow(diameter_max[iD], gamma);
1213 dpsd(iD, 0) = (DMu * exp(-DGamma * Lambda) * gamma * M0 * Lmg *
1214 (-1 - mu + DGamma * gamma * Lambda) / (M1 * beta * c1));
1217 dpsd(iD, 1) = (DMu * exp(-DGamma * Lambda) * gamma * Lmg *
1218 (1 + beta + mu - DGamma * gamma * Lambda) / (beta * c1));
1230 return pow(256.0 * iwc /
PI / rho / n0, 0.25);
1236 return 256.0 * iwc /
PI / rho / pow(dm, 4.0);
1242Numeric
n0_from_t(Numeric t) {
return exp(-0.076586 * (t - 273.15) + 17.948); }
base min(const Array< base > &x)
Min function.
The global header file for ARTS.
Constants of physical expressions as constexpr.
Index nelem() const ARTS_NOEXCEPT
A C++ standards dependent random number generator class.
auto get(Ts &&...x) const
Returns a random number generator of some random distribution.
Internal cloudbox functions.
#define ARTS_ASSERT(condition,...)
Numeric mod_gamma_dist(Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma)
Generalized Modified Gamma Distribution.
void mgd_with_derivatives(VectorView psd, MatrixView jac_data, const Vector &x, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga, const bool &do_n0_jac, const bool &do_mu_jac, const bool &do_la_jac, const bool &do_ga_jac)
Workspace functions for the solution of cloud-box radiative transfer by Monte Carlo methods....
Declarations having to do with the four output streams.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric density_of_ice_at_0c
Global constant, Density of water ice at 0C [kg/m3] source: http://en.wikipedia.org/wiki/Ice.
constexpr Numeric denity_of_water_at_4c
Global constant, Density of liquid water +4C [kg/m3] source: http://en.wikipedia.org/wiki/Water.
This file contains declerations of functions of physical character.
Propagation path structure and functions.
void psd_rain_W16(Vector &psd, const Vector &diameter, const Numeric &rwc)
The Wang16 rain DSD DEPRECATED BY NEW MGD_SMM_COMMON Only included for compatibility with "old" pnd_f...
void psd_snow_F07(Vector &psd, const Vector &diameter, const Numeric &swc, const Numeric &t, const Numeric alpha, const Numeric beta, const String ®ime)
The F07 snow PSD.
void psd_MY05(Vector &psd, Matrix &dpsd, const Vector &diameter_max, const Numeric N_tot, const Numeric WC, const String psd_type)
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.
Numeric n0_from_t(Numeric t)
Sets N0star based on temperature.
void psd_cloudice_MH97(Vector &psd, const Vector &diameter, const Numeric &iwc, const Numeric &t, const bool noisy)
The MH97 cloud ice PSD.
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.
constexpr Numeric DENSITY_OF_ICE
Numeric n0_from_iwc_dm(Numeric iwc, Numeric dm, Numeric rho)
Derives N0star from IWC and Dm.
Numeric dm_from_iwc_n0(Numeric iwc, Numeric n0, Numeric rho)
Derives Dm from IWC and N0star.
void psd_SB06(Vector &psd, Matrix &dpsd, const Vector &mass, const Numeric &N_tot, const Numeric &WC, const String &hydrometeor_type)
constexpr Numeric DENSITY_OF_WATER
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.
Internal functions associated with size distributions.
#define START_OF_PSD_METHODS()
Contains sorting routines.