ARTS 2.5.10 (git: 2f1c442c)
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 */
914 Tensor7& cloudbox_field,
915 const Index& cloudbox_on,
916 const Vector& aa_grid,
917 const Numeric& local_los_azimuth_angle,
918 const Index& aa_interp_order,
919 const Verbosity&) {
920 //--- Check input -----------------------------------------------------------
921 ARTS_USER_ERROR_IF (!cloudbox_on,
922 "No need to use this method with cloudbox=0.");
923 ARTS_USER_ERROR_IF (!(aa_interp_order < aa_grid.nelem()),
924 "Azimuth angle interpolation order *aa_interp_order*"
925 " must be smaller\n"
926 "than number of angles in *aa_grid*.");
927 //---------------------------------------------------------------------------
928
929 if (cloudbox_field.nrows()>1 && aa_grid.nelem()==cloudbox_field.nrows()){
930 Index nf = cloudbox_field.nlibraries();
931 Index np = cloudbox_field.nvitrines();
932 Index nla= cloudbox_field.nshelves();
933 Index nlo= cloudbox_field.nbooks();
934 Index nz = cloudbox_field.npages();
935 Index ns = cloudbox_field.ncols();
936
937 const Tensor7 cloudbox_field_in = std::move(cloudbox_field);
938 Numeric azimuth_los = local_los_azimuth_angle;
939
940 //Convert azimuth from -180,180 to 0,360 convention, as aa_grid is defined in 0,360
941 if (azimuth_los<0) azimuth_los+=360;
942
943 cloudbox_field.resize(nf,np,1,1,nz,1,ns);
944
945 // define interpolations compute cyclic for a azimuth grid [0, 360]
946 const auto lag_aa = LagrangeInterpolation(0,
947 azimuth_los,
948 aa_grid,
949 aa_interp_order,
950 false,
951 Interpolation::GridType::Cyclic,
952 {0, 360});
953
954 // Corresponding interpolation weights
955 const auto itw_aa=interpweights(lag_aa);
956
957 for (Index jf = 0; jf < nf; jf++) {
958 for (Index jp = 0; jp < np; jp++) {
959 for (Index jla = 0; jla < nla; jla++) {
960 for (Index jlo = 0; jlo < nlo; jlo++) {
961 for (Index jz = 0; jz < nz; jz++) {
962 for (Index js = 0; js < ns; js++) {
963 cloudbox_field(jf, jp, jla, jlo, jz, 0, js) =
964 interp(cloudbox_field_in(jf, jp, jla, jlo, jz, joker, js),
965 itw_aa,
966 lag_aa);
967 }
968 }
969 }
970 }
971 }
972 }
973 }
974}
975
976/* Workspace method: Doxygen documentation will be auto-generated */
977void cloudbox_fieldCrop(Tensor7& cloudbox_field,
978 ArrayOfIndex& cloudbox_limits,
979 const Index& atmosphere_dim,
980 const Index& cloudbox_on,
981 const Index& new_limit0,
982 const Index& new_limit1,
983 const Index& new_limit2,
984 const Index& new_limit3,
985 const Index& new_limit4,
986 const Index& new_limit5,
987 const Verbosity&) {
988 ARTS_USER_ERROR_IF (!cloudbox_on,
989 "No need to use this method with cloudbox=0.");
990 ARTS_USER_ERROR_IF (new_limit0 < cloudbox_limits[0],
991 "new_limits0 < cloudbox_limits[0], not allowed!");
992 ARTS_USER_ERROR_IF (new_limit1 > cloudbox_limits[1],
993 "new_limits1 > cloudbox_limits[1], not allowed!");
994
995 Tensor7 fcopy = cloudbox_field;
996
997 if (atmosphere_dim == 1) {
998 cloudbox_field = fcopy(
999 joker,
1000 Range(new_limit0 - cloudbox_limits[0], new_limit1 - new_limit0 + 1),
1001 joker,
1002 joker,
1003 joker,
1004 joker,
1005 joker);
1006 cloudbox_limits[0] = new_limit0;
1007 cloudbox_limits[1] = new_limit1;
1008 } else {
1009 ARTS_USER_ERROR_IF (new_limit2 < cloudbox_limits[2],
1010 "new_limits2 < cloudbox_limits[2], not allowed!");
1011 ARTS_USER_ERROR_IF (new_limit3 > cloudbox_limits[3],
1012 "new_limits3 > cloudbox_limits[3], not allowed!");
1013
1014 if (atmosphere_dim == 2) {
1015 cloudbox_field = fcopy(
1016 joker,
1017 Range(new_limit0 - cloudbox_limits[0], new_limit1 - new_limit0 + 1),
1018 Range(new_limit2 - cloudbox_limits[2], new_limit3 - new_limit2 - 1),
1019 joker,
1020 joker,
1021 joker,
1022 joker);
1023 cloudbox_limits[0] = new_limit0;
1024 cloudbox_limits[1] = new_limit1;
1025 cloudbox_limits[2] = new_limit2;
1026 cloudbox_limits[3] = new_limit3;
1027 } else {
1028 ARTS_USER_ERROR_IF (new_limit4 < cloudbox_limits[4],
1029 "new_limits4 < cloudbox_limits[4], not allowed!");
1030 ARTS_USER_ERROR_IF (new_limit5 > cloudbox_limits[5],
1031 "new_limits5 > cloudbox_limits[5], not allowed!");
1032 cloudbox_field = fcopy(
1033 joker,
1034 Range(new_limit0 - cloudbox_limits[0], new_limit1 - new_limit0 + 1),
1035 Range(new_limit2 - cloudbox_limits[2], new_limit3 - new_limit2 + 1),
1036 Range(new_limit4 - cloudbox_limits[4], new_limit5 - new_limit4 + 1),
1037 joker,
1038 joker,
1039 joker);
1040 cloudbox_limits[0] = new_limit0;
1041 cloudbox_limits[1] = new_limit1;
1042 cloudbox_limits[2] = new_limit2;
1043 cloudbox_limits[3] = new_limit3;
1044 cloudbox_limits[4] = new_limit4;
1045 cloudbox_limits[5] = new_limit5;
1046 }
1047 }
1048}
1049
1050/* Workspace method: Doxygen documentation will be auto-generated */
1051void particle_fieldCleanup( //WS Output:
1052 Tensor4& particle_field_out,
1053 //WS Input:
1054 const Tensor4& particle_field_in,
1055 const Numeric& threshold,
1056 const Verbosity&) {
1057 if (&particle_field_out != &particle_field_in) {
1058 particle_field_out = particle_field_in;
1059 }
1060
1061 // Check that particle_field contains realistic values
1062 //(values smaller than the threshold will be set to 0)
1063 for (Index i = 0; i < particle_field_out.nbooks(); i++) {
1064 for (Index j = 0; j < particle_field_out.npages(); j++) {
1065 for (Index k = 0; k < particle_field_out.nrows(); k++) {
1066 for (Index l = 0; l < particle_field_out.ncols(); l++) {
1067 if (particle_field_out(i, j, k, l) < threshold) {
1068 particle_field_out(i, j, k, l) = 0.0;
1069 }
1070 }
1071 }
1072 }
1073 }
1074}
1075
1076/* Workspace method: Doxygen documentation will be auto-generated */
1077void ScatSpeciesInit( //WS Output:
1078 ArrayOfString& scat_species,
1081 Index& scat_data_checked,
1082 ArrayOfGriddedField3& pnd_field_raw,
1083 const Verbosity&) {
1084 scat_species.resize(0);
1085 scat_data_raw.resize(0);
1086 scat_meta.resize(0);
1087 pnd_field_raw.resize(0);
1088 scat_data_checked = 0;
1089}
1090
1091/* Workspace method: Doxygen documentation will be auto-generated */
1094 ArrayOfGriddedField3& pnd_field_raw,
1095 // WS Input (needed for checking the datafiles):
1096 const Index& atmosphere_dim,
1097 // Keywords:
1098 const ArrayOfString& scat_data_files,
1099 const ArrayOfString& pnd_field_files,
1100 const Verbosity& verbosity) {
1102
1103 //--- Check input ---------------------------------------------------------
1104
1105 // Atmosphere
1106 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1107 //chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
1108
1109 //--- Reading the data ---------------------------------------------------
1110
1111 ARTS_USER_ERROR_IF (scat_data_files.nelem() != pnd_field_files.nelem(),
1112 "Number of elements in scat_data and pnd_field filelists is"
1113 "inconsistent.")
1114
1115 Index last_species = scat_data_raw.nelem() - 1;
1116 if (last_species == -1) {
1117 scat_data_raw.resize(1);
1118 last_species = 0;
1119 }
1120
1121 // Create empty dummy elements to append to *scat_data_raw* and *pnd_field_raw*.
1122 SingleScatteringData scat_data_single;
1123 GriddedField3 pnd_field_data;
1124
1125 for (Index i = 0; i < scat_data_files.nelem(); i++) {
1126 // Append *scat_data_raw* and *pnd_field_raw* with empty Arrays of Tensors.
1127 scat_data_raw[last_species].push_back(scat_data_single);
1128 pnd_field_raw.push_back(pnd_field_data);
1129
1130 out2 << " Read single scattering data file " << scat_data_files[i] << "\n";
1132 scat_data_files[i],
1133 scat_data_raw[last_species][scat_data_raw[last_species].nelem() - 1],
1134 verbosity);
1135
1136 out2 << " Read particle number density field\n";
1137 if (pnd_field_files[i].nelem() < 1) {
1139 out1 << "Warning: No pnd_field_file specified. Ignored here,\n"
1140 << "but user HAS TO add that later on!. \n";
1141 } else {
1142 xml_read_from_file(pnd_field_files[i],
1143 pnd_field_raw[pnd_field_raw.nelem() - 1],
1144 verbosity);
1145
1146 chk_pnd_data(pnd_field_raw[pnd_field_raw.nelem() - 1],
1147 pnd_field_files[i],
1148 atmosphere_dim,
1149 verbosity);
1150 }
1151 }
1152}
1153
1154/* Workspace method: Doxygen documentation will be auto-generated */
1157 ArrayOfGriddedField3& pnd_field_raw,
1158 // WS Input(needed for checking the datafiles):
1159 const Index& atmosphere_dim,
1160 // Keywords:
1161 const ArrayOfString& scat_data_files,
1162 const String& pnd_fieldarray_file,
1163 const Verbosity& verbosity) {
1165
1166 //--- Check input ---------------------------------------------------------
1167
1168 // Atmosphere
1169 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1170 //chk_atm_grids ( atmosphere_dim, p_grid, lat_grid, lon_grid );
1171
1172 //--- Reading the data ---------------------------------------------------
1174 arr_ssd.resize(scat_data_files.nelem());
1175
1176 for (Index i = 0; i < scat_data_files.nelem(); i++) {
1177 out2 << " Read single scattering data file " << scat_data_files[i] << "\n";
1178 xml_read_from_file(scat_data_files[i], arr_ssd[i], verbosity);
1179 }
1180
1181 // append as new scattering species
1182 if (scat_data_raw.nelem() == 0) {
1183 scat_data_raw.resize(1);
1184 scat_data_raw[0] = arr_ssd;
1185 } else
1186 scat_data_raw.push_back(arr_ssd);
1187
1188 out2 << " Read particle number density data \n";
1189 ArrayOfGriddedField3 pnd_tmp;
1190 xml_read_from_file(pnd_fieldarray_file, pnd_tmp, verbosity);
1191
1192 chk_pnd_raw_data(pnd_tmp, pnd_fieldarray_file, atmosphere_dim, verbosity);
1193
1194 // append to pnd_field_raw
1195 for (Index i = 0; i < pnd_tmp.nelem(); ++i)
1196 pnd_field_raw.push_back(pnd_tmp[i]);
1197}
1198
1199/* Workspace method: Doxygen documentation will be auto-generated */
1202 ArrayOfGriddedField3& vmr_field_raw,
1203 ArrayOfArrayOfSpeciesTag& abs_species,
1204 Index& propmat_clearsky_agenda_checked,
1205 // WS Input (needed for checking the datafiles):
1206 const Index& atmosphere_dim,
1207 const Vector& f_grid,
1208 // Keywords:
1209 const ArrayOfString& scat_data_files,
1210 const ArrayOfString& pnd_field_files,
1211 const Verbosity& verbosity) {
1213
1214 //--- Check input ---------------------------------------------------------
1215
1216 // Atmosphere
1217 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1218 //chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
1219
1220 // Frequency grid
1221 ARTS_USER_ERROR_IF (f_grid.empty(), "The frequency grid is empty.");
1222 chk_if_increasing("f_grid", f_grid);
1223
1224 //--- Reading the data ---------------------------------------------------
1225
1226 ARTS_USER_ERROR_IF (scat_data_files.nelem() != pnd_field_files.nelem(),
1227 "Number of elements in scat_data and pnd_field filelists is"
1228 "inconsistent.")
1229
1230 Index last_species = scat_data_raw.nelem() - 1;
1231 if (last_species == -1) {
1232 scat_data_raw.resize(1);
1233 last_species = 0;
1234 }
1235
1236 // Create empty dummy elements to append to *scat_data_raw* and *pnd_field_raw*.
1237 SingleScatteringData scat_data_single;
1238 GriddedField3 pnd_field_data;
1239 ArrayOfString species(1);
1240 species[0] = "particles";
1241
1242 for (Index i = 0; i < scat_data_files.nelem(); i++) {
1243 // Append *scat_data_raw* and *pnd_field_raw* with empty Arrays of Tensors.
1244 scat_data_raw[last_species].push_back(scat_data_single);
1245 vmr_field_raw.push_back(pnd_field_data);
1246
1247 out2 << " Read single scattering data file " << scat_data_files[i] << "\n";
1249 scat_data_files[i],
1250 scat_data_raw[last_species][scat_data_raw[last_species].nelem() - 1],
1251 verbosity);
1252
1253 out2 << " Check single scattering properties\n";
1255 "scat_data_single.f_grid to f_grid",
1256 scat_data_raw[last_species][scat_data_raw[last_species].nelem() - 1]
1257 .f_grid,
1258 f_grid);
1259
1260 out2 << " Read particle number density field\n";
1261 if (pnd_field_files[i].nelem() < 1) {
1263 out1 << "Warning: No pnd_field_file specified. Ignored here,\n"
1264 << "but user HAS TO add that later on!\n";
1265 } else {
1266 try {
1267 xml_read_from_file(pnd_field_files[i],
1268 vmr_field_raw[vmr_field_raw.nelem() - 1],
1269 verbosity);
1270 } catch (...) {
1272 try {
1273 xml_read_from_file(pnd_field_files[i], tmp, verbosity);
1274 if (tmp.nelem() == 1) {
1275 vmr_field_raw[vmr_field_raw.nelem() - 1] = tmp[0];
1276 } else {
1278 "The file ", pnd_field_files[i], "\n"
1279 "is neither GriddedField3 nor a 1-long ArrayOfGriddedField3.\n")
1280 }
1281 } catch (...) {
1283 "The file ", pnd_field_files[i], " does not exist or\n"
1284 "its type is neither GriddedField3 nor a 1-long ArrayOfGriddedField3.\n")
1285 }
1286 }
1287
1288 chk_pnd_data(vmr_field_raw[vmr_field_raw.nelem() - 1],
1289 pnd_field_files[i],
1290 atmosphere_dim,
1291 verbosity);
1292 }
1293
1294 out2 << " Append 'particle' field to abs_species\n";
1295 abs_speciesAdd(abs_species,
1296 propmat_clearsky_agenda_checked,
1297 species,
1298 verbosity);
1299 }
1300 scat_dataCheck(scat_data_raw, "sane", 1e-2, verbosity);
1301}
1302
1303/* Workspace method: Doxygen documentation will be auto-generated */
1307 // Keywords:
1308 const ArrayOfString& scat_data_files,
1309 const Verbosity& verbosity) {
1312
1313 //--- Reading the data ---------------------------------------------------
1316
1317 arr_ssd.resize(scat_data_files.nelem());
1318 arr_smd.resize(scat_data_files.nelem());
1319
1320 Index meta_naming_conv = 0;
1321
1322 for (Index i = 0; i < 1 && i < scat_data_files.nelem(); i++) {
1323 out3 << " Read single scattering data file " << scat_data_files[i] << "\n";
1324 xml_read_from_file(scat_data_files[i], arr_ssd[i], verbosity);
1325
1326 // make meta data name from scat data name
1327 ArrayOfString strarr;
1328 String scat_meta_file;
1329
1330 if (i == 0) {
1331 scat_data_files[i].split(strarr, ".xml");
1332 scat_meta_file = strarr[0] + ".meta.xml";
1333
1334 try {
1335 find_xml_file(scat_meta_file, verbosity);
1336 } catch (const runtime_error&) {
1337 }
1338
1339 if (file_exists(scat_meta_file)) {
1340 out3 << " Read scattering meta data\n";
1341
1342 xml_read_from_file(scat_meta_file, arr_smd[i], verbosity);
1343
1344 meta_naming_conv = 1;
1345 } else {
1346 try {
1347 scat_data_files[i].split(strarr, "scat_data");
1348 ARTS_USER_ERROR_IF (strarr.nelem() < 2,
1349 "Splitting scattering data filename up at 'scat_data' also failed.");
1350 scat_meta_file = strarr[0] + "scat_meta" + strarr[1];
1351
1352 out3 << " Read scattering meta data\n";
1353 xml_read_from_file(scat_meta_file, arr_smd[i], verbosity);
1354
1355 meta_naming_conv = 2;
1356 } catch (const std::runtime_error& e) {
1358 "No meta data file following one of the allowed naming "
1359 "conventions was found.\n"
1360 "Allowed are "
1361 "*.meta.xml from *.xml and "
1362 "*scat_meta* from *scat_data*\n"
1363 "Scattering meta data file not found: ", scat_meta_file,
1364 "\n",
1365 e.what())
1366 }
1367 }
1368 }
1369 }
1370
1371 ArrayOfString fail_msg;
1372
1373#pragma omp parallel for if (!arts_omp_in_parallel() && \
1374 scat_data_files.nelem() > 1) \
1375 num_threads(arts_omp_get_max_threads() > 16 ? 16 \
1376 : arts_omp_get_max_threads()) \
1377 shared(out3, scat_data_files, arr_ssd, arr_smd)
1378 for (Index i = 1; i < scat_data_files.nelem(); i++) {
1379 // make meta data name from scat data name
1380 ArrayOfString strarr;
1381 String scat_meta_file;
1384
1385 try {
1386 out3 << " Read single scattering data file " << scat_data_files[i]
1387 << "\n";
1388 xml_read_from_file(scat_data_files[i], ssd, verbosity);
1389
1390 scat_data_files[i].split(strarr, ".xml");
1391 scat_meta_file = strarr[0] + ".meta.xml";
1392
1393 if (meta_naming_conv == 1) {
1394 scat_data_files[i].split(strarr, ".xml");
1395 scat_meta_file = strarr[0] + ".meta.xml";
1396
1397 out3 << " Read scattering meta data\n";
1398 xml_read_from_file(scat_meta_file, smd, verbosity);
1399 } else {
1400 scat_data_files[i].split(strarr, "scat_data");
1401 scat_meta_file = strarr[0] + "scat_meta" + strarr[1];
1402
1403 out3 << " Read scattering meta data\n";
1404 xml_read_from_file(scat_meta_file, smd, verbosity);
1405 }
1406 } catch (const std::exception& e) {
1407 ostringstream os;
1408 os << "Run-time error reading scattering data : \n" << e.what();
1409#pragma omp critical(ybatchCalc_push_fail_msg)
1410 fail_msg.push_back(os.str());
1411 }
1412
1413 //FIXME: currently nothing is done in chk_scattering_meta_data!
1414 chk_scattering_meta_data(smd, scat_meta_file, verbosity);
1415
1416#pragma omp critical(ScatSpeciesScatAndMetaRead_assign_ssd)
1417 arr_ssd[i] = std::move(ssd);
1418#pragma omp critical(ScatSpeciesScatAndMetaRead_assign_smd)
1419 arr_smd[i] = std::move(smd);
1420 }
1421
1422 if (fail_msg.nelem()) {
1423 std::ostringstream os;
1424 for (auto& msg : fail_msg) os << msg << '\n';
1425
1426 ARTS_USER_ERROR ( os.str());
1427 }
1428
1429 // check if arrays have same size
1430 chk_scattering_data(arr_ssd, arr_smd, verbosity);
1431
1432 // append as new scattering species
1433 scat_data_raw.push_back(std::move(arr_ssd));
1434 scat_meta.push_back(std::move(arr_smd));
1435}
1436
1437/* Workspace method: Doxygen documentation will be auto-generated */
1438void ScatElementsSelect( //WS Output:
1441 // WS Input:
1442 const ArrayOfString& scat_species,
1443 const String& species,
1444 const String& sizeparam,
1445 const Numeric& sizemin,
1446 const Numeric& sizemax,
1447 const Numeric& tolerance,
1448 const String& delim,
1449 const Verbosity&) {
1450 // first check that sizes of scat_species and scat_data_raw/scat_meta agree
1451 Index nspecies = scat_species.nelem();
1452 ARTS_USER_ERROR_IF (nspecies != scat_data_raw.nelem() || nspecies != scat_meta.nelem(),
1453 "Number of scattering species specified by scat_species does\n"
1454 "not agree with number of scattering species in\n"
1455 "scat_data_raw or scat_meta:\n"
1456 "scat_species has ", nspecies,
1457 " entries, while scat_data_raw has ", scat_data_raw.nelem(),
1458 " and scat_meta has ", scat_meta.nelem(), ".")
1459
1460 // create temporary containers for selected elements
1461 ArrayOfSingleScatteringData scat_data_raw_tmp;
1462 ArrayOfScatteringMetaData scat_meta_tmp;
1463
1464 String partfield_name;
1465 //find the species to handle: compare 'species' to 'partfield' part of
1466 //scat_species tags
1467 Index i_ss = -1;
1468 for (Index i = 0; i < scat_species.nelem(); i++) {
1469 parse_partfield_name(partfield_name, scat_species[i], delim);
1470 if (partfield_name == species) i_ss = i;
1471 }
1472 ARTS_USER_ERROR_IF (i_ss < 0,
1473 "Scattering species ", species, " not found among scat_species.")
1474
1475 // choosing the specified SingleScatteringData and ScatteringMetaData
1476 if (sizeparam == "diameter_max")
1477 for (Index i_se = 0; i_se < scat_meta[i_ss].nelem(); i_se++) {
1478 // scattering element diameter is extracted from the
1479 // scattering element's meta data and checked whether it's within size
1480 // selected range (sizemax < 0 check follows from wildcard usage and
1481 // means consider all sizes on the upper end)
1482 if (scat_meta[i_ss][i_se].diameter_max > sizemin - sizemin * tolerance &&
1483 (sizemax + sizemax * tolerance > scat_meta[i_ss][i_se].diameter_max ||
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_volume_equ")
1491 for (Index i_se = 0; i_se < scat_meta[i_ss].nelem(); i_se++) {
1492 if (scat_meta[i_ss][i_se].diameter_volume_equ >
1493 sizemin - sizemin * tolerance &&
1494 (sizemax + sizemax * tolerance >
1495 scat_meta[i_ss][i_se].diameter_volume_equ ||
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 if (sizeparam == "diameter_area_equ_aerodynamical")
1503 for (Index i_se = 0; i_se < scat_meta[i_ss].nelem(); i_se++) {
1504 if (scat_meta[i_ss][i_se].diameter_area_equ_aerodynamical >
1505 sizemin - sizemin * tolerance &&
1506 (sizemax + sizemax * tolerance >
1507 scat_meta[i_ss][i_se].diameter_area_equ_aerodynamical ||
1508 sizemax < 0.)) {
1509 // copy selected scattering element to temp arrays
1510 scat_data_raw_tmp.push_back(scat_data_raw[i_ss][i_se]);
1511 scat_meta_tmp.push_back(scat_meta[i_ss][i_se]);
1512 }
1513 }
1514 else {
1516 "Size parameter ", sizeparam, "is unknown.")
1517 }
1518
1519 // To use a particle species field without associated scattering element
1520 // data poses a high risk of accidentially neglecting these species. That's
1521 // unlikely what the user intends. Hence throw error.
1522 ARTS_USER_ERROR_IF (scat_meta_tmp.nelem() < 1,
1523 "For scattering species ", species, " no scattering "
1524 "element matching the requested size range found.\n"
1525 "Check *scat_data_raw* and *scat_meta* input as well as your size limit "
1526 "selection!")
1527
1528 scat_meta[i_ss] = std::move(scat_meta_tmp);
1529 scat_data_raw[i_ss] = std::move(scat_data_raw_tmp);
1530
1531 // check if array is empty. should never apply (since we checked the re-worked
1532 // data before and that error should also catch cases that are empty from the
1533 // beginning).
1535}
1536
1537/* Workspace method: Doxygen documentation will be auto-generated */
1540 // Keywords:
1541 const ArrayOfString& scat_species,
1542 const String& species,
1543 const String& scat_species_delim,
1544 const Numeric& T_low,
1545 const Numeric& T_high,
1546 const Verbosity&) {
1547 const bool do_tl = (T_low >= 0.);
1548 const bool do_th = (T_high >= 0.);
1549
1550 if (do_tl || do_th) {
1551 Index i_ss = -1;
1552 if (species == "") {
1553 i_ss = scat_data_raw.nelem() - 1;
1554 ARTS_USER_ERROR_IF (i_ss == -1,
1555 "No *scat_data* available. Can not extend temperature range on "
1556 "inexistent data.")
1557 } else {
1558 // first check that sizes of scat_species and scat_data_raw agree
1559 Index nspecies = scat_species.nelem();
1560 ARTS_USER_ERROR_IF (nspecies != scat_data_raw.nelem(),
1561 "Number of scattering species specified by scat_species does\n"
1562 "not agree with number of scattering species in *scat_data*:\n"
1563 "scat_species has ", nspecies,
1564 " entries, while *scat_data* has ", scat_data_raw.nelem(),
1565 ".")
1566 String partfield_name;
1567 //find the species to handle: compare 'species' to 'partfield' part of
1568 //scat_species tags
1569 for (Index i = 0; i < scat_species.nelem(); i++) {
1571 partfield_name, scat_species[i], scat_species_delim);
1572 if (partfield_name == species) i_ss = i;
1573 }
1574 ARTS_USER_ERROR_IF (i_ss < 0,
1575 "Scattering species ", species, " not found among scat_species.")
1576 }
1577
1578 for (Index i_se = 0; i_se < scat_data_raw[i_ss].nelem(); i_se++) {
1579 const SingleScatteringData& ssdo = scat_data_raw[i_ss][i_se];
1580 const Index nTo = ssdo.T_grid.nelem();
1581 Index nTn = nTo;
1582 bool do_htl, do_hth;
1583 if (nTo > 1) {
1584 do_htl = (do_tl && (T_low < ssdo.T_grid[0]));
1585 do_hth = (do_th && (T_high > last(ssdo.T_grid)));
1586 } else {
1587 do_htl = false;
1588 do_hth = false;
1589 }
1590
1591 if (do_htl || do_hth) {
1592 // create new instance of SingleScatteringData
1594 Index iToff;
1595
1596 // determine new temperature grid
1597 if (do_htl) nTn += 1;
1598 if (do_hth) nTn += 1;
1599 Vector T_grid_new(nTn);
1600 if (do_htl) {
1601 T_grid_new[0] = T_low;
1602 iToff = 1;
1603 } else {
1604 iToff = 0;
1605 }
1606 for (Index iT = 0; iT < nTo; iT++)
1607 T_grid_new[iT + iToff] = scat_data_raw[i_ss][i_se].T_grid[iT];
1608 if (do_hth) T_grid_new[nTo + iToff] = T_high;
1609 ssdn.T_grid = std::move(T_grid_new);
1610
1611 // copy grids and other descriptive data that is to remain identical
1612 ssdn.ptype = ssdo.ptype;
1613 ostringstream description;
1614 description << ssdo.description; // here just copy. we append further
1615 // info below if applicable.
1616 ssdn.f_grid = ssdo.f_grid;
1617 ssdn.za_grid = ssdo.za_grid;
1618 ssdn.aa_grid = ssdo.aa_grid;
1619
1620 // determine size of current optical property data
1621 const Index nf = ssdo.f_grid.nelem();
1622 const Index nzas = ssdo.pha_mat_data.nshelves();
1623 const Index naas = ssdo.pha_mat_data.nbooks();
1624 const Index nzai = ssdo.pha_mat_data.npages();
1625 const Index naai = ssdo.pha_mat_data.nrows();
1626 const Index nmep = ssdo.pha_mat_data.ncols();
1627 const Index nmee = ssdo.ext_mat_data.ncols();
1628 const Index nvea = ssdo.abs_vec_data.ncols();
1629
1630 // create containers for extended optical property data
1631 ssdn.pha_mat_data.resize(nf, nTn, nzas, naas, nzai, naai, nmep);
1632 ssdn.ext_mat_data.resize(nf, nTn, nzai, naai, nmee);
1633 ssdn.abs_vec_data.resize(nf, nTn, nzai, naai, nvea);
1634
1635 // copy optical property data
1636 for (Index iT = 0; iT < nTo; iT++) {
1637 ssdn.pha_mat_data(
1638 joker, iT + iToff, joker, joker, joker, joker, joker) =
1639 ssdo.pha_mat_data(joker, iT, joker, joker, joker, joker, joker);
1640 ssdn.ext_mat_data(joker, iT + iToff, joker, joker, joker) =
1641 ssdo.ext_mat_data(joker, iT, joker, joker, joker);
1642 ssdn.abs_vec_data(joker, iT + iToff, joker, joker, joker) =
1643 ssdo.abs_vec_data(joker, iT, joker, joker, joker);
1644 }
1645
1646 // duplicate optical property data on T-edges if applicable
1647 if (do_htl) {
1648 ssdn.pha_mat_data(joker, 0, joker, joker, joker, joker, joker) =
1650 ssdn.ext_mat_data(joker, 0, joker, joker, joker) =
1651 ssdn.ext_mat_data(joker, 1, joker, joker, joker);
1652 ssdn.abs_vec_data(joker, 0, joker, joker, joker) =
1653 ssdn.abs_vec_data(joker, 1, joker, joker, joker);
1654 description << "\n"
1655 << "Low temperature limit extended by"
1656 << " duplicating previous low temperature limit"
1657 << " single scattering properties.";
1658 }
1659 if (do_hth) {
1660 ssdn.pha_mat_data(joker, nTn - 1, joker, joker, joker, joker, joker) =
1661 ssdn.pha_mat_data(
1662 joker, nTn - 2, joker, joker, joker, joker, joker);
1663 ssdn.ext_mat_data(joker, nTn - 1, joker, joker, joker) =
1664 ssdn.ext_mat_data(joker, nTn - 2, joker, joker, joker);
1665 ssdn.abs_vec_data(joker, nTn - 1, joker, joker, joker) =
1666 ssdn.abs_vec_data(joker, nTn - 2, joker, joker, joker);
1667 description << "\n"
1668 << "High temperature limit extended by"
1669 << " duplicating previous high temperature limit"
1670 << " single scattering properties.";
1671 }
1672 ssdn.description = description.str();
1673 scat_data_raw[i_ss][i_se] = std::move(ssdn);
1674 }
1675 }
1676 }
1677}
1678
1679/* Workspace method: Doxygen documentation will be auto-generated */
1681 Tensor4& pnd_field,
1682 ArrayOfTensor4& dpnd_field_dx,
1683 //WS Input
1684 const Vector& p_grid,
1685 const Vector& lat_grid,
1686 const Vector& lon_grid,
1687 const ArrayOfGriddedField3& pnd_field_raw,
1688 const Index& atmosphere_dim,
1689 const ArrayOfIndex& cloudbox_limits,
1690 const ArrayOfRetrievalQuantity& jacobian_quantities,
1691 const Index& zeropadding,
1692 const Verbosity& verbosity) {
1693 // Basic checks of input variables
1694 //
1695 // Particle number density data
1696 //
1697 ARTS_USER_ERROR_IF (pnd_field_raw.empty(),
1698 "No particle number density data given. Please use WSMs\n"
1699 "*ParticleTypeInit* and *ParticleTypeAdd(All)* for reading\n"
1700 "scattering element data.\n")
1701
1702 chk_atm_grids(atmosphere_dim, p_grid, lat_grid, lon_grid);
1703 ArrayOfIndex cloudbox_limits_tmp;
1704 /*if ( cloudbox_limits.empty() )
1705 {
1706 //If no limits set, the cloud(box) is supposed to cover the
1707 //complete atmosphere. This particularly to facilitate use of
1708 //scat_data_single&pnd_fields for non-scatt, greybody cloud calculations.
1709 //Hence, set the limits to first & last elements of the different grids.
1710 //Note: no margin left at lat/lon_grid edges!.
1711 cloudbox_limits_tmp.resize(2*atmosphere_dim);
1712
1713 // Pressure limits
1714 cloudbox_limits_tmp[0] = 0;
1715 cloudbox_limits_tmp[1] = p_grid.nelem() - 1;
1716 // Latitude limits
1717 if( atmosphere_dim >= 2 )
1718 {
1719 cloudbox_limits_tmp[2] = 0;
1720 cloudbox_limits_tmp[3] = lat_grid.nelem() - 1;
1721 }
1722 // Longitude limits
1723 if( atmosphere_dim == 3 )
1724 {
1725 cloudbox_limits_tmp[4] = 0;
1726 cloudbox_limits_tmp[5] = lon_grid.nelem() - 1;
1727 }
1728 }
1729 else */
1730 ARTS_USER_ERROR_IF (cloudbox_limits.nelem() != 2 * atmosphere_dim,
1731 "*cloudbox_limits* is a vector which contains the"
1732 "upper and lower limit of the cloud for all "
1733 "atmospheric dimensions. So its dimension must"
1734 "be 2 x *atmosphere_dim*");
1735 cloudbox_limits_tmp = cloudbox_limits;
1736
1737 // Check that pnd_field_raw has at least 2 grid-points in each dimension.
1738 // Otherwise, interpolation further down will fail with assertion.
1739
1740 for (Index d = 0; d < atmosphere_dim; d++) {
1741 for (Index i = 0; i < pnd_field_raw.nelem(); i++) {
1742 ARTS_USER_ERROR_IF (pnd_field_raw[i].get_grid_size(d) < 2,
1743 "Error in pnd_field_raw data. "
1744 "Dimension ", d, " (name: \"",
1745 pnd_field_raw[i].get_grid_name(d),
1746 "\") has only ",
1747 pnd_field_raw[i].get_grid_size(d),
1748 " element",
1749 ((pnd_field_raw[i].get_grid_size(d) == 1) ? "" : "s"),
1750 ". Must be at least 2.")
1751 }
1752 }
1753 const Index Np_cloud = cloudbox_limits_tmp[1] - cloudbox_limits_tmp[0] + 1;
1754
1755 ConstVectorView p_grid_cloud =
1756 p_grid[Range(cloudbox_limits_tmp[0], Np_cloud)];
1757
1758 // Check that no scatterers exist outside the cloudbox
1760 pnd_field_raw,
1761 p_grid,
1762 lat_grid,
1763 lon_grid,
1764 cloudbox_limits_tmp);
1765
1766 //==========================================================================
1767 if (atmosphere_dim == 1) {
1768 ArrayOfGriddedField3 pnd_field_tmp;
1769
1771 pnd_field_tmp, p_grid_cloud, pnd_field_raw, 1, zeropadding, verbosity);
1772
1773 FieldFromGriddedField(pnd_field,
1774 p_grid_cloud,
1775 pnd_field_tmp[0].get_numeric_grid(1),
1776 pnd_field_tmp[0].get_numeric_grid(2),
1777 pnd_field_tmp,
1778 verbosity);
1779 } else if (atmosphere_dim == 2) {
1780 const Index Nlat_cloud =
1781 cloudbox_limits_tmp[3] - cloudbox_limits_tmp[2] + 1;
1782
1783 ConstVectorView lat_grid_cloud =
1784 lat_grid[Range(cloudbox_limits_tmp[2], Nlat_cloud)];
1785
1786 if (zeropadding) {
1787 // FIXME: OLE: Implement this
1789 out0 << "WARNING: zeropadding currently only supported for 1D.";
1790 }
1791
1792 //Resize variables
1793 pnd_field.resize(pnd_field_raw.nelem(), Np_cloud, Nlat_cloud, 1);
1794
1795 // Gridpositions:
1796 ArrayOfGridPos gp_p(Np_cloud);
1797 ArrayOfGridPos gp_lat(Nlat_cloud);
1798
1799 // Interpolate pnd_field.
1800 // Loop over the scattering element:
1801 for (Index i = 0; i < pnd_field_raw.nelem(); ++i) {
1802 // Calculate grid positions:
1803 p2gridpos(gp_p,
1804 pnd_field_raw[i].get_numeric_grid(GFIELD3_P_GRID),
1805 p_grid_cloud);
1806 gridpos(gp_lat,
1807 pnd_field_raw[i].get_numeric_grid(GFIELD3_LAT_GRID),
1808 lat_grid_cloud);
1809
1810 // Interpolation weights:
1811 Tensor3 itw(Np_cloud, Nlat_cloud, 4);
1812 interpweights(itw, gp_p, gp_lat);
1813
1814 // Interpolate:
1815 interp(pnd_field(i, joker, joker, 0),
1816 itw,
1817 pnd_field_raw[i].data(joker, joker, 0),
1818 gp_p,
1819 gp_lat);
1820 }
1821 } else {
1822 const Index Nlat_cloud =
1823 cloudbox_limits_tmp[3] - cloudbox_limits_tmp[2] + 1;
1824 const Index Nlon_cloud =
1825 cloudbox_limits_tmp[5] - cloudbox_limits_tmp[4] + 1;
1826
1827 if (zeropadding) {
1828 // FIXME: OLE: Implement this
1830 out0 << "WARNING: zeropadding currently only supported for 1D.";
1831 }
1832
1833 ConstVectorView lat_grid_cloud =
1834 lat_grid[Range(cloudbox_limits_tmp[2], Nlat_cloud)];
1835 ConstVectorView lon_grid_cloud =
1836 lon_grid[Range(cloudbox_limits_tmp[4], Nlon_cloud)];
1837
1838 //Resize variables
1839 pnd_field.resize(pnd_field_raw.nelem(), Np_cloud, Nlat_cloud, Nlon_cloud);
1840
1841 // Gridpositions:
1842 ArrayOfGridPos gp_p(Np_cloud);
1843 ArrayOfGridPos gp_lat(Nlat_cloud);
1844 ArrayOfGridPos gp_lon(Nlon_cloud);
1845
1846 // Interpolate pnd_field.
1847 // Loop over the scattering element types:
1848 for (Index i = 0; i < pnd_field_raw.nelem(); ++i) {
1849 // Calculate grid positions:
1850 p2gridpos(gp_p,
1851 pnd_field_raw[i].get_numeric_grid(GFIELD3_P_GRID),
1852 p_grid_cloud);
1853 gridpos(gp_lat,
1854 pnd_field_raw[i].get_numeric_grid(GFIELD3_LAT_GRID),
1855 lat_grid_cloud);
1856 gridpos(gp_lon,
1857 pnd_field_raw[i].get_numeric_grid(GFIELD3_LON_GRID),
1858 lon_grid_cloud);
1859
1860 // Interpolation weights:
1861 Tensor4 itw(Np_cloud, Nlat_cloud, Nlon_cloud, 8);
1862 interpweights(itw, gp_p, gp_lat, gp_lon);
1863
1864 // Interpolate:
1865 interp(pnd_field(i, joker, joker, joker),
1866 itw,
1867 pnd_field_raw[i].data,
1868 gp_p,
1869 gp_lat,
1870 gp_lon);
1871 }
1872 }
1873
1874 // no (cloudy) Jacobians with this WSM, hence no calc.
1875 // but we need to size dpnd_field to be consistent with jacobian_quantities.
1876 dpnd_field_dx.resize(jacobian_quantities.nelem());
1877}
1878
1879/* Workspace method: Doxygen documentation will be auto-generated */
1881 const Index& atmosphere_dim,
1882 const Index& cloudbox_on,
1883 const ArrayOfIndex& cloudbox_limits,
1884 const Index& nzero,
1885 const Verbosity&) {
1886 /* if( !cloudbox_checked )
1887 throw runtime_error( "The cloudbox must be flagged to have passed a "
1888 "consistency check (cloudbox_checked=1)." );
1889*/
1890 ARTS_USER_ERROR_IF (atmosphere_dim == 1,
1891 "No use in calling this method for 1D.");
1892 ARTS_USER_ERROR_IF (!cloudbox_on,
1893 "No use in calling this method with cloudbox off.");
1894 ARTS_USER_ERROR_IF (nzero < 1,
1895 "The argument *nzero* must be > 0.");
1896
1897 // Sizes
1898 const Index npart = pnd_field.nbooks();
1899 const Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
1900 const Index nlat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
1901 Index nlon = 1;
1902 if (atmosphere_dim == 3) {
1903 nlon = cloudbox_limits[5] - cloudbox_limits[4] + 1;
1904 }
1905
1906 ARTS_USER_ERROR_IF (pnd_field.npages() != np ||
1907 pnd_field.nrows() != 1 ||
1908 pnd_field.ncols() != 1,
1909 "The input *pnd_field* is either not 1D or does not "
1910 "match pressure size of cloudbox.");
1911
1912 // Temporary container
1913 Tensor4 pnd_temp = pnd_field;
1914
1915 // Resize and fill
1916 pnd_field.resize(npart, np, nlat, nlon);
1917 pnd_field = 0;
1918 //
1919 for (Index ilon = nzero; ilon < nlon - nzero; ilon++) {
1920 for (Index ilat = nzero; ilat < nlat - nzero; ilat++) {
1921 for (Index ip = 0; ip < np; ip++) {
1922 for (Index is = 0; is < npart; is++) {
1923 pnd_field(is, ip, ilat, ilon) = pnd_temp(is, ip, 0, 0);
1924 }
1925 }
1926 }
1927 }
1928}
1929
1930/* Workspace method: Doxygen documentation will be auto-generated */
1931void pnd_fieldZero( //WS Output:
1932 Tensor4& pnd_field,
1933 ArrayOfTensor4& dpnd_field_dx,
1935 //WS Input:
1936 const Index& atmosphere_dim,
1937 const Vector& f_grid,
1938 const ArrayOfIndex& cloudbox_limits,
1939 const ArrayOfRetrievalQuantity& jacobian_quantities,
1940 const Verbosity&) {
1941 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1942
1943 ARTS_USER_ERROR_IF (cloudbox_limits.nelem() != 2 * atmosphere_dim,
1944 "*cloudbox_limits* is a vector which contains the"
1945 "upper and lower limit of the cloud for all "
1946 "atmospheric dimensions. So its dimension must"
1947 "be 2 x *atmosphere_dim*");
1948
1949 // Resize pnd_field and set it to 0:
1950 Index np = cloudbox_limits[1] - cloudbox_limits[0] + 1;
1951 Index nlat = 1, nlon = 1;
1952 if (atmosphere_dim > 1) {
1953 nlat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
1954 if (atmosphere_dim > 2) {
1955 nlon = cloudbox_limits[5] - cloudbox_limits[4] + 1;
1956 }
1957 }
1958
1959 // no (cloudy) Jacobians with this WSM, hence no setting.
1960 // but we need to size dpnd_field to be consistent with jacobian_quantities.
1961 dpnd_field_dx.resize(jacobian_quantities.nelem());
1962
1963 // Do only reset scat_data if it has not been set yet.
1964 // There's no need otherwise, and it's rather unpractical for testing when
1965 // doing so: we might want to do some actual calcs with the scat_data later
1966 // on. So why throw it away?
1967 const Index N_se = TotalNumberOfElements(scat_data);
1968 if (N_se > 0) {
1969 pnd_field.resize(N_se, np, nlat, nlon);
1970 } else {
1971 pnd_field.resize(1, np, nlat, nlon);
1972
1973 //Resize scat_data and set it to 0:
1974 // Number of scattering elements
1975 scat_data.resize(1);
1976 scat_data[0].resize(1);
1977 scat_data[0][0].ptype = PTYPE_TOTAL_RND;
1978 scat_data[0][0].description = " ";
1979 // Grids which contain full ranges which one wants to calculate
1980 Index nf = f_grid.nelem();
1981 scat_data[0][0].f_grid.resize(nf);
1982 scat_data[0][0].f_grid = f_grid;
1983 Index nT = 1;
1984 scat_data[0][0].T_grid.resize(nT);
1985 scat_data[0][0].T_grid = 270.;
1986 Index nza = 5;
1987 nlinspace(scat_data[0][0].za_grid, 0, 180, nza);
1988 // Resize the data arrays
1989 scat_data[0][0].pha_mat_data.resize(nf, nT, nza, 1, 1, 1, 6);
1990 scat_data[0][0].pha_mat_data = 0.;
1991 scat_data[0][0].ext_mat_data.resize(nf, nT, 1, 1, 1);
1992 scat_data[0][0].ext_mat_data = 0.;
1993 scat_data[0][0].abs_vec_data.resize(nf, nT, 1, 1, 1);
1994 scat_data[0][0].abs_vec_data = 0.;
1995 }
1996
1997 pnd_field = 0.;
1998}
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
base max(const Array< base > &x)
Max function.
Definition: array.h:145
base min(const Array< base > &x)
Min function.
Definition: array.h:161
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:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
Index ncols() const noexcept
Definition: matpackIV.h:146
bool empty() const noexcept
Definition: matpackIV.h:152
Index nrows() const noexcept
Definition: matpackIV.h:145
Index nbooks() const noexcept
Definition: matpackIV.h:143
Index npages() const noexcept
Definition: matpackIV.h:144
Index ncols() const noexcept
Definition: matpackV.h:158
Index ncols() const noexcept
Definition: matpackVII.h:164
Index npages() const noexcept
Definition: matpackVII.h:162
Index nrows() const noexcept
Definition: matpackVII.h:163
Index nlibraries() const noexcept
Definition: matpackVII.h:158
Index nvitrines() const noexcept
Definition: matpackVII.h:159
Index nshelves() const noexcept
Definition: matpackVII.h:160
Index nbooks() const noexcept
Definition: matpackVII.h:161
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
bool empty() const noexcept
Returns true if variable size is zero.
Definition: matpackI.h:536
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
The range class.
Definition: matpackI.h:160
The Tensor3 class.
Definition: matpackIII.h:352
The Tensor4 class.
Definition: matpackIV.h:435
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:2407
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:910
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:71
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
#define ARTS_USER_ERROR(...)
Definition: debug.h:169
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
void find_xml_file(String &filename, const Verbosity &verbosity)
Find an xml file.
Definition: file.cc:355
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
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:977
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:1538
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:1200
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:1155
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:1438
void ScatSpeciesScatAndMetaRead(ArrayOfArrayOfSingleScatteringData &scat_data_raw, ArrayOfArrayOfScatteringMetaData &scat_meta, const ArrayOfString &scat_data_files, const Verbosity &verbosity)
WORKSPACE METHOD: ScatSpeciesScatAndMetaRead.
Definition: m_cloudbox.cc:1304
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:1092
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:1051
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:1880
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:1931
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:1680
constexpr Numeric PI
Definition: m_cloudbox.cc:74
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:913
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:1077
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
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
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.)
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 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.
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