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