ARTS  2.2.66
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 <stdexcept>
27 #include <cmath>
28 #include "messages.h"
29 #include "tmatrix.h"
30 #include "check_input.h"
31 #include "make_array.h"
32 #include "refraction.h"
33 #include "special_interp.h"
34 
35 
36 
37 void TMatrixTest(const Verbosity& verbosity)
38 {
39  tmatrix_tmd_test(verbosity);
40  tmatrix_ampld_test(verbosity);
41  calc_ssp_random_test(verbosity);
42  calc_ssp_fixed_test(verbosity);
43 }
44 
45 //-----------------------------------
46 
47 
48 void scat_meta_arrayInit(// WS Output:
49  ArrayOfScatteringMetaData& scat_meta_array,
50  //WS Input
51  const Verbosity&)
52 {
53  scat_meta_array.resize(0);
54 }
55 
56 //-----------------------------------
57 
58 
59 void scat_meta_arrayAddTmatrix(// WS Output:
60  ArrayOfScatteringMetaData& scat_meta_array,
61  // WS Input:
62  const GriddedField3& complex_refr_index,
63  // WS Generic input:
64  const String& description,
65  const String& material,
66  const String& shape,
67  const String& particle_type,
68  const Numeric& density,
69  const Vector& aspect_ratio_grid,
70  const Vector& diameter_max_grid,
71  const Vector& scat_f_grid,
72  const Vector& scat_T_grid,
73  const Verbosity&)
74 
75 {
76  for(Index k=0; k<diameter_max_grid.nelem(); k++)
77  {
78  for ( Index i=0; i < aspect_ratio_grid.nelem(); i++ )
79  {
80  extern const Numeric PI;
81  Numeric volume;
82 
83  if (shape == "spheroidal")
84  {
85  if (aspect_ratio_grid[i]>1) // If oblate spheroid
86  {
87  Numeric a; //rotational axis (perpendicular to a)
88 
89  a=diameter_max_grid[k]/2.;
90  volume=pow(a,3.)*4.*PI/(3.*aspect_ratio_grid[i]);
91  }
92 
93  else if (aspect_ratio_grid[i]<1) // If prolate spheroid
94  {
95  Numeric b; //non-rotational axis (perpendicular to a)
96 
97  b=diameter_max_grid[k]/2.;
98  volume=pow(b,3.)*4.*PI*pow(aspect_ratio_grid[i],2.)/3.;
99  }
100 
101  else
102  {
103  ostringstream os;
104  os << "Incorrect aspect ratio: " << aspect_ratio_grid[i] << "\n"
105  << "Can not be equal to one";
106  throw std::runtime_error(os.str());
107  }
108  }
109  else if (shape == "cylindrical")
110  {
111  // Formulas to determine cylindrical volume from diamter_max:
112 
113  // D/L=aspect_ratio
114  // diameter_max=sqrt(D²+L²)
115  // volume=D³*pi/(4*aspect_ratio)
116 
117  volume=pow(diameter_max_grid[k]/pow((pow(aspect_ratio_grid[i], 2.)+1.), 1./2.), 3)*pow(aspect_ratio_grid[i], 2.)*PI/4.;
118  }
119 
120  else
121  {
122  ostringstream os;
123  os << "Unknown particle shape: " << shape << "\n"
124  << "Must be spheroidal or cylindrical";
125  throw std::runtime_error(os.str());
126  }
127 
128  ScatteringMetaData smd;
129  if (description=="")
130  {
131  ostringstream os;
132  os << shape<< " "<< material << " particle of type " << particle_type<<
133  ", with volume equivalent diameter "
134  <<diameter_max_grid[k]<<" meters.";
135  smd.description=os.str();
136  }
137  else
138  smd.description = description;
139 
140  smd.material = material;
141  smd.shape = shape;
142  smd.particle_type = ParticleTypeFromString(particle_type);
143  smd.ssd_method = PARTICLE_SSDMETHOD_TMATRIX;
144  smd.density = density;
145  smd.diameter_max =diameter_max_grid[k];
146  smd.volume = volume;
147  smd.area_projected = 0;
148  smd.aspect_ratio = aspect_ratio_grid[i];
149  smd.scat_f_grid = scat_f_grid;
150  smd.scat_T_grid = scat_T_grid;
151  smd.complex_refr_index = complex_refr_index;
152 
153  scat_meta_array.push_back(smd);
154  }
155  }
156 }
157 
158 //-----------------------------------
159 
160 
161 void scat_data_arrayFromMeta(// WS Output:
162  ArrayOfSingleScatteringData& scat_data_array,
163  //WS Input
164  const ArrayOfScatteringMetaData& scat_meta_array,
165  // WS Generic input:
166  const Vector& za_grid,
167  const Vector& aa_grid,
168  const Numeric& precision,
169  const Verbosity&)
170 {
171  for(Index ii=0;ii<scat_meta_array.nelem();ii++)
172  {
173 
174  //Interpolation
175 
176  Tensor3 complex_refr_index;
177 
178 
179  const Index nf = scat_meta_array[ii].scat_f_grid.nelem();
180  const Index nt = scat_meta_array[ii].scat_T_grid.nelem();
181 
182  complex_refr_index.resize(nf,nt,2);
183 
184  complex_n_interp(complex_refr_index(joker, joker,0), complex_refr_index(joker, joker,1),
185  scat_meta_array[ii].complex_refr_index, "complex_refr_index",
186  scat_meta_array[ii].scat_f_grid, scat_meta_array[ii].scat_T_grid);
187 
188 
189  extern const Numeric PI;
190  Index np;
191 
193  sdd.f_grid = scat_meta_array[ii].scat_f_grid;
194  sdd.T_grid = scat_meta_array[ii].scat_T_grid;
195  sdd.za_grid = za_grid;
196  sdd.aa_grid = aa_grid;
197  sdd.particle_type = scat_meta_array[ii].particle_type;
198 
199  if (scat_meta_array[ii].shape == "spheroidal" )
200  np=-1;
201 
202  else if (scat_meta_array[ii].shape == "cylindrical")
203  np=-2;
204  else
205  {
206  ostringstream os;
207  os << "Unknown particle shape: " << scat_meta_array[ii].shape << "\n"
208  << "Must be spheroidal or cylindrical";
209  throw std::runtime_error(os.str());
210  }
211 
213  complex_refr_index(joker,joker,0),
214  complex_refr_index(joker,joker,1),
215  pow(scat_meta_array[ii].volume*1e18*3./(4.*PI),1./3.),
216  np,
217  scat_meta_array[ii].aspect_ratio,
218  precision);
219 
220  scat_data_array.push_back(sdd);
221  }
222 }
223 
224 
225 
226 
227 
228 
229 
230 
231 
232 
233 
235  ArrayOfScatteringMetaData& scat_meta_array,
236  // WS Input:
237  const GriddedField3& complex_refr_index,
238  // WS Generic input:
239  const String& description,
240  const String& material,
241  const String& shape,
242  const String& particle_type,
243  const Numeric& density,
244  const Numeric& aspect_ratio,
245  const Vector& diameter_grid,
246  const Vector& scat_f_grid,
247  const Vector& scat_T_grid,
248  const Verbosity&)
249 
250 {
251  for(Index k=0; k<diameter_grid.nelem(); k++)
252  {
253  extern const Numeric PI;
254 
255  Numeric volume;
256 
257  volume=4./3.*PI*pow(diameter_grid[k]/2.,3);
258 
259 
260  Numeric diameter_max;
261 
262  if (shape == "spheroidal")
263  {
264  if (aspect_ratio>1) // If oblate spheroid
265  {
266  Numeric a; //rotational axis
267 
268  a=pow(3*volume*aspect_ratio/(4*PI), 1./3.);
269  diameter_max=2*a;
270  }
271 
272  else if (aspect_ratio<1) // If prolate spheroid
273  {
274  Numeric b; //non-rotational axis (perpendicular to a)
275 
276  b=pow(3*volume/(4*PI*pow(aspect_ratio, 2)),1./3.);
277  diameter_max=2*b;
278  }
279 
280  else
281  {
282  ostringstream os;
283  os << "Incorrect aspect ratio: " << aspect_ratio << "\n"
284  << "Can not be equal to one";
285  throw std::runtime_error(os.str());
286  }
287  }
288  else if (shape == "cylindrical")
289  {
290  Numeric D;
291  Numeric L;
292 
293  //aspect_ratio=D/L
294 
295  D=pow(volume*4./PI*aspect_ratio, 1./3.);
296  L=D/aspect_ratio;
297  diameter_max=pow((pow(D,2)+pow(L,2)), 1./2.);
298  }
299 
300  else
301  {
302  ostringstream os;
303  os << "Unknown particle shape: " << shape << "\n"
304  << "Must be spheroidal or cylindrical";
305  throw std::runtime_error(os.str());
306  }
307 
308  ScatteringMetaData smd;
309  if (description=="")
310  {
311  ostringstream os;
312  os << shape<< " "<< material << " particle of type " << particle_type<<
313  ", with volume equivalent diameter "
314  <<diameter_grid[k]<<" meters.";
315  smd.description=os.str();
316  }
317  else
318  smd.description = description;
319 
320  smd.material = material;
321  smd.shape = shape;
322  smd.particle_type = ParticleTypeFromString(particle_type);
323  smd.ssd_method = PARTICLE_SSDMETHOD_TMATRIX;
324  smd.density = density;
325  smd.diameter_max =diameter_max;
326  smd.volume = volume;
327  smd.area_projected = 0;
328  smd.aspect_ratio = aspect_ratio;
329  smd.scat_f_grid = scat_f_grid;
330  smd.scat_T_grid = scat_T_grid;
331  smd.complex_refr_index = complex_refr_index;
332 
333  scat_meta_array.push_back(smd);
334  }
335 }
336 
337 //-----------------------------------
338 
339 
void tmatrix_tmd_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1256
Declarations for the T-Matrix interface.
void TMatrixTest(const Verbosity &verbosity)
WORKSPACE METHOD: TMatrixTest.
Definition: m_tmatrix.cc:37
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
#define precision
Definition: logic.cc:45
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)
Calculate SingleScatteringData properties.
Definition: tmatrix.cc:927
Index nelem() const
Number of elements.
Definition: array.h:176
Declarations having to do with the four output streams.
void scat_data_arrayFromMeta(ArrayOfSingleScatteringData &scat_data_array, const ArrayOfScatteringMetaData &scat_meta_array, const Vector &za_grid, const Vector &aa_grid, const Numeric &precision, const Verbosity &)
WORKSPACE METHOD: scat_data_arrayFromMeta.
Definition: m_tmatrix.cc:161
The Vector class.
Definition: matpackI.h:556
void scat_meta_arrayAddTmatrixOldVersion(ArrayOfScatteringMetaData &scat_meta_array, const GriddedField3 &complex_refr_index, const String &description, const String &material, const String &shape, const String &particle_type, const Numeric &density, const Numeric &aspect_ratio, const Vector &diameter_grid, const Vector &scat_f_grid, const Vector &scat_T_grid, const Verbosity &)
WORKSPACE METHOD: scat_meta_arrayAddTmatrixOldVersion.
Definition: m_tmatrix.cc:234
void scat_meta_arrayAddTmatrix(ArrayOfScatteringMetaData &scat_meta_array, const GriddedField3 &complex_refr_index, const String &description, const String &material, const String &shape, const String &particle_type, const Numeric &density, const Vector &aspect_ratio_grid, const Vector &diameter_max_grid, const Vector &scat_f_grid, const Vector &scat_T_grid, const Verbosity &)
WORKSPACE METHOD: scat_meta_arrayAddTmatrix.
Definition: m_tmatrix.cc:59
Structure which describes the single scattering properties of a particle or a particle distribution...
Definition: optproperties.h:84
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
The implementation for String, the ARTS string class.
Definition: mystring.h:63
The Tensor3 class.
Definition: matpackIII.h:348
void calc_ssp_fixed_test(const Verbosity &verbosity)
Single scattering properties calculation for particles with fixed orientation.
Definition: tmatrix.cc:1409
const Numeric PI
void calc_ssp_random_test(const Verbosity &verbosity)
Single scattering properties calculation for randomly oriented particles.
Definition: tmatrix.cc:1351
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
Implements the class MakeArray, which is a derived class of Array, allowing explicit initialization...
Header file for special_interp.cc.
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:862
void scat_meta_arrayInit(ArrayOfScatteringMetaData &scat_meta_array, const Verbosity &)
WORKSPACE METHOD: scat_meta_arrayInit.
Definition: m_tmatrix.cc:48
This can be used to make arrays out of anything.
Definition: array.h:40
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)
complex_n_interp
ParticleType particle_type
Definition: optproperties.h:85
void tmatrix_ampld_test(const Verbosity &verbosity)
T-Matrix validation test.
Definition: tmatrix.cc:1189
ParticleType ParticleTypeFromString(const String &particle_type_string)
Convert particle name to enum value.