ARTS  2.2.66
m_montecarlo.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-2012 Cory Davis <cory@met.ed.ac.uk>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  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, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
18 
19 
20 /*===========================================================================
21  === File description
22  ===========================================================================*/
23 
35 /*===========================================================================
36  === External declarations
37  ===========================================================================*/
38 
39 #include "messages.h"
40 #include "arts.h"
41 #include "auto_md.h"
42 #include "ppath.h"
43 #include "matpackI.h"
44 #include "special_interp.h"
45 #include "check_input.h"
46 #include <stdexcept>
47 #include <cmath>
48 #include "rte.h"
49 #include "lin_alg.h"
50 #include "logic.h"
51 #include "physics_funcs.h"
52 #include "xml_io.h"
53 #include "montecarlo.h"
54 #include "rng.h"
55 #include <ctime>
56 #include <fstream>
57 #include "mc_interp.h"
58 #include "math_funcs.h"
59 
60 extern const Numeric DEG2RAD;
61 extern const Numeric RAD2DEG;
62 extern const Numeric PI;
63 extern const Numeric BOLTZMAN_CONST;
64 extern const Numeric SPEED_OF_LIGHT;
65 
66 
67 /*===========================================================================
68  === The functions (in alphabetical order)
69  ===========================================================================*/
70 
71 
72 /* Workspace method: Doxygen documentation will be auto-generated */
74  //keyword arguments
75  const Numeric& za_sigma,
76  const Numeric& aa_sigma,
77  const Verbosity&)
78 {
79  mc_antenna.set_gaussian(za_sigma,aa_sigma);
80 }
81 
82 
83 /* Workspace method: Doxygen documentation will be auto-generated */
85  //keyword arguments
86  const Numeric& za_fwhm,
87  const Numeric& aa_fwhm,
88  const Verbosity&)
89 {
90  mc_antenna.set_gaussian_fwhm(za_fwhm,aa_fwhm);
91 }
92 
93 
94 /* Workspace method: Doxygen documentation will be auto-generated */
96  const Verbosity&)
97 {
98  mc_antenna.set_pencil_beam();
99 }
100 
101 
102 /* Workspace method: Doxygen documentation will be auto-generated */
104  Vector& y,
105  Index& mc_iteration_count,
106  Vector& mc_error,
107  Tensor3& mc_points,
108  const MCAntenna& mc_antenna,
109  const Vector& f_grid,
110  const Index& f_index,
111  const Matrix& sensor_pos,
112  const Matrix& sensor_los,
113  const Index& stokes_dim,
114  const Index& atmosphere_dim,
115  const Agenda& ppath_step_agenda,
116  const Numeric& ppath_lraytrace,
117  const Agenda& iy_space_agenda,
118  const Agenda& surface_rtprop_agenda,
119  const Agenda& propmat_clearsky_agenda,
120  const Vector& p_grid,
121  const Vector& lat_grid,
122  const Vector& lon_grid,
123  const Tensor3& z_field,
124  const Vector& refellipsoid,
125  const Matrix& z_surface,
126  const Tensor3& t_field,
127  const Tensor4& vmr_field,
128  const Index& cloudbox_on,
129  const ArrayOfIndex& cloudbox_limits,
130  const Tensor4& pnd_field,
131  const ArrayOfSingleScatteringData& scat_data_array_mono,
132  const Index& atmfields_checked,
133  const Index& atmgeom_checked,
134  const Index& cloudbox_checked,
135  const Index& mc_seed,
136  const String& y_unit,
137  const Numeric& std_err,
138  const Index& max_time,
139  const Index& max_iter,
140  const Index& min_iter,
141  const Verbosity& verbosity)
142 {
143  // Basics
144  //
145  chk_if_in_range( "stokes_dim", stokes_dim, 1, 4 );
146  if( atmfields_checked != 1 )
147  throw runtime_error( "The atmospheric fields must be flagged to have "
148  "passed a consistency check (atmfields_checked=1)." );
149  if( atmgeom_checked != 1 )
150  throw runtime_error( "The atmospheric geometry must be flagged to have "
151  "passed a consistency check (atmgeom_checked=1)." );
152  if( cloudbox_checked != 1 )
153  throw runtime_error( "The cloudbox must be flagged to have "
154  "passed a consistency check (cloudbox_checked=1)." );
155  if( !cloudbox_on )
156  throw runtime_error( "The cloudbox must be activated (cloudbox_on=1)." );
157 
158  if( min_iter < 100 )
159  throw runtime_error( "*mc_min_iter* must be >= 100." );
160 
161  if( max_time < 0 && max_iter < 0 && std_err < 0 )
162  throw runtime_error( "At least one of std_err, max_time, and max_iter "
163  "must be positive." );
164 
165  if( f_index < 0 )
166  throw runtime_error( "The option of f_index < 0 is not handled by this "
167  "method." );
168  if( f_index >= f_grid.nelem() )
169  throw runtime_error( "*f_index* is outside the range of *f_grid*." );
170 
171  if( atmosphere_dim != 3 )
172  throw runtime_error( "Only 3D atmospheres are handled. " );
173 
174  if( sensor_pos.ncols() != 3 )
175  {
176  ostringstream os;
177  os << "Expected number of columns in sensor_pos: 3.\n";
178  os << "Found: " << sensor_pos.ncols();
179  throw runtime_error(os.str());
180  }
181 
182  if (! (sensor_los.ncols() == 2))
183  {
184  ostringstream os;
185  os << "Expected number of columns in sensor_los: 2.\n";
186  os << "Found: " << sensor_los.ncols();
187  throw runtime_error(os.str());
188  }
189 
190  Ppath ppath_step;
191  Rng rng; //Random Number generator
192  time_t start_time=time(NULL);
193  Index N_pt = pnd_field.nbooks();//Number of particle types
194  Vector pnd_vec(N_pt); //Vector of particle number densities used at each point
195  Vector Z11maxvector;//Vector holding the maximum phase function for each
196  bool anyptype30 = is_anyptype30(scat_data_array_mono);
197  if (anyptype30)
198  { findZ11max(Z11maxvector,scat_data_array_mono); }
199  rng.seed(mc_seed, verbosity);
200  Numeric g,temperature,albedo,g_los_csc_theta;
201  Matrix A(stokes_dim,stokes_dim), Q(stokes_dim,stokes_dim);
202  Matrix evol_op(stokes_dim,stokes_dim), ext_mat_mono(stokes_dim,stokes_dim);
203  Matrix q(stokes_dim,stokes_dim), newQ(stokes_dim,stokes_dim);
204  Matrix Z(stokes_dim,stokes_dim);
205  q = 0.0;
206  newQ = 0.0;
207  Vector vector1(stokes_dim), abs_vec_mono(stokes_dim), I_i(stokes_dim);
208  Vector Isum(stokes_dim), Isquaredsum(stokes_dim);
209  Index termination_flag = 0;
210  const Numeric f_mono = f_grid[f_index];
211 
212  CREATE_OUT0;
213 
214  y.resize(stokes_dim);
215  y = 0;
216 
217  mc_iteration_count = 0;
218  mc_error.resize(stokes_dim);
219  mc_points.resize( p_grid.nelem(), lat_grid.nelem(), lon_grid.nelem() );
220  mc_points = 0;
221 
222  //local versions of workspace
223  Matrix local_iy(1,stokes_dim), local_surface_emission(1,stokes_dim);
224  Matrix local_surface_los;
225  Tensor4 local_surface_rmatrix;
226  Vector local_rte_pos(2);
227  Vector local_rte_los(2);
228  Vector new_rte_los(2);
229  Index np;
230  Isum=0.0;Isquaredsum=0.0;
231  Numeric std_err_i;
232  bool convert_to_rjbt = false;
233  if ( y_unit == "RJBT" )
234  {
235  std_err_i = f_mono*f_mono*2*BOLTZMAN_CONST/
237  convert_to_rjbt = true;
238  }
239  else if ( y_unit == "1" )
240  { std_err_i = std_err; }
241  else
242  {
243  ostringstream os;
244  os << "Invalid value for *y_unit*:" << y_unit <<".\n"
245  << "This method allows only the options \"RJBT\" and \"1\".";
246  throw runtime_error( os.str() );
247  }
248 
249 
250  //Begin Main Loop
251  //
252  bool keepgoing, oksampling;
253  //
254  while( true )
255  {
256  bool inside_cloud;
257 
258  mc_iteration_count += 1;
259 
260  keepgoing = true; // indicating whether to continue tracing a photon
261  oksampling = true; // gets false if g becomes zero
262 
263  //Sample a FOV direction
264  mc_antenna.draw_los(local_rte_los,rng,sensor_los(0,joker));
265  id_mat(Q);
266  local_rte_pos=sensor_pos(0,joker);
267  I_i=0.0;
268 
269  while( keepgoing )
270  {
271  mcPathTraceGeneral( ws, evol_op, abs_vec_mono, temperature,
272  ext_mat_mono, rng, local_rte_pos, local_rte_los,
273  pnd_vec, g,ppath_step, termination_flag,
274  inside_cloud, ppath_step_agenda, ppath_lraytrace,
275  propmat_clearsky_agenda, stokes_dim, f_mono,
276  p_grid, lat_grid, lon_grid, z_field, refellipsoid,
277  z_surface, t_field, vmr_field,
278  cloudbox_limits, pnd_field, scat_data_array_mono,
279  verbosity );
280 
281  np = ppath_step.np;
282  mc_points(ppath_step.gp_p[np-1].idx,ppath_step.gp_lat[np-1].idx,
283  ppath_step.gp_lon[np-1].idx) += 1;
284 
285  // GH 2011-09-08: if the lowest layer has large
286  // extent and a thick cloud, g may be 0 due to
287  // underflow, but then I_i should be 0 as well.
288  // Don't turn it into nan for no reason.
289  // If reaching underflow, no point in going on;
290  // hence new photon.
291  // GH 2011-09-14: moved this check to outside the different
292  // scenarios, as this goes wrong regardless of the scenario.
293  if( g == 0 )
294  {
295  keepgoing = false;
296  oksampling = false;
297  mc_iteration_count -= 1;
298  out0 << "WARNING: A rejected path sampling (g=0)!\n(if this"
299  << "happens repeatedly, try to decrease *ppath_lmax*)";
300  }
301  else if( termination_flag == 1 )
302  {
303  iy_space_agendaExecute( ws, local_iy, Vector(1,f_mono),
304  local_rte_pos, local_rte_los,
305  iy_space_agenda );
306  mult( vector1, evol_op, local_iy(0,joker) );
307  mult( I_i, Q, vector1 );
308  I_i /= g;
309  keepgoing=false; //stop here. New photon.
310  }
311  else if( termination_flag == 2 )
312  {
313  //Calculate surface properties
314  surface_rtprop_agendaExecute( ws, local_surface_emission,
315  local_surface_los,
316  local_surface_rmatrix,
317  Vector(1,f_mono),
318  local_rte_pos, local_rte_los,
319  surface_rtprop_agenda );
320 
321  if( local_surface_los.nrows() > 1 )
322  throw runtime_error(
323  "The method handles only specular reflections." );
324 
325  //deal with blackbody case
326  if( local_surface_los.nrows() == 0 )
327  {
328  mult( vector1, evol_op, local_surface_emission(0,joker) );
329  mult( I_i, Q, vector1);
330  I_i /= g;
331  keepgoing = false;
332  }
333  else
334  //decide between reflection and emission
335  {
336  Numeric R11 = local_surface_rmatrix(0,0,0,0);
337  if( rng.draw() > R11 )
338  {
339  //then we have emission
340  mult( vector1, evol_op, local_surface_emission(0,joker) );
341  mult( I_i, Q, vector1 );
342  I_i /= g*(1-R11);
343  keepgoing = false;
344  }
345  else
346  {
347  //we have reflection
348  local_rte_los = local_surface_los( 0, joker );
349 
350  mult( q, evol_op, local_surface_rmatrix(0,0,joker,joker));
351  mult( newQ, Q, q );
352  Q = newQ;
353  Q /= g*R11;
354  }
355  }
356  }
357  else if (inside_cloud)
358  {
359  //we have another scattering/emission point
360  //Estimate single scattering albedo
361  albedo = 1 - abs_vec_mono[0]/ext_mat_mono(0,0);
362 
363  //determine whether photon is emitted or scattered
364  if( rng.draw() > albedo )
365  {
366  //Calculate emission
367  Numeric planck_value = planck( f_mono, temperature );
368  Vector emission = abs_vec_mono;
369  emission *= planck_value;
370  Vector emissioncontri(stokes_dim);
371  mult( emissioncontri, evol_op, emission );
372  emissioncontri /= (g*(1-albedo)); //yuck!
373  mult( I_i, Q, emissioncontri );
374  keepgoing = false;
375  }
376  else
377  {
378  //we have a scattering event
379  Sample_los( new_rte_los, g_los_csc_theta, Z, rng,
380  local_rte_los, scat_data_array_mono, stokes_dim,
381  pnd_vec, anyptype30, Z11maxvector,
382  ext_mat_mono(0,0)-abs_vec_mono[0], temperature,
383  verbosity );
384 
385  Z /= g * g_los_csc_theta * albedo;
386 
387  mult( q, evol_op, Z );
388  mult( newQ, Q, q );
389  Q = newQ;
390  //scattering_order+=1;
391  local_rte_los = new_rte_los;
392  }
393  }
394  else
395  {
396  //Must be clear sky emission point
397  //Calculate emission
398  Numeric planck_value = planck( f_mono, temperature );
399  Vector emission = abs_vec_mono;
400  emission *= planck_value;
401  Vector emissioncontri(stokes_dim);
402  mult( emissioncontri, evol_op, emission );
403  emissioncontri /= g;
404  mult( I_i, Q, emissioncontri );
405  keepgoing = false;
406  }
407  } // keepgoing
408 
409  if( oksampling )
410  {
411  Isum += I_i;
412 
413  for( Index j=0; j<stokes_dim; j++ )
414  {
415  assert( !isnan(I_i[j]) );
416  Isquaredsum[j] += I_i[j]*I_i[j];
417  }
418  y = Isum;
419  y /= (Numeric)mc_iteration_count;
420  for(Index j=0; j<stokes_dim; j++)
421  {
422  mc_error[j] = sqrt( ( Isquaredsum[j] /
423  (Numeric)mc_iteration_count - y[j]*y[j] ) /
424  (Numeric)mc_iteration_count );
425  }
426  if( std_err > 0 && mc_iteration_count >= min_iter &&
427  mc_error[0] < std_err_i )
428  { break; }
429  if( max_time > 0 && (Index)(time(NULL)-start_time) >= max_time )
430  { break; }
431  if( max_iter > 0 && mc_iteration_count >= max_iter )
432  { break; }
433  }
434  }
435 
436  if( convert_to_rjbt )
437  {
438  for(Index j=0; j<stokes_dim; j++)
439  {
440  y[j] = invrayjean( y[j], f_mono );
441  mc_error[j] = invrayjean( mc_error[j], f_mono );
442  }
443  }
444 }
445 
446 
447 /* Workspace method: Doxygen documentation will be auto-generated */
448 void MCSetSeedFromTime(Index& mc_seed,
449  const Verbosity&)
450 {
451  mc_seed=(Index)time(NULL);
452 }
453 
454 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
const Numeric DEG2RAD
void set_pencil_beam(void)
makes the antenna pattern a pencil beam
Definition: mc_antenna.cc:88
void MCGeneral(Workspace &ws, Vector &y, Index &mc_iteration_count, Vector &mc_error, Tensor3 &mc_points, const MCAntenna &mc_antenna, const Vector &f_grid, const Index &f_index, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &stokes_dim, const Index &atmosphere_dim, const Agenda &ppath_step_agenda, const Numeric &ppath_lraytrace, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &cloudbox_checked, const Index &mc_seed, const String &y_unit, const Numeric &std_err, const Index &max_time, const Index &max_iter, const Index &min_iter, const Verbosity &verbosity)
WORKSPACE METHOD: MCGeneral.
The Agenda class.
Definition: agenda_class.h:44
Interpolation classes and functions created for use within Monte Carlo scattering simulations...
Declarations having to do with the four output streams.
void set_gaussian(const Numeric &za_sigma, const Numeric &aa_sigma)
makes the antenna pattern a 2D gaussian specified by za and aa standard deviations ...
Definition: mc_antenna.cc:100
#define q
Definition: continua.cc:21469
Numeric invrayjean(const Numeric &i, const Numeric &f)
invrayjean
void findZ11max(Vector &Z11maxvector, const ArrayOfSingleScatteringData &scat_data_array_mono)
findZ11max
Definition: montecarlo.cc:337
The Vector class.
Definition: matpackI.h:556
void draw_los(VectorView &sampled_rte_los, Rng &rng, ConstVectorView bore_sight_los) const
draws a line of sight by sampling the antenna response function
Definition: mc_antenna.cc:167
The Tensor4 class.
Definition: matpackIV.h:383
Linear algebra functions.
This file contains basic functions to handle XML data files.
void mc_antennaSetPencilBeam(MCAntenna &mc_antenna, const Verbosity &)
WORKSPACE METHOD: mc_antennaSetPencilBeam.
Definition: m_montecarlo.cc:95
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
void mc_antennaSetGaussianByFWHM(MCAntenna &mc_antenna, const Numeric &za_fwhm, const Numeric &aa_fwhm, const Verbosity &)
WORKSPACE METHOD: mc_antennaSetGaussianByFWHM.
Definition: m_montecarlo.cc:84
const Numeric SPEED_OF_LIGHT
The implementation for String, the ARTS string class.
Definition: mystring.h:63
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The Tensor3 class.
Definition: matpackIII.h:348
const Numeric PI
The global header file for ARTS.
const Numeric RAD2DEG
void Sample_los(VectorView new_rte_los, Numeric &g_los_csc_theta, MatrixView Z, Rng &rng, ConstVectorView rte_los, const ArrayOfSingleScatteringData &scat_data_array_mono, const Index stokes_dim, ConstVectorView pnd_vec, const bool anyptype30, ConstVectorView Z11maxvector, const Numeric Csca, const Numeric rtp_temperature, const Verbosity &verbosity)
Sample_los.
Definition: montecarlo.cc:1187
void iy_space_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:13886
const Joker joker
Defines the Rng random number generator class.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
const Numeric BOLTZMAN_CONST
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:56
Header file for special_interp.cc.
Propagation path structure and functions.
void resize(Index p, Index r, Index c)
Resize function.
Definition: matpackIII.cc:862
double draw()
Definition: rng.cc:120
Header file for logic.cc.
This can be used to make arrays out of anything.
Definition: array.h:40
void mcPathTraceGeneral(Workspace &ws, MatrixView evol_op, Vector &abs_vec_mono, Numeric &temperature, MatrixView ext_mat_mono, Rng &rng, Vector &rte_pos, Vector &rte_los, Vector &pnd_vec, Numeric &g, Ppath &ppath_step, Index &termination_flag, bool &inside_cloud, const Agenda &ppath_step_agenda, const Agenda &propmat_clearsky_agenda, const Index stokes_dim, const Numeric &f_mono, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Tensor3 &edensity_field, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfSingleScatteringData &scat_data_array_mono, const Verbosity &verbosity)
mcPathTraceGeneral
Definition: mc_NotUsed.cc:150
Definition: rng.h:569
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:101
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
void mc_antennaSetGaussian(MCAntenna &mc_antenna, const Numeric &za_sigma, const Numeric &aa_sigma, const Verbosity &)
WORKSPACE METHOD: mc_antennaSetGaussian.
Definition: m_montecarlo.cc:73
Numeric planck(const Numeric &f, const Numeric &t)
planck
Workspace class.
Definition: workspace_ng.h:47
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
bool is_anyptype30(const ArrayOfSingleScatteringData &scat_data_array_mono)
is_anyptype30
Definition: montecarlo.cc:373
#define CREATE_OUT0
Definition: messages.h:211
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:299
void seed(unsigned long int n, const Verbosity &verbosity)
Definition: rng.cc:72
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
void set_gaussian_fwhm(const Numeric &za_fwhm, const Numeric &aa_fwhm)
makes the antenna pattern a 2D gaussian specified by za and aa FWHM
Definition: mc_antenna.cc:116
void surface_rtprop_agendaExecute(Workspace &ws, Matrix &surface_emission, Matrix &surface_los, Tensor4 &surface_rmatrix, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:15046
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
Declaration of functions in rte.cc.