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