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