ARTS 2.5.9 (git: 825fa5f2)
microphysics.cc
Go to the documentation of this file.
1/* Copyright (C) 2011-2017 Jana Mendrok <jana.mendrok@gmail.com>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA.
17*/
18
27#include "microphysics.h"
28#include "arts_constants.h"
29#include "arts_conversions.h"
30
31inline constexpr Numeric PI=Constant::pi;
35
36
37/*===========================================================================
38 === External declarations
39 ===========================================================================*/
40#include <algorithm>
41#include <cmath>
42#include <ctime>
43#include <limits>
44#include <stdexcept>
45
46#include "arts.h"
47#include "check_input.h"
48#include "cloudbox.h"
49#include "lin_alg.h"
50#include "logic.h"
51#include "math_funcs.h"
52#include "mc_antenna.h"
53#include "messages.h"
54#include "physics_funcs.h"
55#include "ppath.h"
56#include "psd.h"
57#include "rng.h"
58#include "sorting.h"
59
60
62 ConstVectorView pfun)
63{
64 const Index n = sa_grid.nelem();
65
66 ARTS_ASSERT(abs(sa_grid[0]-0.0) < 1.0e-3);
67 ARTS_ASSERT(abs(sa_grid[n-1]-180.0) < 1.0e-3);
68 ARTS_ASSERT(pfun.nelem() == n);
69
70 Vector sa = sa_grid;
71 sa *= DEG2RAD;
72
73 // Sine and cosine of scattering angle
74 Vector sterm = sa;
75 transform(sterm, sin, sterm);
76 Vector cterm = sa;
77 transform(cterm, cos, cterm);
78
79 // Functions to integrate
80 Vector f1(n), f2(n);
81 for (Index i=0; i<n; ++i) {
82 f1[i] = sterm[i] * pfun[i];
83 f2[i] = cterm[i] * f1[i];
84 }
85
86 // We skip some 2, 4 and pi, as they all cancel in the end
87 const Numeric normfac = trapz(sa, f1);
88
89 if (normfac < 1e-12) {
90 return 0.0; // If very small scattering cross-section, set to zero
91 } else { // to avoid numerical issues
92 return trapz(sa, f2) / normfac;
93 }
94}
95
96
98 Numeric& b,
99 const Vector& x,
100 const Vector& mass,
101 const Numeric& x_fit_start,
102 const Numeric& x_fit_end) {
103 const Index nse = x.nelem();
104 ARTS_ASSERT(nse > 1);
105
106 ArrayOfIndex intarr_sort, intarr_unsort(0);
107 Vector x_unsorted(nse), m_unsorted(nse);
108 Vector q;
109 Index nsev = 0;
110
111 for (Index i = 0; i < nse; i++) {
112 if (std::isnan(x[i]))
113 throw runtime_error("NaN found in selected size grid data.");
114 if (std::isnan(mass[i]))
115 throw runtime_error("NaN found among particle mass data.");
116
117 if (x[i] >= x_fit_start && x[i] <= x_fit_end) {
118 x_unsorted[nsev] = x[i];
119 m_unsorted[nsev] = mass[i];
120 nsev += 1;
121 }
122 }
123
124 if (nsev < 2)
125 throw runtime_error(
126 "Less than two size points found in the range "
127 "[x_fit_start,x_fit_end]. It is then not possible "
128 "to determine the a and b parameters.");
129
130 get_sorted_indexes(intarr_sort, x_unsorted[Range(0, nsev)]);
131 Vector log_x(nsev), log_m(nsev);
132
133 for (Index i = 0; i < nsev; i++) {
134 log_x[i] = log(x_unsorted[intarr_sort[i]]);
135 log_m[i] = log(m_unsorted[intarr_sort[i]]);
136 }
137
138 linreg(q, log_x, log_m);
139 a = exp(q[0]);
140 b = q[1];
141}
142
The global header file for ARTS.
Constants of physical expressions as constexpr.
Common ARTS conversions.
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The range class.
Definition: matpackI.h:160
The Vector class.
Definition: matpackI.h:910
Internal cloudbox functions.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:82
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:548
Linear algebra functions.
Header file for logic.cc.
Numeric trapz(ConstVectorView x, ConstVectorView y)
trapz
Definition: math_funcs.cc:293
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1433
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
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:61
constexpr Numeric DEG2RAD
Definition: microphysics.cc:34
constexpr Numeric DENSITY_OF_ICE
Definition: microphysics.cc:32
constexpr Numeric PI
Definition: microphysics.cc:31
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:97
constexpr Numeric DENSITY_OF_WATER
Definition: microphysics.cc:33
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.
member functions of the Rng class and gsl_rng code
Contains sorting routines.
void get_sorted_indexes(ArrayOfIndex &sorted, const T &data)
get_sorted_indexes
Definition: sorting.h:57
#define a
#define b