ARTS 2.5.11 (git: 6827797f)
microphysics.cc
Go to the documentation of this file.
1
9#include "microphysics.h"
10#include "arts_constants.h"
11#include "arts_conversions.h"
12
13inline constexpr Numeric PI=Constant::pi;
16inline constexpr Numeric DEG2RAD=Conversion::deg2rad(1);
17
18
19/*===========================================================================
20 === External declarations
21 ===========================================================================*/
22#include <algorithm>
23#include <cmath>
24#include <ctime>
25#include <limits>
26#include <stdexcept>
27
28#include "arts.h"
29#include "check_input.h"
30#include "cloudbox.h"
31#include "lin_alg.h"
32#include "logic.h"
33#include "math_funcs.h"
34#include "mc_antenna.h"
35#include "messages.h"
36#include "physics_funcs.h"
37#include "ppath.h"
38#include "psd.h"
39#include "rng.h"
40#include "sorting.h"
41
42
43Numeric asymmetry_parameter(ConstVectorView sa_grid,
44 ConstVectorView pfun)
45{
46 const Index n = sa_grid.nelem();
47
48 ARTS_ASSERT(abs(sa_grid[0]-0.0) < 1.0e-3);
49 ARTS_ASSERT(abs(sa_grid[n-1]-180.0) < 1.0e-3);
50 ARTS_ASSERT(pfun.nelem() == n);
51
52 Vector sa{sa_grid};
53 sa *= DEG2RAD;
54
55 // Sine and cosine of scattering angle
56 Vector sterm = sa;
57 transform(sterm, sin, sterm);
58 Vector cterm = sa;
59 transform(cterm, cos, cterm);
60
61 // Functions to integrate
62 Vector f1(n), f2(n);
63 for (Index i=0; i<n; ++i) {
64 f1[i] = sterm[i] * pfun[i];
65 f2[i] = cterm[i] * f1[i];
66 }
67
68 // We skip some 2, 4 and pi, as they all cancel in the end
69 const Numeric normfac = trapz(sa, f1);
70
71 if (normfac < 1e-12) {
72 return 0.0; // If very small scattering cross-section, set to zero
73 } else { // to avoid numerical issues
74 return trapz(sa, f2) / normfac;
75 }
76}
77
78
80 Numeric& b,
81 const Vector& x,
82 const Vector& mass,
83 const Numeric& x_fit_start,
84 const Numeric& x_fit_end) {
85 const Index nse = x.nelem();
86 ARTS_ASSERT(nse > 1);
87
88 ArrayOfIndex intarr_sort, intarr_unsort(0);
89 Vector x_unsorted(nse), m_unsorted(nse);
90 Vector q;
91 Index nsev = 0;
92
93 for (Index i = 0; i < nse; i++) {
94 if (std::isnan(x[i]))
95 throw runtime_error("NaN found in selected size grid data.");
96 if (std::isnan(mass[i]))
97 throw runtime_error("NaN found among particle mass data.");
98
99 if (x[i] >= x_fit_start && x[i] <= x_fit_end) {
100 x_unsorted[nsev] = x[i];
101 m_unsorted[nsev] = mass[i];
102 nsev += 1;
103 }
104 }
105
106 if (nsev < 2)
107 throw runtime_error(
108 "Less than two size points found in the range "
109 "[x_fit_start,x_fit_end]. It is then not possible "
110 "to determine the a and b parameters.");
111
112 get_sorted_indexes(intarr_sort, x_unsorted[Range(0, nsev)]);
113 Vector log_x(nsev), log_m(nsev);
114
115 for (Index i = 0; i < nsev; i++) {
116 log_x[i] = log(x_unsorted[intarr_sort[i]]);
117 log_m[i] = log(m_unsorted[intarr_sort[i]]);
118 }
119
120 linreg(q, log_x, log_m);
121 a = exp(q[0]);
122 b = q[1];
123}
124
The global header file for ARTS.
Constants of physical expressions as constexpr.
Common ARTS conversions.
Internal cloudbox functions.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:84
Numeric trapz(ConstVectorView x, ConstVectorView y)
trapz
Definition: math_funcs.cc:274
Workspace functions for the solution of cloud-box radiative transfer by Monte Carlo methods....
Declarations having to do with the four output streams.
Numeric asymmetry_parameter(ConstVectorView sa_grid, ConstVectorView pfun)
asymmetry_parameter
Definition: microphysics.cc:43
constexpr Numeric DEG2RAD
Definition: microphysics.cc:16
constexpr Numeric DENSITY_OF_ICE
Definition: microphysics.cc:14
constexpr Numeric PI
Definition: microphysics.cc:13
void derive_scat_species_a_and_b(Numeric &a, Numeric &b, const Vector &x, const Vector &mass, const Numeric &x_fit_start, const Numeric &x_fit_end)
Definition: microphysics.cc:79
constexpr Numeric DENSITY_OF_WATER
Definition: microphysics.cc:15
Internal functions for microphysics calculations (size distributions etc.)
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric density_of_ice_at_0c
Global constant, Density of water ice at 0C [kg/m3] source: http://en.wikipedia.org/wiki/Ice.
constexpr Numeric denity_of_water_at_4c
Global constant, Density of liquid water +4C [kg/m3] source: http://en.wikipedia.org/wiki/Water.
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
This file contains declerations of functions of physical character.
Propagation path structure and functions.
Internal functions associated with size distributions.
Contains sorting routines.
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition: sorting.h:39
#define a
#define b