ARTS  2.2.66
m_doit2.cc
Go to the documentation of this file.
1 #include <stdexcept>
2 #include <iostream>
3 #include <cstdlib>
4 #include <cmath>
5 #include "arts.h"
6 #include "array.h"
7 #include "auto_md.h"
8 #include "check_input.h"
9 #include "matpackVII.h"
10 #include "logic.h"
11 #include "ppath.h"
12 #include "agenda_class.h"
13 #include "physics_funcs.h"
14 #include "lin_alg.h"
15 #include "math_funcs.h"
16 #include "messages.h"
17 #include "xml_io.h"
18 #include "rte.h"
19 #include "special_interp.h"
20 #include "doit.h"
21 #include "m_general.h"
22 #include "wsv_aux.h"
23 #include "geodetic.h"
24 
25 
26 
27 /* Workspace method: Doxygen documentation will be auto-generated */
29  Workspace& ws,
30  Tensor7& doit_i_field,
31  const Agenda& iy_main_agenda,
32  const Index& atmosphere_dim,
33  const Vector& lat_grid,
34  const Vector& lon_grid,
35  const Tensor3& z_field,
36  const Tensor3& t_field,
37  const Tensor4& vmr_field,
38  const Matrix& z_surface,
39  const Index& cloudbox_on,
40  const ArrayOfIndex& cloudbox_limits,
41  const Index& basics_checked,
42  const Index& cloudbox_checked,
43  const Index& sensor_checked,
44  const Vector& f_grid,
45  const Index& stokes_dim,
46  const String& iy_unit,
47  const Agenda& blackbody_radiation_agenda,
48  const Vector& scat_za_grid,
49  const Vector& scat_aa_grid,
50  const Index& fill_interior,
51  const Verbosity&)
52 {
53  // Don't do anything if there's no cloudbox defined.
54  if (!cloudbox_on) return;
55 
56  // Basics and cloudbox OK?
57  if( basics_checked<2 )
58  throw runtime_error("The atmosphere, surface and basic control variables "
59  "must be flagged to have passed a consistency check\n"
60  "by basics_checkedCalc (basics_checked=2)!" );
61  if( !cloudbox_checked )
62  throw runtime_error( "The cloudbox must be flagged to have passed a "
63  "consistency check (cloudbox_checked=1)." );
64  if( !sensor_checked )
65  throw runtime_error( "The sensor variables must be flagged to have passed"
66  "a consistency check (sensor_checked=1)." );
67 
68  // Main sizes
69  const Index Nf = f_grid.nelem();
70  const Index Np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
71  const Index Nza = scat_za_grid.nelem();
72  const Index Ni = stokes_dim;
73  Index Nlat=1, Nlon=1, Naa=1;
74 
75  //--- Special checks -------------------------------------------------------
76  if( !(atmosphere_dim == 1 || atmosphere_dim == 3) )
77  throw runtime_error( "The atmospheric dimensionality must be 1 or 3.");
78  // DOIT requires frequency based radiance:
79  if( iy_unit != "1" ||
80  !chk_if_std_blackbody_agenda( ws, blackbody_radiation_agenda ) )
81  {
82  ostringstream os;
83  os << "It is assumed that you use this method together with DOIT.\n"
84  << "Usage of this method then demands that the *iy_main_agenda*\n"
85  << "returns frequency based radiance (ie. [W/m2/Hz/sr]).\n"
86  << "This requires that *iy_unit* is set to \"1\" and that\n"
87  << "*blackbody_radiation_agenda uses *blackbody_radiationPlanck*\n"
88  << "or a corresponding WSM.\n"
89  << "At least one of these requirements is not met.\n";
90  throw runtime_error( os.str() );
91  }
92  if( scat_za_grid[0] != 0. || scat_za_grid[Nza-1] != 180. )
93  throw runtime_error(
94  "*scat_za_grid* must include 0 and 180 degrees as endpoints." );
95  if( atmosphere_dim == 3 )
96  {
97  Naa = scat_aa_grid.nelem();
98  if( scat_aa_grid[0] != 0. || scat_aa_grid[Naa-1] != 360. )
99  throw runtime_error(
100  "*scat_aa_grid* must include 0 and 360 degrees as endpoints." );
101  Nlat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
102  Nlon = cloudbox_limits[5] - cloudbox_limits[4] + 1;
103  }
104  //--------------------------------------------------------------------------
105 
106  // Size doit_i_field and set to dummy value
107  doit_i_field.resize( Nf, Np, Nlat, Nlon, Nza, Naa, Ni );
108  doit_i_field = -999e9;
109 
110  // A matrix to flag if surface done for lat/lon
111  Matrix surface_done( Nlat, Nlon, -1 );
112 
113  // Convert scat_aa_grid to "sensor coordinates"
114  // (-180° < azimuth angle < 180°)
115  Vector aa_grid(Naa);
116  for(Index i = 0; i<Naa; i++)
117  { aa_grid[i] = scat_aa_grid[i] - 180; }
118 
119  // Define the variables for position and direction.
120  Vector pos( atmosphere_dim );
121  Vector los( max( Index(1), atmosphere_dim-1 ) );
122 
123  for( Index o=0; o<Nlon; o++ )
124  {
125  for( Index a=0; a<Nlat; a++ )
126  {
127 
128  if( atmosphere_dim == 3 )
129  {
130  pos[1] = lat_grid[ a + cloudbox_limits[2] ];
131  pos[2] = lon_grid[ o + cloudbox_limits[4] ];
132  }
133 
134  for( Index p=Np-1; p>=0; p-- ) // Note looping downwards
135  {
136 
137  if( fill_interior || p==0 || p==Np-1 || ( atmosphere_dim==3 &&
138  ( a==0 || a==Nlat-1 || o==0 || o==Nlon-1 ) ) )
139  {
140 
141  pos[0] = z_field( p+cloudbox_limits[0], 0, 0 );
142 
143  // Above the surface
144  if( pos[0] > z_surface(a,o) )
145  {
146  for( Index za=0; za<Nza; za++ ) {
147 
148  los[0] = scat_za_grid[za];
149  Matrix iy;
150 
151  for( Index aa=0; aa<Naa; aa++ ) {
152  // For end points of scat_za_grid, we need only to
153  // perform calculations for first aa
154  if( ( za != 0 && za != (Nza-1) ) || aa == 0 )
155  {
156  if( atmosphere_dim == 3 )
157  { los[1] = aa_grid[aa]; }
158 
159  get_iy( ws, iy, t_field, z_field, vmr_field, 0,
160  f_grid, pos, los, Vector(0),
161  iy_main_agenda );
162  }
163 
164  doit_i_field( joker, p, a, o, za, aa, joker ) = iy;
165  }
166  }
167  }
168  // Surface or below
169  else
170  {
171  // If surface already done for lat/lon, just copy
172  // from above
173  if( surface_done( a, o ) > 0 )
174  {
175  for( Index za=0; za<Nza; za++ ) {
176  for( Index aa=0; aa<Naa; aa++ ) {
177  for( Index iv=0; iv<Nf; iv++ ) {
178  for( Index is=0; is<stokes_dim; is++ ) {
179  doit_i_field(iv,p,a,o,za,aa,is) =
180  doit_i_field(iv,p+1,a,o,za,aa,is);
181  } } } }
182  }
183  // Calculate for surface altitude
184  else
185  {
186  pos[0] = z_surface(a,o);
187  for( Index za=0; za<Nza; za++ ) {
188  //
189  los[0] = scat_za_grid[za];
190  Matrix iy;
191  //
192  for( Index aa=0; aa<Naa; aa++ ) {
193  // For end points of scat_za_grid, we need only to
194  // perform calculations for first aa
195  if( ( za != 0 && za != (Nza-1) ) || aa == 0 )
196  {
197  if( atmosphere_dim == 3 )
198  { los[1] = aa_grid[aa]; }
199 
200  get_iy( ws, iy, t_field, z_field, vmr_field,
201  0, f_grid, pos, los, Vector(0),
202  iy_main_agenda );
203  }
204  doit_i_field(joker,p,a,o,za,aa,joker) = iy;
205  }
206  }
207  surface_done( a, o ) = 1;
208  }
209  }
210  }
211  }
212  }
213  }
214 }
215 
216 
217 
218 
219 /* Workspace method: Doxygen documentation will be auto-generated */
221  Matrix& iy,
222  const Tensor7& doit_i_field,
223  const Vector& rte_pos,
224  const Vector& rte_los,
225  const Index& jacobian_do,
226  const Index& cloudbox_on,
227  const ArrayOfIndex& cloudbox_limits,
228  const Index& basics_checked,
229  const Index& cloudbox_checked,
230  const Index& atmosphere_dim,
231  const Vector& p_grid,
232  const Vector& lat_grid,
233  const Vector& lon_grid,
234  const Tensor3& z_field,
235  const Index& stokes_dim,
236  const Vector& scat_za_grid,
237  const Vector& scat_aa_grid,
238  const Vector& f_grid,
239  const Index& za_order,
240  const Verbosity& )
241 {
242  // Basics and cloudbox OK?
243  if( !(atmosphere_dim == 1 || atmosphere_dim == 3) )
244  throw runtime_error( "The atmospheric dimensionality must be 1 or 3.");
245  if( !basics_checked )
246  throw runtime_error( "The atmosphere and basic control variables must be "
247  "flagged to have passed a consistency check (basics_checked=1)." );
248  if( !cloudbox_checked )
249  throw runtime_error( "The cloudbox must be flagged to have passed a "
250  "consistency check (cloudbox_checked=1)." );
251  if( !cloudbox_on )
252  throw runtime_error( "The cloud box is not activated and no outgoing "
253  "field can be returned." );
254 
255  // Main sizes
256  const Index nf = f_grid.nelem();
257  const Index nza = scat_za_grid.nelem();
258  const Index np = cloudbox_limits[1]-cloudbox_limits[0]+1;
259  Index naa = 1;
260  Index nlat = 1;
261  Index nlon = 1;
262  if( atmosphere_dim == 3 )
263  {
264  naa = scat_aa_grid.nelem();
265  nlat = cloudbox_limits[3]-cloudbox_limits[2]+1;
266  nlon = cloudbox_limits[5]-cloudbox_limits[4]+1;
267  }
268 
269  //--- Check input -----------------------------------------------------------
270  // rte_pos checked as part rte_pos2gridpos
271  chk_rte_los( atmosphere_dim, rte_los );
272  if( nza == 0 )
273  throw runtime_error( "The variable *scat_za_grid* is empty. Are dummy "
274  "values from *cloudboxOff used?" );
275  if( za_order < 1 || za_order > 5 )
276  throw runtime_error( "The argument *za_order* must be in the range [1,5].");
277  if( doit_i_field.ncols() != stokes_dim )
278  throw runtime_error( "Inconsistency in size between *stokes_dim* and "
279  "*doit_i_field*." );
280  if( doit_i_field.nrows() != naa )
281  throw runtime_error( "Inconsistency in size between *scat_aa_grid* and "
282  "*doit_i_field*." );
283  if( doit_i_field.npages() != nza )
284  throw runtime_error( "Inconsistency in size between *scat_za_grid* and "
285  "*doit_i_field*." );
286  if( doit_i_field.nshelves() != nlat )
287  throw runtime_error( "Inconsistency in size between *doit_i_field* and "
288  "latitude part of *cloudbox_limits*." );
289  if( doit_i_field.nbooks() != nlon )
290  throw runtime_error( "Inconsistency in size between *doit_i_field* and "
291  "longitude part of *cloudbox_limits*." );
292  if( doit_i_field.nvitrines() != np )
293  throw runtime_error( "Inconsistency in size between *doit_i_field* and "
294  "pressure part of *cloudbox_limits*." );
295  if( doit_i_field.nlibraries() != nf )
296  throw runtime_error( "Inconsistency in size between *f_grid* and "
297  "*doit_i_field*." );
298  if( min( doit_i_field ) <= -999 )
299  throw runtime_error( "A very large negative value found in *doit_i_field*. "
300  "Has the radiation field been calculated properly?" );
301  if( jacobian_do )
302  throw runtime_error(
303  "This method does not provide any jacobians (jacobian_do must be 0)" );
304  //---------------------------------------------------------------------------
305 
306  // The approach is:
307  // 1. if 3D shrink this to 1D.
308  // 2. Interpolate in pressure (surface comes in here)
309  // 3. Interpolate in zenith angle
310 
311  // Convert rte_pos to grid positions
312  GridPos gp_p, gp_lat, gp_lon;
313  rte_pos2gridpos( gp_p, gp_lat, gp_lon, atmosphere_dim, p_grid, lat_grid,
314  lon_grid, z_field, rte_pos );
315 
316  // Remove lat, lon and azimuth angle dimensions
317  //
318  Tensor4 I4;
319  //
320  if( atmosphere_dim == 1 )
321  { I4 = doit_i_field( joker, joker, 0, 0, joker, 0, joker ); }
322  else
323  {
324  // Not ready!!!
325 
326  // Check and shift lat and lon grid positions
327  Numeric fg = fractional_gp( gp_lat );
328  if( fg < Numeric(cloudbox_limits[2]) ||
329  fg > Numeric(cloudbox_limits[3]) )
330  throw runtime_error( "The given *rte_pos* is outside the cloudbox "
331  "(in latitude dimension)!" );
332  gp_lat.idx -= cloudbox_limits[2];
333  gridpos_upperend_check( gp_lat, nlat );
334  //
335  fg = fractional_gp( gp_lon );
336  if( fg < Numeric(cloudbox_limits[4]) ||
337  fg > Numeric(cloudbox_limits[5]) )
338  throw runtime_error( "The given *rte_pos* is outside the cloudbox "
339  "(in longitude dimension)!" );
340  gp_lon.idx -= cloudbox_limits[4];
341  gridpos_upperend_check( gp_lon, nlon );
342  }
343 
344  // Check pressure dimension
345  Numeric fg = fractional_gp( gp_p );
346  if( fg < Numeric(cloudbox_limits[0])-1e-6 || // 1-e6 to have some margin
347  fg > Numeric(cloudbox_limits[1])+1e-6 ) // for numerical issues
348  throw runtime_error( "The given *rte_pos* is outside the cloudbox "
349  "(in pressure dimension)!" );
350 
351  // Pressure (no interpolation if at top or bottom)
352  //
353  Tensor3 I3;
354  //
355  if( fg <= Numeric(cloudbox_limits[0]) )
356  { I3 = I4(joker,0,joker,joker); }
357  else if( fg >= Numeric(cloudbox_limits[1]) )
358  { I3 = I4(joker,np-1,joker,joker); }
359  else
360  { // Interpolation needed!:
361  // Shift grid position
362  gp_p.idx -= cloudbox_limits[0];
363  gridpos_upperend_check( gp_p, np );
364 
365  // Interpolate in pressure
366  Vector itw(2);
367  interpweights( itw, gp_p );
368  //
369  I3.resize(nf,nza,stokes_dim);
370  //
371  for( Index iv=0; iv<nf; iv++ ) {
372  for( Index iza=0; iza<nza; iza++ ) {
373  for( Index is=0; is<stokes_dim; is++ ) {
374  I3(iv,iza,is) = interp( itw, I4(iv,joker,iza,is), gp_p );
375  } } }
376  }
377 
378  iy.resize(nf,stokes_dim);
379 
380  // Interpolate in zenith angle
381  GridPosPoly gp;
382  gridpos_poly( gp, scat_za_grid, rte_los[0], za_order );
383  Vector itw( za_order+1 );
384  interpweights( itw, gp );
385  for( Index iv=0; iv<nf; iv++ )
386  {
387  for( Index is=0; is<stokes_dim; is++ )
388  {
389  iy(iv,is) = interp( itw, I3(iv,joker,is), gp );
390  }
391  }
392 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
Template functions for general supergeneric ws methods.
bool chk_if_std_blackbody_agenda(Workspace &ws, const Agenda &blackbody_radiation_agenda)
Checks if blackbody_radiation_agenda returns frequency based radiance.
The Agenda class.
Definition: agenda_class.h:44
Declarations having to do with the four output streams.
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
The Vector class.
Definition: matpackI.h:556
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
The Tensor4 class.
Definition: matpackIV.h:383
Index nrows() const
Returns the number of rows.
Definition: matpackVII.cc:62
Linear algebra functions.
Numeric fractional_gp(const GridPos &gp)
fractional_gp
void chk_rte_los(const Index &atmosphere_dim, ConstVectorView rte_los)
chk_rte_los
This file contains basic functions to handle XML data files.
The Tensor7 class.
Definition: matpackVII.h:1931
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Structure to store a grid position.
Definition: interpolation.h:74
This file contains the definition of Array.
The implementation for String, the ARTS string class.
Definition: mystring.h:63
The Tensor3 class.
Definition: matpackIII.h:348
The global header file for ARTS.
Declarations for agendas.
#define max(a, b)
Definition: continua.cc:20461
const Joker joker
Index idx
Definition: interpolation.h:75
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
Index nlibraries() const
Returns the number of libraries.
Definition: matpackVII.cc:32
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation.
void iyInterpCloudboxField2(Matrix &iy, const Tensor7 &doit_i_field, const Vector &rte_pos, const Vector &rte_los, const Index &jacobian_do, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &basics_checked, const Index &cloudbox_checked, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Index &stokes_dim, const Vector &scat_za_grid, const Vector &scat_aa_grid, const Vector &f_grid, const Index &za_order, const Verbosity &)
Definition: m_doit2.cc:220
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
Header file for logic.cc.
Radiative transfer in cloudbox.
This can be used to make arrays out of anything.
Definition: array.h:40
void rte_pos2gridpos(GridPos &gp_p, GridPos &gp_lat, GridPos &gp_lon, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView rte_pos)
rte_pos2gridpos
#define min(a, b)
Definition: continua.cc:20460
Index nshelves() const
Returns the number of shelves.
Definition: matpackVII.cc:44
Index npages() const
Returns the number of pages.
Definition: matpackVII.cc:56
Workspace class.
Definition: workspace_ng.h:47
Structure to store a grid position for higher order interpolation.
Index nvitrines() const
Returns the number of vitrines.
Definition: matpackVII.cc:38
Index ncols() const
Returns the number of columns.
Definition: matpackVII.cc:68
void get_iy(Workspace &ws, Matrix &iy, ConstTensor3View t_field, ConstTensor3View z_field, ConstTensor4View vmr_field, const Index &cloudbox_on, ConstVectorView f_grid, ConstVectorView rte_pos, ConstVectorView rte_los, ConstVectorView rte_pos2, const Agenda &iy_main_agenda)
get_iy
Definition: rte.cc:1046
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5379
void CloudboxGetIncoming2(Workspace &ws, Tensor7 &doit_i_field, const Agenda &iy_main_agenda, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &basics_checked, const Index &cloudbox_checked, const Index &sensor_checked, const Vector &f_grid, const Index &stokes_dim, const String &iy_unit, const Agenda &blackbody_radiation_agenda, const Vector &scat_za_grid, const Vector &scat_aa_grid, const Index &fill_interior, const Verbosity &)
Definition: m_doit2.cc:28
Index nbooks() const
Returns the number of books.
Definition: matpackVII.cc:50
Declaration of functions in rte.cc.
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580
Auxiliary header stuff related to workspace variable groups.