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