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