ARTS 2.5.4 (git: 31ce4f0e)
m_tmatrix.cc
Go to the documentation of this file.
1/* Copyright (C) 2013 Oliver Lemke
2 *
3 * This program is free software: you can redistribute it and/or modify
4 * it under the terms of the GNU General Public License as published by
5 * the Free Software Foundation, either version 3 of the License, or
6 * (at your option) any 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, see <http://www.gnu.org/licenses/>.
15 *
16 */
17
26#include <cfloat>
27#include <cmath>
28#include <stdexcept>
29#include "arts_constants.h"
30#include "check_input.h"
31#include "logic.h"
32#include "math_funcs.h"
33#include "messages.h"
34#include "refraction.h"
35#include "special_interp.h"
36#include "tmatrix.h"
37
38inline constexpr Numeric PI=Constant::pi;
39
40/* Workspace method: Doxygen documentation will be auto-generated */
42 Numeric& diameter_aspect_area_max,
43 const String& shape,
44 const Numeric& diameter_volume_equ,
45 const Numeric& aspect_ratio,
46 const Verbosity&) {
47 const Numeric volume = (PI * pow(diameter_volume_equ, 3)) / 6.;
48
49 if (shape == "spheroidal") {
50 if (aspect_ratio < 1) // prolate spheroid
51 {
52 //non-rotational axis (perpendicular to a)
53 const Numeric b =
54 pow((3. * volume) / (4. * PI * pow(aspect_ratio, 2)), 1. / 3.);
55 diameter_max = 2. * b;
56 diameter_aspect_area_max = 2. * b;
57 } else // oblate spheriod
58 {
59 //rotational axis
60 const Numeric a = pow((3. * volume * aspect_ratio) / (4 * PI), 1. / 3.);
61 diameter_max = 2. * a;
62 diameter_aspect_area_max = 2. * a;
63 }
64 }
65
66 else if (shape == "cylindrical") {
67 //aspect_ratio=D/L
68 const Numeric D = pow((volume * 4 * aspect_ratio) / PI, 1. / 3.);
69 const Numeric L = D / aspect_ratio;
70 diameter_max = pow(pow(D, 2) + pow(L, 2), 1. / 2.);
71 diameter_aspect_area_max = max(D, pow(4 / PI * D * L, 1. / 2.));
72 }
73
74 else {
75 ostringstream os;
76 os << "Unknown particle shape: " << shape << "\n"
77 << "Must be spheroidal or cylindrical";
78 throw runtime_error(os.str());
79 }
80}
81
82/* Workspace method: Doxygen documentation will be auto-generated */
84 Numeric& volume,
85 const String& shape,
86 const Numeric& diameter_max,
87 const Numeric& aspect_ratio,
88 const Verbosity&) {
89 if (shape == "spheroidal") {
90 if (aspect_ratio < 1) // prolate spheroid
91 {
92 const Numeric b = diameter_max / 2.;
93 volume = (pow(b, 3.) * 4. * PI * pow(aspect_ratio, 2.)) / 3.;
94 } else // oblate spheriod
95 {
96 const Numeric a = diameter_max / 2.;
97 volume = (pow(a, 3.) * 4. * PI) / (3. * aspect_ratio);
98 }
99 }
100
101 else if (shape == "cylindrical") {
102 volume = pow(diameter_max / pow((pow(aspect_ratio, 2.) + 1.), 1. / 2.), 3) *
103 pow(aspect_ratio, 2.) * PI / 4.;
104 }
105
106 else {
107 ostringstream os;
108 os << "Unknown particle shape: " << shape << "\n"
109 << "Must be spheroidal or cylindrical";
110 throw runtime_error(os.str());
111 }
112
113 diameter_volume_equ = pow((6. * volume) / PI, 1. / 3.);
114}
115
116/* Workspace method: Doxygen documentation will be auto-generated */
118 ScatteringMetaData& scat_meta_single,
119 const GriddedField3& complex_refr_index,
120 const String& shape,
121 const Numeric& diameter_volume_equ,
122 const Numeric& aspect_ratio,
123 const Numeric& mass,
124 const String& ptype,
125 const Vector& data_f_grid,
126 const Vector& data_t_grid,
127 const Vector& data_za_grid,
128 const Vector& data_aa_grid,
129 const Numeric& precision,
130 const String& cri_source,
131 const Index& ndgs,
132 const Index& robust,
133 const Index& quiet,
134 const Verbosity& verbosity) {
135 // Get internal coding for ptype
136 scat_data_single.ptype = PTypeFromString(ptype);
137
138 // Set description
139 {
140 ostringstream os;
141 os << "T-matrix calculation for a " << shape << " particle, with "
142 << "diameter_volume_equ = " << 1e6 * diameter_volume_equ << "um and "
143 << "aspect ratio = " << aspect_ratio << ".";
144 scat_data_single.description = os.str();
145 }
146
147 // Add grids to scat_data_single
148 //
149 scat_data_single.f_grid = data_f_grid;
150 scat_data_single.T_grid = data_t_grid;
151
152 if (scat_data_single.ptype == PTYPE_TOTAL_RND) {
153 // tmatrix random orient requires equidistant angular grid. checking here
154 // that given data_za_grid fulfills this requirement
155 if (!(is_same_within_epsilon(data_za_grid[0], 0., 2 * DBL_EPSILON) &&
156 is_same_within_epsilon(last(data_za_grid), 180., 2 * DBL_EPSILON))) {
157 ostringstream os;
158 os << "Zenith angle (=scattering angle) grid needs to include\n"
159 << "0 deg and 180 deg as first and last grid points, respectively.\n"
160 << "At least one of them does not fit.";
161 throw runtime_error(os.str());
162 }
163 Index nza = data_za_grid.nelem();
164 Numeric dza = 180. / ((Numeric)nza - 1.);
165 for (Index iza = 1; iza < nza; iza++) {
167 data_za_grid[iza], (Numeric)iza * dza, 2 * DBL_EPSILON))) {
168 ostringstream os;
169 os << "Input zenith angle grid *data_za_grid* is required to be\n"
170 << "equidistant for randomly oriented particles, but it is not.";
171 throw runtime_error(os.str());
172 }
173 }
174 }
175 scat_data_single.za_grid = data_za_grid;
176
177 if (scat_data_single.ptype == PTYPE_TOTAL_RND) {
178 // in case of random orientation, azimuth grid should be empty. We just
179 // set that here, ignoring whatever is in data_aa_grid.
180 Vector empty_grid(0);
181 scat_data_single.aa_grid = empty_grid;
182 } else {
183 // For azimuthally-random oriented particles, the azimuth angle grid must cover
184 // 0-180 degrees.
185 if (scat_data_single.ptype == PTYPE_AZIMUTH_RND &&
186 data_aa_grid.nelem() == 0) {
187 ostringstream os;
188 os << "For ptype = \"azimuthally_random\""
189 << " the azimuth angle grid can not be empty.";
190 throw runtime_error(os.str());
191 }
192 if (scat_data_single.ptype == PTYPE_AZIMUTH_RND && data_aa_grid[0] != 0.) {
193 ostringstream os;
194 os << "For ptype = \"azimuthally_random\""
195 << " the first value of the aa grid must be 0.";
196 throw runtime_error(os.str());
197 }
198
199 if (scat_data_single.ptype == PTYPE_AZIMUTH_RND &&
200 last(data_aa_grid) != 180.) {
201 ostringstream os;
202 os << "For ptype = \"azimuthally_random\""
203 << " the last value of the aa grid must be 180.";
204 throw runtime_error(os.str());
205 }
206
207 scat_data_single.aa_grid = data_aa_grid;
208 }
209
210 // Index coding for shape
211 Index np;
212 Numeric ar = aspect_ratio;
213 if (shape == "spheroidal") {
214 np = -1;
215 if (aspect_ratio == 1)
216 /* do not throw error, but slightly increase aspect ratio to value
217 * recommended by original tmatrix code such that numerical issues are
218 * avoided.
219 throw runtime_error( "For spheroidal particles, the aspect ratio "
220 "is not allowed to be exactly 1 (due to "
221 "numerical problems)." );
222*/
223 ar += 1e-6;
224 } else if (shape == "cylindrical") {
225 np = -2;
226 } else {
227 ostringstream os;
228 os << "Unknown particle shape: " << shape << "\n"
229 << "Must be \"spheroidal\" or \"cylindrical\".";
230 throw runtime_error(os.str());
231 }
232
233 // Interpolate refractive index to relevant grids
234 //
235 const Index nf = data_f_grid.nelem();
236 const Index nt = data_t_grid.nelem();
237 //
238 Tensor3 ncomp(nf, nt, 2);
239 complex_n_interp(ncomp(joker, joker, 0),
240 ncomp(joker, joker, 1),
241 complex_refr_index,
242 "complex_refr_index",
243 data_f_grid,
244 data_t_grid);
245
246 // Run T-matrix and we are ready (T-matrix takes size as volume equiv radius(!) )
247 calcSingleScatteringDataProperties(scat_data_single,
248 ncomp(joker, joker, 0),
249 ncomp(joker, joker, 1),
250 0.5 * diameter_volume_equ,
251 np,
252 ar,
253 precision,
254 ndgs,
255 robust,
256 quiet);
257
258 // Meta data
259 scat_meta_single.description =
260 "Meta data for associated file with single scattering data.";
261 scat_meta_single.source =
262 "ARTS interface to T-matrix code by Mishchenko et al.";
263 scat_meta_single.refr_index = cri_source;
264 //
265 Numeric diameter_max, area_max;
267 area_max,
268 shape,
269 diameter_volume_equ,
270 aspect_ratio,
271 verbosity);
272 //
273 scat_meta_single.mass = mass;
274 scat_meta_single.diameter_max = diameter_max;
275 scat_meta_single.diameter_volume_equ = diameter_volume_equ;
276 scat_meta_single.diameter_area_equ_aerodynamical = area_max;
277}
278
279void TMatrixTest(const Verbosity& verbosity) {
280 tmatrix_tmd_test(verbosity);
281 tmatrix_ampld_test(verbosity);
282 calc_ssp_random_test(verbosity);
283 calc_ssp_fixed_test(verbosity);
284}
Constants of physical expressions as constexpr.
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The Tensor3 class.
Definition: matpackIII.h:346
The Vector class.
Definition: matpackI.h:899
#define max(a, b)
#define precision
Definition: logic.cc:46
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:354
Header file for logic.cc.
void diameter_maxFromDiameter_volume_equ(Numeric &diameter_max, Numeric &diameter_aspect_area_max, const String &shape, const Numeric &diameter_volume_equ, const Numeric &aspect_ratio, const Verbosity &)
WORKSPACE METHOD: diameter_maxFromDiameter_volume_equ.
Definition: m_tmatrix.cc:41
void scat_data_singleTmatrix(SingleScatteringData &scat_data_single, ScatteringMetaData &scat_meta_single, const GriddedField3 &complex_refr_index, const String &shape, const Numeric &diameter_volume_equ, const Numeric &aspect_ratio, const Numeric &mass, const String &ptype, const Vector &data_f_grid, const Vector &data_t_grid, const Vector &data_za_grid, const Vector &data_aa_grid, const Numeric &precision, const String &cri_source, const Index &ndgs, const Index &robust, const Index &quiet, const Verbosity &verbosity)
WORKSPACE METHOD: scat_data_singleTmatrix.
Definition: m_tmatrix.cc:117
void diameter_volume_equFromDiameter_max(Numeric &diameter_volume_equ, Numeric &volume, const String &shape, const Numeric &diameter_max, const Numeric &aspect_ratio, const Verbosity &)
WORKSPACE METHOD: diameter_volume_equFromDiameter_max.
Definition: m_tmatrix.cc:83
constexpr Numeric PI
Definition: m_tmatrix.cc:38
void TMatrixTest(const Verbosity &verbosity)
WORKSPACE METHOD: TMatrixTest.
Definition: m_tmatrix.cc:279
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:161
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
const Joker joker
Declarations having to do with the four output streams.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric e
Elementary charge convenience name [C].
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.
PType PTypeFromString(const String &ptype_string)
Convert ptype name to enum value.
@ PTYPE_AZIMUTH_RND
Definition: optproperties.h:54
@ PTYPE_TOTAL_RND
Definition: optproperties.h:55
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
Refraction functions.
void complex_n_interp(MatrixView n_real, MatrixView n_imag, const GriddedField3 &complex_n, const String &varname, ConstVectorView f_grid, ConstVectorView t_grid)
General function for interpolating data of complex n type.
Header file for special_interp.cc.
Numeric diameter_area_equ_aerodynamical
Numeric diameter_volume_equ
void calc_ssp_fixed_test(const Verbosity &verbosity)
Single scattering properties calculation for particles with fixed orientation.
Definition: tmatrix.cc:1508
void calcSingleScatteringDataProperties(SingleScatteringData &ssd, ConstMatrixView ref_index_real, ConstMatrixView ref_index_imag, const Numeric equiv_radius, const Index np, const Numeric aspect_ratio, const Numeric precision, const Index ndgs, const Index robust, const Index quiet)
Calculate SingleScatteringData properties.
Definition: tmatrix.cc:980
void tmatrix_tmd_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1362
void calc_ssp_random_test(const Verbosity &verbosity)
Single scattering properties calculation for randomly oriented particles.
Definition: tmatrix.cc:1455
void tmatrix_ampld_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1298
Declarations for the T-Matrix interface.
#define a
#define b