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