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