ARTS  2.4.0(git:4fb77825)
m_cloudbox.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2017
2  Patrick Eriksson <patrick.eriksson@chalmers.se>
3  Stefan Buehler <sbuehler@ltu.se>
4  Claudia Emde <claudia.emde@dlr.de>
5  Cory Davis <cory.davis@metservice.com>
6  Jana Mendrok <jana.mendrok@gmail.com>
7  Daniel Kreyling <daniel.kreyling@nict.go.jp>
8  Manfred Brath <manfred.brath@uni-hamburg.de>
9 
10  This program is free software; you can redistribute it and/or modify it
11  under the terms of the GNU General Public License as published by the
12  Free Software Foundation; either version 2, or (at your option) any
13  later version.
14 
15  This program is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with this program; if not, write to the Free Software
22  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
23  USA. */
24 
25 /*===========================================================================
26  === File description
27  ===========================================================================*/
28 
39 /*===========================================================================
40  === External declarations
41  ===========================================================================*/
42 #include <cmath>
43 #include <cstdlib>
44 #include <stdexcept>
45 
46 #include "array.h"
47 #include "arts.h"
48 #include "auto_md.h"
49 #include "check_input.h"
50 #include "cloudbox.h"
51 #include "file.h"
52 #include "gridded_fields.h"
53 #include "interpolation.h"
54 #include "lin_alg.h"
55 #include "logic.h"
56 #include "math_funcs.h"
57 #include "messages.h"
58 #include "microphysics.h"
59 #include "optproperties.h"
60 #include "parameters.h"
61 #include "physics_funcs.h"
62 #include "rte.h"
63 #include "sorting.h"
64 #include "special_interp.h"
65 #include "xml_io.h"
66 
67 extern const Index GFIELD3_P_GRID;
68 extern const Index GFIELD3_LAT_GRID;
69 extern const Index GFIELD3_LON_GRID;
70 extern const Numeric PI;
71 extern const Numeric DEG2RAD;
72 extern const Numeric RAD2DEG;
73 extern const Numeric LAT_LON_MIN;
74 extern const Numeric DENSITY_OF_ICE;
75 
76 /*===========================================================================
77  === The functions (in alphabetical order)
78  ===========================================================================*/
79 
80 /* Workspace method: Doxygen documentation will be auto-generated */
93  const Verbosity&) {
94  cloudbox_on = 0;
96  cloudbox_limits.resize(0);
98  iy_cloudbox_agenda.set_name("iy_cloudbox_agenda");
99  pnd_field.resize(0, 0, 0, 0);
100  // we need to size dpnd_field to be consistent with jacobian_quantities.
101  dpnd_field_dx.resize(jacobian_quantities.nelem());
102  scat_data.resize(0);
103  scat_species.resize(0);
104  // remove scat_data_raw resizing once all scat solvers have been convert to
105  // use of (new-type) scat_data
106  scat_data_raw.resize(0);
107  scat_data_checked = 0;
108  particle_masses.resize(0, 0);
109 }
110 
111 /* Workspace method: Doxygen documentation will be auto-generated */
112 void cloudboxSetAutomatically( // WS Output:
113  //Workspace& /* ws */,
116  //Agenda& iy_cloudbox_agenda,
117  // WS Input:
118  const Index& atmosphere_dim,
119  const Vector& p_grid,
120  const Vector& lat_grid,
121  const Vector& lon_grid,
122  const Tensor4& particle_field,
123  // Control Parameters
124  const ArrayOfIndex& cloudbox_limits_old,
125  const Numeric& cloudbox_margin,
126  const Verbosity& verbosity) {
127  // Check existing WSV
128  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
129  // includes p_grid chk_if_decresasing
131  // Set cloudbox_on
132  cloudbox_on = 1;
133 
134  if (atmosphere_dim > 1) {
135  ostringstream os;
136  os << "cloudboxSetAutomatically not yet available for 2D and 3D cases.";
137  throw runtime_error(os.str());
138  }
139 
140  Index np = p_grid.nelem();
141 
142  // Allocate cloudbox_limits
143  cloudbox_limits.resize(atmosphere_dim * 2);
144 
145  bool cb_preset = (min(cloudbox_limits_old) > -1);
146  if (cb_preset) {
147  if (cloudbox_limits_old.nelem() != atmosphere_dim * 2) {
148  ostringstream os;
149  os << "The array *cloudbox_limits_old* has incorrect length.\n"
150  << "For atmospheric dim. = " << atmosphere_dim
151  << " the length shall be " << atmosphere_dim * 2 << " but it is "
152  << cloudbox_limits_old.nelem() << ".";
153  throw runtime_error(os.str());
154  }
155  }
156 
157  // Initialize boundary counters
158  Index p1 = 0, p2 = 0;
159  if (cloudbox_margin != -1) {
160  if (cb_preset)
161  p1 = cloudbox_limits_old[0] + 1;
162  else
163  p1 = np - 1;
164  }
165  if (cb_preset) {
166  p2 = cloudbox_limits_old[1] - 1;
167  }
168 
169  // OLE: Commented out until code that uses it at the end of this function is commented back in
170  /*
171  if ( atmosphere_dim > 1 )
172  {
173  Index lat1 = particle_field.nrows()-1;
174  Index lat2 = 0;
175  }
176  if ( atmosphere_dim > 2 )
177  {
178  Index lon1 = particle_field.ncols()-1;
179  Index lon2 = 0;
180  }
181 */
182 
183  bool any_not_empty = false;
184 
185  if (!particle_field.empty()) {
186  bool one_not_empty = false;
187  Index nss = particle_field.nbooks();
188 
189  //--------- Start loop over fields in particle_field------------------------
190  for (Index l = 0; l < nss; l++) {
191  //not_empty is set to true, if a single value of particle_field
192  //is unequal zero (and not NaN), i.e. if we actually have some amount of
193  //these scattering species in the atmosphere.
194  chk_scat_species_field(one_not_empty,
195  particle_field(l, joker, joker, joker),
196  "particle_field",
198  p_grid,
199  lat_grid,
200  lon_grid);
201 
202  //if particles found, enter detailed search
203  if (one_not_empty) {
204  any_not_empty = true;
205  find_cloudlimits(p1,
206  p2,
207  particle_field(l, joker, joker, joker),
209  cloudbox_margin);
210  }
211  }
212  }
213 
214  if (any_not_empty || cb_preset) {
215  // decrease lower cb limit by one to ensure that linear interpolation of
216  // particle number densities is possible.
217  p1 = max(p1 - 1, Index(0));
218 
219  Numeric p_margin1;
220 
221  // alter lower cloudbox_limit by cloudbox_margin, using barometric
222  // height formula
223  p_margin1 = barometric_heightformula(p_grid[p1], cloudbox_margin);
224  while ((p_grid[p1] < p_margin1) && (p1 > 0)) p1--;
225  cloudbox_limits[0] = p1;
226 
227  // increase upper cb limit by one to ensure that linear interpolation of
228  // particle number densities is possible.
229  p2 = min(p2 + 1, np - 1);
230  // set upper cloudbox_limit
231  // if cloudbox reaches to the upper most pressure level
232  if (p2 >= np - 1) {
233  CREATE_OUT2;
234  out2 << "The cloud reaches to TOA!\n"
235  << "Check your *particle_field* data, if realistic!\n";
236  }
237  cloudbox_limits[1] = p2;
238 
239  // assert keyword arguments
240 
241  // The pressure in *p1* must be bigger than the pressure in *p2*.
242  assert(p_grid[p1] > p_grid[p2]);
243  // The pressure in *p1* must be larger than the last value in *p_grid*.
244  assert(p_grid[p1] > p_grid[p_grid.nelem() - 1]);
245  // The pressure in *p2* must be smaller than the first value in *p_grid*."
246  assert(p_grid[p2] < p_grid[0]);
247 
248  /*
249  if ( atmosphere_dim >= 2 )
250  {
251  // The latitude in *lat2* must be bigger than the latitude in *lat1*.
252  assert ( lat_grid[lat2] > lat_grid[lat1] );
253  // The latitude in *lat1* must be >= the second value in *lat_grid*.
254  assert ( lat_grid[lat1] >= lat_grid[1] );
255  // The latitude in *lat2* must be <= the next to last value in *lat_grid*.
256  assert ( lat_grid[lat2] <= lat_grid[lat_grid.nelem()-2] );
257  }
258  if ( atmosphere_dim == 3 )
259  {
260  // The longitude in *lon2* must be bigger than the longitude in *lon1*.
261  assert ( lon_grid[lon2] > lon_grid[lon1] );
262  // The longitude in *lon1* must be >= the second value in *lon_grid*.
263  assert ( lon_grid[lon1] >= lon_grid[1] );
264  // The longitude in *lon2* must be <= the next to last value in *lon_grid*.
265  assert ( lon_grid[lon2] <= lon_grid[lon_grid.nelem()-2] );
266  }
267  */
268  }
269 
270  else
271  // if all particle fields are zero at each level and cloudbox was not preset,
272  // switch cloudbox off.
273  {
274  CREATE_OUT0;
275  cloudbox_on = 0;
276  cloudbox_limits[1] =
277  -1; // just for consistency with cloudboxSetAutomatically
278  out0 << "Cloudbox is switched off!\n";
279  }
280 }
281 
282 /* Workspace method: Doxygen documentation will be auto-generated */
283 void cloudboxSetFullAtm( //WS Output
286  // WS Input
287  const Index& atmosphere_dim,
288  const Vector& p_grid,
289  const Vector& lat_grid,
290  const Vector& lon_grid,
291  const Verbosity&) {
292  cloudbox_on = 1;
293  cloudbox_limits.resize(2 * atmosphere_dim);
294 
295  cloudbox_limits[0] = 0;
296  cloudbox_limits[1] = p_grid.nelem() - 1;
297 
298  if (atmosphere_dim > 1) {
299  Index last_lat = lat_grid.nelem() - 1;
300 
301  // find minimum lat_grid point i with lat_grid[i]-lat_grid[0]>=LAT_LON_MIN
302  Index i = 1;
303  while ((i < last_lat - 1) && (lat_grid[i] - lat_grid[0] < LAT_LON_MIN)) i++;
304  if (i == last_lat - 1) {
305  ostringstream os;
306  os << "Can not define lower latitude edge of cloudbox:\n"
307  << "Extend of atmosphere too small. Distance to minimum latitude\n"
308  << "has to be at least " << LAT_LON_MIN << "deg, but only "
309  << lat_grid[i - 1] - lat_grid[0] << " available here.";
310  throw runtime_error(os.str());
311  }
312  cloudbox_limits[2] = i;
313 
314  // find maximum lat_grid point j with lat_grid[-1]-lat_grid[j]>=LAT_LON_MIN
315  // and j>i
316  Index j = last_lat - 1;
317  while ((j > i) && (lat_grid[last_lat] - lat_grid[j] < LAT_LON_MIN)) j--;
318  if (j == i) {
319  ostringstream os;
320  os << "Can not define upper latitude edge of cloudbox:\n"
321  << "Extend of atmosphere too small. Distance to maximum latitude\n"
322  << "has to be at least " << LAT_LON_MIN << "deg, but only "
323  << lat_grid[last_lat] - lat_grid[j + 1] << " available here.";
324  throw runtime_error(os.str());
325  }
326  cloudbox_limits[3] = j;
327  }
328 
329  if (atmosphere_dim > 2) {
330  const Numeric latmax = max(abs(lat_grid[cloudbox_limits[2]]),
332  const Numeric lfac = 1 / cos(DEG2RAD * latmax);
333  Index last_lon = lon_grid.nelem() - 1;
334 
335  // find minimum lon_grid point i with lon_grid[i]-lon_grid[0]>=LAT_LON_MIN/lfac
336  Index i = 1;
337  while ((i < last_lon - 1) &&
338  (lon_grid[i] - lon_grid[0] < LAT_LON_MIN / lfac))
339  i++;
340  if (i == last_lon - 1) {
341  ostringstream os;
342  os << "Can not define lower longitude edge of cloudbox:\n"
343  << "Extend of atmosphere too small. Distance to minimum longitude\n"
344  << "has to be at least " << LAT_LON_MIN / lfac << "deg, but only "
345  << lon_grid[i - 1] - lon_grid[0] << " available here.";
346  throw runtime_error(os.str());
347  }
348  cloudbox_limits[4] = i;
349 
350  // find maximum lon_grid point j with lon_grid[-1]-lon_grid[j]>=LAT_LON_MIN/lfac
351  // and j>i
352  Index j = last_lon - 1;
353  while ((j > i) && (lon_grid[last_lon] - lon_grid[j] < LAT_LON_MIN / lfac))
354  j--;
355  if (j == i) {
356  ostringstream os;
357  os << "Can not define upper longitude edge of cloudbox:\n"
358  << "Extend of atmosphere too small. Distance to maximum longitude\n"
359  << "has to be at least " << LAT_LON_MIN / lfac << "deg, but only "
360  << lon_grid[last_lon] - lon_grid[j + 1] << " available here.";
361  throw runtime_error(os.str());
362  }
363  cloudbox_limits[5] = j;
364  }
365 }
366 
367 /* Workspace method: Doxygen documentation will be auto-generated */
368 void cloudboxSetManually( // WS Output:
371  // WS Input:
372  const Index& atmosphere_dim,
373  const Vector& p_grid,
374  const Vector& lat_grid,
375  const Vector& lon_grid,
376  // Control Parameters:
377  const Numeric& p1,
378  const Numeric& p2,
379  const Numeric& lat1,
380  const Numeric& lat2,
381  const Numeric& lon1,
382  const Numeric& lon2,
383  const Verbosity&) {
384  // Check existing WSV
385  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
387 
388  // Check keyword arguments
389  if (p1 <= p2)
390  throw runtime_error(
391  "The pressure in *p1* must be bigger than the "
392  "pressure in *p2*.");
393  if (p1 <= p_grid[p_grid.nelem() - 1])
394  throw runtime_error(
395  "The pressure in *p1* must be larger than the "
396  "last value in *p_grid*.");
397  if (p2 >= p_grid[0])
398  throw runtime_error(
399  "The pressure in *p2* must be smaller than the "
400  "first value in *p_grid*.");
401  if (atmosphere_dim >= 2) {
402  if (lat2 <= lat1)
403  throw runtime_error(
404  "The latitude in *lat2* must be bigger than the "
405  "latitude in *lat1*.");
406  if (lat1 < lat_grid[1])
407  throw runtime_error(
408  "The latitude in *lat1* must be >= the "
409  "second value in *lat_grid*.");
410  if (lat2 > lat_grid[lat_grid.nelem() - 2])
411  throw runtime_error(
412  "The latitude in *lat2* must be <= the "
413  "next to last value in *lat_grid*.");
414  }
415  if (atmosphere_dim == 3) {
416  if (lon2 <= lon1)
417  throw runtime_error(
418  "The longitude in *lon2* must be bigger than the "
419  "longitude in *lon1*.");
420  if (lon1 < lon_grid[1])
421  throw runtime_error(
422  "The longitude in *lon1* must be >= the "
423  "second value in *lon_grid*.");
424  if (lon2 > lon_grid[lon_grid.nelem() - 2])
425  throw runtime_error(
426  "The longitude in *lon2* must be <= the "
427  "next to last value in *lon_grid*.");
428  }
429 
430  // Set cloudbox_on
431  cloudbox_on = 1;
432 
433  // Allocate cloudbox_limits
434  cloudbox_limits.resize(atmosphere_dim * 2);
435 
436  // Pressure limits
437  if (p1 > p_grid[1]) {
438  cloudbox_limits[0] = 0;
439  } else {
440  for (cloudbox_limits[0] = 1; p_grid[cloudbox_limits[0] + 1] >= p1;
441  cloudbox_limits[0]++) {
442  }
443  }
444  if (p2 < p_grid[p_grid.nelem() - 2]) {
445  cloudbox_limits[1] = p_grid.nelem() - 1;
446  } else {
447  for (cloudbox_limits[1] = p_grid.nelem() - 2;
448  p_grid[cloudbox_limits[1] - 1] <= p2;
449  cloudbox_limits[1]--) {
450  }
451  }
452 
453  // Latitude limits
454  if (atmosphere_dim >= 2) {
455  for (cloudbox_limits[2] = 1; lat_grid[cloudbox_limits[2] + 1] <= lat1;
456  cloudbox_limits[2]++) {
457  }
458  for (cloudbox_limits[3] = lat_grid.nelem() - 2;
459  lat_grid[cloudbox_limits[3] - 1] >= lat2;
460  cloudbox_limits[3]--) {
461  }
462  }
463 
464  // Longitude limits
465  if (atmosphere_dim == 3) {
466  for (cloudbox_limits[4] = 1; lon_grid[cloudbox_limits[4] + 1] <= lon1;
467  cloudbox_limits[4]++) {
468  }
469  for (cloudbox_limits[5] = lon_grid.nelem() - 2;
470  lon_grid[cloudbox_limits[5] - 1] >= lon2;
471  cloudbox_limits[5]--) {
472  }
473  }
474 }
475 
476 /* Workspace method: Doxygen documentation will be auto-generated */
477 void cloudboxSetManuallyAltitude( // WS Output:
480  // WS Input:
481  const Index& atmosphere_dim,
482  const Tensor3& z_field,
483  const Vector& lat_grid,
484  const Vector& lon_grid,
485  // Control Parameters:
486  const Numeric& z1,
487  const Numeric& z2,
488  const Numeric& lat1,
489  const Numeric& lat2,
490  const Numeric& lon1,
491  const Numeric& lon2,
492  const Verbosity&) {
493  // Check existing WSV
494  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
495 
496  // Check keyword arguments
497  if (z1 >= z2)
498  throw runtime_error(
499  "The altitude in *z1* must be smaller than the "
500  "altitude in *z2*.");
501  /* These cases are in fact handled
502  if( z1 <= z_field(0, 0, 0) )
503  throw runtime_error( "The altitude in *z1* must be larger than the "
504  "first value in *z_field*." );
505  if( z2 >= z_field(z_field.npages()-1, 0, 0) )
506  throw runtime_error( "The altitude in *z2* must be smaller than the "
507  "last value in *z_field*." );
508  */
509  if (atmosphere_dim == 3) {
510  if (lat2 <= lat1)
511  throw runtime_error(
512  "The latitude in *lat2* must be bigger than the "
513  " latitude in *lat1*.");
514  if (lat1 < lat_grid[1])
515  throw runtime_error(
516  "The latitude in *lat1* must be >= the "
517  "second value in *lat_grid*.");
518  if (lat2 > lat_grid[lat_grid.nelem() - 2])
519  throw runtime_error(
520  "The latitude in *lat2* must be <= the "
521  "next to last value in *lat_grid*.");
522  if (lon2 <= lon1)
523  throw runtime_error(
524  "The longitude in *lon2* must be bigger than the "
525  "longitude in *lon1*.");
526  if (lon1 < lon_grid[1])
527  throw runtime_error(
528  "The longitude in *lon1* must be >= the "
529  "second value in *lon_grid*.");
530  if (lon2 > lon_grid[lon_grid.nelem() - 2])
531  throw runtime_error(
532  "The longitude in *lon2* must be <= the "
533  "next to last value in *lon_grid*.");
534  }
535 
536  // Set cloudbox_on
537  cloudbox_on = 1;
538 
539  // Allocate cloudbox_limits
540  cloudbox_limits.resize(atmosphere_dim * 2);
541 
542  // Pressure/altitude limits
543  if (z1 < z_field(1, 0, 0)) {
544  cloudbox_limits[0] = 0;
545  } else {
546  for (cloudbox_limits[0] = 1; z_field(cloudbox_limits[0] + 1, 0, 0) <= z1;
547  cloudbox_limits[0]++) {
548  }
549  }
550  if (z2 > z_field(z_field.npages() - 2, 0, 0)) {
551  cloudbox_limits[1] = z_field.npages() - 1;
552  } else {
553  for (cloudbox_limits[1] = z_field.npages() - 2;
554  z_field(cloudbox_limits[1] - 1, 0, 0) >= z2;
555  cloudbox_limits[1]--) {
556  }
557  }
558 
559  // Latitude limits
560  if (atmosphere_dim >= 2) {
561  for (cloudbox_limits[2] = 1; lat_grid[cloudbox_limits[2] + 1] <= lat1;
562  cloudbox_limits[2]++) {
563  }
564  for (cloudbox_limits[3] = lat_grid.nelem() - 2;
565  lat_grid[cloudbox_limits[3] - 1] >= lat2;
566  cloudbox_limits[3]--) {
567  }
568  }
569 
570  // Longitude limits
571  if (atmosphere_dim == 3) {
572  for (cloudbox_limits[4] = 1; lon_grid[cloudbox_limits[4] + 1] <= lon1;
573  cloudbox_limits[4]++) {
574  }
575  for (cloudbox_limits[5] = lon_grid.nelem() - 2;
576  lon_grid[cloudbox_limits[5] - 1] >= lon2;
577  cloudbox_limits[5]--) {
578  }
579  }
580 }
581 
582 /* Workspace method: Doxygen documentation will be auto-generated */
584  const Tensor7& cloudbox_field,
585  const Vector& rte_pos,
586  const Vector& rte_los,
587  const Index& jacobian_do,
588  const Index& cloudbox_on,
590  const Index& atmosphere_dim,
591  const Vector& p_grid,
592  const Vector& lat_grid,
593  const Vector& lon_grid,
594  const Tensor3& z_field,
595  const Matrix& z_surface,
596  const Index& stokes_dim,
597  const Vector& za_grid,
598  const Vector& aa_grid,
599  const Vector& f_grid,
600  //const Index& p_interp_order,
601  const Index& za_interp_order,
602  const Index& za_restrict,
603  const Index& cos_za_interp,
604  const Numeric& za_extpolfac,
605  const Index& aa_interp_order,
606  const Verbosity&) {
607  //--- Check input -----------------------------------------------------------
608  if (!(atmosphere_dim == 1 || atmosphere_dim == 3))
609  throw runtime_error("The atmospheric dimensionality must be 1 or 3.");
610  if (jacobian_do)
611  throw runtime_error(
612  "This method does not support jacobians (jacobian_do must be 0)");
613  if (!cloudbox_on)
614  throw runtime_error(
615  "The cloud box is not activated and no outgoing "
616  "field can be returned.");
617  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
618  throw runtime_error(
619  "*cloudbox_limits* is a vector which contains the upper and lower\n"
620  "limit of the cloud for all atmospheric dimensions.\n"
621  "So its length must be 2 x *atmosphere_dim*");
622  if (za_grid.nelem() == 0)
623  throw runtime_error(
624  "The variable *za_grid* is empty. Are dummy "
625  "values from *cloudboxOff used?");
626  if (!(za_interp_order < za_grid.nelem()))
627  throw runtime_error(
628  "Zenith angle interpolation order *za_interp_order*"
629  " must be smaller\n"
630  "than number of angles in *za_grid*.");
631  if (atmosphere_dim > 1 && !(aa_interp_order < aa_grid.nelem()))
632  throw runtime_error(
633  "Azimuth angle interpolation order *aa_interp_order*"
634  " must be smaller\n"
635  "than number of angles in *aa_grid*.");
636  if (cloudbox_field.nlibraries() != f_grid.nelem())
637  throw runtime_error(
638  "Inconsistency in size between f_grid and cloudbox_field! "
639  "(This method does not yet handle dispersion type calculations.)");
640  //---------------------------------------------------------------------------
641 
642  // At each lat/lon, one value of the intensity field is defined at the
643  // surface. Simplest way to handle this is to insert z_surface into z_field
644  Tensor3 z_with_surface = z_field;
645  for (Index ilat = 0; ilat < z_surface.nrows(); ilat++) {
646  for (Index ilon = 0; ilon < z_surface.ncols(); ilon++) {
647  Index ip = 0;
648  while (z_surface(ilat, ilon) >= z_field(ip + 1, ilat, ilon)) {
649  ip++;
650  }
651  z_with_surface(ip, ilat, ilon) = z_surface(ilat, ilon);
652  }
653  }
654 
655  // Convert rte_pos to grid positions
656  GridPos gp_p, gp_lat, gp_lon;
657  rte_pos2gridpos(gp_p,
658  gp_lat,
659  gp_lon,
661  p_grid,
662  lat_grid,
663  lon_grid,
664  z_with_surface,
665  rte_pos);
666 
667  //--- Determine if at border or inside of cloudbox (or outside!)
668  //
669  // Let us introduce a number coding for cloud box borders.
670  // Borders have the same number as position in *cloudbox_limits*.
671  // Inside cloud box is coded as 99, and outside as > 100.
672  Index border = 999;
673  //
674  //- Check if at any border
675  if (is_gridpos_at_index_i(gp_p, cloudbox_limits[0], false)) {
676  border = 0;
677  } else if (is_gridpos_at_index_i(gp_p, cloudbox_limits[1], false)) {
678  border = 1;
679  } else if (atmosphere_dim > 1) {
680  if (is_gridpos_at_index_i(gp_lat, cloudbox_limits[2], false)) {
681  border = 2;
682  } else if (is_gridpos_at_index_i(gp_lat, cloudbox_limits[3], false)) {
683  border = 3;
684  } else if (atmosphere_dim > 2) {
685  if (is_gridpos_at_index_i(gp_lon, cloudbox_limits[4], false)) {
686  border = 4;
687  } else if (is_gridpos_at_index_i(gp_lon, cloudbox_limits[5], false)) {
688  border = 5;
689  }
690  }
691  }
692 
693  //
694  //- Check if inside (when border<100 here, it means we are on a cloudbox border)
695  if (border > 100) {
696  // Assume inside as it is easiest to detect if outside (can be detected
697  // check in one dimension at the time)
698  bool inside = true;
699  Numeric fgp;
700 
701  // Check in pressure dimension
702  fgp = fractional_gp(gp_p);
703  if (fgp < Numeric(cloudbox_limits[0]) ||
704  fgp > Numeric(cloudbox_limits[1])) {
705  inside = false;
706  }
707 
708  // Check in lat and lon dimensions
709  if (atmosphere_dim == 3 && inside) {
710  fgp = fractional_gp(gp_lat);
711  if (fgp < Numeric(cloudbox_limits[2]) ||
712  fgp > Numeric(cloudbox_limits[3])) {
713  inside = false;
714  }
715  fgp = fractional_gp(gp_lon);
716  if (fgp < Numeric(cloudbox_limits[4]) ||
717  fgp > Numeric(cloudbox_limits[5])) {
718  inside = false;
719  }
720  }
721 
722  if (inside) {
723  border = 99;
724  }
725  }
726 
727  // If outside, something is wrong
728  if (border > 100) {
729  throw runtime_error(
730  "Given position has been found to be outside the cloudbox.");
731  }
732 
733  //- Sizes
734  const Index nf = f_grid.nelem();
735  DEBUG_ONLY(const Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1);
736  const Index nza = za_grid.nelem();
737  const Index naa = cloudbox_field.nrows();
738 
739  //- Resize *iy*
740  iy.resize(nf, stokes_dim);
741  //- Make space for spatially interpolated field (we perfrom angle
742  //interpolation separately on this then)
743  Tensor4 i_field_local(nf, nza, naa, stokes_dim);
744 
745  // Index of the p/lat/lon slice (first or last element slice in this
746  // dimension), where rte_pos is located. if located on a border.
747  Index border_index;
748  if (border < 99) {
749  if (border % 2) // odd border id, ie top or north or east
750  border_index = cloudbox_limits[border] - cloudbox_limits[border - 1];
751  else // even border id, ie bottom or south or west
752  border_index = 0;
753  }
754 
755  // Sensor points inside the cloudbox
756  if (border == 99) {
757  if (atmosphere_dim == 3) {
758  throw runtime_error(
759  "Radiation extraction for a position inside cloudbox\n"
760  "is not yet implemented for 3D cases.\n");
761  } else {
762  assert(atmosphere_dim == 1);
763 
764  assert(is_size(cloudbox_field, nf, np, 1, 1, nza, 1, stokes_dim));
765 
766  // Grid position in *p_grid*
767  gp_p.idx = gp_p.idx - cloudbox_limits[0];
769  Vector itw_p(2);
770  interpweights(itw_p, gp_p);
771 
772  for (Index is = 0; is < stokes_dim; is++)
773  for (Index iv = 0; iv < nf; iv++)
774  for (Index i_za = 0; i_za < nza; i_za++)
775  i_field_local(iv, i_za, 0, is) = interp(
776  itw_p, cloudbox_field(iv, joker, 0, 0, i_za, 0, is), gp_p);
777  }
778  }
779 
780  // Sensor outside the cloudbox
781 
782  // --- 1D ------------------------------------------------------------------
783  else if (atmosphere_dim == 1) {
784  i_field_local =
785  cloudbox_field(joker, border_index, 0, 0, joker, joker, joker);
786  }
787 
788  // --- 3D ------------------------------------------------------------------
789  else {
790  assert(is_size(cloudbox_field,
791  nf,
792  cloudbox_field.nvitrines(),
793  cloudbox_field.nshelves(),
794  cloudbox_field.nbooks(),
795  za_grid.nelem(),
796  aa_grid.nelem(),
797  stokes_dim));
798 
799  // Interpolation weights (for 2D "red" interpolation)
800  Vector itw(4);
801 
802  // Outgoing from cloudbox top or bottom, i.e. from a pressure level
803  if (border <= 1) {
804  // Lat and lon grid positions with respect to cloud box
805  GridPos cb_gp_lat, cb_gp_lon;
806  cb_gp_lat = gp_lat;
807  cb_gp_lon = gp_lon;
808  cb_gp_lat.idx -= cloudbox_limits[2];
809  cb_gp_lon.idx -= cloudbox_limits[4];
810  //
811  gridpos_upperend_check(cb_gp_lat,
813  gridpos_upperend_check(cb_gp_lon,
815 
816  interpweights(itw, cb_gp_lat, cb_gp_lon);
817 
818  for (Index is = 0; is < stokes_dim; is++)
819  for (Index iv = 0; iv < nf; iv++)
820  for (Index i_za = 0; i_za < nza; i_za++)
821  for (Index i_aa = 0; i_aa < naa; i_aa++)
822  i_field_local(iv, i_za, i_aa, is) =
823  interp(itw,
825  iv, border_index, joker, joker, i_za, i_aa, is),
826  cb_gp_lat,
827  cb_gp_lon);
828  }
829 
830  // Outgoing from cloudbox north or south border, i.e. from a latitude level
831  else if (border <= 3) {
832  // Pressure and lon grid positions with respect to cloud box
833  GridPos cb_gp_p, cb_gp_lon;
834  cb_gp_p = gp_p;
835  cb_gp_lon = gp_lon;
836  cb_gp_p.idx -= cloudbox_limits[0];
837  cb_gp_lon.idx -= cloudbox_limits[4];
838  //
840  gridpos_upperend_check(cb_gp_lon,
842 
843  interpweights(itw, cb_gp_p, cb_gp_lon);
844 
845  for (Index is = 0; is < stokes_dim; is++)
846  for (Index iv = 0; iv < nf; iv++)
847  for (Index i_za = 0; i_za < nza; i_za++)
848  for (Index i_aa = 0; i_aa < naa; i_aa++)
849  i_field_local(iv, i_za, i_aa, is) =
850  interp(itw,
852  iv, joker, border_index, joker, i_za, i_aa, is),
853  cb_gp_p,
854  cb_gp_lon);
855  }
856 
857  // Outgoing from cloudbox east or west border, i.e. from a longitude level
858  else {
859  // Pressure and lat grid positions with respect to cloud box
860  GridPos cb_gp_p, cb_gp_lat;
861  cb_gp_p = gp_p;
862  cb_gp_lat = gp_lat;
863  cb_gp_p.idx -= cloudbox_limits[0];
864  cb_gp_lat.idx -= cloudbox_limits[2];
865  //
867  gridpos_upperend_check(cb_gp_lat,
869 
870  interpweights(itw, cb_gp_p, cb_gp_lat);
871 
872  for (Index is = 0; is < stokes_dim; is++)
873  for (Index iv = 0; iv < nf; iv++)
874  for (Index i_za = 0; i_za < nza; i_za++)
875  for (Index i_aa = 0; i_aa < naa; i_aa++)
876  i_field_local(iv, i_za, i_aa, is) =
877  interp(itw,
879  iv, joker, joker, border_index, i_za, i_aa, is),
880  cb_gp_p,
881  cb_gp_lat);
882  }
883  }
884 
885  // now, do Nth-oder polynomial interpoaltion of angle(s)
886  //
887  // a bunch of options for testing:
888  // a) interpolation over full zenith angle grid vs over hemispheres
889  // separately.
890  // b) interpolation in plain zenith angles vs. in cosines of zenith angle.
891 
892  // find range of za_grid that we will do interpolation over.
893  Index za_start = 0;
894  Index za_extend = za_grid.nelem();
895  if (za_restrict) {
896  // which hemisphere do we need?
897  if (is_same_within_epsilon(rte_los[0], 90., 1e-6)) {
898  //FIXME: we should allow this though, if za_grid has a grid point
899  //at 90deg.
900  throw runtime_error(
901  "Hemisphere-restricted zenith angle interpolation not allowed\n"
902  "for 90degree views.");
903  } else if (rte_los[0] > 90) {
904  // upwelling, i.e. second part of za_grid. that is, we need to find
905  // the first point in za_grid where za>90. and update za_start
906  // accordingly.
907  while (za_start < za_grid.nelem() && za_grid[za_start] < 90.) {
908  za_start++;
909  }
910  if (za_start == za_grid.nelem())
911  throw runtime_error(
912  "No za_grid grid point found in 90-180deg hemisphere.\n"
913  "No hemispheric interpolation possible.");
914  za_extend -= za_start;
915  } else {
916  // downwelling, i.e. first part of za_grid. that is, we need to
917  // find the last point in za_grid where za<90. and update za_extend
918  // accordingly.
919  while (za_extend > 0 && za_grid[za_extend - 1] > 90.) {
920  za_extend--;
921  }
922  if (za_extend == 0)
923  throw runtime_error(
924  "No za_grid grid point found in 0-90deg hemisphere.\n"
925  "No hemispheric interpolation possible.");
926  }
927  if (!(za_interp_order < za_extend)) {
928  ostringstream os;
929  os << "Zenith angle interpolation order *za_interp_order* ("
930  << za_interp_order << ") must be smaller\n"
931  << "than number of angles in respective hemisphere (" << za_extend
932  << ").";
933  throw runtime_error(os.str());
934  }
935  }
936 
937  // Grid position in *za_grid*
938  GridPosPoly gp_za, gp_aa;
939 
940  if (cos_za_interp) {
941  Vector cosza_grid(za_extend);
942  const Numeric cosza = cos(rte_los[0] * DEG2RAD);
943 
944  for (Index i_za = 0; i_za < za_extend; i_za++)
945  cosza_grid[i_za] = cos(za_grid[i_za + za_start] * DEG2RAD);
946 
947  // OBS: cosza is a decreasing grid!
948  const Numeric cosza_min =
949  cosza_grid[0] + za_extpolfac * (cosza_grid[0] - cosza_grid[1]);
950  if (cosza > cosza_min) {
951  ostringstream os;
952  os << "Zenith angle " << rte_los[0] << "deg is outside the range"
953  << " covered by za_grid.\n"
954  << "Lower limit of allowed range is " << acos(cosza_min) * RAD2DEG
955  << ".\n"
956  << "Increase za_extpolfac (now=" << za_extpolfac << ") or"
957  << " use wider za_grid.\n";
958  throw runtime_error(os.str());
959  }
960  const Numeric cosza_max =
961  cosza_grid[za_extend - 1] -
962  za_extpolfac * (cosza_grid[za_extend - 2] - cosza_grid[za_extend - 1]);
963  if (cosza < cosza_max) {
964  ostringstream os;
965  os << "Zenith angle " << rte_los[0] << "deg is outside the range"
966  << " covered by za_grid.\n"
967  << "Upper limit of allowed range is " << acos(cosza_max) * RAD2DEG
968  << ".\n"
969  << "Increase za_extpolfac (now=" << za_extpolfac << ") or"
970  << " use wider za_grid.\n";
971  throw runtime_error(os.str());
972  }
973 
974  gridpos_poly(gp_za, cosza_grid, cosza, za_interp_order, za_extpolfac);
975  } else {
976  Vector za_g = za_grid[Range(za_start, za_extend)];
977  const Numeric za = rte_los[0];
978 
979  const Numeric za_min = za_g[0] - za_extpolfac * (za_g[1] - za_g[0]);
980  if (za < za_min) {
981  ostringstream os;
982  os << "Zenith angle " << rte_los[0] << "deg is outside the range"
983  << " covered by za_grid.\n"
984  << "Lower limit of allowed range is " << za_min << ".\n"
985  << "Increase za_extpolfac (now=" << za_extpolfac << ") or"
986  << " use wider za_grid.\n";
987  throw runtime_error(os.str());
988  }
989  const Numeric za_max =
990  za_g[za_g.nelem() - 1] +
991  za_extpolfac * (za_g[za_extend - 1] - za_g[za_extend - 2]);
992  if (za > za_max) {
993  ostringstream os;
994  os << "Zenith angle " << za << "deg is outside the range"
995  << " covered by za_grid.\n"
996  << "Upper limit of allowed range is " << za_max << ".\n"
997  << "Increase za_extpolfac (now=" << za_extpolfac << ") or"
998  << " use wider za_grid.\n";
999  throw runtime_error(os.str());
1000  }
1001  gridpos_poly(gp_za, za_g, za, za_interp_order, za_extpolfac);
1002  }
1003 
1004  if (atmosphere_dim > 1) {
1006  gp_aa, aa_grid, rte_los[1], aa_interp_order);
1007  } else {
1008  gp_aa.idx.resize(1);
1009  gp_aa.idx[0] = 0;
1010  gp_aa.w.resize(1);
1011  gp_aa.w[0] = 1;
1012  }
1013 
1014  // Corresponding interpolation weights
1015  Vector itw_angs(gp_za.idx.nelem() * gp_aa.idx.nelem());
1016  interpweights(itw_angs, gp_za, gp_aa);
1017 
1018  for (Index is = 0; is < stokes_dim; is++) {
1019  for (Index iv = 0; iv < nf; iv++) {
1020  iy(iv, is) =
1021  interp(itw_angs,
1022  i_field_local(iv, Range(za_start, za_extend), joker, is),
1023  gp_za,
1024  gp_aa);
1025  }
1026  }
1027 }
1028 
1029 /* Workspace method: Doxygen documentation will be auto-generated */
1032  const Index& atmosphere_dim,
1033  const Index& cloudbox_on,
1034  const Index& new_limit0,
1035  const Index& new_limit1,
1036  const Index& new_limit2,
1037  const Index& new_limit3,
1038  const Index& new_limit4,
1039  const Index& new_limit5,
1040  const Verbosity&) {
1041  if (!cloudbox_on)
1042  throw runtime_error("No need to use this method with cloudbox=0.");
1043  if (new_limit0 < cloudbox_limits[0])
1044  throw runtime_error("new_limits0 < cloudbox_limits[0], not allowed!");
1045  if (new_limit1 > cloudbox_limits[1])
1046  throw runtime_error("new_limits1 > cloudbox_limits[1], not allowed!");
1047 
1048  Tensor7 fcopy = cloudbox_field;
1049 
1050  if (atmosphere_dim == 1) {
1051  cloudbox_field = fcopy(
1052  joker,
1053  Range(new_limit0 - cloudbox_limits[0], new_limit1 - new_limit0 + 1),
1054  joker,
1055  joker,
1056  joker,
1057  joker,
1058  joker);
1059  cloudbox_limits[0] = new_limit0;
1060  cloudbox_limits[1] = new_limit1;
1061  } else {
1062  if (new_limit2 < cloudbox_limits[2])
1063  throw runtime_error("new_limits2 < cloudbox_limits[2], not allowed!");
1064  if (new_limit3 > cloudbox_limits[3])
1065  throw runtime_error("new_limits3 > cloudbox_limits[3], not allowed!");
1066 
1067  if (atmosphere_dim == 2) {
1068  cloudbox_field = fcopy(
1069  joker,
1070  Range(new_limit0 - cloudbox_limits[0], new_limit1 - new_limit0 + 1),
1071  Range(new_limit2 - cloudbox_limits[2], new_limit3 - new_limit2 - 1),
1072  joker,
1073  joker,
1074  joker,
1075  joker);
1076  cloudbox_limits[0] = new_limit0;
1077  cloudbox_limits[1] = new_limit1;
1078  cloudbox_limits[2] = new_limit2;
1079  cloudbox_limits[3] = new_limit3;
1080  } else {
1081  if (new_limit4 < cloudbox_limits[4])
1082  throw runtime_error("new_limits4 < cloudbox_limits[4], not allowed!");
1083  if (new_limit5 > cloudbox_limits[5])
1084  throw runtime_error("new_limits5 > cloudbox_limits[5], not allowed!");
1085  cloudbox_field = fcopy(
1086  joker,
1087  Range(new_limit0 - cloudbox_limits[0], new_limit1 - new_limit0 + 1),
1088  Range(new_limit2 - cloudbox_limits[2], new_limit3 - new_limit2 + 1),
1089  Range(new_limit4 - cloudbox_limits[4], new_limit5 - new_limit4 + 1),
1090  joker,
1091  joker,
1092  joker);
1093  cloudbox_limits[0] = new_limit0;
1094  cloudbox_limits[1] = new_limit1;
1095  cloudbox_limits[2] = new_limit2;
1096  cloudbox_limits[3] = new_limit3;
1097  cloudbox_limits[4] = new_limit4;
1098  cloudbox_limits[5] = new_limit5;
1099  }
1100  }
1101 }
1102 
1103 /* Workspace method: Doxygen documentation will be auto-generated */
1104 void particle_fieldCleanup( //WS Output:
1105  Tensor4& particle_field_out,
1106  //WS Input:
1107  const Tensor4& particle_field_in,
1108  const Numeric& threshold,
1109  const Verbosity&) {
1110  if (&particle_field_out != &particle_field_in) {
1111  particle_field_out = particle_field_in;
1112  }
1113 
1114  // Check that particle_field contains realistic values
1115  //(values smaller than the threshold will be set to 0)
1116  for (Index i = 0; i < particle_field_out.nbooks(); i++) {
1117  for (Index j = 0; j < particle_field_out.npages(); j++) {
1118  for (Index k = 0; k < particle_field_out.nrows(); k++) {
1119  for (Index l = 0; l < particle_field_out.ncols(); l++) {
1120  if (particle_field_out(i, j, k, l) < threshold) {
1121  particle_field_out(i, j, k, l) = 0.0;
1122  }
1123  }
1124  }
1125  }
1126  }
1127 }
1128 
1129 /* Workspace method: Doxygen documentation will be auto-generated */
1130 void ScatSpeciesInit( //WS Output:
1136  const Verbosity&) {
1137  scat_species.resize(0);
1138  scat_data_raw.resize(0);
1139  scat_meta.resize(0);
1140  pnd_field_raw.resize(0);
1141  scat_data_checked = 0;
1142 }
1143 
1144 /* Workspace method: Doxygen documentation will be auto-generated */
1145 void ScatElementsPndAndScatAdd( //WS Output:
1148  // WS Input (needed for checking the datafiles):
1149  const Index& atmosphere_dim,
1150  // Keywords:
1151  const ArrayOfString& scat_data_files,
1152  const ArrayOfString& pnd_field_files,
1153  const Verbosity& verbosity) {
1154  CREATE_OUT2;
1155 
1156  //--- Check input ---------------------------------------------------------
1157 
1158  // Atmosphere
1159  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1160  //chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
1161 
1162  //--- Reading the data ---------------------------------------------------
1163 
1164  if (scat_data_files.nelem() != pnd_field_files.nelem()) {
1165  ostringstream os;
1166  os << "Number of elements in scat_data and pnd_field filelists is"
1167  << "inconsistent.";
1168  throw runtime_error(os.str());
1169  }
1170 
1171  Index last_species = scat_data_raw.nelem() - 1;
1172  if (last_species == -1) {
1173  scat_data_raw.resize(1);
1174  last_species = 0;
1175  }
1176 
1177  // Create empty dummy elements to append to *scat_data_raw* and *pnd_field_raw*.
1179  GriddedField3 pnd_field_data;
1180 
1181  for (Index i = 0; i < scat_data_files.nelem(); i++) {
1182  // Append *scat_data_raw* and *pnd_field_raw* with empty Arrays of Tensors.
1183  scat_data_raw[last_species].push_back(scat_data_single);
1184  pnd_field_raw.push_back(pnd_field_data);
1185 
1186  out2 << " Read single scattering data file " << scat_data_files[i] << "\n";
1188  scat_data_files[i],
1189  scat_data_raw[last_species][scat_data_raw[last_species].nelem() - 1],
1190  verbosity);
1191 
1192  out2 << " Read particle number density field\n";
1193  if (pnd_field_files[i].nelem() < 1) {
1194  CREATE_OUT1;
1195  out1 << "Warning: No pnd_field_file specified. Ignored here,\n"
1196  << "but user HAS TO add that later on!. \n";
1197  } else {
1198  xml_read_from_file(pnd_field_files[i],
1199  pnd_field_raw[pnd_field_raw.nelem() - 1],
1200  verbosity);
1201 
1203  pnd_field_files[i],
1205  verbosity);
1206  }
1207  }
1208 }
1209 
1210 /* Workspace method: Doxygen documentation will be auto-generated */
1211 void ScatSpeciesPndAndScatAdd( //WS Output:
1214  // WS Input(needed for checking the datafiles):
1215  const Index& atmosphere_dim,
1216  // Keywords:
1217  const ArrayOfString& scat_data_files,
1218  const String& pnd_fieldarray_file,
1219  const Verbosity& verbosity) {
1220  CREATE_OUT2;
1221 
1222  //--- Check input ---------------------------------------------------------
1223 
1224  // Atmosphere
1225  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1226  //chk_atm_grids ( atmosphere_dim, p_grid, lat_grid, lon_grid );
1227 
1228  //--- Reading the data ---------------------------------------------------
1230  arr_ssd.resize(scat_data_files.nelem());
1231 
1232  for (Index i = 0; i < scat_data_files.nelem(); i++) {
1233  out2 << " Read single scattering data file " << scat_data_files[i] << "\n";
1234  xml_read_from_file(scat_data_files[i], arr_ssd[i], verbosity);
1235  }
1236 
1237  // append as new scattering species
1238  if (scat_data_raw.nelem() == 0) {
1239  scat_data_raw.resize(1);
1240  scat_data_raw[0] = arr_ssd;
1241  } else
1242  scat_data_raw.push_back(arr_ssd);
1243 
1244  out2 << " Read particle number density data \n";
1245  ArrayOfGriddedField3 pnd_tmp;
1246  xml_read_from_file(pnd_fieldarray_file, pnd_tmp, verbosity);
1247 
1248  chk_pnd_raw_data(pnd_tmp, pnd_fieldarray_file, atmosphere_dim, verbosity);
1249 
1250  // append to pnd_field_raw
1251  for (Index i = 0; i < pnd_tmp.nelem(); ++i)
1252  pnd_field_raw.push_back(pnd_tmp[i]);
1253 }
1254 
1255 /* Workspace method: Doxygen documentation will be auto-generated */
1262  // WS Input (needed for checking the datafiles):
1263  const Index& atmosphere_dim,
1264  const Vector& f_grid,
1265  // Keywords:
1266  const ArrayOfString& scat_data_files,
1267  const ArrayOfString& pnd_field_files,
1268  const Verbosity& verbosity) {
1269  CREATE_OUT2;
1270 
1271  //--- Check input ---------------------------------------------------------
1272 
1273  // Atmosphere
1274  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1275  //chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
1276 
1277  // Frequency grid
1278  if (f_grid.empty()) throw runtime_error("The frequency grid is empty.");
1279  chk_if_increasing("f_grid", f_grid);
1280 
1281  //--- Reading the data ---------------------------------------------------
1282 
1283  if (scat_data_files.nelem() != pnd_field_files.nelem()) {
1284  ostringstream os;
1285  os << "Number of elements in scat_data and pnd_field filelists is"
1286  << "inconsistent.";
1287  throw runtime_error(os.str());
1288  }
1289 
1290  Index last_species = scat_data_raw.nelem() - 1;
1291  if (last_species == -1) {
1292  scat_data_raw.resize(1);
1293  last_species = 0;
1294  }
1295 
1296  // Create empty dummy elements to append to *scat_data_raw* and *pnd_field_raw*.
1298  GriddedField3 pnd_field_data;
1299  ArrayOfString species(1);
1300  species[0] = "particles";
1301 
1302  for (Index i = 0; i < scat_data_files.nelem(); i++) {
1303  // Append *scat_data_raw* and *pnd_field_raw* with empty Arrays of Tensors.
1304  scat_data_raw[last_species].push_back(scat_data_single);
1305  vmr_field_raw.push_back(pnd_field_data);
1306 
1307  out2 << " Read single scattering data file " << scat_data_files[i] << "\n";
1309  scat_data_files[i],
1310  scat_data_raw[last_species][scat_data_raw[last_species].nelem() - 1],
1311  verbosity);
1312 
1313  out2 << " Check single scattering properties\n";
1315  "scat_data_single.f_grid to f_grid",
1316  scat_data_raw[last_species][scat_data_raw[last_species].nelem() - 1]
1317  .f_grid,
1318  f_grid);
1319 
1320  out2 << " Read particle number density field\n";
1321  if (pnd_field_files[i].nelem() < 1) {
1322  CREATE_OUT1;
1323  out1 << "Warning: No pnd_field_file specified. Ignored here,\n"
1324  << "but user HAS TO add that later on!\n";
1325  } else {
1326  try {
1327  xml_read_from_file(pnd_field_files[i],
1328  vmr_field_raw[vmr_field_raw.nelem() - 1],
1329  verbosity);
1330  } catch (...) {
1332  try {
1333  xml_read_from_file(pnd_field_files[i], tmp, verbosity);
1334  if (tmp.nelem() == 1) {
1335  vmr_field_raw[vmr_field_raw.nelem() - 1] = tmp[0];
1336  } else {
1337  std::ostringstream os;
1338  os << "The file " << pnd_field_files[i] << "\n"
1339  << "is neither GriddedField3 nor a 1-long ArrayOfGriddedField3.\n";
1340  throw std::runtime_error(os.str());
1341  }
1342  } catch (...) {
1343  std::ostringstream os;
1344  os << "The file " << pnd_field_files[i] << " does not exist or\n"
1345  << "its type is neither GriddedField3 nor a 1-long ArrayOfGriddedField3.\n";
1346  throw std::runtime_error(os.str());
1347  }
1348  }
1349 
1351  pnd_field_files[i],
1353  verbosity);
1354  }
1355 
1356  out2 << " Append 'particle' field to abs_species\n";
1360  species,
1361  verbosity);
1362  }
1363  scat_dataCheck(scat_data_raw, "sane", 1e-2, verbosity);
1364 }
1365 
1366 /* Workspace method: Doxygen documentation will be auto-generated */
1367 void ScatSpeciesScatAndMetaRead( //WS Output:
1370  // Keywords:
1371  const ArrayOfString& scat_data_files,
1372  const Verbosity& verbosity) {
1373  CREATE_OUT2;
1374  CREATE_OUT3;
1375 
1376  //--- Reading the data ---------------------------------------------------
1378  ArrayOfScatteringMetaData arr_smd;
1379 
1380  arr_ssd.resize(scat_data_files.nelem());
1381  arr_smd.resize(scat_data_files.nelem());
1382 
1383  Index meta_naming_conv = 0;
1384 
1385  for (Index i = 0; i < 1 && i < scat_data_files.nelem(); i++) {
1386  out3 << " Read single scattering data file " << scat_data_files[i] << "\n";
1387  xml_read_from_file(scat_data_files[i], arr_ssd[i], verbosity);
1388 
1389  // make meta data name from scat data name
1390  ArrayOfString strarr;
1391  String scat_meta_file;
1392 
1393  if (i == 0) {
1394  scat_data_files[i].split(strarr, ".xml");
1395  scat_meta_file = strarr[0] + ".meta.xml";
1396 
1397  try {
1398  find_xml_file(scat_meta_file, verbosity);
1399  } catch (const runtime_error&) {
1400  }
1401 
1402  if (file_exists(scat_meta_file)) {
1403  out3 << " Read scattering meta data\n";
1404 
1405  xml_read_from_file(scat_meta_file, arr_smd[i], verbosity);
1406 
1407  meta_naming_conv = 1;
1408  } else {
1409  try {
1410  scat_data_files[i].split(strarr, "scat_data");
1411  if (strarr.nelem() < 2) {
1412  throw std::runtime_error(
1413  "Splitting scattering data filename up at 'scat_data' also failed.");
1414  }
1415  scat_meta_file = strarr[0] + "scat_meta" + strarr[1];
1416 
1417  out3 << " Read scattering meta data\n";
1418  xml_read_from_file(scat_meta_file, arr_smd[i], verbosity);
1419 
1420  meta_naming_conv = 2;
1421  } catch (const std::runtime_error& e) {
1422  ostringstream os;
1423  os << "No meta data file following one of the allowed naming "
1424  << "conventions was found.\n"
1425  << "Allowed are "
1426  << "*.meta.xml from *.xml and "
1427  << "*scat_meta* from *scat_data*\n"
1428  << "Scattering meta data file not found: " << scat_meta_file
1429  << "\n"
1430  << e.what();
1431  throw runtime_error(os.str());
1432  }
1433  }
1434  }
1435  }
1436 
1437  ArrayOfString fail_msg;
1438 
1439 #pragma omp parallel for if (!arts_omp_in_parallel() && \
1440  scat_data_files.nelem() > 1) \
1441  num_threads(arts_omp_get_max_threads() > 16 ? 16 \
1442  : arts_omp_get_max_threads()) \
1443  shared(out3, scat_data_files, arr_ssd, arr_smd)
1444  for (Index i = 1; i < scat_data_files.nelem(); i++) {
1445  // make meta data name from scat data name
1446  ArrayOfString strarr;
1447  String scat_meta_file;
1449  ScatteringMetaData smd;
1450 
1451  try {
1452  out3 << " Read single scattering data file " << scat_data_files[i]
1453  << "\n";
1454  xml_read_from_file(scat_data_files[i], ssd, verbosity);
1455 
1456  scat_data_files[i].split(strarr, ".xml");
1457  scat_meta_file = strarr[0] + ".meta.xml";
1458 
1459  if (meta_naming_conv == 1) {
1460  scat_data_files[i].split(strarr, ".xml");
1461  scat_meta_file = strarr[0] + ".meta.xml";
1462 
1463  out3 << " Read scattering meta data\n";
1464  xml_read_from_file(scat_meta_file, smd, verbosity);
1465  } else {
1466  scat_data_files[i].split(strarr, "scat_data");
1467  scat_meta_file = strarr[0] + "scat_meta" + strarr[1];
1468 
1469  out3 << " Read scattering meta data\n";
1470  xml_read_from_file(scat_meta_file, smd, verbosity);
1471  }
1472  } catch (const std::exception& e) {
1473  ostringstream os;
1474  os << "Run-time error reading scattering data : \n" << e.what();
1475 #pragma omp critical(ybatchCalc_push_fail_msg)
1476  fail_msg.push_back(os.str());
1477  }
1478 
1479  //FIXME: currently nothing is done in chk_scattering_meta_data!
1480  chk_scattering_meta_data(smd, scat_meta_file, verbosity);
1481 
1482 #pragma omp critical(ScatSpeciesScatAndMetaRead_assign_ssd)
1483  arr_ssd[i] = std::move(ssd);
1484 #pragma omp critical(ScatSpeciesScatAndMetaRead_assign_smd)
1485  arr_smd[i] = std::move(smd);
1486  }
1487 
1488  if (fail_msg.nelem()) {
1489  ostringstream os;
1490  for (auto& msg : fail_msg) os << msg << '\n';
1491 
1492  throw runtime_error(os.str());
1493  }
1494 
1495  // check if arrays have same size
1496  chk_scattering_data(arr_ssd, arr_smd, verbosity);
1497 
1498  // append as new scattering species
1499  scat_data_raw.push_back(std::move(arr_ssd));
1500  scat_meta.push_back(std::move(arr_smd));
1501 }
1502 
1503 /* Workspace method: Doxygen documentation will be auto-generated */
1504 void ScatElementsSelect( //WS Output:
1507  // WS Input:
1508  const ArrayOfString& scat_species,
1509  const String& species,
1510  const String& sizeparam,
1511  const Numeric& sizemin,
1512  const Numeric& sizemax,
1513  const Numeric& tolerance,
1514  const String& delim,
1515  const Verbosity&) {
1516  // first check that sizes of scat_species and scat_data_raw/scat_meta agree
1517  Index nspecies = scat_species.nelem();
1518  if (nspecies != scat_data_raw.nelem() || nspecies != scat_meta.nelem()) {
1519  ostringstream os;
1520  os << "Number of scattering species specified by scat_species does\n"
1521  << "not agree with number of scattering species in\n"
1522  << "scat_data_raw or scat_meta:\n"
1523  << "scat_species has " << nspecies
1524  << " entries, while scat_data_raw has " << scat_data_raw.nelem()
1525  << " and scat_meta has " << scat_meta.nelem() << ".";
1526  throw runtime_error(os.str());
1527  }
1528 
1529  // create temporary containers for selected elements
1530  ArrayOfSingleScatteringData scat_data_raw_tmp;
1531  ArrayOfScatteringMetaData scat_meta_tmp;
1532 
1533  String partfield_name;
1534  //find the species to handle: compare 'species' to 'partfield' part of
1535  //scat_species tags
1536  Index i_ss = -1;
1537  for (Index i = 0; i < scat_species.nelem(); i++) {
1538  parse_partfield_name(partfield_name, scat_species[i], delim);
1539  if (partfield_name == species) i_ss = i;
1540  }
1541  if (i_ss < 0) {
1542  ostringstream os;
1543  os << "Scattering species " << species << " not found among scat_species.";
1544  throw runtime_error(os.str());
1545  }
1546 
1547  // choosing the specified SingleScatteringData and ScatteringMetaData
1548  if (sizeparam == "diameter_max")
1549  for (Index i_se = 0; i_se < scat_meta[i_ss].nelem(); i_se++) {
1550  // scattering element diameter is extracted from the
1551  // scattering element's meta data and checked whether it's within size
1552  // selected range (sizemax < 0 check follows from wildcard usage and
1553  // means consider all sizes on the upper end)
1554  if (scat_meta[i_ss][i_se].diameter_max > sizemin - sizemin * tolerance &&
1555  (sizemax + sizemax * tolerance > scat_meta[i_ss][i_se].diameter_max ||
1556  sizemax < 0.)) {
1557  // copy selected scattering element to temp arrays
1558  scat_data_raw_tmp.push_back(scat_data_raw[i_ss][i_se]);
1559  scat_meta_tmp.push_back(scat_meta[i_ss][i_se]);
1560  }
1561  }
1562  else if (sizeparam == "diameter_volume_equ")
1563  for (Index i_se = 0; i_se < scat_meta[i_ss].nelem(); i_se++) {
1564  if (scat_meta[i_ss][i_se].diameter_volume_equ >
1565  sizemin - sizemin * tolerance &&
1566  (sizemax + sizemax * tolerance >
1567  scat_meta[i_ss][i_se].diameter_volume_equ ||
1568  sizemax < 0.)) {
1569  // copy selected scattering element to temp arrays
1570  scat_data_raw_tmp.push_back(scat_data_raw[i_ss][i_se]);
1571  scat_meta_tmp.push_back(scat_meta[i_ss][i_se]);
1572  }
1573  }
1574  else if (sizeparam == "diameter_area_equ_aerodynamical")
1575  for (Index i_se = 0; i_se < scat_meta[i_ss].nelem(); i_se++) {
1576  if (scat_meta[i_ss][i_se].diameter_area_equ_aerodynamical >
1577  sizemin - sizemin * tolerance &&
1578  (sizemax + sizemax * tolerance >
1579  scat_meta[i_ss][i_se].diameter_area_equ_aerodynamical ||
1580  sizemax < 0.)) {
1581  // copy selected scattering element to temp arrays
1582  scat_data_raw_tmp.push_back(scat_data_raw[i_ss][i_se]);
1583  scat_meta_tmp.push_back(scat_meta[i_ss][i_se]);
1584  }
1585  }
1586  else {
1587  ostringstream os;
1588  os << "Size parameter " << sizeparam << "is unknown.";
1589  throw runtime_error(os.str());
1590  }
1591 
1592  // To use a particle species field without associated scattering element
1593  // data poses a high risk of accidentially neglecting these species. That's
1594  // unlikely what the user intends. Hence throw error.
1595  if (scat_meta_tmp.nelem() < 1) {
1596  ostringstream os;
1597  os << "For scattering species " << species << " no scattering "
1598  << "element matching the requested size range found.\n"
1599  << "Check *scat_data_raw* and *scat_meta* input as well as your size limit "
1600  << "selection!";
1601  throw runtime_error(os.str());
1602  }
1603 
1604  scat_meta[i_ss] = std::move(scat_meta_tmp);
1605  scat_data_raw[i_ss] = std::move(scat_data_raw_tmp);
1606 
1607  // check if array is empty. should never apply (since we checked the re-worked
1608  // data before and that error should also catch cases that are empty from the
1609  // beginning).
1611 }
1612 
1613 /* Workspace method: Doxygen documentation will be auto-generated */
1616  // Keywords:
1617  const ArrayOfString& scat_species,
1618  const String& species,
1619  const String& scat_species_delim,
1620  const Numeric& T_low,
1621  const Numeric& T_high,
1622  const Verbosity&) {
1623  const bool do_tl = (T_low >= 0.);
1624  const bool do_th = (T_high >= 0.);
1625 
1626  if (do_tl || do_th) {
1627  Index i_ss = -1;
1628  if (species == "") {
1629  i_ss = scat_data_raw.nelem() - 1;
1630  if (i_ss == -1) {
1631  ostringstream os;
1632  os << "No *scat_data* available. Can not extend temperature range on "
1633  << "inexistent data.";
1634  throw runtime_error(os.str());
1635  }
1636  } else {
1637  // first check that sizes of scat_species and scat_data_raw agree
1638  Index nspecies = scat_species.nelem();
1639  if (nspecies != scat_data_raw.nelem()) {
1640  ostringstream os;
1641  os << "Number of scattering species specified by scat_species does\n"
1642  << "not agree with number of scattering species in *scat_data*:\n"
1643  << "scat_species has " << nspecies
1644  << " entries, while *scat_data* has " << scat_data_raw.nelem()
1645  << ".";
1646  throw runtime_error(os.str());
1647  }
1648  String partfield_name;
1649  //find the species to handle: compare 'species' to 'partfield' part of
1650  //scat_species tags
1651  for (Index i = 0; i < scat_species.nelem(); i++) {
1653  partfield_name, scat_species[i], scat_species_delim);
1654  if (partfield_name == species) i_ss = i;
1655  }
1656  if (i_ss < 0) {
1657  ostringstream os;
1658  os << "Scattering species " << species
1659  << " not found among scat_species.";
1660  throw runtime_error(os.str());
1661  }
1662  }
1663 
1664  for (Index i_se = 0; i_se < scat_data_raw[i_ss].nelem(); i_se++) {
1665  const SingleScatteringData& ssdo = scat_data_raw[i_ss][i_se];
1666  const Index nTo = ssdo.T_grid.nelem();
1667  Index nTn = nTo;
1668  bool do_htl, do_hth;
1669  if (nTo > 1) {
1670  do_htl = (do_tl && (T_low < ssdo.T_grid[0]));
1671  do_hth = (do_th && (T_high > last(ssdo.T_grid)));
1672  } else {
1673  do_htl = false;
1674  do_hth = false;
1675  }
1676 
1677  if (do_htl || do_hth) {
1678  // create new instance of SingleScatteringData
1679  SingleScatteringData ssdn;
1680  Index iToff;
1681 
1682  // determine new temperature grid
1683  if (do_htl) nTn += 1;
1684  if (do_hth) nTn += 1;
1685  Vector T_grid_new(nTn);
1686  if (do_htl) {
1687  T_grid_new[0] = T_low;
1688  iToff = 1;
1689  } else {
1690  iToff = 0;
1691  }
1692  for (Index iT = 0; iT < nTo; iT++)
1693  T_grid_new[iT + iToff] = scat_data_raw[i_ss][i_se].T_grid[iT];
1694  if (do_hth) T_grid_new[nTo + iToff] = T_high;
1695  ssdn.T_grid = std::move(T_grid_new);
1696 
1697  // copy grids and other descriptive data that is to remain identical
1698  ssdn.ptype = ssdo.ptype;
1699  ostringstream description;
1700  description << ssdo.description; // here just copy. we append further
1701  // info below if applicable.
1702  ssdn.f_grid = ssdo.f_grid;
1703  ssdn.za_grid = ssdo.za_grid;
1704  ssdn.aa_grid = ssdo.aa_grid;
1705 
1706  // determine size of current optical property data
1707  const Index nf = ssdo.f_grid.nelem();
1708  const Index nzas = ssdo.pha_mat_data.nshelves();
1709  const Index naas = ssdo.pha_mat_data.nbooks();
1710  const Index nzai = ssdo.pha_mat_data.npages();
1711  const Index naai = ssdo.pha_mat_data.nrows();
1712  const Index nmep = ssdo.pha_mat_data.ncols();
1713  const Index nmee = ssdo.ext_mat_data.ncols();
1714  const Index nvea = ssdo.abs_vec_data.ncols();
1715 
1716  // create containers for extended optical property data
1717  ssdn.pha_mat_data.resize(nf, nTn, nzas, naas, nzai, naai, nmep);
1718  ssdn.ext_mat_data.resize(nf, nTn, nzai, naai, nmee);
1719  ssdn.abs_vec_data.resize(nf, nTn, nzai, naai, nvea);
1720 
1721  // copy optical property data
1722  for (Index iT = 0; iT < nTo; iT++) {
1723  ssdn.pha_mat_data(
1724  joker, iT + iToff, joker, joker, joker, joker, joker) =
1725  ssdo.pha_mat_data(joker, iT, joker, joker, joker, joker, joker);
1726  ssdn.ext_mat_data(joker, iT + iToff, joker, joker, joker) =
1727  ssdo.ext_mat_data(joker, iT, joker, joker, joker);
1728  ssdn.abs_vec_data(joker, iT + iToff, joker, joker, joker) =
1729  ssdo.abs_vec_data(joker, iT, joker, joker, joker);
1730  }
1731 
1732  // duplicate optical property data on T-edges if applicable
1733  if (do_htl) {
1734  ssdn.pha_mat_data(joker, 0, joker, joker, joker, joker, joker) =
1735  ssdn.pha_mat_data(joker, 1, joker, joker, joker, joker, joker);
1736  ssdn.ext_mat_data(joker, 0, joker, joker, joker) =
1737  ssdn.ext_mat_data(joker, 1, joker, joker, joker);
1738  ssdn.abs_vec_data(joker, 0, joker, joker, joker) =
1739  ssdn.abs_vec_data(joker, 1, joker, joker, joker);
1740  description << "\n"
1741  << "Low temperature limit extended by"
1742  << " duplicating previous low temperature limit"
1743  << " single scattering properties.";
1744  }
1745  if (do_hth) {
1746  ssdn.pha_mat_data(joker, nTn - 1, joker, joker, joker, joker, joker) =
1747  ssdn.pha_mat_data(
1748  joker, nTn - 2, joker, joker, joker, joker, joker);
1749  ssdn.ext_mat_data(joker, nTn - 1, joker, joker, joker) =
1750  ssdn.ext_mat_data(joker, nTn - 2, joker, joker, joker);
1751  ssdn.abs_vec_data(joker, nTn - 1, joker, joker, joker) =
1752  ssdn.abs_vec_data(joker, nTn - 2, joker, joker, joker);
1753  description << "\n"
1754  << "High temperature limit extended by"
1755  << " duplicating previous high temperature limit"
1756  << " single scattering properties.";
1757  }
1758  ssdn.description = description.str();
1759  scat_data_raw[i_ss][i_se] = std::move(ssdn);
1760  }
1761  }
1762  }
1763 }
1764 
1765 /* Workspace method: Doxygen documentation will be auto-generated */
1767  Tensor4& pnd_field,
1769  //WS Input
1770  const Vector& p_grid,
1771  const Vector& lat_grid,
1772  const Vector& lon_grid,
1774  const Index& atmosphere_dim,
1777  const Index& zeropadding,
1778  const Verbosity& verbosity) {
1779  // Basic checks of input variables
1780  //
1781  // Particle number density data
1782  //
1783  if (pnd_field_raw.empty()) {
1784  ostringstream os;
1785  os << "No particle number density data given. Please use WSMs\n"
1786  << "*ParticleTypeInit* and *ParticleTypeAdd(All)* for reading\n"
1787  << "scattering element data.\n";
1788  throw runtime_error(os.str());
1789  }
1790 
1792  ArrayOfIndex cloudbox_limits_tmp;
1793  /*if ( cloudbox_limits.empty() )
1794  {
1795  //If no limits set, the cloud(box) is supposed to cover the
1796  //complete atmosphere. This particularly to facilitate use of
1797  //scat_data_single&pnd_fields for non-scatt, greybody cloud calculations.
1798  //Hence, set the limits to first & last elements of the different grids.
1799  //Note: no margin left at lat/lon_grid edges!.
1800  cloudbox_limits_tmp.resize(2*atmosphere_dim);
1801 
1802  // Pressure limits
1803  cloudbox_limits_tmp[0] = 0;
1804  cloudbox_limits_tmp[1] = p_grid.nelem() - 1;
1805  // Latitude limits
1806  if( atmosphere_dim >= 2 )
1807  {
1808  cloudbox_limits_tmp[2] = 0;
1809  cloudbox_limits_tmp[3] = lat_grid.nelem() - 1;
1810  }
1811  // Longitude limits
1812  if( atmosphere_dim == 3 )
1813  {
1814  cloudbox_limits_tmp[4] = 0;
1815  cloudbox_limits_tmp[5] = lon_grid.nelem() - 1;
1816  }
1817  }
1818  else */
1819  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
1820  throw runtime_error(
1821  "*cloudbox_limits* is a vector which contains the"
1822  "upper and lower limit of the cloud for all "
1823  "atmospheric dimensions. So its dimension must"
1824  "be 2 x *atmosphere_dim*");
1825  else
1826  cloudbox_limits_tmp = cloudbox_limits;
1827 
1828  // Check that pnd_field_raw has at least 2 grid-points in each dimension.
1829  // Otherwise, interpolation further down will fail with assertion.
1830 
1831  for (Index d = 0; d < atmosphere_dim; d++) {
1832  for (Index i = 0; i < pnd_field_raw.nelem(); i++) {
1833  if (pnd_field_raw[i].get_grid_size(d) < 2) {
1834  ostringstream os;
1835  os << "Error in pnd_field_raw data. ";
1836  os << "Dimension " << d << " (name: \"";
1837  os << pnd_field_raw[i].get_grid_name(d);
1838  os << "\") has only ";
1839  os << pnd_field_raw[i].get_grid_size(d);
1840  os << " element";
1841  os << ((pnd_field_raw[i].get_grid_size(d) == 1) ? "" : "s");
1842  os << ". Must be at least 2.";
1843  throw runtime_error(os.str());
1844  }
1845  }
1846  }
1847  const Index Np_cloud = cloudbox_limits_tmp[1] - cloudbox_limits_tmp[0] + 1;
1848 
1849  ConstVectorView p_grid_cloud =
1850  p_grid[Range(cloudbox_limits_tmp[0], Np_cloud)];
1851 
1852  // Check that no scatterers exist outside the cloudbox
1854  pnd_field_raw,
1855  p_grid,
1856  lat_grid,
1857  lon_grid,
1858  cloudbox_limits_tmp);
1859 
1860  //==========================================================================
1861  if (atmosphere_dim == 1) {
1862  ArrayOfGriddedField3 pnd_field_tmp;
1863 
1865  pnd_field_tmp, p_grid_cloud, pnd_field_raw, 1, zeropadding, verbosity);
1866 
1868  p_grid_cloud,
1869  pnd_field_tmp[0].get_numeric_grid(1),
1870  pnd_field_tmp[0].get_numeric_grid(2),
1871  pnd_field_tmp,
1872  verbosity);
1873  } else if (atmosphere_dim == 2) {
1874  const Index Nlat_cloud =
1875  cloudbox_limits_tmp[3] - cloudbox_limits_tmp[2] + 1;
1876 
1877  ConstVectorView lat_grid_cloud =
1878  lat_grid[Range(cloudbox_limits_tmp[2], Nlat_cloud)];
1879 
1880  if (zeropadding) {
1881  // FIXME: OLE: Implement this
1882  CREATE_OUT0;
1883  out0 << "WARNING: zeropadding currently only supported for 1D.";
1884  }
1885 
1886  //Resize variables
1887  pnd_field.resize(pnd_field_raw.nelem(), Np_cloud, Nlat_cloud, 1);
1888 
1889  // Gridpositions:
1890  ArrayOfGridPos gp_p(Np_cloud);
1891  ArrayOfGridPos gp_lat(Nlat_cloud);
1892 
1893  // Interpolate pnd_field.
1894  // Loop over the scattering element:
1895  for (Index i = 0; i < pnd_field_raw.nelem(); ++i) {
1896  // Calculate grid positions:
1897  p2gridpos(gp_p,
1898  pnd_field_raw[i].get_numeric_grid(GFIELD3_P_GRID),
1899  p_grid_cloud);
1900  gridpos(gp_lat,
1901  pnd_field_raw[i].get_numeric_grid(GFIELD3_LAT_GRID),
1902  lat_grid_cloud);
1903 
1904  // Interpolation weights:
1905  Tensor3 itw(Np_cloud, Nlat_cloud, 4);
1906  interpweights(itw, gp_p, gp_lat);
1907 
1908  // Interpolate:
1909  interp(pnd_field(i, joker, joker, 0),
1910  itw,
1911  pnd_field_raw[i].data(joker, joker, 0),
1912  gp_p,
1913  gp_lat);
1914  }
1915  } else {
1916  const Index Nlat_cloud =
1917  cloudbox_limits_tmp[3] - cloudbox_limits_tmp[2] + 1;
1918  const Index Nlon_cloud =
1919  cloudbox_limits_tmp[5] - cloudbox_limits_tmp[4] + 1;
1920 
1921  if (zeropadding) {
1922  // FIXME: OLE: Implement this
1923  CREATE_OUT0;
1924  out0 << "WARNING: zeropadding currently only supported for 1D.";
1925  }
1926 
1927  ConstVectorView lat_grid_cloud =
1928  lat_grid[Range(cloudbox_limits_tmp[2], Nlat_cloud)];
1929  ConstVectorView lon_grid_cloud =
1930  lon_grid[Range(cloudbox_limits_tmp[4], Nlon_cloud)];
1931 
1932  //Resize variables
1933  pnd_field.resize(pnd_field_raw.nelem(), Np_cloud, Nlat_cloud, Nlon_cloud);
1934 
1935  // Gridpositions:
1936  ArrayOfGridPos gp_p(Np_cloud);
1937  ArrayOfGridPos gp_lat(Nlat_cloud);
1938  ArrayOfGridPos gp_lon(Nlon_cloud);
1939 
1940  // Interpolate pnd_field.
1941  // Loop over the scattering element types:
1942  for (Index i = 0; i < pnd_field_raw.nelem(); ++i) {
1943  // Calculate grid positions:
1944  p2gridpos(gp_p,
1945  pnd_field_raw[i].get_numeric_grid(GFIELD3_P_GRID),
1946  p_grid_cloud);
1947  gridpos(gp_lat,
1948  pnd_field_raw[i].get_numeric_grid(GFIELD3_LAT_GRID),
1949  lat_grid_cloud);
1950  gridpos(gp_lon,
1951  pnd_field_raw[i].get_numeric_grid(GFIELD3_LON_GRID),
1952  lon_grid_cloud);
1953 
1954  // Interpolation weights:
1955  Tensor4 itw(Np_cloud, Nlat_cloud, Nlon_cloud, 8);
1956  interpweights(itw, gp_p, gp_lat, gp_lon);
1957 
1958  // Interpolate:
1960  itw,
1961  pnd_field_raw[i].data,
1962  gp_p,
1963  gp_lat,
1964  gp_lon);
1965  }
1966  }
1967 
1968  // no (cloudy) Jacobians with this WSM, hence no calc.
1969  // but we need to size dpnd_field to be consistent with jacobian_quantities.
1970  dpnd_field_dx.resize(jacobian_quantities.nelem());
1971 }
1972 
1973 /* Workspace method: Doxygen documentation will be auto-generated */
1975  const Index& atmosphere_dim,
1976  const Index& cloudbox_on,
1978  const Index& nzero,
1979  const Verbosity&) {
1980  /* if( !cloudbox_checked )
1981  throw runtime_error( "The cloudbox must be flagged to have passed a "
1982  "consistency check (cloudbox_checked=1)." );
1983 */
1984  if (atmosphere_dim == 1) {
1985  throw runtime_error("No use in calling this method for 1D.");
1986  }
1987  if (!cloudbox_on) {
1988  throw runtime_error("No use in calling this method with cloudbox off.");
1989  }
1990 
1991  if (nzero < 1) {
1992  throw runtime_error("The argument *nzero* must be > 0.");
1993  }
1994 
1995  // Sizes
1996  const Index npart = pnd_field.nbooks();
1997  const Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
1998  const Index nlat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
1999  Index nlon = 1;
2000  if (atmosphere_dim == 3) {
2001  nlon = cloudbox_limits[5] - cloudbox_limits[4] + 1;
2002  }
2003 
2004  if (pnd_field.npages() != np || pnd_field.nrows() != 1 ||
2005  pnd_field.ncols() != 1) {
2006  throw runtime_error(
2007  "The input *pnd_field* is either not 1D or does not "
2008  "match pressure size of cloudbox.");
2009  }
2010 
2011  // Temporary container
2012  Tensor4 pnd_temp = pnd_field;
2013 
2014  // Resize and fill
2015  pnd_field.resize(npart, np, nlat, nlon);
2016  pnd_field = 0;
2017  //
2018  for (Index ilon = nzero; ilon < nlon - nzero; ilon++) {
2019  for (Index ilat = nzero; ilat < nlat - nzero; ilat++) {
2020  for (Index ip = 0; ip < np; ip++) {
2021  for (Index is = 0; is < npart; is++) {
2022  pnd_field(is, ip, ilat, ilon) = pnd_temp(is, ip, 0, 0);
2023  }
2024  }
2025  }
2026  }
2027 }
2028 
2029 /* Workspace method: Doxygen documentation will be auto-generated */
2030 void pnd_fieldZero( //WS Output:
2031  Tensor4& pnd_field,
2034  //WS Input:
2035  const Index& atmosphere_dim,
2036  const Vector& f_grid,
2039  const Verbosity&) {
2040  chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
2041 
2042  if (cloudbox_limits.nelem() != 2 * atmosphere_dim)
2043  throw runtime_error(
2044  "*cloudbox_limits* is a vector which contains the"
2045  "upper and lower limit of the cloud for all "
2046  "atmospheric dimensions. So its dimension must"
2047  "be 2 x *atmosphere_dim*");
2048 
2049  // Resize pnd_field and set it to 0:
2050  Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2051  Index nlat = 1, nlon = 1;
2052  if (atmosphere_dim > 1) {
2053  nlat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
2054  if (atmosphere_dim > 2) {
2055  nlon = cloudbox_limits[5] - cloudbox_limits[4] + 1;
2056  }
2057  }
2058 
2059  // no (cloudy) Jacobians with this WSM, hence no setting.
2060  // but we need to size dpnd_field to be consistent with jacobian_quantities.
2061  dpnd_field_dx.resize(jacobian_quantities.nelem());
2062 
2063  // Do only reset scat_data if it has not been set yet.
2064  // There's no need otherwise, and it's rather unpractical for testing when
2065  // doing so: we might want to do some actual calcs with the scat_data later
2066  // on. So why throw it away?
2067  const Index N_se = TotalNumberOfElements(scat_data);
2068  if (N_se > 0) {
2069  pnd_field.resize(N_se, np, nlat, nlon);
2070  } else {
2071  pnd_field.resize(1, np, nlat, nlon);
2072 
2073  //Resize scat_data and set it to 0:
2074  // Number of scattering elements
2075  scat_data.resize(1);
2076  scat_data[0].resize(1);
2077  scat_data[0][0].ptype = PTYPE_TOTAL_RND;
2078  scat_data[0][0].description = " ";
2079  // Grids which contain full ranges which one wants to calculate
2080  Index nf = f_grid.nelem();
2081  scat_data[0][0].f_grid.resize(nf);
2082  scat_data[0][0].f_grid = f_grid;
2083  Index nT = 1;
2084  scat_data[0][0].T_grid.resize(nT);
2085  scat_data[0][0].T_grid = 270.;
2086  Index nza = 5;
2087  nlinspace(scat_data[0][0].za_grid, 0, 180, nza);
2088  // Resize the data arrays
2089  scat_data[0][0].pha_mat_data.resize(nf, nT, nza, 1, 1, 1, 6);
2090  scat_data[0][0].pha_mat_data = 0.;
2091  scat_data[0][0].ext_mat_data.resize(nf, nT, 1, 1, 1);
2092  scat_data[0][0].ext_mat_data = 0.;
2093  scat_data[0][0].abs_vec_data.resize(nf, nT, 1, 1, 1);
2094  scat_data[0][0].abs_vec_data = 0.;
2095  }
2096 
2097  pnd_field = 0.;
2098 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
ARTS::Var::scat_data_single
SingleScatteringData scat_data_single(Workspace &ws) noexcept
Definition: autoarts.h:5926
gridded_fields.h
Implementation of gridded fields.
is_same_within_epsilon
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
ARTS::Var::z_field
Tensor3 z_field(Workspace &ws) noexcept
Definition: autoarts.h:7690
ConstTensor7View::nshelves
Index nshelves() const
Returns the number of shelves.
Definition: matpackVII.cc:48
particle_fieldCleanup
void particle_fieldCleanup(Tensor4 &particle_field_out, const Tensor4 &particle_field_in, const Numeric &threshold, const Verbosity &)
WORKSPACE METHOD: particle_fieldCleanup.
Definition: m_cloudbox.cc:1104
SingleScatteringData::za_grid
Vector za_grid
Definition: optproperties.h:85
SingleScatteringData::f_grid
Vector f_grid
Definition: optproperties.h:83
ScatSpeciesScatAndMetaRead
void ScatSpeciesScatAndMetaRead(ArrayOfArrayOfSingleScatteringData &scat_data_raw, ArrayOfArrayOfScatteringMetaData &scat_meta, const ArrayOfString &scat_data_files, const Verbosity &verbosity)
WORKSPACE METHOD: ScatSpeciesScatAndMetaRead.
Definition: m_cloudbox.cc:1367
ARTS::Var::scat_meta
ArrayOfArrayOfScatteringMetaData scat_meta(Workspace &ws) noexcept
Definition: autoarts.h:5983
auto_md.h
ARTS::Var::iy_cloudbox_agenda
Agenda iy_cloudbox_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3742
ARTS::Var::jacobian_quantities
ArrayOfRetrievalQuantity jacobian_quantities(Workspace &ws) noexcept
Definition: autoarts.h:3924
abs_speciesAdd
void abs_speciesAdd(ArrayOfArrayOfSpeciesTag &abs_species, Index &propmat_clearsky_agenda_checked, Index &abs_xsec_agenda_checked, const ArrayOfString &names, const Verbosity &verbosity)
WORKSPACE METHOD: abs_speciesAdd.
Definition: m_abs_lookup.cc:1870
chk_if_increasing
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:117
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
ConstTensor5View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackV.cc:56
cloudboxSetManually
void cloudboxSetManually(Index &cloudbox_on, ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Numeric &p1, const Numeric &p2, const Numeric &lat1, const Numeric &lat2, const Numeric &lon1, const Numeric &lon2, const Verbosity &)
WORKSPACE METHOD: cloudboxSetManually.
Definition: m_cloudbox.cc:368
p2gridpos
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
Calculates grid positions for pressure values.
Definition: special_interp.cc:759
ConstTensor7View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackVII.cc:60
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
joker
const Joker joker
interpolation.h
Header file for interpolation.cc.
GFIELD3_LAT_GRID
const Index GFIELD3_LAT_GRID
ARTS::Var::pnd_field_raw
ArrayOfGriddedField3 pnd_field_raw(Workspace &ws) noexcept
Definition: autoarts.h:5105
ARTS::Var::pnd_field
Tensor4 pnd_field(Workspace &ws) noexcept
Definition: autoarts.h:5081
SingleScatteringData::description
String description
Definition: optproperties.h:82
ARTS::Var::ppath_inside_cloudbox_do
Index ppath_inside_cloudbox_do(Workspace &ws) noexcept
Definition: autoarts.h:5172
last
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:165
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
chk_interpolation_grids
void chk_interpolation_grids(const String &which_interpolation, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac, const bool islog)
Check interpolation grids.
Definition: check_input.cc:989
ARTS::Var::cloudbox_limits
ArrayOfIndex cloudbox_limits(Workspace &ws) noexcept
Definition: autoarts.h:2762
ARTS::Var::jacobian_do
Index jacobian_do(Workspace &ws) noexcept
Definition: autoarts.h:3912
SingleScatteringData::aa_grid
Vector aa_grid
Definition: optproperties.h:86
DEBUG_ONLY
#define DEBUG_ONLY(...)
Definition: debug.h:36
cloudbox_fieldCrop
void cloudbox_fieldCrop(Tensor7 &cloudbox_field, ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Index &cloudbox_on, const Index &new_limit0, const Index &new_limit1, const Index &new_limit2, const Index &new_limit3, const Index &new_limit4, const Index &new_limit5, const Verbosity &)
WORKSPACE METHOD: cloudbox_fieldCrop.
Definition: m_cloudbox.cc:1030
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
ScatElementsPndAndScatAdd
void ScatElementsPndAndScatAdd(ArrayOfArrayOfSingleScatteringData &scat_data_raw, ArrayOfGriddedField3 &pnd_field_raw, const Index &atmosphere_dim, const ArrayOfString &scat_data_files, const ArrayOfString &pnd_field_files, const Verbosity &verbosity)
WORKSPACE METHOD: ScatElementsPndAndScatAdd.
Definition: m_cloudbox.cc:1145
data
G0 G2 FVC Y DV Numeric Numeric Numeric Zeeman LowerQuantumNumbers void * data
Definition: arts_api_classes.cc:232
fractional_gp
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Definition: interpolation.cc:492
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
cloudboxOff
void cloudboxOff(Index &cloudbox_on, Index &ppath_inside_cloudbox_do, ArrayOfIndex &cloudbox_limits, Agenda &iy_cloudbox_agenda, Tensor4 &pnd_field, ArrayOfTensor4 &dpnd_field_dx, ArrayOfString &scat_species, ArrayOfArrayOfSingleScatteringData &scat_data, ArrayOfArrayOfSingleScatteringData &scat_data_raw, Index &scat_data_checked, Matrix &particle_masses, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
WORKSPACE METHOD: cloudboxOff.
Definition: m_cloudbox.cc:81
chk_pnd_field_raw_only_in_cloudbox
void chk_pnd_field_raw_only_in_cloudbox(const Index &dim, const ArrayOfGriddedField3 &pnd_field_raw, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, const ArrayOfIndex &cloudbox_limits)
chk_pnd_field_raw_only_in_cloudbox
Definition: cloudbox.cc:146
DEG2RAD
const Numeric DEG2RAD
CREATE_OUT2
#define CREATE_OUT2
Definition: messages.h:206
ARTS::Var::iy
Matrix iy(Workspace &ws) noexcept
Definition: autoarts.h:3690
TotalNumberOfElements
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
array.h
This file contains the definition of Array.
is_size
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
is_gridpos_at_index_i
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
Definition: interpolation.cc:657
CREATE_OUT1
#define CREATE_OUT1
Definition: messages.h:205
Agenda
The Agenda class.
Definition: agenda_class.h:44
GriddedField3
Definition: gridded_fields.h:284
SingleScatteringData
Definition: optproperties.h:80
GFIELD3_LON_GRID
const Index GFIELD3_LON_GRID
chk_pnd_raw_data
void chk_pnd_raw_data(const ArrayOfGriddedField3 &pnd_field_raw, const String &pnd_field_file, const Index &atmosphere_dim, const Verbosity &verbosity)
Check particle number density files (pnd_field_raw)
Definition: cloudbox.cc:116
microphysics.h
Internal functions for microphysics calculations (size distributions etc.)
Array
This can be used to make arrays out of anything.
Definition: array.h:108
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
parse_partfield_name
void parse_partfield_name(String &partfield_name, const String &part_string, const String &delim)
Definition: cloudbox.cc:944
min
#define min(a, b)
Definition: legacy_continua.cc:20628
xml_read_from_file
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.cc:901
Tensor7::resize
void resize(Index l, Index v, Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackVII.cc:5484
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:207
SingleScatteringData::ptype
PType ptype
Definition: optproperties.h:81
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:204
Tensor5::resize
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1743
PI
const Numeric PI
gridpos_poly_cyclic_longitudinal
void gridpos_poly_cyclic_longitudinal(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
Set up grid positions for higher order interpolation.
Definition: interpolation_poly.cc:326
messages.h
Declarations having to do with the four output streams.
ARTS::Var::rte_pos
Vector rte_pos(Workspace &ws) noexcept
Definition: autoarts.h:5672
pnd_fieldExpand1D
void pnd_fieldExpand1D(Tensor4 &pnd_field, const Index &atmosphere_dim, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &nzero, const Verbosity &)
WORKSPACE METHOD: pnd_fieldExpand1D.
Definition: m_cloudbox.cc:1974
ARTS::Var::particle_masses
Matrix particle_masses(Workspace &ws) noexcept
Definition: autoarts.h:4850
gridpos_upperend_check
void gridpos_upperend_check(GridPos &gp, const Index &ie)
gridpos_upperend_check
Definition: interpolation.cc:589
GridPosPoly
Structure to store a grid position for higher order interpolation.
Definition: interpolation_poly.h:64
my_basic_string< char >
ScatSpeciesExtendTemperature
void ScatSpeciesExtendTemperature(ArrayOfArrayOfSingleScatteringData &scat_data_raw, const ArrayOfString &scat_species, const String &species, const String &scat_species_delim, const Numeric &T_low, const Numeric &T_high, const Verbosity &)
WORKSPACE METHOD: ScatSpeciesExtendTemperature.
Definition: m_cloudbox.cc:1614
ConstTensor4View::empty
bool empty() const
Check if variable is empty.
Definition: matpackIV.cc:52
ScatteringMetaData
Definition: optproperties.h:109
SingleScatteringData::ext_mat_data
Tensor5 ext_mat_data
Definition: optproperties.h:88
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
ScatElementsToabs_speciesAdd
void ScatElementsToabs_speciesAdd(ArrayOfArrayOfSingleScatteringData &scat_data_raw, ArrayOfGriddedField3 &vmr_field_raw, ArrayOfArrayOfSpeciesTag &abs_species, Index &propmat_clearsky_agenda_checked, Index &abs_xsec_agenda_checked, const Index &atmosphere_dim, const Vector &f_grid, const ArrayOfString &scat_data_files, const ArrayOfString &pnd_field_files, const Verbosity &verbosity)
WORKSPACE METHOD: ScatElementsToabs_speciesAdd.
Definition: m_cloudbox.cc:1256
ARTS::Var::rte_los
Vector rte_los(Workspace &ws) noexcept
Definition: autoarts.h:5651
physics_funcs.h
This file contains declerations of functions of physical character.
ARTS::Var::abs_species
ArrayOfArrayOfSpeciesTag abs_species(Workspace &ws) noexcept
Definition: autoarts.h:2157
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
ConstTensor4View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:66
optproperties.h
Scattering database structure and functions.
GridPosPoly::idx
ArrayOfIndex idx
Definition: interpolation_poly.h:67
ConstTensor7View::npages
Index npages() const
Returns the number of pages.
Definition: matpackVII.cc:54
find_xml_file
void find_xml_file(String &filename, const Verbosity &verbosity)
Find an xml file.
Definition: file.cc:414
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
Verbosity
Definition: messages.h:49
ConstTensor4View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:60
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:57
math_funcs.h
chk_scattering_meta_data
void chk_scattering_meta_data(const ScatteringMetaData &scat_meta_single, const String &scat_meta_file, const Verbosity &verbosity)
Check scattering data meta.
Definition: cloudbox.cc:281
find_cloudlimits
void find_cloudlimits(Index &lower, Index &upper, const Tensor3 &scat_species_field, const Index &atmosphere_dim, const Numeric &cloudbox_margin)
Adjust uppermost and lowermost cloudy level for one scat_species_*_*_field.
Definition: cloudbox.cc:714
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
max
#define max(a, b)
Definition: legacy_continua.cc:20629
file_exists
bool file_exists(const String &filename)
Checks if the given file exists.
Definition: file.cc:303
ARTS::Var::scat_data
ArrayOfArrayOfSingleScatteringData scat_data(Workspace &ws) noexcept
Definition: autoarts.h:5840
parameters.h
This file contains header information for the dealing with command line parameters.
ScatSpeciesPndAndScatAdd
void ScatSpeciesPndAndScatAdd(ArrayOfArrayOfSingleScatteringData &scat_data_raw, ArrayOfGriddedField3 &pnd_field_raw, const Index &atmosphere_dim, const ArrayOfString &scat_data_files, const String &pnd_fieldarray_file, const Verbosity &verbosity)
WORKSPACE METHOD: ScatSpeciesPndAndScatAdd.
Definition: m_cloudbox.cc:1211
ARTS::Var::za_grid
Vector za_grid(Workspace &ws) noexcept
Definition: autoarts.h:7771
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:231
ARTS::Var::dpnd_field_dx
ArrayOfTensor4 dpnd_field_dx(Workspace &ws) noexcept
Definition: autoarts.h:3251
cloudbox.h
Internal cloudbox functions.
chk_pnd_data
void chk_pnd_data(const GriddedField3 &pnd_field_raw, const String &pnd_field_file, const Index &atmosphere_dim, const Verbosity &verbosity)
Check particle number density files.
Definition: cloudbox.cc:67
my_basic_string::split
void split(Array< my_basic_string< charT > > &aos, const my_basic_string< charT > &delim) const
Split string into substrings.
Definition: mystring.h:206
ConstTensor4View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:63
pnd_fieldCalcFrompnd_field_raw
void pnd_fieldCalcFrompnd_field_raw(Tensor4 &pnd_field, ArrayOfTensor4 &dpnd_field_dx, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfGriddedField3 &pnd_field_raw, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &zeropadding, const Verbosity &verbosity)
WORKSPACE METHOD: pnd_fieldCalcFrompnd_field_raw.
Definition: m_cloudbox.cc:1766
DENSITY_OF_ICE
const Numeric DENSITY_OF_ICE
ARTS::Var::scat_data_checked
Index scat_data_checked(Workspace &ws) noexcept
Definition: autoarts.h:5850
Range
The range class.
Definition: matpackI.h:160
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
ConstTensor7View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackVII.cc:51
ARTS::Var::propmat_clearsky_agenda_checked
Index propmat_clearsky_agenda_checked(Workspace &ws) noexcept
Definition: autoarts.h:5414
logic.h
Header file for logic.cc.
iyInterpCloudboxField
void iyInterpCloudboxField(Matrix &iy, const Tensor7 &cloudbox_field, const Vector &rte_pos, const Vector &rte_los, const Index &jacobian_do, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Matrix &z_surface, const Index &stokes_dim, const Vector &za_grid, const Vector &aa_grid, const Vector &f_grid, const Index &za_interp_order, const Index &za_restrict, const Index &cos_za_interp, const Numeric &za_extpolfac, const Index &aa_interp_order, const Verbosity &)
WORKSPACE METHOD: iyInterpCloudboxField.
Definition: m_cloudbox.cc:583
chk_scattering_data
void chk_scattering_data(const ArrayOfSingleScatteringData &scat_data, const ArrayOfScatteringMetaData &scat_meta, const Verbosity &)
Check scattering data general.
Definition: cloudbox.cc:259
ScatElementsSelect
void ScatElementsSelect(ArrayOfArrayOfSingleScatteringData &scat_data_raw, ArrayOfArrayOfScatteringMetaData &scat_meta, const ArrayOfString &scat_species, const String &species, const String &sizeparam, const Numeric &sizemin, const Numeric &sizemax, const Numeric &tolerance, const String &delim, const Verbosity &)
WORKSPACE METHOD: ScatElementsSelect.
Definition: m_cloudbox.cc:1504
cloudboxSetManuallyAltitude
void cloudboxSetManuallyAltitude(Index &cloudbox_on, ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Tensor3 &z_field, const Vector &lat_grid, const Vector &lon_grid, const Numeric &z1, const Numeric &z2, const Numeric &lat1, const Numeric &lat2, const Numeric &lon1, const Numeric &lon2, const Verbosity &)
WORKSPACE METHOD: cloudboxSetManuallyAltitude.
Definition: m_cloudbox.cc:477
chk_atm_grids
void chk_atm_grids(const Index &dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid)
chk_atm_grids
Definition: check_input.cc:1201
pnd_fieldZero
void pnd_fieldZero(Tensor4 &pnd_field, ArrayOfTensor4 &dpnd_field_dx, ArrayOfArrayOfSingleScatteringData &scat_data, const Index &atmosphere_dim, const Vector &f_grid, const ArrayOfIndex &cloudbox_limits, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &)
WORKSPACE METHOD: pnd_fieldZero.
Definition: m_cloudbox.cc:2030
GriddedFieldPRegrid
void GriddedFieldPRegrid(GriddedField3 &gfraw_out, const Vector &p_grid, const GriddedField3 &gfraw_in_orig, const Index &interp_order, const Index &zeropadding, const Verbosity &verbosity)
WORKSPACE METHOD: GriddedFieldPRegrid.
Definition: m_atmosphere.cc:703
rte.h
Declaration of functions in rte.cc.
GFIELD3_P_GRID
const Index GFIELD3_P_GRID
ARTS::Var::abs_xsec_agenda_checked
Index abs_xsec_agenda_checked(Workspace &ws) noexcept
Definition: autoarts.h:2283
special_interp.h
Header file for special_interp.cc.
scat_dataCheck
void scat_dataCheck(const ArrayOfArrayOfSingleScatteringData &scat_data, const String &check_type, const Numeric &threshold, const Verbosity &verbosity)
WORKSPACE METHOD: scat_dataCheck.
Definition: m_optproperties.cc:1094
SingleScatteringData::abs_vec_data
Tensor5 abs_vec_data
Definition: optproperties.h:89
RAD2DEG
const Numeric RAD2DEG
LAT_LON_MIN
const Numeric LAT_LON_MIN
GridPosPoly::w
Vector w
Definition: interpolation_poly.h:70
ARTS::Var::cloudbox_field
Tensor7 cloudbox_field(Workspace &ws) noexcept
Definition: autoarts.h:2676
ARTS::Var::scat_data_raw
ArrayOfArrayOfSingleScatteringData scat_data_raw(Workspace &ws) noexcept
Definition: autoarts.h:5898
cloudboxSetAutomatically
void cloudboxSetAutomatically(Index &cloudbox_on, ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor4 &particle_field, const ArrayOfIndex &cloudbox_limits_old, const Numeric &cloudbox_margin, const Verbosity &verbosity)
WORKSPACE METHOD: cloudboxSetAutomatically.
Definition: m_cloudbox.cc:112
chk_scat_species_field
void chk_scat_species_field(bool &empty_flag, const Tensor3 &scat_species_field, const String &fieldname, const Index &dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Check whether field of a specific scattering species zero everywhere.
Definition: cloudbox.cc:646
file.h
This file contains basic functions to handle ASCII files.
rte_pos2gridpos
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)
Converts a geographical position (rte_pos) to grid positions for p, lat and lon.
Definition: special_interp.cc:788
GridPos::idx
Index idx
Definition: interpolation.h:74
cloudboxSetFullAtm
void cloudboxSetFullAtm(Index &cloudbox_on, ArrayOfIndex &cloudbox_limits, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Verbosity &)
WORKSPACE METHOD: cloudboxSetFullAtm.
Definition: m_cloudbox.cc:283
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
chk_if_in_range
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:89
ConstTensor7View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackVII.cc:57
ARTS::Var::vmr_field_raw
ArrayOfGriddedField3 vmr_field_raw(Workspace &ws) noexcept
Definition: autoarts.h:7158
ARTS::Var::scat_species
ArrayOfString scat_species(Workspace &ws) noexcept
Definition: autoarts.h:6086
check_input.h
gridpos_poly
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
The maximum difference from 1 that we allow for a sum check.
Definition: interpolation_poly.cc:120
Vector
The Vector class.
Definition: matpackI.h:860
SingleScatteringData::pha_mat_data
Tensor7 pha_mat_data
Definition: optproperties.h:87
ARTS::Var::aa_grid
Vector aa_grid(Workspace &ws) noexcept
Definition: autoarts.h:1717
PTYPE_TOTAL_RND
@ PTYPE_TOTAL_RND
Definition: optproperties.h:56
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
SingleScatteringData::T_grid
Vector T_grid
Definition: optproperties.h:84
sorting.h
Contains sorting routines.
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
barometric_heightformula
Numeric barometric_heightformula(const Numeric &p, const Numeric &dh)
barometric_heightformula
Definition: physics_funcs.cc:62
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
ARTS::Var::z_surface
Matrix z_surface(Workspace &ws) noexcept
Definition: autoarts.h:7754
ScatSpeciesInit
void ScatSpeciesInit(ArrayOfString &scat_species, ArrayOfArrayOfSingleScatteringData &scat_data_raw, ArrayOfArrayOfScatteringMetaData &scat_meta, Index &scat_data_checked, ArrayOfGriddedField3 &pnd_field_raw, const Verbosity &)
WORKSPACE METHOD: ScatSpeciesInit.
Definition: m_cloudbox.cc:1130
Tensor7
The Tensor7 class.
Definition: matpackVII.h:2382
arts.h
The global header file for ARTS.
FieldFromGriddedField
void FieldFromGriddedField(Matrix &field_out, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const GriddedField2 &gfraw_in, const Verbosity &)
WORKSPACE METHOD: FieldFromGriddedField.
Definition: m_atmosphere.cc:339
ARTS::Var::cloudbox_on
Index cloudbox_on(Workspace &ws) noexcept
Definition: autoarts.h:2782
xml_io.h
This file contains basic functions to handle XML data files.
ARTS::Group::Agenda
Agenda Agenda
Definition: autoarts.h:28