ARTS  2.3.1285(git:92a29ea9-dirty)
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 
37 extern 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 */
82 void diameter_volume_equFromDiameter_max(Numeric& diameter_volume_equ,
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 
278 void 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 }
void tmatrix_tmd_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1361
Declarations for the T-Matrix interface.
void TMatrixTest(const Verbosity &verbosity)
WORKSPACE METHOD: TMatrixTest.
Definition: m_tmatrix.cc:278
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
#define precision
Definition: logic.cc:43
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
Declarations having to do with the four output streams.
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
The Vector class.
Definition: matpackI.h:860
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:165
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
Numeric diameter_volume_equ
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
The Tensor3 class.
Definition: matpackIII.h:339
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
_CS_string_type str() const
Definition: sstream.h:491
void calc_ssp_fixed_test(const Verbosity &verbosity)
Single scattering properties calculation for particles with fixed orientation.
Definition: tmatrix.cc:1507
PType PTypeFromString(const String &ptype_string)
Convert ptype name to enum value.
const Numeric PI
void calc_ssp_random_test(const Verbosity &verbosity)
Single scattering properties calculation for randomly oriented particles.
Definition: tmatrix.cc:1454
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Header file for special_interp.cc.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
Header file for logic.cc.
Numeric diameter_area_equ_aerodynamical
#define max(a, b)
Refraction functions.
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
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.
void tmatrix_ampld_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1297