ARTS 2.5.11 (git: 725533f0)
m_psd.cc
Go to the documentation of this file.
1
9/*===========================================================================
10 === External declarations
11 ===========================================================================*/
12#include <cmath>
13#include <cstdlib>
14#include <stdexcept>
15#include "array.h"
16#include "arts.h"
17#include "auto_md.h"
18#include "check_input.h"
19#include "lin_alg.h"
20#include "math_funcs.h"
21#include "messages.h"
22#include "physics_funcs.h"
23#include "psd.h"
24
25/*===========================================================================
26 === PSDs of Mono type
27 ===========================================================================*/
28
29/* Workspace method: Doxygen documentation will be auto-generated */
30void psdMonoDispersive(Matrix& psd_data,
31 Tensor3& dpsd_data_dx,
32 const Vector& pnd_agenda_input_t,
33 const Matrix& pnd_agenda_input,
34 const ArrayOfString& pnd_agenda_input_names,
35 const ArrayOfString& dpnd_data_dx_names,
36 const ArrayOfArrayOfScatteringMetaData& scat_meta,
37 const Index& species_index,
38 const Numeric& t_min,
39 const Numeric& t_max,
40 const Index& picky,
41 const Verbosity& verbosity) {
42 psd_mono_common(psd_data,
43 dpsd_data_dx,
44 "ntot",
45 pnd_agenda_input_t,
46 pnd_agenda_input,
47 pnd_agenda_input_names,
48 dpnd_data_dx_names,
49 scat_meta,
50 species_index,
51 t_min,
52 t_max,
53 picky,
54 verbosity);
55}
56
57/* Workspace method: Doxygen documentation will be auto-generated */
58void psdMonoMass(Matrix& psd_data,
59 Tensor3& dpsd_data_dx,
60 const Vector& pnd_agenda_input_t,
61 const Matrix& pnd_agenda_input,
62 const ArrayOfString& pnd_agenda_input_names,
63 const ArrayOfString& dpnd_data_dx_names,
64 const ArrayOfArrayOfScatteringMetaData& scat_meta,
65 const Index& species_index,
66 const Numeric& t_min,
67 const Numeric& t_max,
68 const Index& picky,
69 const Verbosity& verbosity) {
70 psd_mono_common(psd_data,
71 dpsd_data_dx,
72 "mass",
73 pnd_agenda_input_t,
74 pnd_agenda_input,
75 pnd_agenda_input_names,
76 dpnd_data_dx_names,
77 scat_meta,
78 species_index,
79 t_min,
80 t_max,
81 picky,
82 verbosity);
83}
84
85/*===========================================================================
86 === PSDs of MGD type
87 ===========================================================================*/
88
89/* Workspace method: Doxygen documentation will be auto-generated */
90void psdModifiedGamma(Matrix& psd_data,
91 Tensor3& dpsd_data_dx,
92 const Vector& psd_size_grid,
93 const Vector& pnd_agenda_input_t,
94 const Matrix& pnd_agenda_input,
95 const ArrayOfString& pnd_agenda_input_names,
96 const ArrayOfString& dpnd_data_dx_names,
97 const Numeric& n0,
98 const Numeric& mu,
99 const Numeric& la,
100 const Numeric& ga,
101 const Numeric& t_min,
102 const Numeric& t_max,
103 const Index& picky,
104 const Verbosity&) {
105 // Standard checks
106 START_OF_PSD_METHODS();
107
108 // Additional (basic) checks
109 ARTS_USER_ERROR_IF (nin > 4,
110 "The number of columns in *pnd_agenda_input* must "
111 "be 0, 1, 2, 3 or 4.");
112
113 // Check fixed parameters
114 const Index n0_fixed = (Index) !(std::isnan(n0));
115 const Index mu_fixed = (Index) !(std::isnan(mu));
116 const Index la_fixed = (Index) !(std::isnan(la));
117 const Index ga_fixed = (Index) !(std::isnan(ga));
118 //
119 ARTS_USER_ERROR_IF (nin + n0_fixed + mu_fixed + la_fixed + ga_fixed != 4,
120 "This PSD has four free parameters. This means that "
121 "the number\nof columns in *pnd_agenda_input* and the "
122 "number of numerics\n(i.e. non-NaN) and among "
123 "the GIN arguments n0, mu, la and\nga must add up to "
124 "four. And this was found not to be the case.");
125
126 // Create vectors to hold the four MGD and the "extra" parameters
127 Vector mgd_pars(4);
128 ArrayOfIndex mgd_i_pai = {-1, -1, -1, -1}; // Position in pnd_agenda_input
129 {
130 Index nhit = 0;
131 if (n0_fixed) {
132 mgd_pars[0] = n0;
133 } else {
134 mgd_i_pai[0] = nhit++;
135 }
136 if (mu_fixed) {
137 mgd_pars[1] = mu;
138 } else {
139 mgd_i_pai[1] = nhit++;
140 }
141 if (la_fixed) {
142 mgd_pars[2] = la;
143 } else {
144 mgd_i_pai[2] = nhit++;
145 }
146 if (ga_fixed) {
147 mgd_pars[3] = ga;
148 } else {
149 mgd_i_pai[3] = nhit++;
150 }
151 }
152
153 // Determine what derivatives to do and their positions
154 ArrayOfIndex mgd_do_jac = {
155 0,
156 0,
157 0,
158 0,
159 };
160 ArrayOfIndex mgd_i_jac = {
161 -1, -1, -1, -1}; // Position among jacobian quantities
162 //
163 for (Index i = 0; i < ndx; i++) {
164 for (Index j = 0; j < 4; j++) {
165 if (dx2in[i] == mgd_i_pai[j]) {
166 mgd_do_jac[j] = 1;
167 mgd_i_jac[j] = i;
168 break;
169 }
170 }
171 }
172
173 // Loop input data and calculate PSDs
174 for (Index ip = 0; ip < np; ip++) {
175 // Extract MGD parameters
176 for (Index i = 0; i < 4; i++) {
177 if (mgd_i_pai[i] >= 0) {
178 mgd_pars[i] = pnd_agenda_input(ip, mgd_i_pai[i]);
179 }
180 }
181 Numeric t = pnd_agenda_input_t[ip];
182
183 // No calc needed if n0==0 and no jacobians requested.
184 if ((mgd_pars[0] == 0.) && (!ndx)) {
185 continue;
186 } // If here, we are ready with this point!
187
188 // Outside of [t_min,tmax]?
189 if (t < t_min || t > t_max) {
190 ARTS_USER_ERROR_IF (picky,
191 "Method called with a temperature of ", t, " K.\n"
192 "This is outside the specified allowed range: [ max(0.,", t_min,
193 "), ", t_max, " ]")
194 continue;
195 // If here, we are ready with this point!
196 }
197
198 // Check that la and ga are OK
199 ARTS_USER_ERROR_IF (mgd_pars[2] <= 0,
200 "Bad MGD parameter detected: la <= 0");
201 ARTS_USER_ERROR_IF (mgd_pars[3] <= 0,
202 "Bad MGD parameter detected: ga <= 0");
203
204 // Calculate PSD and derivatives
205 Matrix jac_data(4, nsi);
206 //
207 mgd_with_derivatives(psd_data(ip, joker),
208 jac_data,
209 psd_size_grid,
210 mgd_pars[0],
211 mgd_pars[1],
212 mgd_pars[2],
213 mgd_pars[3],
214 mgd_do_jac[0],
215 mgd_do_jac[1],
216 mgd_do_jac[2],
217 mgd_do_jac[3]);
218 //
219 for (Index i = 0; i < 4; i++) {
220 if (mgd_do_jac[i]) {
221 dpsd_data_dx(mgd_i_jac[i], ip, joker) = jac_data(i, joker);
222 }
223 }
224 }
225}
226
227/* Workspace method: Doxygen documentation will be auto-generated */
228void psdModifiedGammaMass(Matrix& psd_data,
229 Tensor3& dpsd_data_dx,
230 const Vector& psd_size_grid,
231 const Vector& pnd_agenda_input_t,
232 const Matrix& pnd_agenda_input,
233 const ArrayOfString& pnd_agenda_input_names,
234 const ArrayOfString& dpnd_data_dx_names,
235 const Numeric& scat_species_a,
236 const Numeric& scat_species_b,
237 const Numeric& n0,
238 const Numeric& mu,
239 const Numeric& la,
240 const Numeric& ga,
241 const Numeric& t_min,
242 const Numeric& t_max,
243 const Index& picky,
244 const Verbosity&) {
245 // Standard checks
246 START_OF_PSD_METHODS();
247
248 // Additional (basic) checks
249 ARTS_USER_ERROR_IF (nin < 1 || nin > 4,
250 "The number of columns in *pnd_agenda_input* must "
251 "be 1, 2, 3 or 4.");
252 ARTS_USER_ERROR_IF (scat_species_a <= 0,
253 "*scat_species_a* should be > 0.");
254 ARTS_USER_ERROR_IF (scat_species_b <= 0 || scat_species_b >= 5,
255 "*scat_species_b* should be > 0 and < 5.");
256
257 // Check and determine dependent and fixed parameters
258 const Index n0_depend = (Index)n0 == -999;
259 const Index mu_depend = (Index)mu == -999;
260 const Index la_depend = (Index)la == -999;
261 const Index ga_depend = (Index)ga == -999;
262 //
263 ARTS_USER_ERROR_IF (n0_depend + mu_depend + la_depend + ga_depend != 1,
264 "One (but only one) of n0, mu, la and ga must be NaN, "
265 "to flag that this parameter is the one dependent of "
266 "mass content.");
267 ARTS_USER_ERROR_IF (mu_depend || ga_depend,
268 "Sorry, mu and la are not yet allowed to be the "
269 "dependent parameter.");
270 //
271 const Index n0_fixed = (Index) !(n0_depend || std::isnan(n0));
272 const Index mu_fixed = (Index) !(mu_depend || std::isnan(mu));
273 const Index la_fixed = (Index) !(la_depend || std::isnan(la));
274 const Index ga_fixed = (Index) !(ga_depend || std::isnan(ga));
275 //
276 ARTS_USER_ERROR_IF (nin + n0_fixed + mu_fixed + la_fixed + ga_fixed != 4,
277 "This PSD has four free parameters. This means that "
278 "the number\nof columns in *pnd_agenda_input* and the "
279 "number of numerics\n(i.e. not -999 or NaN) and among "
280 "the GIN arguments n0, mu, la and\nga must add up to "
281 "four. And this was found not to be the case.");
282
283 // Create vectors to hold the four MGD and the "extra" parameters
284 Vector mgd_pars(4), ext_pars(1);
285 ArrayOfIndex mgd_i_pai = {-1, -1, -1, -1}; // Position in pnd_agenda_input
286 const ArrayOfIndex ext_i_pai = {0}; // Position in pnd_agenda_input
287 {
288 Index nhit = 1; // As mass always occupies first position
289 if (n0_fixed) {
290 mgd_pars[0] = n0;
291 } else if (!n0_depend) {
292 mgd_i_pai[0] = nhit++;
293 }
294 if (mu_fixed) {
295 mgd_pars[1] = mu;
296 } else if (!mu_depend) {
297 mgd_i_pai[1] = nhit++;
298 }
299 if (la_fixed) {
300 mgd_pars[2] = la;
301 } else if (!la_depend) {
302 mgd_i_pai[2] = nhit++;
303 }
304 if (ga_fixed) {
305 mgd_pars[3] = ga;
306 } else if (!ga_depend) {
307 mgd_i_pai[3] = nhit++;
308 }
309 }
310
311 // Determine what derivatives to do and their positions
312 ArrayOfIndex mgd_do_jac = {
313 0,
314 0,
315 0,
316 0,
317 };
318 ArrayOfIndex ext_do_jac = {0};
319 ArrayOfIndex mgd_i_jac = {
320 -1, -1, -1, -1}; // Position among jacobian quantities
321 ArrayOfIndex ext_i_jac = {-1}; // Position among jacobian quantities
322 //
323 for (Index i = 0; i < ndx; i++) {
324 if (dx2in[i] == 0) // That is, mass is a derivative
325 {
326 ext_do_jac[0] = 1;
327 ext_i_jac[0] = i;
328 } else // Otherwise, either n0, mu, la or ga
329 {
330 for (Index j = 0; j < 4; j++) {
331 if (dx2in[i] == mgd_i_pai[j]) {
332 mgd_do_jac[j] = 1;
333 mgd_i_jac[j] = i;
334 break;
335 }
336 }
337 }
338 }
339
340 // Loop input data and calculate PSDs
341 for (Index ip = 0; ip < np; ip++) {
342 // Extract mass
343 ext_pars[0] = pnd_agenda_input(ip, ext_i_pai[0]);
344 // Extract core MGD parameters
345 for (Index i = 0; i < 4; i++) {
346 if (mgd_i_pai[i] >= 0) {
347 mgd_pars[i] = pnd_agenda_input(ip, mgd_i_pai[i]);
348 }
349 }
350 Numeric t = pnd_agenda_input_t[ip];
351
352 // No calc needed if mass==0 and no jacobians requested.
353 if ((ext_pars[0] == 0.) && (!ndx)) {
354 continue;
355 } // If here, we are ready with this point!
356
357 // Outside of [t_min,tmax]?
358 if (t < t_min || t > t_max) {
359 ARTS_USER_ERROR_IF (picky,
360 "Method called with a temperature of ", t, " K.\n"
361 "This is outside the specified allowed range: [ max(0.,", t_min,
362 "), ", t_max, " ]")
363 continue;
364 // If here, we are ready with this point!
365 }
366
367 // Derive the dependent parameter
368 //
369 Numeric mub1 = 0, eterm = 0, scfac = 0;
370 //
371 if (n0_depend) {
372 mub1 = mgd_pars[1] + scat_species_b + 1;
373 eterm = mub1 / mgd_pars[3];
374 scfac = (mgd_pars[3] * pow(mgd_pars[2], eterm)) /
375 (scat_species_a * tgamma(eterm));
376 mgd_pars[0] = scfac * ext_pars[0];
377 } else if (la_depend) {
378 ARTS_USER_ERROR_IF (ext_pars[0] <= 0,
379 "The mass content must be > 0 when la is "
380 "the dependent parameter.");
381 mub1 = mgd_pars[1] + scat_species_b + 1;
382 eterm = mub1 / mgd_pars[3];
383 scfac = mgd_pars[3] / (scat_species_a * mgd_pars[0] * tgamma(eterm));
384 scfac = pow(scfac, -1 / eterm);
385 mgd_pars[2] = scfac * pow(ext_pars[0], -1 / eterm);
386 } else {
387 ARTS_ASSERT(0);
388 }
389
390 // Now when all four MGD parameters are set, check that they were OK from
391 // start, or became OK if set
392 ARTS_USER_ERROR_IF (mub1 <= 0,
393 "Bad MGD parameter detected: mu + b + 1 <= 0");
394 ARTS_USER_ERROR_IF (mgd_pars[2] <= 0,
395 "Bad MGD parameter detected: la <= 0");
396 ARTS_USER_ERROR_IF (mgd_pars[3] <= 0,
397 "Bad MGD parameter detected: ga <= 0");
398
399 // Calculate PSS
400 Matrix jac_data(4, nsi);
401 mgd_with_derivatives(psd_data(ip, joker),
402 jac_data,
403 psd_size_grid,
404 mgd_pars[0],
405 mgd_pars[1],
406 mgd_pars[2],
407 mgd_pars[3],
408 (bool)mgd_do_jac[0] || n0_depend,
409 (bool)mgd_do_jac[1] || mu_depend,
410 (bool)mgd_do_jac[2] || la_depend,
411 (bool)mgd_do_jac[3] || ga_depend);
412
413 // Derivative with respect to mass
414 if (ext_do_jac[0]) {
415 if (n0_depend) {
416 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(0, joker);
417 dpsd_data_dx(ext_i_jac[0], ip, joker) *= scfac;
418 } else if (la_depend) {
419 dpsd_data_dx(ext_i_jac[0], ip, joker) = jac_data(2, joker);
420 dpsd_data_dx(ext_i_jac[0], ip, joker) *=
421 scfac * (-1 / eterm) * pow(ext_pars[0], -(1 / eterm + 1));
422 } else {
423 ARTS_ASSERT(0);
424 }
425 }
426
427 // Derivatives for non-dependent native parameters
428 for (Index i = 0; i < 4; i++) {
429 if (mgd_do_jac[i]) {
430 dpsd_data_dx(mgd_i_jac[i], ip, joker) = jac_data(i, joker);
431 }
432 }
433 }
434}
435
436/* Workspace method: Doxygen documentation will be auto-generated */
437void psdModifiedGammaMassNtot(Matrix& psd_data,
438 Tensor3& dpsd_data_dx,
439 const Vector& psd_size_grid,
440 const Vector& pnd_agenda_input_t,
441 const Matrix& pnd_agenda_input,
442 const ArrayOfString& pnd_agenda_input_names,
443 const ArrayOfString& dpnd_data_dx_names,
444 const Numeric& scat_species_a,
445 const Numeric& scat_species_b,
446 const Numeric& n0,
447 const Numeric& mu,
448 const Numeric& la,
449 const Numeric& ga,
450 const Numeric& t_min,
451 const Numeric& t_max,
452 const Index& picky,
453 const Verbosity& verbosity) {
454 psd_mgd_mass_and_something(psd_data,
455 dpsd_data_dx,
456 "Ntot",
457 psd_size_grid,
458 pnd_agenda_input_t,
459 pnd_agenda_input,
460 pnd_agenda_input_names,
461 dpnd_data_dx_names,
462 scat_species_a,
463 scat_species_b,
464 n0,
465 mu,
466 la,
467 ga,
468 t_min,
469 t_max,
470 picky,
471 verbosity);
472}
473
474/* Workspace method: Doxygen documentation will be auto-generated */
475void psdModifiedGammaMassMeanParticleMass(Matrix& psd_data,
476 Tensor3& dpsd_data_dx,
477 const Vector& psd_size_grid,
478 const Vector& pnd_agenda_input_t,
479 const Matrix& pnd_agenda_input,
480 const ArrayOfString& pnd_agenda_input_names,
481 const ArrayOfString& dpnd_data_dx_names,
482 const Numeric& scat_species_a,
483 const Numeric& scat_species_b,
484 const Numeric& n0,
485 const Numeric& mu,
486 const Numeric& la,
487 const Numeric& ga,
488 const Numeric& t_min,
489 const Numeric& t_max,
490 const Index& picky,
491 const Verbosity& verbosity) {
492 psd_mgd_mass_and_something(psd_data,
493 dpsd_data_dx,
494 "mean particle mass",
495 psd_size_grid,
496 pnd_agenda_input_t,
497 pnd_agenda_input,
498 pnd_agenda_input_names,
499 dpnd_data_dx_names,
500 scat_species_a,
501 scat_species_b,
502 n0,
503 mu,
504 la,
505 ga,
506 t_min,
507 t_max,
508 picky,
509 verbosity);
510}
511
512/* Workspace method: Doxygen documentation will be auto-generated */
513void psdModifiedGammaMassXmean(Matrix& psd_data,
514 Tensor3& dpsd_data_dx,
515 const Vector& psd_size_grid,
516 const Vector& pnd_agenda_input_t,
517 const Matrix& pnd_agenda_input,
518 const ArrayOfString& pnd_agenda_input_names,
519 const ArrayOfString& dpnd_data_dx_names,
520 const Numeric& scat_species_a,
521 const Numeric& scat_species_b,
522 const Numeric& n0,
523 const Numeric& mu,
524 const Numeric& la,
525 const Numeric& ga,
526 const Numeric& t_min,
527 const Numeric& t_max,
528 const Index& picky,
529 const Verbosity& verbosity) {
530 psd_mgd_mass_and_something(psd_data,
531 dpsd_data_dx,
532 "mean size",
533 psd_size_grid,
534 pnd_agenda_input_t,
535 pnd_agenda_input,
536 pnd_agenda_input_names,
537 dpnd_data_dx_names,
538 scat_species_a,
539 scat_species_b,
540 n0,
541 mu,
542 la,
543 ga,
544 t_min,
545 t_max,
546 picky,
547 verbosity);
548}
549
550/* Workspace method: Doxygen documentation will be auto-generated */
551void psdModifiedGammaMassXmedian(Matrix& psd_data,
552 Tensor3& dpsd_data_dx,
553 const Vector& psd_size_grid,
554 const Vector& pnd_agenda_input_t,
555 const Matrix& pnd_agenda_input,
556 const ArrayOfString& pnd_agenda_input_names,
557 const ArrayOfString& dpnd_data_dx_names,
558 const Numeric& scat_species_a,
559 const Numeric& scat_species_b,
560 const Numeric& n0,
561 const Numeric& mu,
562 const Numeric& la,
563 const Numeric& ga,
564 const Numeric& t_min,
565 const Numeric& t_max,
566 const Index& picky,
567 const Verbosity& verbosity) {
568 psd_mgd_mass_and_something(psd_data,
569 dpsd_data_dx,
570 "median size",
571 psd_size_grid,
572 pnd_agenda_input_t,
573 pnd_agenda_input,
574 pnd_agenda_input_names,
575 dpnd_data_dx_names,
576 scat_species_a,
577 scat_species_b,
578 n0,
579 mu,
580 la,
581 ga,
582 t_min,
583 t_max,
584 picky,
585 verbosity);
586}
587
588/* Workspace method: Doxygen documentation will be auto-generated */
589void psdModifiedGammaMassSingleMoment(
590 Matrix& psd_data,
591 Tensor3& dpsd_data_dx,
592 const Vector& psd_size_grid,
593 const Vector& pnd_agenda_input_t,
594 const Matrix& pnd_agenda_input,
595 const ArrayOfString& pnd_agenda_input_names,
596 const ArrayOfString& dpnd_data_dx_names,
597 const Numeric& scat_species_a,
598 const Numeric& scat_species_b,
599 const Numeric& n_alpha,
600 const Numeric& n_b,
601 const Numeric& mu,
602 const Numeric& gamma,
603 const Numeric& t_min,
604 const Numeric& t_max,
605 const Index& picky,
606 const Verbosity& verbosity) {
607 psd_mgd_smm_common(psd_data,
608 dpsd_data_dx,
609 "generic",
610 psd_size_grid,
611 pnd_agenda_input_t,
612 pnd_agenda_input,
613 pnd_agenda_input_names,
614 dpnd_data_dx_names,
615 scat_species_a,
616 scat_species_b,
617 n_alpha,
618 n_b,
619 mu,
620 gamma,
621 t_min,
622 t_max,
623 picky,
624 verbosity);
625}
626
627/*===========================================================================
628 === Input: IWC and T
629 ===========================================================================*/
630
631/* Workspace method: Doxygen documentation will be auto-generated */
632void psdDelanoeEtAl14(Matrix& psd_data,
633 Tensor3& dpsd_data_dx,
634 const Vector& psd_size_grid,
635 const Vector& pnd_agenda_input_t,
636 const Matrix& pnd_agenda_input,
637 const ArrayOfString& pnd_agenda_input_names,
638 const ArrayOfString& dpnd_data_dx_names,
639 const Numeric& iwc,
640 const Numeric& n0,
641 const Numeric& dm,
642 const Numeric& rho,
643 const Numeric& alpha,
644 const Numeric& beta,
645 const Numeric& t_min,
646 const Numeric& t_max,
647 const Numeric& dm_min,
648 const Index& picky,
649 const Verbosity&) {
650 // Standard checks
651 START_OF_PSD_METHODS();
652
653 // Additional (basic) checks
654 ARTS_USER_ERROR_IF (nin > 2,
655 "The number of columns in *pnd_agenda_input* must "
656 "be 0, 1 or 2");
657
658 // Check and determine dependent and fixed parameters
659 const bool n0_depend = (Index)n0 == -999;
660 const bool dm_depend = (Index)dm == -999;
661 const bool iwc_depend = (Index)iwc == -999;
662
663 // Check fixed parameters
664 const bool iwc_fixed = !(std::isnan(iwc)) && !iwc_depend;
665 const bool n0_fixed = !(std::isnan(n0)) && !n0_depend;
666 const bool dm_fixed = !(std::isnan(dm)) && !dm_depend;
667
668 ARTS_USER_ERROR_IF (!((nin + iwc_fixed + n0_fixed + dm_fixed == 2) ||
669 (nin + iwc_fixed + n0_fixed + dm_fixed == 1)),
670 "This PSD can have one or two independent parameters, that is \n"
671 "the sum of the number of columns in pnd_agenda_input and\n"
672 "non-NAN, non-dependent values in iwc, n0, dm must be equal to\n"
673 "one or two.");
674
675 ArrayOfIndex i_pai = {-1, -1, -1}; // Position in pnd_agenda_input
676
677 Index nhit = 0;
678
679 if ((n0_depend || dm_depend) && (!iwc_fixed)) {
680 i_pai[0] = nhit++;
681 }
682
683 if ((!n0_depend) && (!n0_fixed)) {
684 i_pai[1] = nhit++;
685 }
686
687 if ((!dm_depend) && (!dm_fixed)) {
688 i_pai[2] = nhit++;
689 }
690
691 // Determine what derivatives to do and their positions
692 ArrayOfIndex do_jac = {0, 0, 0};
693 ArrayOfIndex i_jac = {-1, -1, -1}; // Position among jacobian quantities
694
695 for (Index i = 0; i < ndx; ++i) {
696 for (Index j = 0; j < 3; ++j) {
697 if (dx2in[i] == i_pai[j]) {
698 do_jac[j] = 1;
699 i_jac[j] = i;
700 break;
701 }
702 }
703 }
704
705 if (psd_size_grid[0] < std::numeric_limits<Numeric>::epsilon()) {
706 ARTS_USER_ERROR_IF (psd_size_grid.nelem() < 2,
707 "psd_size_grid has only one element which is 0. This is not allowed.");
708 }
709
710 Numeric iwc_p(0.0), n0_p(0.0), dm_p(0.0);
711 // Loop input data and calculate PSDs
712 for (Index ip = 0; ip < np; ip++) {
713 Numeric t = pnd_agenda_input_t[ip];
714
715 // Extract MGD parameters
716 if (i_pai[0] >= 0) {
717 iwc_p = pnd_agenda_input(ip, i_pai[0]);
718 }
719 if (i_pai[1] >= 0) {
720 n0_p = pnd_agenda_input(ip, i_pai[1]);
721 }
722 if (i_pai[2] >= 0) {
723 dm_p = pnd_agenda_input(ip, i_pai[2]);
724 }
725
726 if (n0_depend && dm_depend) {
727 n0_p = n0_from_t(t);
728 dm_p = dm_from_iwc_n0(iwc_p, n0_p, rho);
729 } else if (n0_depend) {
730 n0_p = n0_from_iwc_dm(iwc_p, dm_p, rho);
731 } else if (dm_depend) {
732 dm_p = dm_from_iwc_n0(iwc_p, n0_p, rho);
733 }
734
735 // Outside of [t_min,tmax]?
736 if ((t < t_min) || (t > t_max)) {
737 ARTS_USER_ERROR_IF (picky,
738 "Method called with a temperature of ", t, " K.\n"
739 "This is outside the specified allowed range: [ max(0.,", t_min,
740 "), ", t_max, " ]")
741 continue;
742 }
743
744 if (iwc > 0.0) {
745 ARTS_USER_ERROR_IF ((iwc > 0.0) && ((dm_p <= 0.0) || (dm_p < dm_min)),
746 "The provided or inferred value of *Dm* (", dm_p, ") is "
747 " less than zero or *Dm_min* and this is not allowed. "
748 "This means that you have very small or zero values "
749 "in *pnd_agenda_input* which is not supported "
750 "by this PSD.\n")
751 }
752
753 // Calculate PSD and derivatives
754 Matrix jac_data(1, nsi);
755 Vector x_grid(psd_size_grid);
756 x_grid *= 1.0 / dm_p;
757
758 if (x_grid[0] < std::numeric_limits<Numeric>::epsilon()) {
759 x_grid[0] = 0.1 * psd_size_grid[1];
760 }
761
762 delanoe_shape_with_derivative(
763 psd_data(ip, joker), jac_data, x_grid, alpha, beta);
764 psd_data(ip, joker) *= n0_p;
765 jac_data(0, joker) *= n0_p;
766
767 Vector dndx{jac_data(0, joker)};
768 Vector dndn0{psd_data(ip, joker)};
769 dndn0 *= (1.0 / n0_p);
770
771 Vector dxddm = x_grid;
772 dxddm *= (-1.0 / dm_p);
773 Numeric dn0diwc = n0_p / iwc_p;
774
775 if (do_jac[0]) {
776 dpsd_data_dx(i_jac[0], ip, joker) = 0.0;
777
778 if (dm_depend) {
779 Numeric ddmdiwc = 0.25 * dm_p / iwc_p;
780 Vector dndiwc = dxddm;
781 dndiwc *= dndx;
782 dndiwc *= ddmdiwc;
783 dpsd_data_dx(i_jac[0], ip, joker) += dndiwc;
784 } else if (n0_depend) {
785 Vector dndiwc = dndn0;
786 dndiwc *= dn0diwc;
787 dpsd_data_dx(i_jac[0], ip, joker) += dndiwc;
788 }
789 }
790
791 if (do_jac[1]) {
792 dpsd_data_dx(i_jac[1], ip, joker) = dndn0;
793 if (dm_depend) {
794 Vector dndn02 = dndx;
795 dndn02 *= dxddm;
796 Numeric ddmdn0 = -0.25 / n0_p * dm_p;
797 dndn02 *= ddmdn0;
798 dpsd_data_dx(i_jac[1], ip, joker) += dndn02;
799 }
800 }
801
802 if (do_jac[2]) {
803 dpsd_data_dx(i_jac[2], ip, joker) = dxddm;
804 dpsd_data_dx(i_jac[2], ip, joker) *= dndx;
805 if (n0_depend) {
806 Vector dndn02 = dndn0;
807 Numeric dn0ddm = -4.0 * n0_p / dm_p;
808 dndn02 *= dn0ddm;
809 dpsd_data_dx(i_jac[2], ip, joker) += dndn02;
810 }
811 }
812
813 // Ensure that results are zero when IWC is zero.
814 if ((!iwc_depend) && (iwc_p == 0.0)) {
815 psd_data(0, joker) = 0.0;
816 for (size_t i = 0; i < 2; ++i) {
817 if (do_jac[i]) {
818 dpsd_data_dx(i_jac[i], ip, joker) = 0.0;
819 }
820 }
821 }
822 }
823}
824
825/* Workspace method: Doxygen documentation will be auto-generated */
826void psdFieldEtAl07(Matrix& psd_data,
827 Tensor3& dpsd_data_dx,
828 const Vector& psd_size_grid,
829 const Vector& pnd_agenda_input_t,
830 const Matrix& pnd_agenda_input,
831 const ArrayOfString& pnd_agenda_input_names,
832 const ArrayOfString& dpnd_data_dx_names,
833 const Numeric& scat_species_a,
834 const Numeric& scat_species_b,
835 const String& regime,
836 const Numeric& t_min,
837 const Numeric& t_max,
838 const Numeric& t_min_psd,
839 const Numeric& t_max_psd,
840 const Numeric& b_min,
841 const Numeric& b_max,
842 const Index& picky,
843 const Verbosity&) {
844 // Standard checcks
845 START_OF_PSD_METHODS();
846
847 // Additional (basic) checks
848 ARTS_USER_ERROR_IF (pnd_agenda_input.ncols() != 1,
849 "*pnd_agenda_input* must have one column.");
850 ARTS_USER_ERROR_IF (regime != "TR" && regime != "ML",
851 "regime must either be \"TR\" or \"ML\".");
852 ARTS_USER_ERROR_IF (scat_species_a <= 0,
853 "*scat_species_a* should be > 0.");
854 ARTS_USER_ERROR_IF (scat_species_b < b_min || scat_species_b > b_max,
855 "Method called with a mass-dimension-relation exponent b of ",
856 scat_species_b, ".\n"
857 "This is outside the specified allowed range: [", b_min, ",",
858 b_max, "]")
859
860 for (Index ip = 0; ip < np; ip++) {
861 // Extract the input variables
862 Numeric swc = pnd_agenda_input(ip, 0);
863 Numeric t = pnd_agenda_input_t[ip];
864
865 // NaN can be generated for extremly small SWC (such as 1e-78)
866 // As a solution set a limit, and consider all below as zero
867 if (abs(swc) < 1e-15) {
868 swc = 0.0;
869 }
870
871 // No calc needed if swc==0 and no jacobians requested.
872 if ((swc == 0.) && (!ndx)) {
873 continue;
874 } // If here, we are ready with this point!
875
876 // Outside of [t_min,tmax]?
877 if (t < t_min || t > t_max) {
878 ARTS_USER_ERROR_IF (picky,
879 "Method called with a temperature of ", t, " K.\n"
880 "This is outside the specified allowed range: [ max(0.,", t_min,
881 "), ", t_max, " ]")
882 continue;
883 // If here, we are ready with this point!
884 }
885
886 // PSD assumed to be constant outside [*t_min_psd*,*t_max_psd*]
887 if (t < t_min_psd) {
888 t = t_min_psd;
889 } else if (t > t_max_psd) {
890 t = t_max_psd;
891 }
892
893 // Negative swc?
894 Numeric psd_weight = 1.0;
895 if (swc < 0) {
896 psd_weight = -1.0;
897 swc *= -1.0;
898 }
899
900 // Calculate PSD
901 Vector psd_1p(nsi);
902 if (swc != 0) {
903 psd_snow_F07(psd_1p,
904 psd_size_grid,
905 swc,
906 t,
907 scat_species_a,
908 scat_species_b,
909 regime);
910 for (Index i = 0; i < nsi; i++) {
911 psd_data(ip, i) = psd_weight * psd_1p[i];
912 }
913 }
914
915 // Calculate derivative with respect to IWC
916 if (ndx) {
917 //const Numeric dswc = max( 0.001*swc, 1e-7 );
918 const Numeric dswc = 1e-9;
919 const Numeric swcp = swc + dswc;
920 psd_snow_F07(psd_1p,
921 psd_size_grid,
922 swcp,
923 t,
924 scat_species_a,
925 scat_species_b,
926 regime);
927 for (Index i = 0; i < nsi; i++) {
928 dpsd_data_dx(0, ip, i) = (psd_1p[i] - psd_data(ip, i)) / dswc;
929 }
930 }
931 }
932}
933
934/* Workspace method: Doxygen documentation will be auto-generated */
935void psdMcFarquaharHeymsfield97(Matrix& psd_data,
936 Tensor3& dpsd_data_dx,
937 const Vector& psd_size_grid,
938 const Vector& pnd_agenda_input_t,
939 const Matrix& pnd_agenda_input,
940 const ArrayOfString& pnd_agenda_input_names,
941 const ArrayOfString& dpnd_data_dx_names,
942 const Numeric& scat_species_a,
943 const Numeric& scat_species_b,
944 const Numeric& t_min,
945 const Numeric& t_max,
946 const Numeric& t_min_psd,
947 const Numeric& t_max_psd,
948 const Index& picky,
949 const Index& noisy,
950 const Verbosity&) {
951 // Standard checcks
952 START_OF_PSD_METHODS();
953
954 // Extra checks for this PSD
955 ARTS_USER_ERROR_IF (pnd_agenda_input.ncols() != 1,
956 "*pnd_agenda_input* must have one column.");
957 ARTS_USER_ERROR_IF (noisy && ndx,
958 "Jacobian calculations and \"noisy\" can not be "
959 "combined.");
960 ARTS_USER_ERROR_IF (scat_species_b < 2.9 || scat_species_b > 3.1,
961 "This PSD treats pure ice, using Dveq as size grid.\n"
962 "This means that *scat_species_b* should be close to 3,\n"
963 "but it is outside of the tolerated range of [2.9,3.1].\n"
964 "Your value of *scat_species_b* is: ", scat_species_b)
965 ARTS_USER_ERROR_IF (scat_species_a < 460 || scat_species_a > 500,
966 "This PSD treats pure ice, using Dveq as size grid.\n"
967 "This means that *scat_species_a* should be close to 480,\n"
968 "but it is outside of the tolerated range of [460,500].\n"
969 "Your value of *scat_species_a* is: ", scat_species_a)
970
971 for (Index ip = 0; ip < np; ip++) {
972 // Extract the input variables
973 Numeric iwc = pnd_agenda_input(ip, 0);
974 Numeric t = pnd_agenda_input_t[ip];
975
976 // No calc needed if iwc==0 and no jacobians requested.
977 if ((iwc == 0.) && (!ndx)) {
978 continue;
979 } // If here, we are ready with this point!
980
981 // Outside of [t_min,tmax]?
982 if (t < t_min || t > t_max) {
983 ARTS_USER_ERROR_IF (picky,
984 "Method called with a temperature of ", t, " K.\n"
985 "This is outside the specified allowed range: [ max(0.,", t_min,
986 "), ", t_max, " ]")
987 continue;
988 // If here, we are ready with this point!
989 }
990
991 // PSD assumed to be constant outside [*t_min_psd*,*t_max_psd*]
992 if (t < t_min_psd) {
993 t = t_min_psd;
994 } else if (t > t_max_psd) {
995 t = t_max_psd;
996 }
997
998 // Negative iwc?
999 Numeric psd_weight = 1.0;
1000 if (iwc < 0) {
1001 psd_weight = -1.0;
1002 iwc *= -1.0;
1003 }
1004
1005 // Calculate PSD
1006 Vector psd_1p(nsi);
1007 if (iwc != 0) {
1008 psd_cloudice_MH97(psd_1p, psd_size_grid, iwc, t, noisy);
1009 for (Index i = 0; i < nsi; i++) {
1010 psd_data(ip, i) = psd_weight * psd_1p[i];
1011 }
1012 }
1013
1014 // Calculate derivative with respect to IWC
1015 if (ndx) {
1016 //const Numeric diwc = max( 0.001*iwc, 1e-9 );
1017 const Numeric diwc = 1e-9;
1018 const Numeric iwcp = iwc + diwc;
1019 psd_cloudice_MH97(psd_1p, psd_size_grid, iwcp, t, noisy);
1020 for (Index i = 0; i < nsi; i++) {
1021 dpsd_data_dx(0, ip, i) = (psd_1p[i] - psd_data(ip, i)) / diwc;
1022 }
1023 }
1024 }
1025}
1026
1027/*===========================================================================
1028 === Input: RWC
1029 ===========================================================================*/
1030
1031/* Workspace method: Doxygen documentation will be auto-generated */
1032void psdAbelBoutle12(Matrix& psd_data,
1033 Tensor3& dpsd_data_dx,
1034 const Vector& psd_size_grid,
1035 const Vector& pnd_agenda_input_t,
1036 const Matrix& pnd_agenda_input,
1037 const ArrayOfString& pnd_agenda_input_names,
1038 const ArrayOfString& dpnd_data_dx_names,
1039 const Numeric& scat_species_a,
1040 const Numeric& scat_species_b,
1041 const Numeric& t_min,
1042 const Numeric& t_max,
1043 const Index& picky,
1044 const Verbosity& verbosity) {
1045 psd_mgd_smm_common(psd_data,
1046 dpsd_data_dx,
1047 "Abel12",
1048 psd_size_grid,
1049 pnd_agenda_input_t,
1050 pnd_agenda_input,
1051 pnd_agenda_input_names,
1052 dpnd_data_dx_names,
1053 scat_species_a,
1054 scat_species_b,
1055 0,
1056 0,
1057 0,
1058 0,
1059 t_min,
1060 t_max,
1061 picky,
1062 verbosity);
1063}
1064
1065/* Workspace method: Doxygen documentation will be auto-generated */
1066void psdWangEtAl16(Matrix& psd_data,
1067 Tensor3& dpsd_data_dx,
1068 const Vector& psd_size_grid,
1069 const Vector& pnd_agenda_input_t,
1070 const Matrix& pnd_agenda_input,
1071 const ArrayOfString& pnd_agenda_input_names,
1072 const ArrayOfString& dpnd_data_dx_names,
1073 const Numeric& scat_species_a,
1074 const Numeric& scat_species_b,
1075 const Numeric& t_min,
1076 const Numeric& t_max,
1077 const Index& picky,
1078 const Verbosity& verbosity) {
1079 psd_mgd_smm_common(psd_data,
1080 dpsd_data_dx,
1081 "Wang16",
1082 psd_size_grid,
1083 pnd_agenda_input_t,
1084 pnd_agenda_input,
1085 pnd_agenda_input_names,
1086 dpnd_data_dx_names,
1087 scat_species_a,
1088 scat_species_b,
1089 0,
1090 0,
1091 0,
1092 0,
1093 t_min,
1094 t_max,
1095 picky,
1096 verbosity);
1097}
1098
1099/*===========================================================================
1100 === Input: GWC
1101 ===========================================================================*/
1102
1103/* Workspace method: Doxygen documentation will be auto-generated */
1104void psdFieldEtAl19(Matrix& psd_data,
1105 Tensor3& dpsd_data_dx,
1106 const Vector& psd_size_grid,
1107 const Vector& pnd_agenda_input_t,
1108 const Matrix& pnd_agenda_input,
1109 const ArrayOfString& pnd_agenda_input_names,
1110 const ArrayOfString& dpnd_data_dx_names,
1111 const Numeric& scat_species_a,
1112 const Numeric& scat_species_b,
1113 const Numeric& t_min,
1114 const Numeric& t_max,
1115 const Index& picky,
1116 const Verbosity& verbosity) {
1117 psd_mgd_smm_common(psd_data,
1118 dpsd_data_dx,
1119 "Field19",
1120 psd_size_grid,
1121 pnd_agenda_input_t,
1122 pnd_agenda_input,
1123 pnd_agenda_input_names,
1124 dpnd_data_dx_names,
1125 scat_species_a,
1126 scat_species_b,
1127 0,
1128 0,
1129 0,
1130 0,
1131 t_min,
1132 t_max,
1133 picky,
1134 verbosity);
1135}
1136
1137/*===========================================================================
1138 === Input: Atmospheric model PSDs
1139 ===========================================================================*/
1140
1141/* Workspace method: Doxygen documentation will be auto-generated */
1142void psdSeifertBeheng06(Matrix& psd_data,
1143 Tensor3& dpsd_data_dx,
1144 const Vector& psd_size_grid,
1145 const Vector& pnd_agenda_input_t,
1146 const Matrix& pnd_agenda_input,
1147 const ArrayOfString& pnd_agenda_input_names,
1148 const ArrayOfString& dpnd_data_dx_names,
1149 const String& hydrometeor_type,
1150 const Numeric& t_min,
1151 const Numeric& t_max,
1152 const Index& picky,
1153 const Verbosity&) {
1154 // Some sizes
1155 const Index nin = pnd_agenda_input_names.nelem();
1156 const Index ndx = dpnd_data_dx_names.nelem();
1157 const Index np = pnd_agenda_input.nrows();
1158 const Index nsi = psd_size_grid.nelem();
1159
1160 // Checks
1161 ARTS_USER_ERROR_IF (pnd_agenda_input.ncols() != nin,
1162 "Length of *pnd_agenda_input_names* and number of "
1163 "columns in *pnd_agenda_input* must be equal.");
1164 ARTS_USER_ERROR_IF (pnd_agenda_input.ncols() != 2,
1165 "*pnd_agenda_input* must have two columns"
1166 "(mass density and number density).");
1167
1168 ARTS_USER_ERROR_IF (ndx > 2,
1169 "*dpnd_data_dx_names* must have length <=2.");
1170
1171 // check name tags
1172 ArrayOfIndex input_idx = {-1, -1};
1173
1174 for (Index i = 0; i < nin; i++) {
1175 if ((Index)pnd_agenda_input_names[i].find("mass_density") != String::npos) {
1176 input_idx[0] = i; //mass density index
1177 } else if ((Index)pnd_agenda_input_names[i].find("number_density") !=
1178 String::npos) {
1179 input_idx[1] = i; //number density index
1180 }
1181 }
1182
1183 ARTS_USER_ERROR_IF (input_idx[0] == -1,
1184 "mass_density-tag not found ");
1185 ARTS_USER_ERROR_IF (input_idx[1] == -1,
1186 "number_density-tag not found ");
1187
1188 // look after jacobian tags
1189 ArrayOfIndex dpnd_data_dx_idx = {-1, -1};
1190
1191 for (Index i = 0; i < ndx; i++) {
1192 if ((Index)dpnd_data_dx_names[i].find("mass_density") != String::npos) {
1193 dpnd_data_dx_idx[0] = i; //mass density index
1194 } else if ((Index)dpnd_data_dx_names[i].find("number_density") !=
1195 String::npos) {
1196 dpnd_data_dx_idx[1] = i; //number density index
1197 }
1198 }
1199
1200 // Init psd_data and dpsd_data_dx with zeros
1201 psd_data.resize(np, nsi);
1202 psd_data = 0.0;
1203 if (ndx != 0) {
1204 dpsd_data_dx.resize(ndx, np, nsi);
1205 dpsd_data_dx = 0.0;
1206 } else {
1207 dpsd_data_dx.resize(0, 0, 0);
1208 }
1209
1210 for (Index ip = 0; ip < np; ip++) {
1211 // Extract the input variables
1212 Numeric WC = pnd_agenda_input(ip, input_idx[0]);
1213 Numeric N_tot = pnd_agenda_input(ip, input_idx[1]);
1214 Numeric t = pnd_agenda_input_t[ip];
1215
1216 // No calc needed if swc==0 and no jacobians requested.
1217 if ((WC == 0.) && (!ndx)) {
1218 continue;
1219 } // If here, we are ready with this point!
1220
1221 // Outside of [t_min,tmax]?
1222 if (t < t_min || t > t_max) {
1223 ARTS_USER_ERROR_IF (picky,
1224 "Method called with a temperature of ", t, " K.\n"
1225 "This is outside the specified allowed range: [ max(0.,", t_min,
1226 "), ", t_max, " ]")
1227 continue;
1228 // If here, we are ready with this point!
1229 }
1230
1231 // Negative swc?
1232 Numeric psd_weight = 1.0;
1233 if (WC < 0) {
1234 psd_weight = -1.0;
1235 WC *= -1.0;
1236 }
1237
1238 // Calculate PSD and derivatives
1239 Vector psd_1p(nsi);
1240 Matrix dpsd_1p(nsi, 2);
1241 if (WC > 0) {
1242 psd_SB06(psd_1p, dpsd_1p, psd_size_grid, N_tot, WC, hydrometeor_type);
1243
1244 for (Index i = 0; i < nsi; i++) {
1245 psd_data(ip, i) = psd_weight * psd_1p[i];
1246
1247 for (Index idx = 0; idx < dpnd_data_dx_idx.nelem(); idx++) {
1248 // with respect to WC
1249
1250 if (dpnd_data_dx_idx[idx] != -1) {
1251 dpsd_data_dx(dpnd_data_dx_idx[idx], ip, i) =
1252 psd_weight * dpsd_1p(i, idx);
1253 }
1254 }
1255 }
1256 }
1257 }
1258}
1259
1260/* Workspace method: Doxygen documentation will be auto-generated */
1261void psdMilbrandtYau05(Matrix& psd_data,
1262 Tensor3& dpsd_data_dx,
1263 const Vector& psd_size_grid,
1264 const Vector& pnd_agenda_input_t,
1265 const Matrix& pnd_agenda_input,
1266 const ArrayOfString& pnd_agenda_input_names,
1267 const ArrayOfString& dpnd_data_dx_names,
1268 const String& hydrometeor_type,
1269 const Numeric& t_min,
1270 const Numeric& t_max,
1271 const Index& picky,
1272 const Verbosity&) {
1273 // Some sizes
1274 const Index nin = pnd_agenda_input_names.nelem();
1275 const Index ndx = dpnd_data_dx_names.nelem();
1276 const Index np = pnd_agenda_input.nrows();
1277 const Index nsi = psd_size_grid.nelem();
1278
1279 // Checks
1280 ARTS_USER_ERROR_IF (pnd_agenda_input.ncols() != nin,
1281 "Length of *pnd_agenda_input_names* and number of "
1282 "columns in *pnd_agenda_input* must be equal.");
1283 ARTS_USER_ERROR_IF (pnd_agenda_input.ncols() != 2,
1284 "*pnd_agenda_input* must have two columns"
1285 "(mass density and number density).");
1286
1287 ARTS_USER_ERROR_IF (ndx > 2,
1288 "*dpnd_data_dx_names* must have length <=2.");
1289
1290 // check name tags
1291 ArrayOfIndex input_idx = {-1, -1};
1292
1293 for (Index i = 0; i < nin; i++) {
1294 if ((Index)pnd_agenda_input_names[i].find("mass_density") != String::npos) {
1295 input_idx[0] = i; //mass density index
1296 } else if ((Index)pnd_agenda_input_names[i].find("number_density") !=
1297 String::npos) {
1298 input_idx[1] = i; //number density index
1299 }
1300 }
1301
1302 ARTS_USER_ERROR_IF (input_idx[0] == -1,
1303 "mass_density-tag not found ");
1304 ARTS_USER_ERROR_IF (input_idx[1] == -1,
1305 "number_density-tag not found ");
1306
1307 // look after jacobian tags
1308 ArrayOfIndex dpnd_data_dx_idx = {-1, -1};
1309
1310 for (Index i = 0; i < ndx; i++) {
1311 if ((Index)dpnd_data_dx_names[i].find("mass_density") != String::npos) {
1312 dpnd_data_dx_idx[0] = i; //mass density index
1313 } else if ((Index)dpnd_data_dx_names[i].find("number_density") !=
1314 String::npos) {
1315 dpnd_data_dx_idx[1] = i; //number density index
1316 }
1317 }
1318
1319 // Init psd_data and dpsd_data_dx with zeros
1320 psd_data.resize(np, nsi);
1321 psd_data = 0.0;
1322 if (ndx != 0) {
1323 dpsd_data_dx.resize(ndx, np, nsi);
1324 dpsd_data_dx = 0.0;
1325 } else {
1326 dpsd_data_dx.resize(0, 0, 0);
1327 }
1328
1329 for (Index ip = 0; ip < np; ip++) {
1330 // Extract the input variables
1331 Numeric WC = pnd_agenda_input(ip, input_idx[0]);
1332 Numeric N_tot = pnd_agenda_input(ip, input_idx[1]);
1333 Numeric t = pnd_agenda_input_t[ip];
1334
1335 // No calc needed if wc==0 and no jacobians requested.
1336 if ((WC == 0.) && (!ndx)) {
1337 continue;
1338 } // If here, we are ready with this point!
1339
1340 // Outside of [t_min,tmax]?
1341 if (t < t_min || t > t_max) {
1342 ARTS_USER_ERROR_IF (picky,
1343 "Method called with a temperature of ", t, " K.\n"
1344 "This is outside the specified allowed range: [ max(0.,", t_min,
1345 "), ", t_max, " ]")
1346 continue;
1347 // If here, we are ready with this point!
1348 }
1349
1350 // Negative wc?
1351 Numeric psd_weight = 1.0;
1352 if (WC < 0) {
1353 psd_weight = -1.0;
1354 WC *= -1.0;
1355 }
1356
1357 // Calculate PSD and derivatives
1358 Vector psd_1p(nsi);
1359 Matrix dpsd_1p(nsi, 2);
1360 if (WC > 0) {
1361 psd_MY05(psd_1p, dpsd_1p, psd_size_grid, N_tot, WC, hydrometeor_type);
1362
1363 for (Index i = 0; i < nsi; i++) {
1364 psd_data(ip, i) = psd_weight * psd_1p[i];
1365
1366 for (Index idx = 0; idx < dpnd_data_dx_idx.nelem(); idx++) {
1367 // with respect to WC
1368
1369 if (dpnd_data_dx_idx[idx] != -1) {
1370 dpsd_data_dx(dpnd_data_dx_idx[idx], ip, i) =
1371 psd_weight * dpsd_1p(i, idx);
1372 }
1373 }
1374 }
1375 }
1376 }
1377}
array.h
This file contains the definition of Array.
arts.h
The global header file for ARTS.
auto_md.h
check_input.h
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
my_basic_string< char >
my_basic_string< char >::npos
static const Index npos
Define npos:
Definition mystring.h:192
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
psdModifiedGammaMassXmean
void psdModifiedGammaMassXmean(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &verbosity)
WORKSPACE METHOD: psdModifiedGammaMassXmean.
Definition m_psd.cc:513
psdMonoMass
void psdMonoMass(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &verbosity)
WORKSPACE METHOD: psdMonoMass.
Definition m_psd.cc:58
psdFieldEtAl07
void psdFieldEtAl07(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 String &regime, const Numeric &t_min, const Numeric &t_max, const Numeric &t_min_psd, const Numeric &t_max_psd, const Numeric &b_min, const Numeric &b_max, const Index &picky, const Verbosity &)
WORKSPACE METHOD: psdFieldEtAl07.
Definition m_psd.cc:826
psdModifiedGammaMassNtot
void psdModifiedGammaMassNtot(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &verbosity)
WORKSPACE METHOD: psdModifiedGammaMassNtot.
Definition m_psd.cc:437
psdModifiedGammaMass
void psdModifiedGammaMass(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &)
WORKSPACE METHOD: psdModifiedGammaMass.
Definition m_psd.cc:228
psdMcFarquaharHeymsfield97
void psdMcFarquaharHeymsfield97(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &t_min, const Numeric &t_max, const Numeric &t_min_psd, const Numeric &t_max_psd, const Index &picky, const Index &noisy, const Verbosity &)
WORKSPACE METHOD: psdMcFarquaharHeymsfield97.
Definition m_psd.cc:935
psdDelanoeEtAl14
void psdDelanoeEtAl14(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &iwc, const Numeric &n0, const Numeric &dm, const Numeric &rho, const Numeric &alpha, const Numeric &beta, const Numeric &t_min, const Numeric &t_max, const Numeric &dm_min, const Index &picky, const Verbosity &)
WORKSPACE METHOD: psdDelanoeEtAl14.
Definition m_psd.cc:632
psdModifiedGammaMassSingleMoment
void psdModifiedGammaMassSingleMoment(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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, const Numeric &n_b, const Numeric &mu, const Numeric &gamma, const Numeric &t_min, const Numeric &t_max, const Index &picky, const Verbosity &verbosity)
WORKSPACE METHOD: psdModifiedGammaMassSingleMoment.
Definition m_psd.cc:589
psdMonoDispersive
void psdMonoDispersive(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &verbosity)
WORKSPACE METHOD: psdMonoDispersive.
Definition m_psd.cc:30
psdModifiedGammaMassMeanParticleMass
void psdModifiedGammaMassMeanParticleMass(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &verbosity)
WORKSPACE METHOD: psdModifiedGammaMassMeanParticleMass.
Definition m_psd.cc:475
psdSeifertBeheng06
void psdSeifertBeheng06(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 String &hydrometeor_type, const Numeric &t_min, const Numeric &t_max, const Index &picky, const Verbosity &)
WORKSPACE METHOD: psdSeifertBeheng06.
Definition m_psd.cc:1142
psdModifiedGammaMassXmedian
void psdModifiedGammaMassXmedian(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &verbosity)
WORKSPACE METHOD: psdModifiedGammaMassXmedian.
Definition m_psd.cc:551
psdAbelBoutle12
void psdAbelBoutle12(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &t_min, const Numeric &t_max, const Index &picky, const Verbosity &verbosity)
WORKSPACE METHOD: psdAbelBoutle12.
Definition m_psd.cc:1032
psdWangEtAl16
void psdWangEtAl16(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &t_min, const Numeric &t_max, const Index &picky, const Verbosity &verbosity)
WORKSPACE METHOD: psdWangEtAl16.
Definition m_psd.cc:1066
psdFieldEtAl19
void psdFieldEtAl19(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &t_min, const Numeric &t_max, const Index &picky, const Verbosity &verbosity)
WORKSPACE METHOD: psdFieldEtAl19.
Definition m_psd.cc:1104
psdMilbrandtYau05
void psdMilbrandtYau05(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 String &hydrometeor_type, const Numeric &t_min, const Numeric &t_max, const Index &picky, const Verbosity &)
WORKSPACE METHOD: psdMilbrandtYau05.
Definition m_psd.cc:1261
psdModifiedGamma
void psdModifiedGamma(Matrix &psd_data, Tensor3 &dpsd_data_dx, 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 &n0, const Numeric &mu, const Numeric &la, const Numeric &ga, const Numeric &t_min, const Numeric &t_max, const Index &picky, const Verbosity &)
WORKSPACE METHOD: psdModifiedGamma.
Definition m_psd.cc:90
delanoe_shape_with_derivative
void delanoe_shape_with_derivative(VectorView psd, MatrixView jac_data, const Vector &x, const Numeric &alpha, const Numeric &beta)
! Shape functions for normalized PSD.
Definition math_funcs.cc:623
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:543
math_funcs.h
messages.h
Declarations having to do with the four output streams.
physics_funcs.h
This file contains declerations of functions of physical character.
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:871
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:1103
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:173
n0_from_t
Numeric n0_from_t(Numeric t)
Sets N0star based on temperature.
Definition psd.cc:1242
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:39
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:592
n0_from_iwc_dm
Numeric n0_from_iwc_dm(Numeric iwc, Numeric dm, Numeric rho)
Derives N0star from IWC and Dm.
Definition psd.cc:1234
dm_from_iwc_n0
Numeric dm_from_iwc_n0(Numeric iwc, Numeric n0, Numeric rho)
Derives Dm from IWC and N0star.
Definition psd.cc:1226
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:962
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:712
psd.h
Internal functions associated with size distributions.
START_OF_PSD_METHODS
#define START_OF_PSD_METHODS()
Definition psd.h:24