ARTS 2.5.4 (git: 7d04b88e)
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
29extern const Numeric PI;
30extern const Numeric DENSITY_OF_ICE;
31extern const Numeric DENSITY_OF_WATER;
32extern const Numeric DEG2RAD;
33
34
35/*===========================================================================
36 === External declarations
37 ===========================================================================*/
38#include <algorithm>
39#include <cmath>
40#include <ctime>
41#include <limits>
42#include <stdexcept>
43
44#include "arts.h"
45#include "check_input.h"
46#include "cloudbox.h"
47#include "lin_alg.h"
48#include "logic.h"
49#include "math_funcs.h"
50#include "mc_antenna.h"
51#include "messages.h"
52#include "physics_funcs.h"
53#include "ppath.h"
54#include "psd.h"
55#include "rng.h"
56#include "sorting.h"
57
58
60 ConstVectorView pfun)
61{
62 const Index n = sa_grid.nelem();
63
64 ARTS_ASSERT(abs(sa_grid[0]-0.0) < 1.0e-3);
65 ARTS_ASSERT(abs(sa_grid[n-1]-180.0) < 1.0e-3);
66 ARTS_ASSERT(pfun.nelem() == n);
67
68 Vector sa = sa_grid;
69 sa *= DEG2RAD;
70
71 // Sine and cosine of scattering angle
72 Vector sterm = sa;
73 transform(sterm, sin, sterm);
74 Vector cterm = sa;
75 transform(cterm, cos, cterm);
76
77 // Functions to integrate
78 Vector f1(n), f2(n);
79 for (Index i=0; i<n; ++i) {
80 f1[i] = sterm[i] * pfun[i];
81 f2[i] = cterm[i] * f1[i];
82 }
83
84 // We skip some 2, 4 and pi, as they all cancel in the end
85 const Numeric normfac = trapz(sa, f1);
86
87 if (normfac < 1e-12) {
88 return 0.0; // If very small scattering cross-section, set to zero
89 } else { // to avoid numerical issues
90 return trapz(sa, f2) / normfac;
91 }
92}
93
94
96 Numeric& b,
97 const Vector& x,
98 const Vector& mass,
99 const Numeric& x_fit_start,
100 const Numeric& x_fit_end) {
101 const Index nse = x.nelem();
102 ARTS_ASSERT(nse > 1);
103
104 ArrayOfIndex intarr_sort, intarr_unsort(0);
105 Vector x_unsorted(nse), m_unsorted(nse);
106 Vector q;
107 Index nsev = 0;
108
109 for (Index i = 0; i < nse; i++) {
110 if (std::isnan(x[i]))
111 throw runtime_error("NaN found in selected size grid data.");
112 if (std::isnan(mass[i]))
113 throw runtime_error("NaN found among particle mass data.");
114
115 if (x[i] >= x_fit_start && x[i] <= x_fit_end) {
116 x_unsorted[nsev] = x[i];
117 m_unsorted[nsev] = mass[i];
118 nsev += 1;
119 }
120 }
121
122 if (nsev < 2)
123 throw runtime_error(
124 "Less than two size points found in the range "
125 "[x_fit_start,x_fit_end]. It is then not possible "
126 "to determine the a and b parameters.");
127
128 get_sorted_indexes(intarr_sort, x_unsorted[Range(0, nsev)]);
129 Vector log_x(nsev), log_m(nsev);
130
131 for (Index i = 0; i < nsev; i++) {
132 log_x[i] = log(x_unsorted[intarr_sort[i]]);
133 log_m[i] = log(m_unsorted[intarr_sort[i]]);
134 }
135
136 linreg(q, log_x, log_m);
137 a = exp(q[0]);
138 b = q[1];
139}
140
The global header file for ARTS.
A constant view of a Vector.
Definition: matpackI.h:530
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:554
The range class.
Definition: matpackI.h:177
The Vector class.
Definition: matpackI.h:922
Internal cloudbox functions.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define abs(x)
#define q
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:1087
Linear algebra functions.
Header file for logic.cc.
Numeric trapz(ConstVectorView x, ConstVectorView y)
trapz
Definition: math_funcs.cc:291
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:1484
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:59
const Numeric DENSITY_OF_WATER
const Numeric PI
const Numeric DEG2RAD
const Numeric DENSITY_OF_ICE
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:95
Internal functions for microphysics calculations (size distributions etc.)
Vector mass(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species, const SpeciesIsotopologueRatios &ir) ARTS_NOEXCEPT
Returns a mass vector for this model's main calculations.
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:53
This file contains declerations of functions of physical character.
Propagation path structure and functions.
#define a
#define b
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