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