ARTS 2.5.0 (git: 9ee3ac6c)
m_ppath.cc
Go to the documentation of this file.
1/* Copyright (C) 2002-2012 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
29/*===========================================================================
30 === External declarations
31 ===========================================================================*/
32
33#include <cmath>
34#include "arts.h"
35#include "auto_md.h"
36#include "check_input.h"
37#include "geodetic.h"
38#include "lin_alg.h"
39#include "m_general.h"
40#include "m_xml.h"
41#include "math_funcs.h"
42#include "messages.h"
43#include "ppath.h"
44#include "refraction.h"
45#include "special_interp.h"
46#include "xml_io.h"
47
48extern const Numeric RAD2DEG;
49extern const Numeric DEG2RAD;
50
51/*===========================================================================
52 === The functions (in alphabetical order)
53 ===========================================================================*/
54
55/* Workspace method: Doxygen documentation will be auto-generated */
56void AddZaAa(Matrix& new_los,
57 const Vector& ref_los,
58 const Matrix& dlos,
59 const Verbosity&) {
60 ARTS_USER_ERROR_IF (ref_los.nelem() != 2,
61 "*ref_los* must have two columns.");
62 ARTS_USER_ERROR_IF (dlos.ncols() != 2,
63 "*dlos* must have two columns.");
64
65 const Index nlos = dlos.nrows();
66
67 new_los.resize(nlos, 2);
68
69 for (Index i = 0; i < nlos; i++) {
70 add_za_aa(new_los(i, 0),
71 new_los(i, 1),
72 ref_los[0],
73 ref_los[1],
74 dlos(i, 0),
75 dlos(i, 1));
76 }
77}
78
79/* Workspace method: Doxygen documentation will be auto-generated */
80void DiffZaAa(Matrix& dlos,
81 const Vector& ref_los,
82 const Matrix& other_los,
83 const Verbosity&) {
84 ARTS_USER_ERROR_IF (ref_los.nelem() != 2,
85 "*ref_los* must have two columns.");
86 ARTS_USER_ERROR_IF (other_los.ncols() != 2,
87 "*other_los* must have two columns.");
88
89 const Index nlos = other_los.nrows();
90
91 dlos.resize(nlos, 2);
92
93 for (Index i = 0; i < nlos; i++) {
94 diff_za_aa(dlos(i, 0),
95 dlos(i, 1),
96 ref_los[0],
97 ref_los[1],
98 other_los(i, 0),
99 other_los(i, 1));
100 }
101}
102
103/* Workspace method: Doxygen documentation will be auto-generated */
105 const Ppath& ppath,
106 const Verbosity& verbosity) {
107 geo_pos.resize(5);
108 geo_pos = NAN;
109
110 geo_pos[Range(0, ppath.pos.ncols())] =
111 ppath.pos(ppath.np - 1, Range(0, ppath.pos.ncols()));
112 geo_pos[Range(3, ppath.los.ncols())] =
113 ppath.los(ppath.np - 1, Range(0, ppath.los.ncols()));
114
116 out2 << " Sets geo-position to:\n" << geo_pos;
117}
118
119/* Workspace method: Doxygen documentation will be auto-generated */
121 const Ppath& ppath,
122 const Verbosity& verbosity) {
123 geo_pos.resize(5);
124 geo_pos = NAN;
125
126 // Take first point of ppath as first guess
127 geo_pos[Range(0, ppath.pos.ncols())] =
128 ppath.pos(0, Range(0, ppath.pos.ncols()));
129 geo_pos[Range(3, ppath.los.ncols())] =
130 ppath.los(0, Range(0, ppath.los.ncols()));
131
132 for (Index i = 1; i < ppath.np; i++) {
133 if (ppath.pos(i, 0) < geo_pos[0]) {
134 geo_pos[Range(0, ppath.pos.ncols())] =
135 ppath.pos(i, Range(0, ppath.pos.ncols()));
136 geo_pos[Range(3, ppath.los.ncols())] =
137 ppath.los(i, Range(0, ppath.los.ncols()));
138 }
139 }
140
142 out2 << " Sets geo-position to:\n" << geo_pos;
143}
144
145/* Workspace method: Doxygen documentation will be auto-generated */
147 const Ppath& ppath,
148 const Numeric& z_ref,
149 const Verbosity& verbosity) {
150 geo_pos.resize(5);
151 geo_pos = NAN;
152
153 bool found = false;
154 Index ihit = 0;
155 bool above = false;
156
157 if (ppath.pos(0, 0) >= z_ref) {
158 above = true;
159 }
160
161 while (!found && ihit < ppath.np - 1) {
162 ihit += 1;
163 if (above && ppath.pos(ihit, 0) < z_ref) {
164 found = true;
165 } else if (!above && ppath.pos(ihit, 0) >= z_ref) {
166 found = true;
167 }
168 }
169
170 if (found) {
171 geo_pos[0] = z_ref;
172
173 // Make a simple linear interpolation to determine lat and lon
174 const Numeric w = (z_ref - ppath.pos(ihit - 1, 0)) /
175 (ppath.pos(ihit, 0) - ppath.pos(ihit - 1, 0));
176
177 geo_pos[3] = w * ppath.los(ihit, 0) + (1 - w) * ppath.los(ihit - 1, 0);
178
179 if (ppath.pos.ncols() > 1) {
180 geo_pos[1] = w * ppath.pos(ihit, 1) + (1 - w) * ppath.pos(ihit - 1, 1);
181
182 if (ppath.pos.ncols() > 2) {
183 geo_pos[2] = w * ppath.pos(ihit, 2) + (1 - w) * ppath.pos(ihit - 1, 2);
184 geo_pos[4] = w * ppath.los(ihit, 1) + (1 - w) * ppath.los(ihit - 1, 1);
185 }
186 }
187 }
188
190 out2 << " Sets geo-position to:\n" << geo_pos;
191}
192
193/* Workspace method: Doxygen documentation will be auto-generated */
195 Ppath& ppath,
196 const Agenda& ppath_agenda,
197 const Numeric& ppath_lmax,
198 const Numeric& ppath_lraytrace,
199 const Index& atmgeom_checked,
200 const Vector& f_grid,
201 const Index& cloudbox_on,
202 const Index& cloudbox_checked,
203 const Index& ppath_inside_cloudbox_do,
204 const Vector& rte_pos,
205 const Vector& rte_los,
206 const Vector& rte_pos2,
207 const Verbosity&) {
208 // Basics
209 //
210 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
211 "The atmospheric geometry must be flagged to have "
212 "passed a consistency check (atmgeom_checked=1).");
213 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
214 "The cloudbox must be flagged to have "
215 "passed a consistency check (cloudbox_checked=1).");
216
218 ppath,
219 ppath_lmax,
220 ppath_lraytrace,
221 rte_pos,
222 rte_los,
223 rte_pos2,
224 cloudbox_on,
225 ppath_inside_cloudbox_do,
226 f_grid,
227 ppath_agenda);
228}
229
230/* Workspace method: Doxygen documentation will be auto-generated */
232 Ppath& ppath,
233 const Agenda& ppath_agenda,
234 const Numeric& ppath_lmax,
235 const Numeric& ppath_lraytrace,
236 const Index& atmgeom_checked,
237 const Vector& f_grid,
238 const Index& cloudbox_on,
239 const Index& cloudbox_checked,
240 const Index& ppath_inside_cloudbox_do,
241 const Vector& rte_pos,
242 const Vector& rte_los,
243 const Vector& rte_pos2,
244 const Numeric& altitude,
245 const Numeric& accuracy,
246 const Verbosity& verbosity) {
247 ppathCalc(ws,
248 ppath,
249 ppath_agenda,
250 ppath_lmax,
251 ppath_lraytrace,
252 atmgeom_checked,
253 f_grid,
254 cloudbox_on,
255 cloudbox_checked,
256 ppath_inside_cloudbox_do,
257 rte_pos,
258 rte_los,
259 rte_pos2,
260 verbosity);
261
262 // Iterate until converging at altitude of interest
263 Index pos = first_pos_before_altitude(ppath, altitude);
264 Numeric lmax = ppath_lmax;
265 while (true) {
266 lmax *= 0.5;
267
268 ppathCalc(ws,
269 ppath,
270 ppath_agenda,
271 lmax,
272 ppath_lraytrace,
273 atmgeom_checked,
274 f_grid,
275 cloudbox_on,
276 cloudbox_checked,
277 ppath_inside_cloudbox_do,
278 ppath.dim == 1 ? ppath.pos(pos, 0) : ppath.pos(pos, joker),
279 ppath.dim == 1 ? ppath.los(pos, 0) : ppath.los(pos, joker),
280 rte_pos2,
281 verbosity);
282 pos = first_pos_before_altitude(ppath, altitude);
283
284 if (std::abs(ppath.pos(pos, 0) - altitude) < accuracy) {
285 ppathCalc(ws,
286 ppath,
287 ppath_agenda,
288 ppath_lmax,
289 ppath_lraytrace,
290 atmgeom_checked,
291 f_grid,
292 cloudbox_on,
293 cloudbox_checked,
294 ppath_inside_cloudbox_do,
295 ppath.dim == 1 ? ppath.pos(pos, 0) : ppath.pos(pos, joker),
296 ppath.dim == 1 ? ppath.los(pos, 0) : ppath.los(pos, joker),
297 rte_pos2,
298 verbosity);
299 break;
300 }
301 }
302}
303
304
305/* A help function for *ppathFixedLstep*. Determines the radius of the path,
306 raius of ellipsoid and surface altitiude for (lat,lon) at a distance l from
307 the sensor.
308 */
310 Numeric& r_e,
311 Numeric& z_surf,
312 const Index& atmosphere_dim,
313 const Vector& lat_grid,
314 const Vector& lon_grid,
315 const Vector& refellipsoid,
316 const Matrix& z_surface,
317 const Numeric& x0,
318 const Numeric& y0,
319 const Numeric& z0,
320 const Numeric& dx,
321 const Numeric& dy,
322 const Numeric& dz,
323 const Numeric& l,
324 const Numeric& lat0,
325 const Numeric& lon0,
326 const Numeric& za0,
327 const Numeric& aa0) {
328 if (atmosphere_dim == 1) { // 1D
329 const Numeric x = x0+l*dx, y = y0+l*dy, z = z0+l*dz;
330 r = sqrt( x*x + y*y + z*z);
331 r_e = refellipsoid[0];
332 z_surf = z_surface(0,0);
333
334 } else { // 2D and 3D
335
336 // Go to spherical coordinates
337 Numeric lat, lon;
338 cart2sph(r, lat, lon, x0+l*dx, y0+l*dy, z0+l*dz, lat0, lon0, za0, aa0);
339
340 // Latitude grid position
341 ARTS_USER_ERROR_IF (lat < lat_grid[0],
342 "Search of surface intersection ended up ", lat_grid[0]-lat,
343 " degrees below start of *lat_grid*. You need to expand the grid.")
344 ARTS_USER_ERROR_IF (lat > last(lat_grid),
345 "Search of surface intersection ended up ", lat-last(lat_grid),
346 " degrees above end of *lat_grid*. You need to expand the grid.")
347
348 GridPos gp_lat, gp_lon;
349 gridpos(gp_lat, lat_grid, lat);
350
351 // Interpolate to get ellipsoid radius and surface altitude
352 // If 3D we also need to calculate lon grid position
353 r_e = refell2d(refellipsoid, lat_grid, gp_lat);
354 if (atmosphere_dim==2) {
355 Vector itw(2);
356 interpweights(itw, gp_lat);
357 z_surf = interp(itw, z_surface(joker, 0), gp_lat);
358 } else {
359 const Numeric lonmax = last(lon_grid);
360 resolve_lon(lon, lon_grid[0], lonmax);
361 ARTS_USER_ERROR_IF (lon < lon_grid[0],
362 "Search of surface intersection ended up ", lon_grid[0]-lon,
363 " degrees below start of *lon_grid*. You need to expand the grid.")
364 ARTS_USER_ERROR_IF (lon > lonmax,
365 "Search of surface intersection ended up ", lon-lonmax,
366 " degrees above end of *lon_grid*. You need to expand the grid.")
367
368 gridpos(gp_lon, lon_grid, lon);
369 Vector itw(4);
370 interpweights(itw, gp_lat, gp_lon);
371 z_surf = interp(itw, z_surface, gp_lat, gp_lon);
372 }
373 }
374}
375
376
377/* Workspace method: Doxygen documentation will be auto-generated */
379 const Index& atmfields_checked,
380 const Index& atmgeom_checked,
381 const Index& atmosphere_dim,
382 const Vector& lat_grid,
383 const Vector& lon_grid,
384 const Tensor3& z_field,
385 const Vector& refellipsoid,
386 const Matrix& z_surface,
387 const Index& cloudbox_on,
388 const Vector& rte_pos,
389 const Vector& rte_los,
390 const Numeric& ppath_lmax,
391 const Index& za_scale,
392 const Numeric& z_coarse,
393 const Numeric& l_coarse,
394 const Verbosity&) {
395 // Basics checks of input
396 ARTS_USER_ERROR_IF (atmfields_checked != 1,
397 "The atmospheric fields must be flagged to have "
398 "passed a consistency check (atmfields_checked=1).");
399 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
400 "The atmospheric geometry must be flagged to have "
401 "passed a consistency check (atmgeom_checked=1).");
402 chk_rte_pos(atmosphere_dim, rte_pos);
403 chk_rte_los(atmosphere_dim, rte_los);
404 ARTS_USER_ERROR_IF (cloudbox_on,
405 "This method does not yet handle an active cloudbox.");
406 ARTS_USER_ERROR_IF (ppath_lmax <= 0,
407 "This method requires that *ppath_lmax* > 0.");
408 ARTS_USER_ERROR_IF (atmosphere_dim == 1 && refellipsoid[1] > 1e-7,
409 "For 1D, second element of *refellipsoid* must be 0.");
410
411 // Set lat etc according to atmosphere_dim
412 const Numeric lon0 = atmosphere_dim == 3 ? rte_pos[2] : 0;
413 const Numeric lat0 = atmosphere_dim >= 2 ? rte_pos[1] : 0;
414 const Numeric r0 = refell2r(refellipsoid, lat0) + rte_pos[0];
415 const Numeric za0 = abs(rte_los[0]);
416 Numeric aa0 = atmosphere_dim == 3 ? rte_los[1] : 0;
417 if (atmosphere_dim == 2 && rte_los[0]<0) { aa0 = 180; } // 2D special
418
419 // Obtain a top of the atmosphere altitude
420 // To make sure that it represents a point inside the atmosphere, we take the lowest
421 // altitude of top pressure level
422 const Numeric z_toa = min( z_field(z_field.npages()-1,joker,joker) );
423
424 // Set step lengths to use
425 const Numeric lstep = ppath_lmax * (za_scale ? 1/abs(cos(DEG2RAD*za0)) : 1);
426 const Numeric lcoarse = l_coarse * (za_scale ? 1/abs(cos(DEG2RAD*za0)) : 1);
427
428 // ECEF pos and los
429 Numeric x0, y0, z0, dx, dy, dz;
430 poslos2cart(x0, y0, z0, dx, dy, dz, r0, lat0, lon0, za0, aa0);
431
432 // Length to z_coarse
433 //
434 // We do this, here and below, by adding the search altitude to the major
435 // axis of the ellipsoid. This is approxinative. In theory, the eccentricity
436 // should be adopted. More important is that in ARTS the radius varies
437 // lineraly between points of the latitude grid, while
438 // line_refellipsoid_intersect operates with a fully analytical description
439 // of the ellipsoid and this causes some inconsistency. That is, the found
440 // length is not exact from ARTS's perspective.
441 //
442 Vector rell = refellipsoid;
443 Numeric l2coarse = -1;
444 if (z_coarse >= 0) {
445 rell[0] += z_coarse;
446 line_refellipsoid_intersect(l2coarse, rell, x0, y0, z0, dx, dy, dz);
447 }
448
449 // Create a vector with the distance from the sensor for each ppath point
450 //
451 Vector lvec(1); lvec[0] = 0;
452 Index background = 1; // Index of radiative background. 1 = space
453 //
454 // Upward
455 if (za0 < 90) {
456
457 if (rte_pos[0] >= z_toa) {
458 // If looking up we get an empty path
459 // lvec initiated to represent this
460
461 } else {
462 // We are inside looking up. Lengths go from 0 to distance to TOA
463 Numeric l2toa;
464 rell[0] = refellipsoid[0] + z_toa;
465 line_refellipsoid_intersect(l2toa, rell, x0, y0, z0, dx, dy, dz);
466
467 // Check that sensor actually is above the surface
468 Numeric r, r_e, z_surf;
469 surf_radius_at_l(r, r_e, z_surf, atmosphere_dim, lat_grid, lon_grid,
470 refellipsoid, z_surface, x0, y0, z0, dx, dy, dz,
471 0, lat0, lon0, za0, aa0);
472 ARTS_USER_ERROR_IF (r < r_e+z_surf,
473 "The sensor is ", r_e+z_surf - r, " m below the surface")
474
475 // Create vector with lengths, considering z_coarse
476 if (z_coarse < 0) {
477 linspace(lvec, 0, l2toa, lstep);
478 } else {
479 if (rte_pos[0] > z_coarse) {
480 linspace(lvec, 0, l2toa, lcoarse);
481 } else {
482 Vector l1, l2;
483 linspace(l1, 0, l2coarse, lstep);
484 linspace(l2, last(l1)+lstep, l2toa, lcoarse);
485 const Index n1=l1.nelem(), n2=l2.nelem();
486 lvec.resize(n1+n2);
487 lvec[Range(0,n1)] = l1;
488 lvec[Range(n1,n2)] = l2;
489 }
490 }
491 }
492 }
493
494 // Downward
495 else {
496 // Determine lowest surface altitude
497 // We can directly find length to the surface if its altitude is constant
498 // so we also look for this. Always true for 1D!
499 Numeric z_surf_min = atmosphere_dim == 1 ? z_surface(0,0) : 1e99;
500 bool surface_found = true;
501 if (atmosphere_dim == 2) {
502 for (Index i=1; i<lat_grid.nelem(); i++) {
503 if (z_surface(i,0) < z_surf_min)
504 z_surf_min = z_surface(i,0);
505 if (z_surface(i,0) != z_surface(0,0))
506 surface_found = false;
507 }
508 } else if (atmosphere_dim == 3) {
509 for (Index i=1; i<lat_grid.nelem(); i++) {
510 for (Index j=1; j<lon_grid.nelem(); j++) {
511 if (z_surface(i,j) < z_surf_min)
512 z_surf_min = z_surface(i,j);
513 if (z_surface(i,j) != z_surface(0,0))
514 surface_found = false;
515 }
516 }
517 }
518
519 // Determine length to z_surf_min (l2s)
520 Numeric l2s;
521 rell[0] = refellipsoid[0] + z_surf_min;
522 line_refellipsoid_intersect(l2s, rell, x0, y0, z0, dx, dy, dz);
523
524 // No intersection with surface
525 if (l2s < 0) {
526 if (!surface_found) {
528 "Cases of limb sounding type are only allowed "
529 "together with constant surface altitude.");
530 } else {
531 // Make function to find proper (geodetic) tangent point to finish this
532 ARTS_USER_ERROR ( "Limb sounding not yet handled.");
533 }
534
535 // Ensured intersection with surface
536 } else {
537
538 background = 2;
539
540 // If we need to search for the surface, we do this in two steps
541 // 1. Move upward until we have a point at or above surface
542 // 2. Use intersection to find length to surface
543 bool inphase1 = true;
544 bool underground = false; // If first l2s happen to be above surface
545 const Numeric lmoveup = 1e3;
546 const Numeric dr = 0.01; // When we stop
547 Numeric llong = l2s;
548 Numeric lshort = -1;
549 while (!surface_found) {
550 // Determine where we have the surface at l2s
551 Numeric r, r_e, z_surf;
552 surf_radius_at_l(r, r_e, z_surf, atmosphere_dim, lat_grid, lon_grid,
553 refellipsoid, z_surface, x0, y0, z0, dx, dy, dz,
554 l2s, lat0, lon0, za0, aa0);
555
556 // Compare radii and update iteration variables
557 const Numeric r_s = r_e + z_surf;
558 if (abs(r-r_s) <= dr) {
559 surface_found = true;
560 } else if (inphase1) {
561 if (r < r_s) {
562 underground = true;
563 l2s -= lmoveup;
564 } else {
565 if (underground) {
566 inphase1 = false;
567 lshort = l2s;
568 l2s = (lshort+llong)/2;
569 } else {
570 l2s += 0.1*lmoveup; // If we end up here the start l2s was too
571 llong = l2s; // short. Can happen due to inconsistencies,
572 } // as explained above.
573 }
574 } else {
575 if (r > r_s) {
576 lshort = l2s;
577 } else {
578 llong = l2s;
579 }
580 l2s = (lshort+llong)/2;
581 }
582 } // while
583
584 // Distance to TOA (same code as above for upward)
585 Numeric l2toa;
586 if (rte_pos[0] <= z_toa) {
587 l2toa = 0;
588 } else {
589 rell[0] = refellipsoid[0] + z_toa;
590 line_refellipsoid_intersect(l2toa, rell, x0, y0, z0, dx, dy, dz);
591 }
592
593 // Create vector with lengths, considering z_coarse
594 // lvec shall start exactly at l2s
595 if (l2s < lstep) {
596 linspace(lvec, 0, l2s, l2s);
597 } else if (z_coarse < 0) {
598 linspace(lvec, l2toa, l2s, lstep);
599 lvec += l2s - last(lvec);
600 } else {
601 if (rte_pos[0] < z_coarse) {
602 linspace(lvec, l2toa, l2s, lstep);
603 lvec += l2s - last(lvec);
604 } else {
605 Vector l1, l2;
606 linspace(l1, l2coarse, l2s, lstep);
607 l1 += l2s - last(l1);
608 linspace(l2, l2toa, l1[0]-lstep, lcoarse);
609 l2 += l1[0]-lstep - last(l2);
610 const Index n1=l1.nelem(), n2=l2.nelem();
611 lvec.resize(n1+n2);
612 lvec[Range(0,n2)] = l2;
613 lvec[Range(n2,n1)] = l1;
614 }
615 }
616
617
618 } // Surface intersection
619 } // Up/down
620
621 // Create ppath structure
622 ppath_init_structure(ppath, atmosphere_dim, lvec.nelem());
623 ppath.constant = geometrical_ppc(r0, abs(rte_los[0]));
624 ppath_set_background(ppath, background);
625 ppath.end_los[joker] = rte_los[joker]; // end_pos done below
626 ppath.end_lstep = lvec[0];
627 ppath.nreal = 1;
628 ppath.ngroup = 1;
629
630 // Empty ppath
631 if (ppath.np == 1) {
632 ppath.r = r0;
633 ppath.pos(0,0) = rte_pos[0];
634 if (atmosphere_dim == 1)
635 ppath.end_pos[1] = 0;
636 ppath.los(0,joker) = rte_los[joker];
637 }
638 // Otherwise loop lvec (split in atm. dim. to make code more efficient)
639 else {
640 // 1D
641 if (atmosphere_dim == 1) {
642 ppath.end_pos[0] = rte_pos[0];
643 ppath.end_pos[1] = 0;
644 for (Index i=0; i<ppath.np; i++) {
645 if (i > 0)
646 ppath.lstep[i-1] = lvec[i]-lvec[i-1];;
647 const Numeric l = lvec[i], x = x0+l*dx, y = y0+l*dy, z = z0+l*dz;
648 ppath.r[i] = sqrt( x*x + y*y + z*z);;
649 ppath.pos(i,0) = ppath.r[i] - refellipsoid[0];
650 ppath.los(i,0) = geompath_za_at_r(ppath.constant,
651 za0, // Will not work for limb sounding!
652 ppath.r[i]);
653 ppath.pos(i,1) = geompath_lat_at_za(za0, 0, ppath.los(i,0));
654 }
655 gridpos(ppath.gp_p, z_field(joker,0,0), ppath.pos(joker,0));
656
657 // 2D
658 } else if (atmosphere_dim == 2) {
659 ppath.end_pos[joker] = rte_pos[joker];
660 const Index nz = z_field.npages();
661 ArrayOfGridPos gp_z, gp_lat(1), gp_lon(0);
662 gridpos_1to1(gp_z, Vector(nz));
663 Tensor3 z_grid(nz,1,1); // The altitudes at one (lat,lon)
664 const Numeric lat1=lat_grid[0], lat2=last(lat_grid);
665 ARTS_ASSERT( abs(dy) < 1e-9 ); // 2D happens strictly inside plane y=0
666 for (Index i=0; i<ppath.np; i++) {
667 if (i > 0)
668 ppath.lstep[i-1] = lvec[i]-lvec[i-1];;
669 const Numeric l = lvec[i], x = x0+l*dx, z = z0+l*dz;
670 cart2poslos(ppath.r[i], ppath.pos(i,1), ppath.los(i,0),
671 x, z, dx, dz, ppath.constant, rte_pos[1], rte_los[0]);
672 ARTS_USER_ERROR_IF (ppath.pos(i,1) < lat1,
673 "The latitude grid must be extended downwards with at "
674 "least ", lat1-ppath.pos(i,1), " degrees to allow "
675 "the ppath to fully be inside of the model atmosphere.")
676 ARTS_USER_ERROR_IF (ppath.pos(i,1) > lat2,
677 "The latitude grid must be extended upwards with at "
678 "least ", ppath.pos(i,1)-lat2, " degrees to allow "
679 "the ppath to fully be inside of the model atmosphere.")
680 gridpos(ppath.gp_lat[i], lat_grid, ppath.pos(i,1));
681 //
682 gridpos_copy(gp_lat[0], ppath.gp_lat[i]);
683 regrid_atmfield_by_gp(z_grid, atmosphere_dim, z_field, gp_z, gp_lat, gp_lon);
684 const Numeric r_e = refell2d(refellipsoid, lat_grid, ppath.gp_lat[i]);
685 ppath.pos(i,0) = ppath.r[i] - r_e;
686 gridpos(ppath.gp_p[i], z_grid(joker,0,0), ppath.pos(i,0));
687 }
688
689 // 3D
690 } else {
691 ppath.end_pos[joker] = rte_pos[joker];
692 const Index nz = z_field.npages();
693 ArrayOfGridPos gp_z, gp_lat(1), gp_lon(1);
694 gridpos_1to1(gp_z, Vector(nz));
695 Tensor3 z_grid(nz,1,1); // The altitudes at one (lat,lon)
696 const Numeric lat1=lat_grid[0], lat2=last(lat_grid);
697 const Numeric lon1=lon_grid[0], lon2=last(lon_grid);
698 for (Index i=0; i<ppath.np; i++) {
699 if (i > 0)
700 ppath.lstep[i-1] = lvec[i]-lvec[i-1];;
701 const Numeric l = lvec[i], x = x0+l*dx, y = y0+l*dy, z = z0+l*dz;
702 cart2poslos(ppath.r[i], ppath.pos(i,1), ppath.pos(i,2),
703 ppath.los(i,0), ppath.los(i,1), x, y, z, dx, dy, dz,
704 ppath.constant, x0, y0, z0, rte_pos[1], rte_pos[2],
705 rte_los[0], rte_los[1]);
706 ARTS_USER_ERROR_IF (ppath.pos(i,1) < lat1,
707 "The latitude grid must be extended downwards with at "
708 "least ", lat1-ppath.pos(i,1), " degrees to allow "
709 "the ppath to fully be inside of the model atmosphere.")
710 ARTS_USER_ERROR_IF (ppath.pos(i,1) > lat2,
711 "The latitude grid must be extended upwards with at "
712 "least ", ppath.pos(i,1)-lat2, " degrees to allow "
713 "the ppath to fully be inside of the model atmosphere.")
714 ARTS_USER_ERROR_IF (ppath.pos(i,2) < lon1,
715 "The latitude grid must be extended downwards with at "
716 "least ", lon1-ppath.pos(i,2), " degrees to allow "
717 "the ppath to fully be inside of the model atmosphere.")
718 ARTS_USER_ERROR_IF (ppath.pos(i,2) > lon2,
719 "The latitude grid must be extended upwards with at "
720 "least ", ppath.pos(i,2)-lon2, " degrees to allow "
721 "the ppath to fully be inside of the model atmosphere.")
722 gridpos(ppath.gp_lat[i], lat_grid, ppath.pos(i,1));
723 gridpos(ppath.gp_lon[i], lon_grid, ppath.pos(i,2));
724 //
725 gridpos_copy(gp_lat[0], ppath.gp_lat[i]);
726 gridpos_copy(gp_lon[0], ppath.gp_lon[i]);
727 regrid_atmfield_by_gp(z_grid, atmosphere_dim, z_field, gp_z, gp_lat, gp_lon);
728 const Numeric r_e = refell2d(refellipsoid, lat_grid, ppath.gp_lat[i]);
729 ppath.pos(i,0) = ppath.r[i] - r_e;
730 gridpos(ppath.gp_p[i], z_grid(joker,0,0), ppath.pos(i,0));
731 }
732 }
733 }
734}
735
736/* Workspace method: Doxygen documentation will be auto-generated */
738 Ppath& ppath,
739 Vector& rte_los,
740 Numeric& ppath_lraytrace,
741 const Agenda& ppath_step_agenda,
742 const Index& atmosphere_dim,
743 const Vector& p_grid,
744 const Vector& lat_grid,
745 const Vector& lon_grid,
746 const Tensor3& z_field,
747 const Vector& f_grid,
748 const Vector& refellipsoid,
749 const Matrix& z_surface,
750 const Vector& rte_pos,
751 const Vector& rte_pos2,
752 const Numeric& ppath_lmax,
753 const Numeric& za_accuracy,
754 const Numeric& pplrt_factor,
755 const Numeric& pplrt_lowest,
756 const Verbosity& verbosity) {
757 //--- Check input -----------------------------------------------------------
758 ARTS_USER_ERROR_IF (atmosphere_dim == 2,
759 "2D atmospheres not yet handled. Support for negative"
760 " zenith angles needed. Remind me (Patrick) to fix this.");
761 //---------------------------------------------------------------------------
762
763 // Geometric LOS from rte_pos to rte_pos2
764 Vector rte_los_geom;
766 atmosphere_dim,
767 lat_grid,
768 lon_grid,
769 refellipsoid,
770 rte_pos,
771 rte_pos2,
772 verbosity);
773
774 // Radius of rte_pos and rte_pos2
775 const Numeric r1 =
776 pos2refell_r(atmosphere_dim, refellipsoid, lat_grid, lon_grid, rte_pos) +
777 rte_pos[0];
778 const Numeric r2 =
779 pos2refell_r(atmosphere_dim, refellipsoid, lat_grid, lon_grid, rte_pos2) +
780 rte_pos2[0];
781
782 // Geometric distance between rte_pos and rte_pos2, effective 2D-lat for
783 // rte_pos and and Cartesian coordinates of rte_pos:
784 Numeric l12, lat1 = 0, x1, y1 = 0, z1;
785 if (atmosphere_dim <= 2) {
786 if (atmosphere_dim == 2) {
787 lat1 = rte_pos[1];
788 }
789 distance2D(l12, r1, lat1, r2, rte_pos2[1]);
790 pol2cart(x1, z1, r1, lat1);
791 } else {
792 distance3D(l12, r1, rte_pos[1], rte_pos[2], r2, rte_pos2[1], rte_pos2[2]);
793 sph2cart(x1, y1, z1, r1, rte_pos[1], rte_pos[2]);
794 }
795
796 // Define remaining variables used in the while-loop below
797 //
798 // Basic bookkeeping variables
799 Numeric za_upp_limit = 180;
800 Numeric za_low_limit = 0;
801 //
802 // Various variables associated with the ppath, and the point of the path
803 // closest to the transmitter
804 Ppath ppt; // "Test ppath"
805 Index ip = -999; // Index of closest ppath point
806 Numeric xip, yip = 0, zip; // Cartesian coords. of the closest ppath point
807 Numeric dxip, dyip = 0, dzip; // Cartesian LOS of the closest ppath point
808 //
809 // Data for the intersection of the l12-sphere
810 Vector posc(max(Index(2), atmosphere_dim));
811 Numeric rc, xc, yc = 0, zc;
812
815
816 const Index maxiter = 99;
817 Vector t_za(maxiter, -999), t_dza(maxiter, -999);
818 Index it = -1;
819
820 // Keep trying until ready or ground intersetion determined
821 //
822 bool ground = false;
823 bool failed = false;
824 Index ntries = 0;
825 //
826 while (true) {
827 // Path for present rte_los (no cloudbox!)
828 ppath_calc(ws,
829 ppt,
830 ppath_step_agenda,
831 atmosphere_dim,
832 p_grid,
833 lat_grid,
834 lon_grid,
835 z_field,
836 f_grid,
837 refellipsoid,
838 z_surface,
839 0,
840 ArrayOfIndex(0),
841 rte_pos,
842 rte_los,
843 ppath_lmax,
844 ppath_lraytrace,
845 0,
846 verbosity);
847
848 // Find the point closest to rte_pos2, on the side towards rte_pos.
849 // We do this by looking at the distance to rte_pos, that should be
850 // as close to l12 as possible, but not exceed it.
851 Numeric lip = 99e99;
852 ip = ppt.np;
853 //
854 while (lip >= l12 && ip > 0) {
855 ip--;
856 if (atmosphere_dim <= 2) {
857 distance2D(lip, r1, lat1, ppt.r[ip], ppt.pos(ip, 1));
858 } else {
859 distance3D(lip,
860 r1,
861 rte_pos[1],
862 rte_pos[2],
863 ppt.r[ip],
864 ppt.pos(ip, 1),
865 ppt.pos(ip, 2));
866 }
867 }
868
869 Numeric za_new, daa = 0;
870
871 // Surface intersection:
872 // Not OK if the ground position is too far from rte_pos2.
873 // (30 km selected to allow misses of smaller size when rte_pos2 is at
874 // surface level, but surface interference never OK if rte_pos above TOA)
875 if (ppath_what_background(ppt) == 2 && ip == ppt.np - 1 &&
876 l12 - lip > 30e3) {
877 za_new = rte_los[0] - 1;
878 za_upp_limit = rte_los[0];
879 }
880
881 // Ppath OK
882 else {
883 // Estimate ppath at the distance of l12, and calculate size
884 // of "miss" (measured in diffference in geometric angles)
885 Vector los;
886 Numeric dza;
887 if (atmosphere_dim <= 2) {
888 // Convert pos and los for point ip to cartesian coordinates
890 xip, zip, dxip, dzip, ppt.r[ip], ppt.pos(ip, 1), ppt.los(ip, 0));
891 // Find where the extension from point ip crosses the l12
892 // sphere: point c
893 Numeric latc;
894 line_circle_intersect(xc, zc, xip, zip, dxip, dzip, x1, z1, l12);
895 cart2pol(rc, latc, xc, zc, ppt.pos(ip, 1), ppt.los(ip, 0));
896 posc[1] = latc;
897 posc[0] =
898 rc - pos2refell_r(
899 atmosphere_dim, refellipsoid, lat_grid, lon_grid, posc);
900 } else {
901 // Convert pos and los for point ip to cartesian coordinates
902 poslos2cart(xip,
903 yip,
904 zip,
905 dxip,
906 dyip,
907 dzip,
908 ppt.r[ip],
909 ppt.pos(ip, 1),
910 ppt.pos(ip, 2),
911 ppt.los(ip, 0),
912 ppt.los(ip, 1));
913 // Find where the extension from point ip crosses the l12
914 // sphere: point c
915 Numeric latc, lonc;
917 xc, yc, zc, xip, yip, zip, dxip, dyip, dzip, x1, y1, z1, l12);
918 cart2sph(rc,
919 latc,
920 lonc,
921 xc,
922 yc,
923 zc,
924 ppt.pos(ip, 1),
925 ppt.pos(ip, 2),
926 ppt.los(ip, 0),
927 ppt.los(ip, 1));
928 posc[1] = latc;
929 posc[2] = lonc;
930 posc[0] =
931 rc - pos2refell_r(
932 atmosphere_dim, refellipsoid, lat_grid, lon_grid, posc);
933 }
934 //
936 atmosphere_dim,
937 lat_grid,
938 lon_grid,
939 refellipsoid,
940 rte_pos,
941 posc,
942 verbosity);
943 //
944 dza = los[0] - rte_los_geom[0];
945
946 // Update bookkeeping variables
947 it++;
948 t_za[it] = rte_los[0];
949 t_dza[it] = dza;
950 //
951 if (dza > 0 && rte_los[0] < za_upp_limit) {
952 za_upp_limit = rte_los[0];
953 } else if (dza < 0 && rte_los[0] > za_low_limit) {
954 za_low_limit = rte_los[0];
955 }
956
957 // Ready ?
958 if (abs(dza) <= za_accuracy) {
959 break;
960 } else if (za_upp_limit - za_low_limit <= za_accuracy / 10) {
961 if (max(t_dza) < -10 * za_accuracy) {
962 ground = true;
963 out3 << " Ground intersection determined !!!\n";
964 break;
965 } else {
966 failed = true;
967 out3 << " Zenith angle search range closed !!!\n";
968 break;
969 }
970 }
971 // Catch non-convergence (just for extra safety, za-range should be
972 // closed quicker than this)
973 ntries += 1;
974 if (ntries >= maxiter) {
975 failed = true;
976 out3 << " Too many iterations !!!\n";
977 break;
978 }
979
980 // Estimate new angle
981 if (it < 1) {
982 za_new = rte_los[0] - dza;
983 } else {
984 // Estimate new angle by linear regression over some of the
985 // last calculations
986 const Index nfit = min(it + 1, (Index)3);
987 const Index i0 = it - nfit + 1;
988 Vector p;
989 linreg(p, t_za[Range(i0, nfit)], t_dza[Range(i0, nfit)]);
990 za_new = -p[0] / p[1];
991 }
992 //
993 if (atmosphere_dim == 3) {
994 daa = los[1] - rte_los_geom[1];
995 }
996 }
997
998 // Update rte_los. Use bisection of za_new is basically
999 // identical to old angle, or is outside lower or upper
1000 // limit. Otherwise use reult of linear reg.
1001 if (std::isinf(za_new) || std::isnan(za_new) ||
1002 abs(za_new - rte_los[0]) < 0.99 * za_accuracy ||
1003 za_new <= za_low_limit || za_new >= za_upp_limit) {
1004 rte_los[0] = (za_low_limit + za_upp_limit) / 2;
1005 } else {
1006 rte_los[0] = za_new;
1007 if (atmosphere_dim == 3) {
1008 rte_los[1] -= daa;
1009 if (rte_los[1] < -180) {
1010 rte_los[1] += 360;
1011 } else if (rte_los[1] > 180) {
1012 rte_los[1] -= 360;
1013 }
1014 }
1015 }
1016 } // while
1017 //--------------------------------------------------------------------------
1018
1019 // If failed re-try with a shorter ppath_lraytrace, if not ending up with
1020 // a too small value.
1021 if (failed) {
1022 ppath_lraytrace /= pplrt_factor;
1023
1024 if (ppath_lraytrace >= pplrt_lowest) {
1025 out2 << " Re-start with ppath_lraytrace = " << ppath_lraytrace;
1027 ppath,
1028 rte_los,
1029 ppath_lraytrace,
1030 ppath_step_agenda,
1031 atmosphere_dim,
1032 p_grid,
1033 lat_grid,
1034 lon_grid,
1035 z_field,
1036 f_grid,
1037 refellipsoid,
1038 z_surface,
1039 rte_pos,
1040 rte_pos2,
1041 ppath_lmax,
1042 za_accuracy,
1043 pplrt_factor,
1044 pplrt_lowest,
1045 verbosity);
1046 } else {
1047 ppath_init_structure(ppath, atmosphere_dim, 1);
1048 ppath_set_background(ppath, 0);
1049 }
1050 return; // --->
1051 }
1052
1053 // Create final ppath.
1054 // If ground intersection: Set to length 1 and ground background,
1055 // to flag non-OK path
1056 // Otherwise: Fill path and set background to transmitter
1057
1058 if (ground) {
1059 ppath_init_structure(ppath, atmosphere_dim, 1);
1060 ppath_set_background(ppath, 2);
1061 }
1062
1063 else {
1064 // Distance between point ip of ppt and posc
1065 Numeric ll;
1066 if (atmosphere_dim <= 2) {
1067 distance2D(ll, rc, posc[1], ppt.r[ip], ppt.pos(ip, 1));
1068 } else {
1069 distance3D(
1070 ll, rc, posc[1], posc[2], ppt.r[ip], ppt.pos(ip, 1), ppt.pos(ip, 2));
1071 }
1072
1073 // Last point of ppt closest to rte_pos2. No point to add, maybe
1074 // calculate start_lstep and start_los:
1075 if (ip == ppt.np - 1) {
1076 ppath_init_structure(ppath, atmosphere_dim, ppt.np);
1077 ppath_copy(ppath, ppt, -1);
1078 if (ppath_what_background(ppath) == 1) {
1079 ppath.start_lstep = ll;
1080 Numeric d1, d2 = 0, d3;
1081 if (atmosphere_dim <= 2) {
1082 cart2poslos(d1,
1083 d3,
1084 ppath.start_los[0],
1085 xc,
1086 zc,
1087 dxip,
1088 dzip,
1089 ppt.r[ip] * sin(DEG2RAD * ppt.los(ip, 0)),
1090 ppt.pos(ip, 1),
1091 ppt.los(ip, 0));
1092 } else {
1093 cart2poslos(d1,
1094 d2,
1095 d3,
1096 ppath.start_los[0],
1097 ppath.start_los[1],
1098 xc,
1099 yc,
1100 zc,
1101 dxip,
1102 dyip,
1103 dzip,
1104 ppt.r[ip] * sin(DEG2RAD * ppt.los(ip, 0)),
1105 xip,
1106 yip,
1107 zip, // Added 161027,PE
1108 ppt.pos(ip, 1),
1109 ppt.pos(ip, 2),
1110 ppt.los(ip, 0),
1111 ppt.los(ip, 1));
1112 }
1113 }
1114 }
1115 // rte_pos2 inside the atmosphere (posc entered as end point)
1116 else {
1117 ppath_init_structure(ppath, atmosphere_dim, ip + 2);
1118 ppath_copy(ppath, ppt, ip + 1);
1119 //
1120 const Index i = ip + 1;
1121 if (atmosphere_dim <= 2) {
1122 cart2poslos(ppath.r[i],
1123 ppath.pos(i, 1),
1124 ppath.los(i, 0),
1125 xc,
1126 zc,
1127 dxip,
1128 dzip,
1129 ppt.r[ip] * sin(DEG2RAD * ppt.los(ip, 0)),
1130 ppt.pos(ip, 1),
1131 ppt.los(ip, 0));
1132 } else {
1133 cart2poslos(ppath.r[i],
1134 ppath.pos(i, 1),
1135 ppath.pos(i, 2),
1136 ppath.los(i, 0),
1137 ppath.los(i, 1),
1138 xc,
1139 yc,
1140 zc,
1141 dxip,
1142 dyip,
1143 dzip,
1144 ppt.r[ip] * sin(DEG2RAD * ppt.los(ip, 0)),
1145 xip,
1146 yip,
1147 zip, // Added 161027,PE
1148 ppt.pos(ip, 1),
1149 ppt.pos(ip, 2),
1150 ppt.los(ip, 0),
1151 ppt.los(ip, 1));
1152 }
1153 //
1154 ppath.pos(i, joker) = posc;
1155 ppath.lstep[i - 1] = ll;
1156 ppath.start_los = ppath.los(i, joker);
1157
1158 // n by linear interpolation
1159 // Gets tripped when ll is very close to (slightly greater than) lstep (ISA)
1160 ARTS_ASSERT(ll < ppt.lstep[i - 1]);
1161 const Numeric w = ll / ppt.lstep[i - 1];
1162 ppath.nreal[i] = (1 - w) * ppt.nreal[i - 1] + w * ppt.nreal[i];
1163 ppath.ngroup[i] = (1 - w) * ppt.ngroup[i - 1] + w * ppt.ngroup[i];
1164
1165 // Grid positions
1166 GridPos gp_lat, gp_lon;
1167 rte_pos2gridpos(ppath.gp_p[i],
1168 gp_lat,
1169 gp_lon,
1170 atmosphere_dim,
1171 p_grid,
1172 lat_grid,
1173 lon_grid,
1174 z_field,
1175 ppath.pos(i, Range(0, atmosphere_dim)));
1176 if (atmosphere_dim >= 2) {
1177 gridpos_copy(ppath.gp_lat[i], gp_lat);
1178 if (atmosphere_dim == 3) {
1179 gridpos_copy(ppath.gp_lon[i], gp_lon);
1180 }
1181 }
1182 }
1183
1184 // Common stuff
1185 ppath_set_background(ppath, 9);
1186 ppath.start_pos = rte_pos2;
1187 }
1188}
1189
1190/* Workspace method: Doxygen documentation will be auto-generated */
1192 const Index& atmosphere_dim,
1193 const Tensor3& z_field,
1194 const Matrix& z_surface,
1195 const Index& cloudbox_on,
1196 const ArrayOfIndex& cloudbox_limits,
1197 const Index& ppath_inside_cloudbox_do,
1198 const Vector& rte_pos,
1199 const Vector& rte_los,
1200 const Numeric& ppath_lmax,
1201 const Verbosity&) {
1202 // This function is a WSM but it is normally only called from yCalc.
1203 // For that reason, this function does not repeat input checks that are
1204 // performed in yCalc, it only performs checks regarding the sensor
1205 // position and LOS.
1206
1207 const Numeric z_sensor = rte_pos[0];
1208 const Numeric za_sensor = rte_los[0];
1209 const Index nz = z_field.npages();
1210 const Numeric z_toa = z_field(nz - 1, 0, 0);
1211 const bool above_toa = z_sensor > z_toa ? true : false;
1212 const Numeric z_end = above_toa ? z_toa : z_sensor;
1213 const Numeric dz2dl = abs(1 / cos(DEG2RAD * za_sensor));
1214 Index background = -99;
1215
1216 // Basics checks of input
1217 ARTS_USER_ERROR_IF (atmosphere_dim != 1,
1218 "The function can only be used for 1D atmospheres.");
1219 chk_rte_pos(atmosphere_dim, rte_pos);
1220 chk_rte_los(atmosphere_dim, rte_los);
1221 ARTS_USER_ERROR_IF (ppath_inside_cloudbox_do && !cloudbox_on,
1222 "The WSV *ppath_inside_cloudbox_do* can only be set "
1223 "to 1 if also *cloudbox_on* is 1.");
1224 ARTS_USER_ERROR_IF (z_sensor < z_surface(0, 0),
1225 "The sensor is below the surface."
1226 " altitude of sensor : ", z_sensor, "\n"
1227 " altitude of surface : ", z_surface(0, 0))
1228 ARTS_USER_ERROR_IF (abs(za_sensor - 90) < 0.1,
1229 "The zenith angle is ", za_sensor, "\n"
1230 "The method does not allow this. The zenith angle must deviate\n"
1231 "from 90 deg with at least 0.1 deg. That is, to be outside [89.9,90.1].")
1232
1233 // Find end grid position
1234 GridPos gp_end;
1235 // To avoid compiler warnings, start to assuming above_toa
1236 gp_end.idx = nz - 2;
1237 gp_end.fd[0] = 1;
1238 gp_end.fd[1] = 0;
1239 if (!above_toa) {
1240 for (Index i = 0; i < nz - 1; i++) {
1241 if (z_sensor < z_field(i + 1, 0, 0)) {
1242 gp_end.idx = i;
1243 gp_end.fd[0] = (z_sensor - z_field(i, 0, 0)) /
1244 (z_field(i + 1, 0, 0) - z_field(i, 0, 0));
1245 gp_end.fd[1] = 1 - gp_end.fd[0];
1246 break;
1247 }
1248 }
1249 }
1250
1251 // Catch cases resulting in a ppath with 1 point
1252 bool path_to_follow = true;
1253 if (above_toa && za_sensor < 90) {
1254 // Path fully in space
1255 ppath_init_structure(ppath, atmosphere_dim, 1);
1256 background = 1;
1257 path_to_follow = false;
1258 } else if (z_sensor == z_surface(0, 0) && za_sensor > 90) {
1259 // On ground, looking down
1260 ppath_init_structure(ppath, atmosphere_dim, 1);
1261 background = 2;
1262 path_to_follow = false;
1263 } else if (cloudbox_on) {
1264 if (!ppath_inside_cloudbox_do &&
1265 z_sensor > z_field(cloudbox_limits[0], 0, 0) &&
1266 z_sensor < z_field(cloudbox_limits[1], 0, 0)) {
1267 // Inside cloud box
1268 ppath_init_structure(ppath, atmosphere_dim, 1);
1269 background = 4;
1270 path_to_follow = false;
1271 } else if ((z_sensor == z_field(cloudbox_limits[0], 0, 0) &&
1272 za_sensor > 90) ||
1273 (z_sensor == z_field(cloudbox_limits[1], 0, 0) &&
1274 za_sensor < 90)) {
1275 // Cloud box boundary
1276 ppath_init_structure(ppath, atmosphere_dim, 1);
1277 background = 3;
1278 path_to_follow = false;
1279 } else if (above_toa && cloudbox_limits[1] == nz - 1) {
1280 // Cloud box boundary is at TOA
1281 ppath_init_structure(ppath, atmosphere_dim, 1);
1282 background = 3;
1283 path_to_follow = false;
1284 }
1285 }
1286
1287 // Determine ppath
1288 if (path_to_follow) {
1289 const Numeric max_dz = ppath_lmax > 0 ? ppath_lmax / dz2dl : 9e99;
1290
1291 // Variables to describe each "break-point" of ppath. Point 0 is the end
1292 // point. Not all nz points are necessarily passed.
1293 ArrayOfIndex l_idx(nz);
1294 ArrayOfVector l_fd0(nz);
1295 ArrayOfVector l_z(nz);
1296 Index nptot = 0;
1297
1298 // Determine number of ppath points in each layer
1299 {
1300 Numeric z = z_end;
1301 Index iout = -1;
1302
1303 // Code similar, but for simplicity, we handle down- and
1304 // up-ward separately:
1305 if (za_sensor > 90) // Downward-looking
1306 {
1307 // Here we go down to next pressure level (or the surface) in each
1308 // step. That is, if above surface, last point of step has fd[0]=0.
1309
1310 // Put in end point
1311 iout++;
1312 nptot++;
1313 l_fd0[0].resize(1);
1314 l_z[0].resize(1);
1315 l_idx[0] = gp_end.idx;
1316 l_fd0[0][0] = gp_end.fd[0];
1317 l_z[0][0] = z_end;
1318
1319 for (Index i = gp_end.idx; i >= 0 && background < 0; i--) {
1320 // Surface inside layer?
1321 Numeric dz_step;
1322 if (z_field(i, 0, 0) > z_surface(0, 0)) {
1323 dz_step = z - z_field(i, 0, 0);
1324 } else {
1325 dz_step = z - z_surface(0, 0);
1326 background = 2;
1327 }
1328
1329 const Index np =
1330 dz_step <= max_dz ? 1 : Index(ceil(dz_step / max_dz));
1331 const Numeric dz = dz_step / Numeric(np);
1332 const Numeric dz_layer = z_field(i + 1, 0, 0) - z_field(i, 0, 0);
1333
1334 // Update counters and resize
1335 iout++;
1336 nptot += np;
1337 l_fd0[iout].resize(np);
1338 l_z[iout].resize(np);
1339
1340 // Intermediate points
1341 for (Index j = 0; j < np - 1; j++) {
1342 l_z[iout][j] = z - (Numeric(j) + 1) * dz;
1343 l_fd0[iout][j] = (l_z[iout][j] - z_field(i, 0, 0)) / dz_layer;
1344 }
1345
1346 // End points handled seperately to avoid numerical problems
1347 l_idx[iout] = i;
1348 if (background == 2) // Surface is reached
1349 {
1350 l_z[iout][np - 1] = z_surface(0, 0);
1351 l_fd0[iout][np - 1] =
1352 (l_z[iout][np - 1] - z_field(i, 0, 0)) / dz_layer;
1353 } else {
1354 l_z[iout][np - 1] = z_field(i, 0, 0);
1355 l_fd0[iout][np - 1] = 0;
1356 //
1357 if (cloudbox_on &&
1358 (i == cloudbox_limits[1] || i == cloudbox_limits[0])) {
1359 background = 3;
1360 }
1361 }
1362
1363 // Update z
1364 z = z_field(i, 0, 0);
1365 }
1366 } else // Upward-looking
1367 {
1368 // Here we have that first point of step has fd[0]=0, if not at
1369 // sensor
1370 for (Index i = gp_end.idx; i < nz && background < 0; i++) {
1371 Numeric dz_layer;
1372 Numeric dz_step;
1373 if (cloudbox_on && i != gp_end.idx &&
1374 (i == cloudbox_limits[0] ||
1375 i == cloudbox_limits[1])) { // At an active cloudbox boundary
1376 dz_step = 0;
1377 dz_layer = 1;
1378 background = 3;
1379 } else if (i == nz - 1) { // At TOA
1380 dz_step = 0;
1381 dz_layer = 1;
1382 background = 1;
1383 } else {
1384 dz_step = z_field(i + 1, 0, 0) - z;
1385 dz_layer = z_field(i + 1, 0, 0) - z_field(i, 0, 0);
1386 }
1387
1388 const Index np =
1389 dz_step <= max_dz ? 1 : Index(ceil(dz_step / max_dz));
1390 const Numeric dz = dz_step / Numeric(np);
1391
1392 // Update counters and resize
1393 iout++;
1394 nptot += np;
1395 l_fd0[iout].resize(np);
1396 l_z[iout].resize(np);
1397
1398 // Start points handled seperately to avoid numerical problems
1399 if (i == gp_end.idx) { // At sensor
1400 l_idx[iout] = i;
1401 l_z[iout][0] = z_sensor;
1402 l_fd0[iout][0] = gp_end.fd[0];
1403 } else if (i == nz - 1) { // At TOA
1404 l_idx[iout] = i - 1;
1405 l_z[iout][0] = z_field(i, 0, 0);
1406 l_fd0[iout][0] = 1;
1407 } else {
1408 l_idx[iout] = i;
1409 l_z[iout][0] = z_field(i, 0, 0);
1410 l_fd0[iout][0] = 0;
1411 }
1412
1413 // Intermediate points
1414 for (Index j = 1; j < np; j++) {
1415 l_z[iout][j] = z + Numeric(j) * dz;
1416 l_fd0[iout][j] = (l_z[iout][j] - z_field(i, 0, 0)) / dz_layer;
1417 }
1418
1419 // Update z
1420 if (background < 0) {
1421 z = z_field(i + 1, 0, 0);
1422 }
1423 }
1424 }
1425 }
1426
1427 ppath_init_structure(ppath, atmosphere_dim, nptot);
1428
1429 // Fill ppath.pos(joker,0), ppath.gp_p and ppath.lstep
1430 Index iout = -1;
1431 Numeric z_last = -999;
1432 for (Index i = 0; i < nz; i++) {
1433 for (Index j = 0; j < l_z[i].nelem(); j++) {
1434 iout++;
1435 ppath.pos(iout, 0) = l_z[i][j];
1436 ppath.gp_p[iout].idx = l_idx[i];
1437 ppath.gp_p[iout].fd[0] = l_fd0[i][j];
1438 ppath.gp_p[iout].fd[1] = 1 - l_fd0[i][j];
1439 if (iout == 0) {
1440 z_last = ppath.pos(iout, 0);
1441 } else {
1442 ppath.lstep[iout - 1] = dz2dl * abs(z_last - l_z[i][j]);
1443 z_last = l_z[i][j];
1444 }
1445 }
1446 }
1447 }
1448
1449 // Remaining data
1450 ppath_set_background(ppath, background);
1451 if (ppath.np == 1) {
1452 ppath.pos(0, 0) = z_end;
1453 ppath.gp_p[0] = gp_end;
1454 }
1455 ppath.pos(joker, 1) = 0;
1456 ppath.los(joker, 0) = za_sensor;
1457 ppath.constant = INFINITY; // Not defined here as r = Inf
1458 ppath.r = INFINITY;
1459 ppath.start_pos[0] = ppath.pos(ppath.np - 1, 0);
1460 ppath.start_pos[1] = 0;
1461 ppath.start_los[0] = za_sensor;
1462 ppath.end_pos[0] = z_sensor;
1463 ppath.end_pos[1] = 0;
1464 ppath.end_los[0] = za_sensor;
1465 if (above_toa) {
1466 ppath.end_lstep = dz2dl * (z_sensor - z_toa);
1467 }
1468 ppath.nreal = 1;
1469 ppath.ngroup = 1;
1470}
1471
1472/* Workspace method: Doxygen documentation will be auto-generated */
1474 Ppath& ppath,
1475 const Agenda& ppath_step_agenda,
1476 const Index& ppath_inside_cloudbox_do,
1477 const Index& atmosphere_dim,
1478 const Vector& p_grid,
1479 const Vector& lat_grid,
1480 const Vector& lon_grid,
1481 const Tensor3& z_field,
1482 const Vector& f_grid,
1483 const Vector& refellipsoid,
1484 const Matrix& z_surface,
1485 const Index& cloudbox_on,
1486 const ArrayOfIndex& cloudbox_limits,
1487 const Vector& rte_pos,
1488 const Vector& rte_los,
1489 const Numeric& ppath_lmax,
1490 const Numeric& ppath_lraytrace,
1491 const Verbosity& verbosity) {
1492 ppath_calc(ws,
1493 ppath,
1494 ppath_step_agenda,
1495 atmosphere_dim,
1496 p_grid,
1497 lat_grid,
1498 lon_grid,
1499 z_field,
1500 f_grid,
1501 refellipsoid,
1502 z_surface,
1503 cloudbox_on,
1504 cloudbox_limits,
1505 rte_pos,
1506 rte_los,
1507 ppath_lmax,
1508 ppath_lraytrace,
1509 ppath_inside_cloudbox_do,
1510 verbosity);
1511}
1512
1513/* Workspace method: Doxygen documentation will be auto-generated */
1514void ppathWriteXMLPartial( //WS Input:
1515 const String& file_format,
1516 const Ppath& ppath,
1517 // WS Generic Input:
1518 const String& f,
1519 const Index& file_index,
1520 const Verbosity& verbosity) {
1521 String filename = f;
1522 Ppath ppath_partial = ppath;
1523 ArrayOfGridPos empty_gp;
1524 //Vector empty_v;
1525
1526 ppath_partial.gp_p = empty_gp;
1527 ppath_partial.gp_lat = empty_gp;
1528 ppath_partial.gp_lon = empty_gp;
1529 //ppath_partial.nreal = empty_v;
1530 //ppath_partial.ngroup = empty_v;
1531
1532 if (file_index >= 0) {
1533 // Create default filename if empty
1534 filename_xml_with_index(filename, file_index, "ppath");
1535 }
1536
1537 WriteXML(file_format, ppath_partial, filename, 0, "ppath", "", "", verbosity);
1538}
1539
1540// FIXMEDOC@Richard TRy to describe the meaning of ppath_field
1541
1542/* Workspace method: Doxygen documentation will be auto-generated */
1544 ArrayOfPpath& ppath_field,
1545 const Agenda& ppath_agenda,
1546 const Numeric& ppath_lmax,
1547 const Numeric& ppath_lraytrace,
1548 const Index& atmgeom_checked,
1549 const Tensor3& z_field,
1550 const Vector& f_grid,
1551 const Index& cloudbox_on,
1552 const Index& cloudbox_checked,
1553 const Index& ppath_inside_cloudbox_do,
1554 const Vector& rte_pos,
1555 const Vector& rte_los,
1556 const Vector& rte_pos2,
1557 const Vector& refellipsoid,
1558 const Index& atmosphere_dim,
1559 const Index& zenith_angles_per_position,
1560 const Verbosity& verbosity) {
1561 ARTS_USER_ERROR_IF (atmosphere_dim not_eq 1,
1562 "Only for 1D atmospheres");
1563 ARTS_USER_ERROR_IF (refellipsoid[1] not_eq 0.0,
1564 "Not allowed for non-spherical planets");
1565 ARTS_USER_ERROR_IF (ppath_lmax >= 0,
1566 "Only allowed for long paths (ppath_lmax < 0)");
1567
1568 // Positions and angles of interest
1569 const Numeric zmin = z_field(0, 0, 0);
1570 const Numeric zmax = z_field(z_field.npages() - 1, 0, 0);
1571 const Numeric r = refellipsoid[0];
1572 const Numeric above_surface_tangent =
1573 90 - RAD2DEG * std::acos((r) / (r + zmax)) + 1e-4;
1574 const Numeric below_surface_tangent =
1575 90 - RAD2DEG * std::acos((r) / (r + zmax)) - 1e-4;
1576 const Numeric top_tangent = 90 - 1e-4;
1577
1578 ppath_field.resize(3 * zenith_angles_per_position);
1579 Index ppath_field_pos = 0;
1580
1581 Vector zenith_angles(zenith_angles_per_position);
1582
1583 // Upwards:
1584 nlinspace(zenith_angles, 0, 90, zenith_angles_per_position);
1585 Vector rte_pos_true = rte_pos;
1586 rte_pos_true[0] = zmin;
1587 Vector rte_los_true = rte_los;
1588 for (Index iz = 0; iz < zenith_angles_per_position; iz++) {
1589 rte_los_true[0] = zenith_angles[iz];
1590
1591 ppathCalc(ws,
1592 ppath_field[ppath_field_pos],
1593 ppath_agenda,
1594 ppath_lmax,
1595 ppath_lraytrace,
1596 atmgeom_checked,
1597 f_grid,
1598 cloudbox_on,
1599 cloudbox_checked,
1600 ppath_inside_cloudbox_do,
1601 rte_pos_true,
1602 rte_los_true,
1603 rte_pos2,
1604 verbosity);
1605
1606 ppath_field_pos++;
1607 }
1608
1609 // Limb:
1610 nlinspace(zenith_angles,
1611 above_surface_tangent,
1612 top_tangent,
1613 zenith_angles_per_position);
1614 rte_pos_true[0] = zmax;
1615 for (Index iz = 0; iz < zenith_angles_per_position; iz++) {
1616 rte_los_true[0] = 180 - zenith_angles[iz];
1617
1618 ppathCalc(ws,
1619 ppath_field[ppath_field_pos],
1620 ppath_agenda,
1621 ppath_lmax,
1622 ppath_lraytrace,
1623 atmgeom_checked,
1624 f_grid,
1625 cloudbox_on,
1626 cloudbox_checked,
1627 ppath_inside_cloudbox_do,
1628 rte_pos_true,
1629 rte_los_true,
1630 rte_pos2,
1631 verbosity);
1632
1633 ppath_field_pos++;
1634 }
1635
1636 // Downwards:
1637 nlinspace(
1638 zenith_angles, 0, below_surface_tangent, zenith_angles_per_position);
1639 for (Index iz = 0; iz < zenith_angles_per_position; iz++) {
1640 rte_los_true[0] = 180 - zenith_angles[iz];
1641
1642 ppathCalc(ws,
1643 ppath_field[ppath_field_pos],
1644 ppath_agenda,
1645 ppath_lmax,
1646 ppath_lraytrace,
1647 atmgeom_checked,
1648 f_grid,
1649 cloudbox_on,
1650 cloudbox_checked,
1651 ppath_inside_cloudbox_do,
1652 rte_pos_true,
1653 rte_los_true,
1654 rte_pos2,
1655 verbosity);
1656
1657 ppath_field_pos++;
1658 }
1659}
1660
1661/* Workspace method: Doxygen documentation will be auto-generated */
1663 ArrayOfPpath& ppath_field,
1664 const Agenda& ppath_agenda,
1665 const Numeric& ppath_lmax,
1666 const Numeric& ppath_lraytrace,
1667 const Index& atmgeom_checked,
1668 const Vector& f_grid,
1669 const Index& cloudbox_on,
1670 const Index& cloudbox_checked,
1671 const Index& ppath_inside_cloudbox_do,
1672 const Matrix& sensor_pos,
1673 const Matrix& sensor_los,
1674 const Vector& rte_pos2,
1675 const Verbosity& verbosity) {
1676 auto n = sensor_pos.nrows();
1677 ppath_field.resize(n);
1678
1679 ARTS_USER_ERROR_IF (sensor_los.nrows() not_eq n,
1680 "Your sensor position matrix and sensor line of sight matrix do not match in size.\n");
1681
1682 for (auto i = 0; i < n; i++)
1683 ppathCalc(ws,
1684 ppath_field[i],
1685 ppath_agenda,
1686 ppath_lmax,
1687 ppath_lraytrace,
1688 atmgeom_checked,
1689 f_grid,
1690 cloudbox_on,
1691 cloudbox_checked,
1692 ppath_inside_cloudbox_do,
1693 sensor_pos(i, joker),
1694 sensor_los(i, joker),
1695 rte_pos2,
1696 verbosity);
1697}
1698
1699/* Workspace method: Doxygen documentation will be auto-generated */
1700void ppath_stepGeometric( // WS Output:
1701 Ppath& ppath_step,
1702 // WS Input:
1703 const Index& atmosphere_dim,
1704 const Vector& lat_grid,
1705 const Vector& lon_grid,
1706 const Tensor3& z_field,
1707 const Vector& refellipsoid,
1708 const Matrix& z_surface,
1709 const Numeric& ppath_lmax,
1710 const Verbosity&) {
1711 // Input checks here would be rather costly as this function is called
1712 // many times. So we perform asserts in the sub-functions, but no checks
1713 // here.
1714
1715 // A call with background set, just wants to obtain the refractive index for
1716 // complete ppaths consistent of a single point.
1717 if (!ppath_what_background(ppath_step)) {
1718 if (atmosphere_dim == 1) {
1719 ppath_step_geom_1d(ppath_step,
1720 z_field(joker, 0, 0),
1721 refellipsoid,
1722 z_surface(0, 0),
1723 ppath_lmax);
1724 }
1725
1726 else if (atmosphere_dim == 2) {
1727 ppath_step_geom_2d(ppath_step,
1728 lat_grid,
1729 z_field(joker, joker, 0),
1730 refellipsoid,
1731 z_surface(joker, 0),
1732 ppath_lmax);
1733 }
1734
1735 else if (atmosphere_dim == 3) {
1736 ppath_step_geom_3d(ppath_step,
1737 lat_grid,
1738 lon_grid,
1739 z_field,
1740 refellipsoid,
1741 z_surface,
1742 ppath_lmax);
1743 }
1744
1745 else {
1746 ARTS_USER_ERROR ( "The atmospheric dimensionality must be 1-3.");
1747 }
1748 }
1749
1750 else {
1751 ARTS_ASSERT(ppath_step.np == 1);
1752 ppath_step.nreal[0] = 1;
1753 ppath_step.ngroup[0] = 1;
1754 }
1755}
1756
1757/* Workspace method: Doxygen documentation will be auto-generated */
1759 Ppath& ppath_step,
1760 const Agenda& refr_index_air_agenda,
1761 const Index& atmosphere_dim,
1762 const Vector& p_grid,
1763 const Vector& lat_grid,
1764 const Vector& lon_grid,
1765 const Tensor3& z_field,
1766 const Tensor3& t_field,
1767 const Tensor4& vmr_field,
1768 const Vector& refellipsoid,
1769 const Matrix& z_surface,
1770 const Vector& f_grid,
1771 const Numeric& ppath_lmax,
1772 const Numeric& ppath_lraytrace,
1773 const Verbosity&) {
1774 // Input checks here would be rather costly as this function is called
1775 // many times.
1776 ARTS_ASSERT(ppath_lraytrace > 0);
1777
1778 // A call with background set, just wants to obtain the refractive index for
1779 // complete ppaths consistent of a single point.
1780 if (!ppath_what_background(ppath_step)) {
1781 if (atmosphere_dim == 1) {
1783 ppath_step,
1784 p_grid,
1785 z_field,
1786 t_field,
1787 vmr_field,
1788 f_grid,
1789 refellipsoid,
1790 z_surface(0, 0),
1791 ppath_lmax,
1792 refr_index_air_agenda,
1793 "linear_basic",
1794 ppath_lraytrace);
1795 } else if (atmosphere_dim == 2) {
1797 ppath_step,
1798 p_grid,
1799 lat_grid,
1800 z_field,
1801 t_field,
1802 vmr_field,
1803 f_grid,
1804 refellipsoid,
1805 z_surface(joker, 0),
1806 ppath_lmax,
1807 refr_index_air_agenda,
1808 "linear_basic",
1809 ppath_lraytrace);
1810 } else if (atmosphere_dim == 3) {
1812 ppath_step,
1813 p_grid,
1814 lat_grid,
1815 lon_grid,
1816 z_field,
1817 t_field,
1818 vmr_field,
1819 f_grid,
1820 refellipsoid,
1821 z_surface,
1822 ppath_lmax,
1823 refr_index_air_agenda,
1824 "linear_basic",
1825 ppath_lraytrace);
1826 } else {
1827 ARTS_USER_ERROR ( "The atmospheric dimensionality must be 1-3.");
1828 }
1829 }
1830
1831 else {
1832 ARTS_ASSERT(ppath_step.np == 1);
1833 if (atmosphere_dim == 1) {
1835 ppath_step.nreal[0],
1836 ppath_step.ngroup[0],
1837 refr_index_air_agenda,
1838 p_grid,
1839 refellipsoid,
1840 z_field,
1841 t_field,
1842 vmr_field,
1843 f_grid,
1844 ppath_step.r[0]);
1845 } else if (atmosphere_dim == 2) {
1847 ppath_step.nreal[0],
1848 ppath_step.ngroup[0],
1849 refr_index_air_agenda,
1850 p_grid,
1851 lat_grid,
1852 refellipsoid,
1853 z_field,
1854 t_field,
1855 vmr_field,
1856 f_grid,
1857 ppath_step.r[0],
1858 ppath_step.pos(0, 1));
1859 } else {
1861 ppath_step.nreal[0],
1862 ppath_step.ngroup[0],
1863 refr_index_air_agenda,
1864 p_grid,
1865 lat_grid,
1866 lon_grid,
1867 refellipsoid,
1868 z_field,
1869 t_field,
1870 vmr_field,
1871 f_grid,
1872 ppath_step.r[0],
1873 ppath_step.pos(0, 1),
1874 ppath_step.pos(0, 2));
1875 }
1876 }
1877}
1878
1879/* Workspace method: Doxygen documentation will be auto-generated */
1880void rte_losSet(Vector& rte_los,
1881 const Index& atmosphere_dim,
1882 const Numeric& za,
1883 const Numeric& aa,
1884 const Verbosity&) {
1885 // Check input
1886 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1887
1888 if (atmosphere_dim == 1) {
1889 rte_los.resize(1);
1890 } else {
1891 rte_los.resize(2);
1892 rte_los[1] = aa;
1893 }
1894 rte_los[0] = za;
1895}
1896
1897/* Workspace method: Doxygen documentation will be auto-generated */
1899 const Index& atmosphere_dim,
1900 const Vector& lat_grid,
1901 const Vector& lon_grid,
1902 const Vector& refellipsoid,
1903 const Vector& rte_pos,
1904 const Vector& rte_pos2,
1905 const Verbosity&) {
1906 // Check input
1907 chk_rte_pos(atmosphere_dim, rte_pos);
1908 chk_rte_pos(atmosphere_dim, rte_pos2, true);
1909
1910 // Radius of rte_pos and rte_pos2
1911 const Numeric r1 =
1912 pos2refell_r(atmosphere_dim, refellipsoid, lat_grid, lon_grid, rte_pos) +
1913 rte_pos[0];
1914 const Numeric r2 =
1915 pos2refell_r(atmosphere_dim, refellipsoid, lat_grid, lon_grid, rte_pos2) +
1916 rte_pos2[0];
1917
1918 // Remaining polar and cartesian coordinates of rte_pos
1919 Numeric lat1, lon1 = 0, x1, y1 = 0, z1;
1920 // Cartesian coordinates of rte_pos2
1921 Numeric x2, y2 = 0, z2;
1922 //
1923 if (atmosphere_dim == 1) {
1924 // Latitude distance implicitly checked by chk_rte_pos
1925 lat1 = 0;
1926 pol2cart(x1, z1, r1, lat1);
1927 pol2cart(x2, z2, r2, rte_pos2[1]);
1928 } else if (atmosphere_dim == 2) {
1929 lat1 = rte_pos[1];
1930 pol2cart(x1, z1, r1, lat1);
1931 pol2cart(x2, z2, r2, rte_pos2[1]);
1932 } else {
1933 lat1 = rte_pos[1];
1934 lon1 = rte_pos[2];
1935 sph2cart(x1, y1, z1, r1, lat1, lon1);
1936 sph2cart(x2, y2, z2, r2, rte_pos2[1], rte_pos2[2]);
1937 }
1938
1939 // Geometrical LOS to transmitter
1940 Numeric za, aa;
1941 //
1942 los2xyz(za, aa, r1, lat1, lon1, x1, y1, z1, x2, y2, z2);
1943 //
1944 if (atmosphere_dim == 3) {
1945 rte_los.resize(2);
1946 rte_los[0] = za;
1947 rte_los[1] = aa;
1948 } else {
1949 rte_los.resize(1);
1950 rte_los[0] = za;
1951 if (atmosphere_dim == 2 && aa < 0) // Should 2D-za be negative?
1952 {
1953 rte_los[0] = -za;
1954 }
1955 }
1956}
1957
1958/* Workspace method: Doxygen documentation will be auto-generated */
1959void rte_posSet(Vector& rte_pos,
1960 const Index& atmosphere_dim,
1961 const Numeric& z,
1962 const Numeric& lat,
1963 const Numeric& lon,
1964 const Verbosity&) {
1965 // Check input
1966 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1967
1968 rte_pos.resize(atmosphere_dim);
1969 rte_pos[0] = z;
1970 if (atmosphere_dim >= 2) {
1971 rte_pos[1] = lat;
1972 }
1973 if (atmosphere_dim == 3) {
1974 rte_pos[2] = lon;
1975 }
1976}
1977
1978/* Workspace method: Doxygen documentation will be auto-generated */
1980 Vector& rte_los,
1981 const Index& atmosphere_dim,
1982 const Ppath& ppath,
1983 const Verbosity&) {
1984 const Index np = ppath.np;
1985
1986 // Check input
1987 chk_if_in_range("atmosphere_dim", atmosphere_dim, 1, 3);
1988 ARTS_USER_ERROR_IF (np == 0, "The input *ppath* is empty.");
1989 ARTS_USER_ERROR_IF (ppath.pos.nrows() != np,
1990 "Internal inconsistency in *ppath* (size of data "
1991 "does not match np).");
1992
1993 rte_pos = ppath.pos(np - 1, Range(0, atmosphere_dim));
1994 if (atmosphere_dim < 3) {
1995 rte_los = ppath.los(np - 1, Range(0, 1));
1996 } else {
1997 rte_los = ppath.los(np - 1, Range(0, 2));
1998 }
1999}
2000
2001/* Workspace method: Doxygen documentation will be auto-generated */
2003 Matrix& sensor_los,
2004 const Index& atmosphere_dim,
2005 const Vector& lat_grid,
2006 const Vector& lon_grid,
2007 const Vector& refellipsoid,
2008 const Matrix& sensor_pos,
2009 const Matrix& target_pos,
2010 const Verbosity& verbosity) {
2011 const Index n = sensor_pos.nrows();
2012
2013 ARTS_USER_ERROR_IF (sensor_pos.ncols() != atmosphere_dim,
2014 "The number of columns of sensor_pos must be "
2015 "equal to the atmospheric dimensionality.");
2016 ARTS_USER_ERROR_IF ((atmosphere_dim == 1 && target_pos.ncols() != 2) ||
2017 (atmosphere_dim >= 2 && target_pos.ncols() != atmosphere_dim),
2018 "The number of columns of targe_pos must be equal to "
2019 "the atmospheric dimensionality, except for 1D where "
2020 "two columns are demended (as for *rte_pos2*).");
2021 ARTS_USER_ERROR_IF (target_pos.nrows() != n,
2022 "*sensor_pos* and *target_pos* must have the same "
2023 "number of rows.");
2024
2025 atmosphere_dim < 3 ? sensor_los.resize(n, 1) : sensor_los.resize(n, 2);
2026 Vector rte_los;
2027 for (Index i = 0; i < n; i++) {
2029 atmosphere_dim,
2030 lat_grid,
2031 lon_grid,
2032 refellipsoid,
2033 sensor_pos(i, joker),
2034 target_pos(i, joker),
2035 verbosity);
2036 sensor_los(i, joker) = rte_los;
2037 }
2038}
2039
2040/* Workspace method: Doxygen documentation will be auto-generated */
2042 const Ppath& ppath,
2043 const Verbosity&) {
2044 Index it;
2045 find_tanpoint(it, ppath);
2046
2047 tan_pos.resize(ppath.pos.ncols());
2048
2049 if (it < 0) {
2050 tan_pos = std::numeric_limits<Numeric>::quiet_NaN();
2051 } else {
2052 tan_pos[0] = ppath.pos(it, 0);
2053 tan_pos[1] = ppath.pos(it, 1);
2054 if (ppath.pos.ncols() == 3) {
2055 tan_pos[2] = ppath.pos(it, 2);
2056 }
2057 }
2058}
2059
2060/* Workspace method: Doxygen documentation will be auto-generated */
2061void TangentPointPrint(const Ppath& ppath,
2062 const Index& level,
2063 const Verbosity& verbosity) {
2064 Index it;
2065 find_tanpoint(it, ppath);
2066
2067 ostringstream os;
2068
2069 if (it < 0) {
2070 os << "Lowest altitude found at the end of the propagation path.\n"
2071 << "This indicates that the tangent point is either above the\n"
2072 << "top-of-the-atmosphere or below the planet's surface.";
2073 } else {
2074 os << "Tangent point position:\n-----------------------\n"
2075 << " z = " << ppath.pos(it, 0) / 1e3 << " km\n"
2076 << " lat = " << ppath.pos(it, 1) << " deg";
2077 if (ppath.pos.ncols() == 3)
2078 os << "\n lon: " << ppath.pos(it, 2) << " deg";
2079 }
2080
2082 SWITCH_OUTPUT(level, os.str());
2083}
2084
2085/* Workspace method: Doxygen documentation will be auto-generated */
2087 Vector& za_vector,
2088 const Agenda& refr_index_air_agenda,
2089 const Matrix& sensor_pos,
2090 const Vector& p_grid,
2091 const Tensor3& t_field,
2092 const Tensor3& z_field,
2093 const Tensor4& vmr_field,
2094 const Vector& refellipsoid,
2095 const Index& atmosphere_dim,
2096 const Vector& f_grid,
2097 const Vector& ztan_vector,
2098 const Verbosity&) {
2099 ARTS_USER_ERROR_IF (atmosphere_dim != 1,
2100 "The function can only be used for 1D atmospheres.");
2101
2102 ARTS_USER_ERROR_IF (ztan_vector.nelem() != sensor_pos.nrows(),
2103 "The number of altitudes in true tangent altitude vector must\n"
2104 "match the number of positions in *sensor_pos*.")
2105
2106 // Set za_vector's size equal to ztan_vector
2107 za_vector.resize(ztan_vector.nelem());
2108
2109 // Define refraction variables
2110 Numeric refr_index_air, refr_index_air_group;
2111
2112 // Calculate refractive index for the tangential altitudes
2113 for (Index i = 0; i < ztan_vector.nelem(); i++) {
2114 ARTS_USER_ERROR_IF (ztan_vector[i] > sensor_pos(i, 0),
2115 "Invalid observation geometry: sensor (at z=", sensor_pos(i, 0),
2116 "m) is located below the requested tangent altitude (tanh=",
2117 ztan_vector[i], "m)")
2118
2120 refr_index_air,
2121 refr_index_air_group,
2122 refr_index_air_agenda,
2123 p_grid,
2124 refellipsoid[0],
2125 z_field,
2126 t_field,
2127 vmr_field,
2128 f_grid,
2129 ztan_vector[i] + refellipsoid[0]);
2130
2131 // Calculate zenith angle
2132 za_vector[i] = 180 - RAD2DEG * asin(refr_index_air *
2133 (refellipsoid[0] + ztan_vector[i]) /
2134 (refellipsoid[0] + sensor_pos(i, 0)));
2135 }
2136}
2137
2138/* Workspace method: Doxygen documentation will be auto-generated */
2139void VectorZtanToZa1D(Vector& za_vector,
2140 const Matrix& sensor_pos,
2141 const Vector& refellipsoid,
2142 const Index& atmosphere_dim,
2143 const Vector& ztan_vector,
2144 const Verbosity&) {
2145 ARTS_USER_ERROR_IF (atmosphere_dim != 1,
2146 "The function can only be used for 1D atmospheres.");
2147
2148 const Index npos = sensor_pos.nrows();
2149
2150 ARTS_USER_ERROR_IF (ztan_vector.nelem() != npos,
2151 "The number of altitudes in the geometric tangent altitude vector\n"
2152 "must match the number of positions in *sensor_pos*.")
2153
2154 za_vector.resize(npos);
2155
2156 for (Index i = 0; i < npos; i++) {
2157 ARTS_USER_ERROR_IF (ztan_vector[i] > sensor_pos(i, 0),
2158 "Invalid observation geometry: sensor (at z=", sensor_pos(i, 0),
2159 "m) is located below the requested tangent altitude (tanh=",
2160 ztan_vector[i], "m)")
2161
2162 za_vector[i] = geompath_za_at_r(refellipsoid[0] + ztan_vector[i],
2163 100,
2164 refellipsoid[0] + sensor_pos(i, 0));
2165 }
2166}
2167
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:39
The global header file for ARTS.
void ppath_agendaExecute(Workspace &ws, Ppath &ppath, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index cloudbox_on, const Index ppath_inside_cloudbox_do, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:24134
void WriteXML(const String &output_file_format, const T &in, const String &filename, const Index &no_clobber, const String &in_wsvname, const String &filename_wsvname, const String &no_clobber_wsvname, const Verbosity &verbosity)
WORKSPACE METHOD: WriteXML.
Definition: m_xml.h:118
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
void chk_rte_los(const Index &atmosphere_dim, ConstVectorView rte_los)
chk_rte_los
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
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The Matrix class.
Definition: matpackI.h:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
The range class.
Definition: matpackI.h:165
The Tensor3 class.
Definition: matpackIII.h:339
The Tensor4 class.
Definition: matpackIV.h:421
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
Workspace class.
Definition: workspace_ng.h:40
#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 distance3D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &r2, const Numeric &lat2, const Numeric &lon2)
distance3D
Definition: geodetic.cc:669
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
Definition: geodetic.cc:182
void sph2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &r, const Numeric &lat, const Numeric &lon)
sph2cart
Definition: geodetic.cc:1366
void los2xyz(Numeric &za, Numeric &aa, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &x1, const Numeric &y1, const Numeric &z1, const Numeric &x2, const Numeric &y2, const Numeric &z2)
los2xyz
Definition: geodetic.cc:1059
void cart2sph(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &lat0, const Numeric &lon0, const Numeric &za0, const Numeric &aa0)
cart2sph
Definition: geodetic.cc:602
void pol2cart(Numeric &x, Numeric &z, const Numeric &r, const Numeric &lat)
pol2cart
Definition: geodetic.cc:331
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1265
Numeric pos2refell_r(const Index &atmosphere_dim, ConstVectorView refellipsoid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView rte_pos)
pos2refell_r
Definition: geodetic.cc:1211
void line_circle_intersect(Numeric &x, Numeric &z, const Numeric &xl, const Numeric &zl, const Numeric &dx, const Numeric &dz, const Numeric &xc, const Numeric &zc, const Numeric &r)
geomtanpoint2d
Definition: geodetic.cc:278
void line_refellipsoid_intersect(Numeric &l, const Vector &refellipsoid, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
geomtanpoint
Definition: geodetic.cc:874
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:355
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1304
void cart2poslos(Numeric &r, Numeric &lat, Numeric &za, const Numeric &x, const Numeric &z, const Numeric &dx, const Numeric &dz, const Numeric &ppc, const Numeric &lat0, const Numeric &za0)
cart2poslos
Definition: geodetic.cc:113
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition: geodetic.cc:74
void line_sphere_intersect(Numeric &x, Numeric &y, Numeric &z, const Numeric &xl, const Numeric &yl, const Numeric &zl, const Numeric &dx, const Numeric &dy, const Numeric &dz, const Numeric &xc, const Numeric &yc, const Numeric &zc, const Numeric &r)
line_sphere_intersect
Definition: geodetic.cc:953
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
void gridpos_1to1(ArrayOfGridPos &gp, ConstVectorView grid)
gridpos_1to1
#define abs(x)
#define min(a, b)
#define dx
#define ll
#define max(a, b)
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:2297
Linear algebra functions.
std::vector< T > linspace(T s, T e, typename std::vector< T >::size_type count) noexcept
Template functions for general supergeneric ws methods.
#define SWITCH_OUTPUT(x, y)
Definition: m_general.h:46
void DiffZaAa(Matrix &dlos, const Vector &ref_los, const Matrix &other_los, const Verbosity &)
WORKSPACE METHOD: DiffZaAa.
Definition: m_ppath.cc:80
void rte_losGeometricFromRtePosToRtePos2(Vector &rte_los, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Vector &rte_pos, const Vector &rte_pos2, const Verbosity &)
WORKSPACE METHOD: rte_losGeometricFromRtePosToRtePos2.
Definition: m_ppath.cc:1898
void ppath_fieldFromDownUpLimbGeoms(Workspace &ws, ArrayOfPpath &ppath_field, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Tensor3 &z_field, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Vector &refellipsoid, const Index &atmosphere_dim, const Index &zenith_angles_per_position, const Verbosity &verbosity)
WORKSPACE METHOD: ppath_fieldFromDownUpLimbGeoms.
Definition: m_ppath.cc:1543
void ppath_fieldCalc(Workspace &ws, ArrayOfPpath &ppath_field, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Matrix &sensor_pos, const Matrix &sensor_los, const Vector &rte_pos2, const Verbosity &verbosity)
WORKSPACE METHOD: ppath_fieldCalc.
Definition: m_ppath.cc:1662
void ppathStepByStep(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &ppath_inside_cloudbox_do, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Verbosity &verbosity)
WORKSPACE METHOD: ppathStepByStep.
Definition: m_ppath.cc:1473
void surf_radius_at_l(Numeric &r, Numeric &r_e, Numeric &z_surf, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Matrix &z_surface, const Numeric &x0, const Numeric &y0, const Numeric &z0, const Numeric &dx, const Numeric &dy, const Numeric &dz, const Numeric &l, const Numeric &lat0, const Numeric &lon0, const Numeric &za0, const Numeric &aa0)
Definition: m_ppath.cc:309
void rte_losSet(Vector &rte_los, const Index &atmosphere_dim, const Numeric &za, const Numeric &aa, const Verbosity &)
WORKSPACE METHOD: rte_losSet.
Definition: m_ppath.cc:1880
void sensor_losGeometricFromSensorPosToOtherPositions(Matrix &sensor_los, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &refellipsoid, const Matrix &sensor_pos, const Matrix &target_pos, const Verbosity &verbosity)
WORKSPACE METHOD: sensor_losGeometricFromSensorPosToOtherPositions.
Definition: m_ppath.cc:2002
void AddZaAa(Matrix &new_los, const Vector &ref_los, const Matrix &dlos, const Verbosity &)
WORKSPACE METHOD: AddZaAa.
Definition: m_ppath.cc:56
void geo_posLowestAltitudeOfPpath(Vector &geo_pos, const Ppath &ppath, const Verbosity &verbosity)
WORKSPACE METHOD: geo_posLowestAltitudeOfPpath.
Definition: m_ppath.cc:120
void ppathCalc(Workspace &ws, Ppath &ppath, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Verbosity &)
WORKSPACE METHOD: ppathCalc.
Definition: m_ppath.cc:194
void rte_pos_losMoveToStartOfPpath(Vector &rte_pos, Vector &rte_los, const Index &atmosphere_dim, const Ppath &ppath, const Verbosity &)
WORKSPACE METHOD: rte_pos_losMoveToStartOfPpath.
Definition: m_ppath.cc:1979
void ppathWriteXMLPartial(const String &file_format, const Ppath &ppath, const String &f, const Index &file_index, const Verbosity &verbosity)
WORKSPACE METHOD: ppathWriteXMLPartial.
Definition: m_ppath.cc:1514
void ppathPlaneParallel(Ppath &ppath, const Index &atmosphere_dim, const Tensor3 &z_field, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Verbosity &)
WORKSPACE METHOD: ppathPlaneParallel.
Definition: m_ppath.cc:1191
void geo_posWherePpathPassesZref(Vector &geo_pos, const Ppath &ppath, const Numeric &z_ref, const Verbosity &verbosity)
WORKSPACE METHOD: geo_posWherePpathPassesZref.
Definition: m_ppath.cc:146
void TangentPointPrint(const Ppath &ppath, const Index &level, const Verbosity &verbosity)
WORKSPACE METHOD: TangentPointPrint.
Definition: m_ppath.cc:2061
void VectorZtanToZa1D(Vector &za_vector, const Matrix &sensor_pos, const Vector &refellipsoid, const Index &atmosphere_dim, const Vector &ztan_vector, const Verbosity &)
WORKSPACE METHOD: VectorZtanToZa1D.
Definition: m_ppath.cc:2139
const Numeric RAD2DEG
void VectorZtanToZaRefr1D(Workspace &ws, Vector &za_vector, const Agenda &refr_index_air_agenda, const Matrix &sensor_pos, const Vector &p_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &refellipsoid, const Index &atmosphere_dim, const Vector &f_grid, const Vector &ztan_vector, const Verbosity &)
WORKSPACE METHOD: VectorZtanToZaRefr1D.
Definition: m_ppath.cc:2086
void TangentPointExtract(Vector &tan_pos, const Ppath &ppath, const Verbosity &)
WORKSPACE METHOD: TangentPointExtract.
Definition: m_ppath.cc:2041
void rte_posSet(Vector &rte_pos, const Index &atmosphere_dim, const Numeric &z, const Numeric &lat, const Numeric &lon, const Verbosity &)
WORKSPACE METHOD: rte_posSet.
Definition: m_ppath.cc:1959
void ppathFromRtePos2(Workspace &ws, Ppath &ppath, Vector &rte_los, Numeric &ppath_lraytrace, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Vector &rte_pos, const Vector &rte_pos2, const Numeric &ppath_lmax, const Numeric &za_accuracy, const Numeric &pplrt_factor, const Numeric &pplrt_lowest, const Verbosity &verbosity)
WORKSPACE METHOD: ppathFromRtePos2.
Definition: m_ppath.cc:737
void geo_posEndOfPpath(Vector &geo_pos, const Ppath &ppath, const Verbosity &verbosity)
WORKSPACE METHOD: geo_posEndOfPpath.
Definition: m_ppath.cc:104
void ppath_stepGeometric(Ppath &ppath_step, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Numeric &ppath_lmax, const Verbosity &)
WORKSPACE METHOD: ppath_stepGeometric.
Definition: m_ppath.cc:1700
void ppath_stepRefractionBasic(Workspace &ws, Ppath &ppath_step, const Agenda &refr_index_air_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Vector &refellipsoid, const Matrix &z_surface, const Vector &f_grid, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Verbosity &)
WORKSPACE METHOD: ppath_stepRefractionBasic.
Definition: m_ppath.cc:1758
const Numeric DEG2RAD
void ppathFixedLstep(Ppath &ppath, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Index &za_scale, const Numeric &z_coarse, const Numeric &l_coarse, const Verbosity &)
WORKSPACE METHOD: ppathFixedLstep.
Definition: m_ppath.cc:378
void ppathCalcFromAltitude(Workspace &ws, Ppath &ppath, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &atmgeom_checked, const Vector &f_grid, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &ppath_inside_cloudbox_do, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Numeric &altitude, const Numeric &accuracy, const Verbosity &verbosity)
WORKSPACE METHOD: ppathCalcFromAltitude.
Definition: m_ppath.cc:231
Workspace methods and template functions for supergeneric XML IO.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:225
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:159
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
Declarations having to do with the four output streams.
#define CREATE_OUT3
Definition: messages.h:207
#define CREATE_OUT2
Definition: messages.h:206
#define CREATE_OUTS
Definition: messages.h:209
constexpr Numeric r0
The reference radius in IGRF13.
Definition: igrf13.cc:203
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
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:39
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
void ppath_step_refr_2d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstVectorView lat_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstVectorView z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
Calculates 2D propagation path steps, with refraction, using a simple and fast ray tracing scheme.
Definition: ppath.cc:3901
void ppath_step_refr_3d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
Calculates 3D propagation path steps, with refraction, using a simple and fast ray tracing scheme.
Definition: ppath.cc:4316
void diff_za_aa(Numeric &dza, Numeric &daa, const Numeric &za0, const Numeric &aa0, const Numeric &za, const Numeric &aa)
Takes the difference of zenith and azimuth angles.
Definition: ppath.cc:429
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1476
Index first_pos_before_altitude(const Ppath &p, const Numeric &alt)
Determines ppath position just below an altitude.
Definition: ppath.cc:522
void ppath_step_geom_3d(Ppath &ppath, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Numeric &lmax)
Calculates 3D geometrical propagation path steps.
Definition: ppath.cc:3251
Numeric geometrical_ppc(const Numeric &r, const Numeric &za)
Calculates the propagation path constant for pure geometrical calculations.
Definition: ppath.cc:96
void add_za_aa(Numeric &za, Numeric &aa, const Numeric &za0, const Numeric &aa0, const Numeric &dza, const Numeric &daa)
Adds up zenith and azimuth angles.
Definition: ppath.cc:391
void ppath_copy(Ppath &ppath1, const Ppath &ppath2, const Index &ncopy)
Copy the content in ppath2 to ppath1.
Definition: ppath.cc:1496
void ppath_step_geom_2d(Ppath &ppath, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, ConstVectorView z_surface, const Numeric &lmax)
Calculates 2D geometrical propagation path steps.
Definition: ppath.cc:2717
Numeric geompath_za_at_r(const Numeric &ppc, const Numeric &a_za, const Numeric &r)
Calculates the zenith angle for a given radius along a geometrical propagation path.
Definition: ppath.cc:103
void find_tanpoint(Index &it, const Ppath &ppath)
Identifies the tangent point of a propagation path.
Definition: ppath.cc:510
void ppath_set_background(Ppath &ppath, const Index &case_nr)
Sets the background field of a Ppath structure.
Definition: ppath.cc:1450
void resolve_lon(Numeric &lon, const Numeric &lon5, const Numeric &lon6)
Resolves which longitude angle that shall be used.
Definition: ppath.cc:500
Numeric geompath_lat_at_za(const Numeric &za0, const Numeric &lat0, const Numeric &za)
Calculates the latitude for a given zenith angle along a geometrical propagation path.
Definition: ppath.cc:148
void ppath_calc(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const bool &ppath_inside_cloudbox_do, const Verbosity &verbosity)
This is the core for the WSM ppathStepByStep.
Definition: ppath.cc:5175
void ppath_step_geom_1d(Ppath &ppath, ConstVectorView z_field, ConstVectorView refellipsoid, const Numeric &z_surface, const Numeric &lmax)
Calculates 1D geometrical propagation path steps.
Definition: ppath.cc:2353
void ppath_init_structure(Ppath &ppath, const Index &atmosphere_dim, const Index &np)
Initiates a Ppath structure to hold the given number of points.
Definition: ppath.cc:1409
void ppath_step_refr_1d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, const Numeric &z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
Calculates 1D propagation path steps including effects of refraction.
Definition: ppath.cc:3557
Propagation path structure and functions.
#define w
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
void get_refr_index_3d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat, const Numeric &lon)
Definition: refraction.cc:357
void get_refr_index_1d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r)
get_refr_index_1d
Definition: refraction.cc:184
void get_refr_index_2d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat)
get_refr_index_2d
Definition: refraction.cc:264
Refraction functions.
void regrid_atmfield_by_gp(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regrids an atmospheric field, for precalculated grid positions.
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
Numeric fd[2]
Definition: interpolation.h:75
Index idx
Definition: interpolation.h:74
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
Matrix los
Line-of-sight at each ppath point.
Definition: ppath.h:66
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Definition: ppath.h:86
Index np
Number of points describing the ppath.
Definition: ppath.h:52
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath.h:64
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Definition: ppath.h:84
Numeric end_lstep
The distance between end pos and the first position in pos.
Definition: ppath.h:76
Vector start_pos
Start position.
Definition: ppath.h:58
Vector lstep
The length between ppath points.
Definition: ppath.h:70
Numeric start_lstep
Length between sensor and atmospheric boundary.
Definition: ppath.h:62
Numeric constant
The propagation path constant (only used for 1D)
Definition: ppath.h:54
Vector r
Radius of each ppath point.
Definition: ppath.h:68
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Definition: ppath.h:82
Vector ngroup
The group index of refraction.
Definition: ppath.h:80
Index dim
Atmospheric dimensionality.
Definition: ppath.h:50
Vector end_pos
End position.
Definition: ppath.h:72
Vector start_los
Start line-of-sight.
Definition: ppath.h:60
Vector nreal
The real part of the refractive index at each path position.
Definition: ppath.h:78
Vector end_los
End line-of-sight.
Definition: ppath.h:74
void filename_xml_with_index(String &filename, const Index &file_index, const String &varname, const Index &digits)
Gives the default filename, with file index, for the XML formats.
Definition: xml_io.cc:234
This file contains basic functions to handle XML data files.