ARTS 2.5.11 (git: 725533f0)
m_tmatrix.cc
Go to the documentation of this file.
1
9#include <cfloat>
10#include <cmath>
11#include <stdexcept>
12#include "arts_constants.h"
13#include "check_input.h"
14#include "logic.h"
15#include "math_funcs.h"
16#include "messages.h"
17#include "refraction.h"
18#include "special_interp.h"
19#include "tmatrix.h"
20
21inline constexpr Numeric PI=Constant::pi;
22
23/* Workspace method: Doxygen documentation will be auto-generated */
24void diameter_maxFromDiameter_volume_equ(Numeric& diameter_max,
25 Numeric& diameter_aspect_area_max,
26 const String& shape,
27 const Numeric& diameter_volume_equ,
28 const Numeric& aspect_ratio,
29 const Verbosity&) {
30 const Numeric volume = (PI * pow(diameter_volume_equ, 3)) / 6.;
31
32 if (shape == "spheroidal") {
33 if (aspect_ratio < 1) // prolate spheroid
34 {
35 //non-rotational axis (perpendicular to a)
36 const Numeric b =
37 pow((3. * volume) / (4. * PI * pow(aspect_ratio, 2)), 1. / 3.);
38 diameter_max = 2. * b;
39 diameter_aspect_area_max = 2. * b;
40 } else // oblate spheriod
41 {
42 //rotational axis
43 const Numeric a = pow((3. * volume * aspect_ratio) / (4 * PI), 1. / 3.);
44 diameter_max = 2. * a;
45 diameter_aspect_area_max = 2. * a;
46 }
47 }
48
49 else if (shape == "cylindrical") {
50 //aspect_ratio=D/L
51 const Numeric D = pow((volume * 4 * aspect_ratio) / PI, 1. / 3.);
52 const Numeric L = D / aspect_ratio;
53 diameter_max = pow(pow(D, 2) + pow(L, 2), 1. / 2.);
54 diameter_aspect_area_max = max(D, pow(4 / PI * D * L, 1. / 2.));
55 }
56
57 else {
58 ostringstream os;
59 os << "Unknown particle shape: " << shape << "\n"
60 << "Must be spheroidal or cylindrical";
61 throw runtime_error(os.str());
62 }
63}
64
65/* Workspace method: Doxygen documentation will be auto-generated */
66void diameter_volume_equFromDiameter_max(Numeric& diameter_volume_equ,
67 Numeric& volume,
68 const String& shape,
69 const Numeric& diameter_max,
70 const Numeric& aspect_ratio,
71 const Verbosity&) {
72 if (shape == "spheroidal") {
73 if (aspect_ratio < 1) // prolate spheroid
74 {
75 const Numeric b = diameter_max / 2.;
76 volume = (pow(b, 3.) * 4. * PI * pow(aspect_ratio, 2.)) / 3.;
77 } else // oblate spheriod
78 {
79 const Numeric a = diameter_max / 2.;
80 volume = (pow(a, 3.) * 4. * PI) / (3. * aspect_ratio);
81 }
82 }
83
84 else if (shape == "cylindrical") {
85 volume = pow(diameter_max / pow((pow(aspect_ratio, 2.) + 1.), 1. / 2.), 3) *
86 pow(aspect_ratio, 2.) * PI / 4.;
87 }
88
89 else {
90 ostringstream os;
91 os << "Unknown particle shape: " << shape << "\n"
92 << "Must be spheroidal or cylindrical";
93 throw runtime_error(os.str());
94 }
95
96 diameter_volume_equ = pow((6. * volume) / PI, 1. / 3.);
97}
98
99/* Workspace method: Doxygen documentation will be auto-generated */
101 ScatteringMetaData& scat_meta_single,
102 const GriddedField3& complex_refr_index,
103 const String& shape,
104 const Numeric& diameter_volume_equ,
105 const Numeric& aspect_ratio,
106 const Numeric& mass,
107 const String& ptype,
108 const Vector& data_f_grid,
109 const Vector& data_t_grid,
110 const Vector& data_za_grid,
111 const Vector& data_aa_grid,
112 const Numeric& precision,
113 const String& cri_source,
114 const Index& ndgs,
115 const Index& robust,
116 const Index& quiet,
117 const Verbosity& verbosity) {
118 // Get internal coding for ptype
119 scat_data_single.ptype = PTypeFromString(ptype);
120
121 // Set description
122 {
123 ostringstream os;
124 os << "T-matrix calculation for a " << shape << " particle, with "
125 << "diameter_volume_equ = " << 1e6 * diameter_volume_equ << "um and "
126 << "aspect ratio = " << aspect_ratio << ".";
127 scat_data_single.description = os.str();
128 }
129
130 // Add grids to scat_data_single
131 //
132 scat_data_single.f_grid = data_f_grid;
133 scat_data_single.T_grid = data_t_grid;
134
135 if (scat_data_single.ptype == PTYPE_TOTAL_RND) {
136 // tmatrix random orient requires equidistant angular grid. checking here
137 // that given data_za_grid fulfills this requirement
138 if (!(is_same_within_epsilon(data_za_grid[0], 0., 2 * DBL_EPSILON) &&
139 is_same_within_epsilon(last(data_za_grid), 180., 2 * DBL_EPSILON))) {
140 ostringstream os;
141 os << "Zenith angle (=scattering angle) grid needs to include\n"
142 << "0 deg and 180 deg as first and last grid points, respectively.\n"
143 << "At least one of them does not fit.";
144 throw runtime_error(os.str());
145 }
146 Index nza = data_za_grid.nelem();
147 Numeric dza = 180. / ((Numeric)nza - 1.);
148 for (Index iza = 1; iza < nza; iza++) {
149 if (!(is_same_within_epsilon(
150 data_za_grid[iza], (Numeric)iza * dza, 2 * DBL_EPSILON))) {
151 ostringstream os;
152 os << "Input zenith angle grid *data_za_grid* is required to be\n"
153 << "equidistant for randomly oriented particles, but it is not.";
154 throw runtime_error(os.str());
155 }
156 }
157 }
158 scat_data_single.za_grid = data_za_grid;
159
160 if (scat_data_single.ptype == PTYPE_TOTAL_RND) {
161 // in case of random orientation, azimuth grid should be empty. We just
162 // set that here, ignoring whatever is in data_aa_grid.
163 Vector empty_grid(0);
164 scat_data_single.aa_grid = empty_grid;
165 } else {
166 // For azimuthally-random oriented particles, the azimuth angle grid must cover
167 // 0-180 degrees.
168 if (scat_data_single.ptype == PTYPE_AZIMUTH_RND &&
169 data_aa_grid.nelem() == 0) {
170 ostringstream os;
171 os << "For ptype = \"azimuthally_random\""
172 << " the azimuth angle grid can not be empty.";
173 throw runtime_error(os.str());
174 }
175 if (scat_data_single.ptype == PTYPE_AZIMUTH_RND && data_aa_grid[0] != 0.) {
176 ostringstream os;
177 os << "For ptype = \"azimuthally_random\""
178 << " the first value of the aa grid must be 0.";
179 throw runtime_error(os.str());
180 }
181
182 if (scat_data_single.ptype == PTYPE_AZIMUTH_RND &&
183 last(data_aa_grid) != 180.) {
184 ostringstream os;
185 os << "For ptype = \"azimuthally_random\""
186 << " the last value of the aa grid must be 180.";
187 throw runtime_error(os.str());
188 }
189
190 scat_data_single.aa_grid = data_aa_grid;
191 }
192
193 // Index coding for shape
194 Index np;
195 Numeric ar = aspect_ratio;
196 if (shape == "spheroidal") {
197 np = -1;
198 if (aspect_ratio == 1)
199 /* do not throw error, but slightly increase aspect ratio to value
200 * recommended by original tmatrix code such that numerical issues are
201 * avoided.
202 throw runtime_error( "For spheroidal particles, the aspect ratio "
203 "is not allowed to be exactly 1 (due to "
204 "numerical problems)." );
205*/
206 ar += 1e-6;
207 } else if (shape == "cylindrical") {
208 np = -2;
209 } else {
210 ostringstream os;
211 os << "Unknown particle shape: " << shape << "\n"
212 << "Must be \"spheroidal\" or \"cylindrical\".";
213 throw runtime_error(os.str());
214 }
215
216 // Interpolate refractive index to relevant grids
217 //
218 const Index nf = data_f_grid.nelem();
219 const Index nt = data_t_grid.nelem();
220 //
221 Tensor3 ncomp(nf, nt, 2);
222 complex_n_interp(ncomp(joker, joker, 0),
223 ncomp(joker, joker, 1),
224 complex_refr_index,
225 "complex_refr_index",
226 data_f_grid,
227 data_t_grid);
228
229 // Run T-matrix and we are ready (T-matrix takes size as volume equiv radius(!) )
230 calcSingleScatteringDataProperties(scat_data_single,
231 ncomp(joker, joker, 0),
232 ncomp(joker, joker, 1),
233 0.5 * diameter_volume_equ,
234 np,
235 ar,
236 precision,
237 ndgs,
238 robust,
239 quiet);
240
241 // Meta data
242 scat_meta_single.description =
243 "Meta data for associated file with single scattering data.";
244 scat_meta_single.source =
245 "ARTS interface to T-matrix code by Mishchenko et al.";
246 scat_meta_single.refr_index = cri_source;
247 //
248 Numeric diameter_max, area_max;
250 area_max,
251 shape,
252 diameter_volume_equ,
253 aspect_ratio,
254 verbosity);
255 //
256 scat_meta_single.mass = mass;
257 scat_meta_single.diameter_max = diameter_max;
258 scat_meta_single.diameter_volume_equ = diameter_volume_equ;
259 scat_meta_single.diameter_area_equ_aerodynamical = area_max;
260}
261
262void TMatrixTest(const Verbosity& verbosity) {
263 tmatrix_tmd_test(verbosity);
264 tmatrix_ampld_test(verbosity);
265 calc_ssp_random_test(verbosity);
266 calc_ssp_fixed_test(verbosity);
267}
base max(const Array< base > &x)
Max function.
Definition array.h:128
Constants of physical expressions as constexpr.
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:24
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:100
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:66
constexpr Numeric PI
Definition m_tmatrix.cc:21
void TMatrixTest(const Verbosity &verbosity)
WORKSPACE METHOD: TMatrixTest.
Definition m_tmatrix.cc:262
Numeric last(ConstVectorView x)
last
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:
PType PTypeFromString(const String &ptype_string)
Convert ptype name to enum value.
@ PTYPE_AZIMUTH_RND
@ PTYPE_TOTAL_RND
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:1491
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:963
void tmatrix_tmd_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition tmatrix.cc:1345
void calc_ssp_random_test(const Verbosity &verbosity)
Single scattering properties calculation for randomly oriented particles.
Definition tmatrix.cc:1438
void tmatrix_ampld_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition tmatrix.cc:1281
Declarations for the T-Matrix interface.
#define a
#define b