84 return r * sin(
DEG2RAD * abs(za));
95 Numeric za =
RAD2DEG * asin(ppc / r);
129 return ppc / sin(
DEG2RAD * abs(za));
137 ARTS_ASSERT((za0 >= 0 && za >= 0) || (za0 < 0 && za < 0));
139 return lat0 + za0 - za;
147 return sqrt(r * r - ppc * ppc);
156 return sqrt(l * l + ppc * ppc);
174 const Numeric& lat) {
177 ARTS_ASSERT((za0 >= 0 && lat >= lat0) || (za0 <= 0 && lat <= lat0));
180 const Numeric za = za0 + lat0 - lat;
216 const bool& tanpoint,
217 const Numeric& lmax) {
237 n =
max(Index(1), Index(ceil(abs(l2 - l1) / lmax)));
243 lstep = (l2 - l1) / (Numeric)n;
254 for (Index i = 1; i < n; i++) {
255 const Numeric l = l1 + lstep * (Numeric)i;
271 for (Index i = 1; i <= n; i++) {
288 const Numeric r = sqrt(dx * dx + dy * dy + dz * dz);
301 const Numeric zarad =
DEG2RAD * za;
302 const Numeric aarad =
DEG2RAD * aa;
306 dy = sin(aarad) * dx;
307 dx = cos(aarad) * dx;
316 aa =
RAD2DEG * atan2( de, dn );
324 const Numeric st = sin(
DEG2RAD*za);
352 sqrt(vrot[0] * vrot[0] + vrot[1] * vrot[1] + vrot[2] * vrot[2]);
354 const Numeric
u = vrot[0] / l;
355 const Numeric
v = vrot[1] / l;
356 const Numeric
w = vrot[2] / l;
357 const Numeric u2 =
u *
u;
358 const Numeric v2 =
v *
v;
359 const Numeric w2 =
w *
w;
364 R(0, 0) = u2 + (v2 + w2) *
c;
365 R(0, 1) =
u *
v * (1 -
c) -
w * s;
366 R(0, 2) =
u *
w * (1 -
c) +
v * s;
367 R(1, 0) =
u *
v * (1 -
c) +
w * s;
368 R(1, 1) = v2 + (u2 + w2) *
c;
369 R(1, 2) =
v *
w * (1 -
c) -
u * s;
370 R(2, 0) =
u *
w * (1 -
c) -
v * s;
371 R(2, 1) =
v *
w * (1 -
c) +
u * s;
372 R(2, 2) = w2 + (u2 + v2) *
c;
380 const Numeric& daa) {
387 zaaa2cart(xyz[0], xyz[1], xyz[2], 90, aa0);
400 zaaa2cart(xyz[0], xyz[1], xyz[2], 90 + dza, aa0 + daa);
427 zaaa2cart(xyz[0], xyz[1], xyz[2], za0, aa0);
440 zaaa2cart(xyz[0], xyz[1], xyz[2], za, aa);
450 Numeric za_tmp, aa_tmp;
477 const Numeric& refr_index_air) {
481 return r * refr_index_air * sin(
DEG2RAD * abs(za));
484void resolve_lon(Numeric& lon,
const Numeric& lon5,
const Numeric& lon6) {
487 if (lon < lon5 && lon + 180 <= lon6) {
489 }
else if (lon > lon6 && lon - 180 >= lon5) {
495 Numeric zmin = 99e99;
497 while (it < ppath.
np - 1 && ppath.
pos(it + 1, 0) < zmin) {
499 zmin = ppath.
pos(it, 0);
501 if (it == 0 || it == ppath.
np - 1) {
508 bool below =
false, above =
false;
511 for (Index i = 0; i < p.
np; i++) {
512 if (p.
pos(i, 0) < alt) {
513 if (above)
return i - 1;
516 if (below)
return i - 1;
526 Numeric signfac =
sign(ppath.
pos(1, 0) - ppath.
pos(0, 0));
527 for (Index i = 2; i < ppath.
np; i++) {
529 "A propagation path of limb character found. Such "
530 "viewing geometries are not supported by this method. "
531 "Propagation paths must result in strictly "
532 "increasing or decreasing altitudes.");
559 const Numeric& lat) {
560 return r1 + (lat - lat1) * (r3 - r1) / (lat3 - lat1);
564 ConstVectorView lat_grid,
565 ConstVectorView refellipsoid,
566 ConstVectorView z_surf,
570 const Numeric r1 =
refell2r(refellipsoid, lat_grid[i1]) + z_surf[i1];
571 const Numeric r2 =
refell2r(refellipsoid, lat_grid[i1 + 1]) + z_surf[i1 + 1];
573 c1 = (r2 - r1) / (lat_grid[i1 + 1] - lat_grid[i1]);
597 c1 = (r2 - r1) / (lat2 - lat1);
610 if (za > (90 - tilt) || za < (-90 - tilt)) {
640 const Numeric& r_hit,
641 const Numeric& r_start,
642 const Numeric& lat_start,
643 const Numeric& za_start,
644 const Numeric& ppc) {
648 const Numeric absza = abs(za_start);
652 if ((r_start >= r_hit && absza <= 90) || ppc > r_hit) {
659 if (absza > 90 && r_start <= r_hit) {
713 const Numeric zaabs = abs(za);
726 const Numeric beta =
DEG2RAD * (180 - zaabs);
727 const Numeric cv = cos(beta);
728 const Numeric sv = sin(beta);
731 const Numeric r0s = r0 * sv;
732 const Numeric r0c = r0 * cv;
733 const Numeric cs = c1 * sv;
734 const Numeric cc = c1 * cv;
741 p0[0] = r0s - rp * sv;
743 p0[2] = -r0s / 2 + cc;
744 p0[3] = -r0c / 6 - cs / 2;
745 p0[4] = r0s / 24 - cc / 6;
746 p0[5] = r0c / 120 + cs / 24;
747 p0[6] = -r0s / 720 + cc / 120;
754 if (abs(90 - zaabs) > 89.9) {
756 }
else if (abs(90 - zaabs) > 75) {
762 int solutionfailure = 1;
764 while (solutionfailure) {
767 p = p0[Range(0, n + 1)];
769 if (solutionfailure) {
778 if (abs(r0 - rp) < 1e-9)
784 Numeric dlat = 1.571;
786 for (Index i = 0; i < n; i++) {
787 if (roots(i, 1) == 0 && roots(i, 0) > dmin && roots(i, 0) < dlat) {
841 const Numeric& r_start0,
842 const Numeric& lat_start,
843 const Numeric& za_start,
850 const Numeric absza = abs(za_start);
853 ARTS_ASSERT(lat_start >= lat1 && lat_start <= lat3);
864 l =
max(1e-9, r - r_start0);
869 else if (absza > 180 -
ANGTOL) {
873 l =
max(1e-9, r_start0 - r);
883 const Numeric rmin =
min(r1, r3);
884 const Numeric rmax =
max(r1, r3);
887 if (rmax - rmin < 1e-6) {
889 Numeric r_start = r_start0;
892 if (r_start < rmax) {
896 if (r_start > rmin) {
904 if (lat > lat3 || lat < lat1) {
913 Numeric r_start = r_start0;
915 if (r_start < rmin) {
919 if (r_start > rmax) {
927 if (r_start > rmax) {
930 }
else if (r_start < rmin) {
941 if (lat < lat1 || lat > lat3) {
950 const Numeric rpl = r1 + cpl * (lat - lat1);
966 za = lat_start + za_start - lat;
975 if (lat < lat1 || lat > lat3) {
982 r = rpl + cpl * dlat;
985 za = lat_start + za_start - lat;
988 if (absza > 90 && abs(za) < 90) {
1023 const Numeric& lat3,
1024 const Numeric& lon5,
1025 const Numeric& lon6,
1031 const Numeric& lon) {
1036 return r15 + (lon - lon5) * (r16 - r15) / (lon6 - lon5);
1037 }
else if (lat == lat3) {
1038 return r35 + (lon - lon5) * (r36 - r35) / (lon6 - lon5);
1039 }
else if (lon == lon5) {
1040 return r15 + (lat - lat1) * (r35 - r15) / (lat3 - lat1);
1041 }
else if (lon == lon6) {
1042 return r16 + (lat - lat1) * (r36 - r16) / (lat3 - lat1);
1044 const Numeric fdlat = (lat - lat1) / (lat3 - lat1);
1045 const Numeric fdlon = (lon - lon5) / (lon6 - lon5);
1046 return (1 - fdlat) * (1 - fdlon) * r15 + fdlat * (1 - fdlon) * r35 +
1047 (1 - fdlat) * fdlon * r16 + fdlat * fdlon * r36;
1053 const Numeric& lat1,
1054 const Numeric& lat3,
1055 const Numeric& lon5,
1056 const Numeric& lon6,
1063 const Numeric& aa) {
1065 if (r15 == r35 && r15 == r36 && r15 == r16 && r35 == r36 && r35 == r16 &&
1074 rsurf_at_latlon(lat1, lat3, lon5, lon6, r15, r35, r36, r16, lat, lon);
1080 const Numeric dang = 1e-4;
1087 lat1, lat3, lon5, lon6, r15, r35, r36, r16, lat2, lon2) -
1094 lat1, lat3, lon5, lon6, r15, r35, r36, r16, lat2, lon2) -
1098 c1 = 0.5 * (4 * dr1 - dr2);
1099 c2 = (dr1 - c1) / (dang * dang);
1131 ConstVectorView lat_grid,
1132 ConstVectorView lon_grid,
1133 ConstVectorView refellipsoid,
1134 ConstMatrixView z_surf,
1137 const Numeric& aa) {
1142 const Index llat = lat_grid.nelem() - 1;
1145 if (lat_grid[llat] >
POLELAT) {
1149 "The upper latitude end of the atmosphere "
1150 "reached, that is not allowed.");
1154 if (ilon >= lon_grid.nelem() - 1) {
1155 if (is_lon_cyclic(lon_grid)) {
1159 "The upper longitude end of the atmosphere "
1160 "reached, that is not allowed.");
1168 lat =
interp(itw, lat_grid, gp_lat);
1170 lon =
interp(itw, lon_grid, gp_lon);
1173 const Numeric lat1 = lat_grid[ilat];
1174 const Numeric lat3 = lat_grid[ilat + 1];
1175 const Numeric lon5 = lon_grid[ilon];
1176 const Numeric lon6 = lon_grid[ilon + 1];
1177 const Numeric re1 =
refell2r(refellipsoid, lat1);
1178 const Numeric re3 =
refell2r(refellipsoid, lat3);
1179 const Numeric r15 = re1 + z_surf(ilat, ilon);
1180 const Numeric r35 = re3 + z_surf(ilat + 1, ilon);
1181 const Numeric r36 = re3 + z_surf(ilat + 1, ilon + 1);
1182 const Numeric r16 = re1 + z_surf(ilat, ilon + 1);
1185 c1, c2, lat1, lat3, lon5, lon6, r15, r35, r36, r16, lat, lon, aa);
1216 const Numeric beta =
DEG2RAD * (180 - za);
1217 const Numeric cv = cos(beta);
1218 const Numeric sv = sin(beta);
1221 const Numeric r0s = r0 * sv;
1222 const Numeric r0c = r0 * cv;
1223 const Numeric c1s = c1 * sv;
1224 const Numeric c1c = c1 * cv;
1225 const Numeric c2s = c2 * sv;
1226 const Numeric c2c = c2 * cv;
1233 p0[0] = r0s - rp * sv;
1235 p0[2] = -r0s / 2 + c1c + c2s;
1236 p0[3] = -r0c / 6 - c1s / 2 + c2c;
1237 p0[4] = r0s / 24 - c1c / 6 - c2s / 2;
1238 p0[5] = r0c / 120 + c1s / 24 - c2c / 6;
1239 p0[6] = -r0s / 720 + c1c / 120 + c2s / 24;
1247 if (abs(90 - za) > 89.9) {
1249 }
else if (abs(90 - za) > 75) {
1255 int solutionfailure = 1;
1257 while (solutionfailure) {
1260 p = p0[Range(0, n + 1)];
1262 if (solutionfailure) {
1277 Numeric dlat = 1.571;
1279 for (Index i = 0; i < n; i++) {
1280 if (roots(i, 1) == 0 && roots(i, 0) > dmin && roots(i, 0) < dlat) {
1330 const Numeric& r_hit,
1331 const Numeric& r_start,
1332 const Numeric& lat_start,
1333 const Numeric& lon_start,
1334 const Numeric& za_start,
1341 const Numeric& dz) {
1347 if ((r_start >= r_hit && za_start <= 90) || ppc > r_hit) {
1355 if (za_start < ANGTOL || za_start > 180 -
ANGTOL) {
1356 l = abs(r_hit - r_start);
1362 const Numeric p = x * dx + y * dy + z * dz;
1363 const Numeric pp = p * p;
1364 const Numeric q = x * x + y * y + z * z - r_hit * r_hit;
1365 const Numeric sq = sqrt(pp - q);
1366 const Numeric l1 = -p + sq;
1367 const Numeric l2 = -p - sq;
1369 const Numeric lmin =
min(l1, l2);
1370 const Numeric lmax =
max(l1, l2);
1383 lat =
RAD2DEG * asin((z + dz * l) / r_hit);
1384 lon =
RAD2DEG * atan2(y + dy * l, x + dx * l);
1394 const Index& atmosphere_dim,
1400 ppath.
dim = atmosphere_dim;
1404 const Index npos =
max(Index(2), atmosphere_dim);
1405 const Index nlos =
max(Index(1), atmosphere_dim - 1);
1416 ppath.
pos.resize(np, npos);
1417 ppath.
los.resize(np, nlos);
1419 ppath.
lstep.resize(np - 1);
1421 ppath.
gp_p.resize(np);
1422 if (atmosphere_dim >= 2) {
1424 if (atmosphere_dim == 3) {
1430 ppath.
nreal.resize(np);
1456 "Case number ", case_nr,
" is not defined.")
1467 }
else if (ppath.
background ==
"cloud box level") {
1469 }
else if (ppath.
background ==
"cloud box interior") {
1471 }
else if (ppath.
background ==
"transmitter") {
1476 " is not a valid background case.")
1502 if (n == ppath1.
np) {
1508 ppath1.
pos(Range(0, n), joker) = ppath2.
pos(Range(0, n), joker);
1509 ppath1.
los(Range(0, n), joker) = ppath2.
los(Range(0, n), joker);
1510 ppath1.
r[Range(0, n)] = ppath2.
r[Range(0, n)];
1511 ppath1.
nreal[Range(0, n)] = ppath2.
nreal[Range(0, n)];
1512 ppath1.
ngroup[Range(0, n)] = ppath2.
ngroup[Range(0, n)];
1514 ppath1.
lstep[Range(0, n - 1)] = ppath2.
lstep[Range(0, n - 1)];
1517 for (Index i = 0; i < n; i++) {
1520 if (ppath1.
dim >= 2) {
1524 if (ppath1.
dim == 3) {
1547 const Index n1 = ppath1.
np;
1548 const Index n2 = ppath2.
np;
1559 for (Index i = 1; i < n2; i++) {
1562 ppath1.
pos(i1, 0) = ppath2.
pos(i, 0);
1563 ppath1.
pos(i1, 1) = ppath2.
pos(i, 1);
1564 ppath1.
los(i1, 0) = ppath2.
los(i, 0);
1565 ppath1.
r[i1] = ppath2.
r[i];
1570 if (ppath1.
dim >= 2) {
1574 if (ppath1.
dim == 3) {
1575 ppath1.
pos(i1, 2) = ppath2.
pos(i, 2);
1576 ppath1.
los(i1, 1) = ppath2.
los(i, 1);
1610 const Ppath& ppath) {
1612 const Index imax = ppath.
np - 1;
1615 r_start = ppath.
r[imax];
1616 lat_start = ppath.
pos(imax, 1);
1617 za_start = ppath.
los(imax, 0);
1636 ConstVectorView r_v,
1637 ConstVectorView lat_v,
1638 ConstVectorView za_v,
1639 ConstVectorView lstep,
1640 ConstVectorView n_v,
1641 ConstVectorView ng_v,
1642 ConstVectorView z_field,
1643 ConstVectorView refellipsoid,
1645 const Index& endface,
1646 const Numeric& ppc) {
1648 const Index np = r_v.nelem();
1656 const Numeric r1 = refellipsoid[0] + z_field[ip];
1657 const Numeric dr = z_field[ip + 1] - z_field[ip];
1659 for (Index i = 0; i < np; i++) {
1660 ppath.
r[i] = r_v[i];
1661 ppath.
pos(i, 0) = r_v[i] - refellipsoid[0];
1662 ppath.
pos(i, 1) = lat_v[i];
1663 ppath.
los(i, 0) = za_v[i];
1664 ppath.
nreal[i] = n_v[i];
1665 ppath.
ngroup[i] = ng_v[i];
1667 ppath.
gp_p[i].idx = ip;
1668 ppath.
gp_p[i].fd[0] = (r_v[i] - r1) / dr;
1669 ppath.
gp_p[i].fd[1] = 1 - ppath.
gp_p[i].fd[0];
1673 ppath.
lstep[i - 1] = lstep[i - 1];
1684 else if (endface <= 4) {
1713 ConstVectorView lat_grid,
1714 ConstMatrixView z_field,
1715 ConstVectorView refellipsoid,
1716 ConstVectorView z_surface) {
1718 const Index imax = ppath.
np - 1;
1721 r_start = ppath.
r[imax];
1722 lat_start = ppath.
pos(imax, 1);
1723 za_start = ppath.
los(imax, 0);
1730 lat1 = lat_grid[ilat];
1731 lat3 = lat_grid[ilat + 1];
1741 const Numeric re1 =
refell2r(refellipsoid, lat_grid[ilat]);
1742 const Numeric re3 =
refell2r(refellipsoid, lat_grid[ilat + 1]);
1744 r1a = re1 + z_field(ip, ilat);
1745 r3a = re3 + z_field(ip, ilat + 1);
1746 r3b = re3 + z_field(ip + 1, ilat + 1);
1747 r1b = re1 + z_field(ip + 1, ilat);
1753 const Numeric rlow =
rsurf_at_lat(lat1, lat3, r1a, r3a, lat_start);
1754 const Numeric rupp =
rsurf_at_lat(lat1, lat3, r1b, r3b, lat_start);
1755 if (abs(r_start - rlow) <
RTOL || abs(r_start - rupp) <
RTOL) {
1777 r1a = re1 + z_field(ip, ilat);
1778 r3a = re3 + z_field(ip, ilat + 1);
1788 r3b = re3 + z_field(ip + 1, ilat + 1);
1789 r1b = re1 + z_field(ip + 1, ilat);
1795 rsurface1 = re1 + z_surface[ilat];
1796 rsurface3 = re3 + z_surface[ilat + 1];
1810 ConstVectorView r_v,
1811 ConstVectorView lat_v,
1812 ConstVectorView za_v,
1813 ConstVectorView lstep,
1814 ConstVectorView n_v,
1815 ConstVectorView ng_v,
1816 ConstVectorView lat_grid,
1817 ConstMatrixView z_field,
1818 ConstVectorView refellipsoid,
1821 const Index& endface,
1822 const Numeric& ppc) {
1824 const Index np = r_v.nelem();
1825 const Index imax = np - 1;
1834 const Numeric dlat = lat_grid[ilat + 1] - lat_grid[ilat];
1835 const Numeric z1low = z_field(ip, ilat);
1836 const Numeric z1upp = z_field(ip + 1, ilat);
1837 const Numeric dzlow = z_field(ip, ilat + 1) - z1low;
1838 const Numeric dzupp = z_field(ip + 1, ilat + 1) - z1upp;
1839 Numeric re =
refell2r(refellipsoid, lat_grid[ilat]);
1840 const Numeric r1low = re + z1low;
1841 const Numeric r1upp = re + z1upp;
1842 re =
refell2r(refellipsoid, lat_grid[ilat + 1]);
1843 const Numeric drlow = re + z_field(ip, ilat + 1) - r1low;
1844 const Numeric drupp = re + z_field(ip + 1, ilat + 1) - r1upp;
1846 for (Index i = 0; i < np; i++) {
1847 ppath.
r[i] = r_v[i];
1848 ppath.
pos(i, 1) = lat_v[i];
1849 ppath.
los(i, 0) = za_v[i];
1850 ppath.
nreal[i] = n_v[i];
1851 ppath.
ngroup[i] = ng_v[i];
1854 Numeric
w = (lat_v[i] - lat_grid[ilat]) / dlat;
1857 const Numeric rlow = r1low +
w * drlow;
1858 const Numeric rupp = r1upp +
w * drupp;
1861 const Numeric zlow = z1low +
w * dzlow;
1862 const Numeric zupp = z1upp +
w * dzupp;
1864 ppath.
gp_p[i].idx = ip;
1865 ppath.
gp_p[i].fd[0] = (r_v[i] - rlow) / (rupp - rlow);
1866 ppath.
gp_p[i].fd[1] = 1 - ppath.
gp_p[i].fd[0];
1869 ppath.
pos(i, 0) = zlow + ppath.
gp_p[i].fd[0] * (zupp - zlow);
1871 ppath.
gp_lat[i].idx = ilat;
1872 ppath.
gp_lat[i].fd[0] = (lat_v[i] - lat_grid[ilat]) / dlat;
1877 ppath.
lstep[i - 1] = lstep[i - 1];
1890 if (endface == 1 || endface == 3) {
1892 }
else if (endface == 2 || endface == 4) {
1898 if (ppath.
gp_p[imax].fd[0] < 0 || ppath.
gp_p[imax].fd[1] < 0) {
1901 if (ppath.
gp_lat[imax].fd[0] < 0 || ppath.
gp_lat[imax].fd[1] < 0) {
1936 Numeric& rsurface15,
1937 Numeric& rsurface35,
1938 Numeric& rsurface36,
1939 Numeric& rsurface16,
1941 ConstVectorView lat_grid,
1942 ConstVectorView lon_grid,
1943 ConstTensor3View z_field,
1944 ConstVectorView refellipsoid,
1945 ConstMatrixView z_surface) {
1947 const Index imax = ppath.
np - 1;
1950 r_start = ppath.
r[imax];
1951 lat_start = ppath.
pos(imax, 1);
1952 lon_start = ppath.
pos(imax, 2);
1953 za_start = ppath.
los(imax, 0);
1954 aa_start = ppath.
los(imax, 1);
1957 const Index nlat = lat_grid.nelem();
1958 const Index nlon = lon_grid.nelem();
1965 if (lat_start == 90) {
1968 gridpos(gp_tmp, lon_grid, aa_start);
1969 if (aa_start < 180) {
1974 }
else if (lat_start == -90) {
1977 gridpos(gp_tmp, lon_grid, aa_start);
1978 if (aa_start < 180) {
1984 if (lat_start > 0) {
1989 if (lon_start < lon_grid[nlon - 1]) {
1996 lat1 = lat_grid[ilat];
1997 lat3 = lat_grid[ilat + 1];
1998 lon5 = lon_grid[ilon];
1999 lon6 = lon_grid[ilon + 1];
2009 const Numeric re1 =
refell2r(refellipsoid, lat_grid[ilat]);
2010 const Numeric re3 =
refell2r(refellipsoid, lat_grid[ilat + 1]);
2012 r15a = re1 + z_field(ip, ilat, ilon);
2013 r35a = re3 + z_field(ip, ilat + 1, ilon);
2014 r36a = re3 + z_field(ip, ilat + 1, ilon + 1);
2015 r16a = re1 + z_field(ip, ilat, ilon + 1);
2016 r15b = re1 + z_field(ip + 1, ilat, ilon);
2017 r35b = re3 + z_field(ip + 1, ilat + 1, ilon);
2018 r36b = re3 + z_field(ip + 1, ilat + 1, ilon + 1);
2019 r16b = re1 + z_field(ip + 1, ilat, ilon + 1);
2026 if (fabs(za_start - 90) <= 10)
2052 r15a = re1 + z_field(ip, ilat, ilon);
2053 r35a = re3 + z_field(ip, ilat + 1, ilon);
2054 r36a = re3 + z_field(ip, ilat + 1, ilon + 1);
2055 r16a = re1 + z_field(ip, ilat, ilon + 1);
2081 r15b = re1 + z_field(ip + 1, ilat, ilon);
2082 r35b = re3 + z_field(ip + 1, ilat + 1, ilon);
2083 r36b = re3 + z_field(ip + 1, ilat + 1, ilon + 1);
2084 r16b = re1 + z_field(ip + 1, ilat, ilon + 1);
2090 rsurface15 = re1 + z_surface(ilat, ilon);
2091 rsurface35 = re3 + z_surface(ilat + 1, ilon);
2092 rsurface36 = re3 + z_surface(ilat + 1, ilon + 1);
2093 rsurface16 = re1 + z_surface(ilat, ilon + 1);
2107 ConstVectorView r_v,
2108 ConstVectorView lat_v,
2109 ConstVectorView lon_v,
2110 ConstVectorView za_v,
2111 ConstVectorView aa_v,
2112 ConstVectorView lstep,
2113 ConstVectorView n_v,
2114 ConstVectorView ng_v,
2115 ConstVectorView lat_grid,
2116 ConstVectorView lon_grid,
2117 ConstTensor3View z_field,
2118 ConstVectorView refellipsoid,
2122 const Index& endface,
2123 const Numeric& ppc) {
2125 const Index np = r_v.nelem();
2126 const Index imax = np - 1;
2135 const Numeric lat1 = lat_grid[ilat];
2136 const Numeric lat3 = lat_grid[ilat + 1];
2137 const Numeric lon5 = lon_grid[ilon];
2138 const Numeric lon6 = lon_grid[ilon + 1];
2139 const Numeric re1 =
refell2r(refellipsoid, lat_grid[ilat]);
2140 const Numeric re3 =
refell2r(refellipsoid, lat_grid[ilat + 1]);
2141 const Numeric r15a = re1 + z_field(ip, ilat, ilon);
2142 const Numeric r35a = re3 + z_field(ip, ilat + 1, ilon);
2143 const Numeric r36a = re3 + z_field(ip, ilat + 1, ilon + 1);
2144 const Numeric r16a = re1 + z_field(ip, ilat, ilon + 1);
2145 const Numeric r15b = re1 + z_field(ip + 1, ilat, ilon);
2146 const Numeric r35b = re3 + z_field(ip + 1, ilat + 1, ilon);
2147 const Numeric r36b = re3 + z_field(ip + 1, ilat + 1, ilon + 1);
2148 const Numeric r16b = re1 + z_field(ip + 1, ilat, ilon + 1);
2149 const Numeric dlat = lat3 - lat1;
2150 const Numeric dlon = lon6 - lon5;
2152 for (Index i = 0; i < np; i++) {
2155 lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, lat_v[i], lon_v[i]);
2157 lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b, lat_v[i], lon_v[i]);
2160 ppath.
r[i] = r_v[i];
2161 ppath.
pos(i, 1) = lat_v[i];
2162 ppath.
pos(i, 2) = lon_v[i];
2163 ppath.
los(i, 0) = za_v[i];
2164 ppath.
los(i, 1) = aa_v[i];
2165 ppath.
nreal[i] = n_v[i];
2166 ppath.
ngroup[i] = ng_v[i];
2169 ppath.
gp_p[i].idx = ip;
2170 ppath.
gp_p[i].fd[0] = (r_v[i] - rlow) / (rupp - rlow);
2171 ppath.
gp_p[i].fd[1] = 1 - ppath.
gp_p[i].fd[0];
2176 lat1, lat3, lon5, lon6, re1, re3, re3, re1, lat_v[i], lon_v[i]);
2177 const Numeric zlow = rlow - re;
2178 const Numeric zupp = rupp - re;
2180 ppath.
pos(i, 0) = zlow + ppath.
gp_p[i].fd[0] * (zupp - zlow);
2183 ppath.
gp_lat[i].idx = ilat;
2184 ppath.
gp_lat[i].fd[0] = (lat_v[i] - lat1) / dlat;
2193 if (abs(lat_v[i]) <
POLELAT) {
2194 ppath.
gp_lon[i].idx = ilon;
2195 ppath.
gp_lon[i].fd[0] = (lon_v[i] - lon5) / dlon;
2200 ppath.
gp_lon[i].fd[0] = 0;
2201 ppath.
gp_lon[i].fd[1] = 1;
2205 ppath.
lstep[i - 1] = lstep[i - 1];
2216 if (endface == 1 || endface == 3) {
2218 }
else if (endface == 2 || endface == 4) {
2220 }
else if (endface == 5 || endface == 6) {
2226 if (ppath.
gp_p[imax].fd[0] < 0 || ppath.
gp_p[imax].fd[1] < 0) {
2229 if (ppath.
gp_lat[imax].fd[0] < 0 || ppath.
gp_lat[imax].fd[1] < 0) {
2232 if (ppath.
gp_lon[imax].fd[0] < 0 || ppath.
gp_lon[imax].fd[1] < 0) {
2270 const Numeric& r_start0,
2271 const Numeric& lat_start,
2272 const Numeric& za_start,
2274 const Numeric& lmax,
2277 const Numeric& rsurface) {
2278 Numeric r_start = r_start0;
2287 }
else if (r_start > rb) {
2292 bool tanpoint =
false;
2296 if (za_start <= 90) {
2303 if (ra > rsurface && ra > ppc) {
2309 else if (rsurface > ppc) {
2338 ConstVectorView z_field,
2339 ConstVectorView refellipsoid,
2340 const Numeric& z_surface,
2341 const Numeric& lmax) {
2343 Numeric r_start, lat_start, za_start;
2365 Vector r_v, lat_v, za_v;
2379 refellipsoid[0] + z_field[ip],
2380 refellipsoid[0] + z_field[ip + 1],
2381 refellipsoid[0] + z_surface);
2384 const Index np = r_v.nelem();
2389 Vector(np - 1, lstep),
2409 const Numeric& r_start0,
2410 const Numeric& lat_start0,
2411 const Numeric& za_start,
2412 const Numeric& l_start,
2415 const Numeric& lmax,
2416 const Numeric& lat1,
2417 const Numeric& lat3,
2422 const Numeric& rsurface1,
2423 const Numeric& rsurface3) {
2424 Numeric r_start = r_start0;
2425 Numeric lat_start = lat_start0;
2434 if (lat_start < lat1) {
2436 }
else if (lat_start > lat3) {
2441 Numeric rlow =
rsurf_at_lat(lat1, lat3, r1a, r3a, lat_start);
2442 Numeric rupp =
rsurf_at_lat(lat1, lat3, r1b, r3b, lat_start);
2449 if (r_start < rlow) {
2451 }
else if (r_start > rupp) {
2456 Numeric x, z, dx, dz;
2457 poslos2cart(x, z, dx, dz, r_start, lat_start, za_start);
2460 bool unsafe =
false;
2461 bool do_surface =
false;
2467 Numeric r_end, lat_end, l_end;
2470 const Numeric absza = abs(za_start);
2474 l_end = rupp - r_start;
2479 else if (absza > 180 -
ANGTOL) {
2480 const Numeric rsurface =
2481 rsurf_at_lat(lat1, lat3, rsurface1, rsurface3, lat_start);
2483 if (rlow > rsurface) {
2484 l_end = r_start - rlow;
2487 l_end = r_start - rsurface;
2498 Numeric r_corr, lat_corr;
2500 cart2pol(r_corr, lat_corr, x, z, lat_start, za_start);
2503 lat_corr -= lat_start;
2518 l_end = 2 * (rupp - rlow);
2521 Numeric l_in = 0, l_out = l_end;
2522 bool ready =
false, startup =
true;
2524 if (rsurface1 +
RTOL >= r1a || rsurface3 +
RTOL >= r3a) {
2530 r_end, lat_end, x + dx * l_end, z + dz * l_end, lat_start, za_start);
2532 lat_end -= lat_corr;
2539 const Numeric r_surface =
2540 rsurf_at_lat(lat1, lat3, rsurface1, rsurface3, lat_end);
2541 if (r_surface >= rlow && r_end <= r_surface) {
2548 if (lat_end < lat1) {
2551 }
else if (lat_end > lat3) {
2554 }
else if (r_end < rlow) {
2572 l_end = (l_out + l_in) / 2;
2582 if ((l_out - l_in) <
LACC) {
2585 l_end = (l_out + l_in) / 2;
2596 n = Index(ceil(abs(l_end / lmax)));
2603 lat_v.resize(n + 1);
2607 lat_v[0] = lat_start;
2610 lstep = l_end / (Numeric)n;
2614 for (Index j = 1; j <= n; j++) {
2615 l = lstep * (Numeric)j;
2639 const Numeric r_surface =
2640 rsurf_at_lat(lat1, lat3, rsurface1, rsurface3, lat_v[j]);
2641 const Numeric r_test =
max(r_surface, rlow);
2642 if (r_v[j] < r_test) {
2646 }
else if (r_v[j] < rlow) {
2652 if (r_v[j] > rupp) {
2664 }
else if (endface == 2) {
2666 }
else if (endface == 3) {
2668 }
else if (endface == 4) {
2670 }
else if (endface == 7) {
2671 r_v[n] =
rsurf_at_lat(lat1, lat3, rsurface1, rsurface3, lat_v[n]);
2702 ConstVectorView lat_grid,
2703 ConstMatrixView z_field,
2704 ConstVectorView refellipsoid,
2705 ConstVectorView z_surface,
2706 const Numeric& lmax) {
2708 Numeric r_start, lat_start, za_start;
2714 Numeric lat1, lat3, r1a, r3a, r3b, r1b, rsurface1, rsurface3;
2746 Vector r_v, lat_v, za_v;
2771 const Index np = r_v.nelem();
2776 Vector(np - 1, lstep),
2800 const Numeric& r_start0,
2801 const Numeric& lat_start0,
2802 const Numeric& lon_start0,
2803 const Numeric& za_start,
2804 const Numeric& aa_start,
2805 const Numeric& l_start,
2808 const Numeric& lmax,
2809 const Numeric& lat1,
2810 const Numeric& lat3,
2811 const Numeric& lon5,
2812 const Numeric& lon6,
2813 const Numeric& r15a,
2814 const Numeric& r35a,
2815 const Numeric& r36a,
2816 const Numeric& r16a,
2817 const Numeric& r15b,
2818 const Numeric& r35b,
2819 const Numeric& r36b,
2820 const Numeric& r16b,
2821 const Numeric& rsurface15,
2822 const Numeric& rsurface35,
2823 const Numeric& rsurface36,
2824 const Numeric& rsurface16) {
2825 Numeric r_start = r_start0;
2826 Numeric lat_start = lat_start0;
2827 Numeric lon_start = lon_start0;
2841 if (lat_start < lat1) {
2843 }
else if (lat_start > lat3) {
2846 if (lon_start < lon5) {
2848 }
else if (lon_start > lon6) {
2854 lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, lat_start, lon_start);
2856 lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b, lat_start, lon_start);
2863 if (r_start < rlow) {
2865 }
else if (r_start > rupp) {
2870 Numeric x, y, z, dx, dy, dz;
2872 x, y, z, dx, dy, dz, r_start, lat_start, lon_start, za_start, aa_start);
2875 bool unsafe =
false;
2876 bool do_surface =
false;
2882 Numeric r_end, lat_end, lon_end, l_end;
2888 l_end = rupp - r_start;
2893 else if (za_start > 180 -
ANGTOL) {
2905 if (rlow > rsurface) {
2906 l_end = r_start - rlow;
2909 l_end = r_start - rsurface;
2920 Numeric r_corr, lat_corr, lon_corr;
2934 lat_corr -= lat_start;
2935 lon_corr -= lon_start;
2950 l_end = 2 * (rupp - rlow);
2953 Numeric l_in = 0, l_out = l_end;
2954 bool ready =
false, startup =
true;
2956 if (rsurface15 +
RTOL >= r15a || rsurface35 +
RTOL >= r35a ||
2957 rsurface36 +
RTOL >= r36a || rsurface16 +
RTOL >= r16a) {
2973 lat_end -= lat_corr;
2974 lon_end -= lon_corr;
2979 (abs(aa_start) <
ANGTOL || abs(aa_start) > 180 -
ANGTOL)) {
2980 lon_end = lon_start;
2986 lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, lat_end, lon_end);
2999 if (r_surface >= rlow && r_end <= r_surface) {
3006 if (lat_end < lat1) {
3009 }
else if (lat_end > lat3) {
3012 }
else if (lon_end < lon5) {
3015 }
else if (lon_end > lon6) {
3018 }
else if (r_end < rlow) {
3023 lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b, lat_end, lon_end);
3037 l_end = (l_out + l_in) / 2;
3047 if ((l_out - l_in) <
LACC) {
3050 l_end = (l_out + l_in) / 2;
3061 n = Index(ceil(abs(l_end / lmax)));
3068 lat_v.resize(n + 1);
3069 lon_v.resize(n + 1);
3074 lat_v[0] = lat_start;
3075 lon_v[0] = lon_start;
3079 lstep = l_end / (Numeric)n;
3083 for (Index j = 1; j <= n; j++) {
3084 l = lstep * (Numeric)j;
3119 lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, lat_v[j], lon_v[j]);
3131 const Numeric r_test =
max(r_surface, rlow);
3132 if (r_v[j] < r_test) {
3136 }
else if (r_v[j] < rlow) {
3142 lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b, lat_v[j], lon_v[j]);
3143 if (r_v[j] > rupp) {
3155 }
else if (endface == 2) {
3166 }
else if (endface == 3) {
3168 }
else if (endface == 4) {
3179 }
else if (endface == 5) {
3181 }
else if (endface == 6) {
3183 }
else if (endface == 7) {
3236 ConstVectorView lat_grid,
3237 ConstVectorView lon_grid,
3238 ConstTensor3View z_field,
3239 ConstVectorView refellipsoid,
3240 ConstMatrixView z_surface,
3241 const Numeric& lmax) {
3243 Numeric r_start, lat_start, lon_start, za_start, aa_start;
3246 Index ip, ilat, ilon;
3250 Numeric lat1, lat3, lon5, lon6;
3251 Numeric r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b;
3252 Numeric rsurface15, rsurface35, rsurface36, rsurface16;
3296 Vector r_v, lat_v, lon_v, za_v, aa_v;
3334 const Index np = r_v.nelem();
3341 Vector(np - 1, lstep),
3405 ConstVectorView p_grid,
3406 ConstVectorView refellipsoid,
3407 ConstTensor3View z_field,
3408 ConstTensor3View t_field,
3409 ConstTensor4View vmr_field,
3410 ConstVectorView f_grid,
3411 const Numeric& lmax,
3412 const Agenda& refr_index_air_agenda,
3413 const Numeric& lraytrace,
3414 const Numeric& rsurface,
3424 Numeric refr_index_air, refr_index_air_group;
3427 refr_index_air_group,
3428 refr_index_air_agenda,
3436 r_array.push_back(r);
3437 lat_array.push_back(lat);
3438 za_array.push_back(za);
3439 n_array.push_back(refr_index_air);
3440 ng_array.push_back(refr_index_air_group);
3443 Vector r_v, lat_v, za_v;
3444 Numeric lstep, lcum = 0, dlat;
3470 Numeric za_flagside = za;
3472 if (lstep <= lraytrace) {
3474 dlat = lat_v[1] - lat;
3485 za_flagside = 180 - za_flagside;
3493 dlat = lat_new - lat;
3503 refr_index_air_group,
3505 refr_index_air_agenda,
3515 const Numeric za_rad =
DEG2RAD * za;
3519 za += (
RAD2DEG * lstep / refr_index_air) * (-sin(za_rad) * dndr);
3524 }
else if (za > 180) {
3529 if (ready || (lmax > 0 && lcum + lraytrace > lmax)) {
3530 r_array.push_back(r);
3531 lat_array.push_back(lat);
3532 za_array.push_back(za);
3533 n_array.push_back(refr_index_air);
3534 ng_array.push_back(refr_index_air_group);
3535 l_array.push_back(lcum);
3543 ConstVectorView p_grid,
3544 ConstTensor3View z_field,
3545 ConstTensor3View t_field,
3546 ConstTensor4View vmr_field,
3547 ConstVectorView f_grid,
3548 ConstVectorView refellipsoid,
3549 const Numeric& z_surface,
3550 const Numeric& lmax,
3551 const Agenda& refr_index_air_agenda,
3552 const String& rtrace_method,
3553 const Numeric& lraytrace) {
3555 Numeric r_start, lat_start, za_start;
3570 Numeric refr_index_air, refr_index_air_group;
3573 refr_index_air_group,
3574 refr_index_air_agenda,
3591 Array<Numeric> r_array, lat_array, za_array, l_array, n_array, ng_array;
3594 if (rtrace_method ==
"linear_basic") {
3618 refr_index_air_agenda,
3620 refellipsoid[0] + z_surface,
3621 refellipsoid[0] + z_field(ip, 0, 0),
3622 refellipsoid[0] + z_field(ip + 1, 0, 0),
3633 const Index np = r_array.
nelem();
3634 Vector r_v(np), lat_v(np), za_v(np), l_v(np - 1), n_v(np), ng_v(np);
3635 for (Index i = 0; i < np; i++) {
3636 r_v[i] = r_array[i];
3637 lat_v[i] = lat_array[i];
3638 za_v[i] = za_array[i];
3639 n_v[i] = n_array[i];
3640 ng_v[i] = ng_array[i];
3642 l_v[i] = l_array[i];
3653 z_field(joker, 0, 0),
3711 ConstVectorView p_grid,
3712 ConstVectorView lat_grid,
3713 ConstVectorView refellipsoid,
3714 ConstTensor3View z_field,
3715 ConstTensor3View t_field,
3716 ConstTensor4View vmr_field,
3717 ConstVectorView f_grid,
3718 const Numeric& lmax,
3719 const Agenda& refr_index_air_agenda,
3720 const Numeric& lraytrace,
3721 const Numeric& lat1,
3722 const Numeric& lat3,
3723 const Numeric& rsurface1,
3724 const Numeric& rsurface3,
3736 Numeric refr_index_air, refr_index_air_group;
3739 refr_index_air_group,
3740 refr_index_air_agenda,
3750 r_array.push_back(r);
3751 lat_array.push_back(lat);
3752 za_array.push_back(za);
3753 n_array.push_back(refr_index_air);
3754 ng_array.push_back(refr_index_air_group);
3757 Vector r_v, lat_v, za_v;
3758 Numeric lstep, lcum = 0, dlat;
3791 Numeric za_flagside = za;
3793 if (lstep <= lraytrace) {
3795 dlat = lat_v[1] - lat;
3801 if (abs(za) <= 90) {
3807 za_flagside =
sign(za) * 180 - za_flagside;
3815 dlat = lat_new - lat;
3824 }
else if (lat > lat3) {
3830 Numeric dndr, dndlat;
3833 refr_index_air_group,
3836 refr_index_air_agenda,
3848 const Numeric za_rad =
DEG2RAD * za;
3852 za += (
RAD2DEG * lstep / refr_index_air) *
3853 (-sin(za_rad) * dndr + cos(za_rad) * dndlat);
3858 }
else if (za > 180) {
3864 if (lat == lat1 && za < 0) {
3867 }
else if (lat == lat3 && za > 0) {
3873 if (ready || (lmax > 0 && lcum + lraytrace > lmax)) {
3874 r_array.push_back(r);
3875 lat_array.push_back(lat);
3876 za_array.push_back(za);
3877 n_array.push_back(refr_index_air);
3878 ng_array.push_back(refr_index_air_group);
3879 l_array.push_back(lcum);
3887 ConstVectorView p_grid,
3888 ConstVectorView lat_grid,
3889 ConstTensor3View z_field,
3890 ConstTensor3View t_field,
3891 ConstTensor4View vmr_field,
3892 ConstVectorView f_grid,
3893 ConstVectorView refellipsoid,
3894 ConstVectorView z_surface,
3895 const Numeric& lmax,
3896 const Agenda& refr_index_air_agenda,
3897 const String& rtrace_method,
3898 const Numeric& lraytrace) {
3900 Numeric r_start, lat_start, za_start;
3906 Numeric lat1, lat3, r1a, r3a, r3b, r1b, rsurface1, rsurface3;
3924 z_field(joker, joker, 0),
3934 Array<Numeric> r_array, lat_array, za_array, l_array, n_array, ng_array;
3937 if (rtrace_method ==
"linear_basic") {
3954 refr_index_air_agenda,
3974 const Index np = r_array.
nelem();
3975 Vector r_v(np), lat_v(np), za_v(np), l_v(np - 1), n_v(np), ng_v(np);
3976 for (Index i = 0; i < np; i++) {
3977 r_v[i] = r_array[i];
3978 lat_v[i] = lat_array[i];
3979 za_v[i] = za_array[i];
3980 n_v[i] = n_array[i];
3981 ng_v[i] = ng_array[i];
3983 l_v[i] = l_array[i];
3995 z_field(joker, joker, 0),
4070 ConstVectorView refellipsoid,
4071 ConstVectorView p_grid,
4072 ConstVectorView lat_grid,
4073 ConstVectorView lon_grid,
4074 ConstTensor3View z_field,
4075 ConstTensor3View t_field,
4076 ConstTensor4View vmr_field,
4077 ConstVectorView f_grid,
4078 const Numeric& lmax,
4079 const Agenda& refr_index_air_agenda,
4080 const Numeric& lraytrace,
4081 const Numeric& lat1,
4082 const Numeric& lat3,
4083 const Numeric& lon5,
4084 const Numeric& lon6,
4085 const Numeric& rsurface15,
4086 const Numeric& rsurface35,
4087 const Numeric& rsurface36,
4088 const Numeric& rsurface16,
4089 const Numeric& r15a,
4090 const Numeric& r35a,
4091 const Numeric& r36a,
4092 const Numeric& r16a,
4093 const Numeric& r15b,
4094 const Numeric& r35b,
4095 const Numeric& r36b,
4096 const Numeric& r16b,
4106 Numeric refr_index_air, refr_index_air_group;
4109 refr_index_air_group,
4110 refr_index_air_agenda,
4122 r_array.push_back(r);
4123 lat_array.push_back(lat);
4124 lon_array.push_back(lon);
4125 za_array.push_back(za);
4126 aa_array.push_back(aa);
4127 n_array.push_back(refr_index_air);
4128 ng_array.push_back(refr_index_air_group);
4131 Vector r_v, lat_v, lon_v, za_v, aa_v;
4132 Numeric lstep, lcum = 0;
4133 Numeric za_new, aa_new;
4177 if (lstep <= lraytrace) {
4187 Numeric x, y, z, dx, dy, dz, lat_new, lon_new;
4189 poslos2cart(x, y, z, dx, dy, dz, r, lat, lon, za, aa);
4220 Numeric dndr, dndlat, dndlon;
4223 refr_index_air_group,
4227 refr_index_air_agenda,
4241 const Numeric aterm =
RAD2DEG * lstep / refr_index_air;
4242 const Numeric za_rad =
DEG2RAD * za;
4243 const Numeric aa_rad =
DEG2RAD * aa;
4244 const Numeric sinza = sin(za_rad);
4245 const Numeric sinaa = sin(aa_rad);
4246 const Numeric cosaa = cos(aa_rad);
4252 if (za < ANGTOL || za > 180 -
ANGTOL) {
4253 los[0] += aterm * (cos(za_rad) * (cosaa * dndlat + sinaa * dndlon));
4254 los[1] =
RAD2DEG * atan2(dndlon, dndlat);
4256 los[0] += aterm * (-sinza * dndr +
4257 cos(za_rad) * (cosaa * dndlat + sinaa * dndlon));
4258 los[1] += aterm * sinza * (cosaa * dndlon - sinaa * dndlat);
4269 if (za > 0 && za < 180) {
4270 if (lon == lon5 && aa < 0) {
4273 }
else if (lon == lon6 && aa > 0) {
4276 }
else if (lat == lat1 && lat != -90 && abs(aa) > 90) {
4279 }
else if (lat == lat3 && lat != 90 && abs(aa) < 90) {
4286 if (ready || (lmax > 0 && lcum + lraytrace > lmax)) {
4287 r_array.push_back(r);
4288 lat_array.push_back(lat);
4289 lon_array.push_back(lon);
4290 za_array.push_back(za);
4291 aa_array.push_back(aa);
4292 n_array.push_back(refr_index_air);
4293 ng_array.push_back(refr_index_air_group);
4294 l_array.push_back(lcum);
4302 ConstVectorView p_grid,
4303 ConstVectorView lat_grid,
4304 ConstVectorView lon_grid,
4305 ConstTensor3View z_field,
4306 ConstTensor3View t_field,
4307 ConstTensor4View vmr_field,
4308 ConstVectorView f_grid,
4309 ConstVectorView refellipsoid,
4310 ConstMatrixView z_surface,
4311 const Numeric& lmax,
4312 const Agenda& refr_index_air_agenda,
4313 const String& rtrace_method,
4314 const Numeric& lraytrace) {
4316 Numeric r_start, lat_start, lon_start, za_start, aa_start;
4319 Index ip, ilat, ilon;
4323 Numeric lat1, lat3, lon5, lon6;
4324 Numeric r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b;
4325 Numeric rsurface15, rsurface35, rsurface36, rsurface16;
4365 Array<Numeric> r_array, lat_array, lon_array, za_array, aa_array;
4369 if (rtrace_method ==
"linear_basic") {
4389 refr_index_air_agenda,
4419 const Index np = r_array.
nelem();
4420 Vector r_v(np), lat_v(np), lon_v(np), za_v(np), aa_v(np), l_v(np - 1);
4421 Vector n_v(np), ng_v(np);
4422 for (Index i = 0; i < np; i++) {
4423 r_v[i] = r_array[i];
4424 lat_v[i] = lat_array[i];
4425 lon_v[i] = lon_array[i];
4426 za_v[i] = za_array[i];
4427 aa_v[i] = aa_array[i];
4428 n_v[i] = n_array[i];
4429 ng_v[i] = ng_array[i];
4431 l_v[i] = l_array[i];
4461 const Index& atmosphere_dim,
4462 ConstVectorView p_grid,
4463 ConstVectorView lat_grid,
4464 ConstVectorView lon_grid,
4465 ConstTensor3View z_field,
4466 ConstVectorView refellipsoid,
4467 ConstMatrixView z_surface,
4468 const Index& cloudbox_on,
4470 const bool& ppath_inside_cloudbox_do,
4471 ConstVectorView rte_pos,
4472 ConstVectorView rte_los,
4482 const Index lp = p_grid.nelem() - 1;
4487 if (atmosphere_dim == 1) {
4489 ppath.
end_pos[0] = rte_pos[0];
4491 ppath.
end_los[0] = rte_los[0];
4494 if (rte_pos[0] < z_field(lp, 0, 0)) {
4497 "The ppath starting point is placed ",
4498 (z_surface(0, 0) - rte_pos[0]) / 1e3,
" km below the surface.")
4502 ppath.
r[0] = refellipsoid[0] + rte_pos[0];
4505 gridpos(ppath.
gp_p, z_field(joker, 0, 0), ExhaustiveConstVectorView{ppath.
pos(0, 0)});
4511 if (ppath.
pos(0, 0) <= z_surface(0, 0) && ppath.
los(0, 0) > 90) {
4517 if (cloudbox_on && !ppath_inside_cloudbox_do) {
4519 if (fgp >= (Numeric)cloudbox_limits[0] &&
4520 fgp <= (Numeric)cloudbox_limits[1]) {
4529 fgp <= (Numeric)cloudbox_limits[1]);
4542 if (rte_los[0] <= 90 ||
4543 ppath.
constant >= refellipsoid[0] + z_field(lp, 0, 0)) {
4544 ppath.
pos(0, 0) = rte_pos[0];
4545 ppath.
pos(0, 1) = 0;
4546 ppath.
r[0] = refellipsoid[0] + rte_pos[0];
4547 ppath.
los(0, 0) = rte_los[0];
4550 out1 <<
" --- WARNING ---, path is totally outside of the "
4551 <<
"model atmosphere\n";
4556 ppath.
r[0] = refellipsoid[0] + z_field(lp, 0, 0);
4557 ppath.
pos(0, 0) = z_field(lp, 0, 0);
4566 ppath.
gp_p[0].idx = lp - 1;
4567 ppath.
gp_p[0].fd[0] = 1;
4568 ppath.
gp_p[0].fd[1] = 0;
4571 if (cloudbox_on && cloudbox_limits[1] == lp) {
4579 else if (atmosphere_dim == 2) {
4581 ppath.
end_pos[0] = rte_pos[0];
4582 ppath.
end_pos[1] = rte_pos[1];
4583 ppath.
end_los[0] = rte_los[0];
4586 const Index llat = lat_grid.nelem() - 1;
4592 bool islatin =
false;
4594 Numeric z_toa = -99e99;
4595 if (rte_pos[1] > lat_grid[0] && rte_pos[1] < lat_grid[llat]) {
4597 gridpos(gp_lat, lat_grid, rte_pos[1]);
4599 z_toa =
interp(itw, z_field(lp, joker, 0), gp_lat);
4600 r_e =
refell2d(refellipsoid, lat_grid, gp_lat);
4602 r_e =
refell2r(refellipsoid, rte_pos[1]);
4606 if (islatin && rte_pos[0] < z_toa) {
4607 const Numeric z_s =
interp(itw, z_surface(joker, 0), gp_lat);
4611 "The ppath starting point is placed ", (z_s - rte_pos[0]) / 1e3,
4612 " km below the surface.")
4616 ppath.
r[0] = r_e + rte_pos[0];
4636 if (ppath.
pos(0, 0) <= z_s) {
4642 z_surface(joker, 0),
4659 if (cloudbox_on && !ppath_inside_cloudbox_do) {
4662 if (fgp >= (Numeric)cloudbox_limits[0] &&
4663 fgp <= (Numeric)cloudbox_limits[1] &&
4664 fgl >= (Numeric)cloudbox_limits[2] &&
4665 fgl <= (Numeric)cloudbox_limits[3]) {
4675 fgp <= (Numeric)cloudbox_limits[1] &&
4676 fgl >= (Numeric)cloudbox_limits[2] &&
4677 fgl <= (Numeric)cloudbox_limits[3]);
4685 (rte_pos[1] >= lat_grid[llat] && rte_los[0] >= 0),
4686 "The sensor is outside (or at the limit) of the model "
4687 "atmosphere but\nlooks in the wrong direction (wrong sign "
4688 "for the zenith angle?).\nThis case includes nadir "
4689 "looking exactly at the latitude end points.")
4694 const Numeric r_p = r_e + rte_pos[0];
4698 Vector r_toa(llat + 1);
4699 Numeric r_toa_min = 99e99, r_toa_max = -1;
4700 for (Index ilat = 0; ilat <= llat; ilat++) {
4702 refell2r(refellipsoid, lat_grid[ilat]) + z_field(lp, ilat, 0);
4703 if (r_toa[ilat] < r_toa_min) {
4704 r_toa_min = r_toa[ilat];
4706 if (r_toa[ilat] > r_toa_max) {
4707 r_toa_max = r_toa[ilat];
4711 "The sensor is horizontally outside (or at the limit) of "
4712 "the model\natmosphere, but is at a radius smaller than "
4713 "the maximum value of\nthe top-of-the-atmosphere radii. "
4714 "This is not allowed. Make the\nmodel atmosphere larger "
4715 "to also cover the sensor position?")
4718 if (abs(rte_los[0]) <= 90) {
4719 ppath.
pos(0, 0) = rte_pos[0];
4720 ppath.
pos(0, 1) = rte_pos[1];
4721 ppath.
r[0] = r_e + rte_pos[0];
4722 ppath.
los(0, 0) = rte_los[0];
4725 out1 <<
" ------- WARNING -------: path is totally outside of "
4726 <<
"the model atmosphere\n";
4731 bool above =
false, ready =
false, failed =
false;
4741 if (islatin || ppath.
constant > r_toa_min) {
4752 while (!ready && !failed) {
4762 latt, lt, rt, r_p, rte_pos[1], rte_los[0], ppath.
constant);
4766 if (latt < lat_grid[0] || latt > lat_grid[llat]) {
4773 if (abs(lt - lt_old) <
LACC || ntries == 10000) {
4779 gridpos(gp_latt, lat_grid, latt);
4781 rt =
interp(itwt, r_toa, gp_latt);
4788 if (ntries == 10000) {
4794 "The path does not enter the model atmosphere. It "
4795 "reaches the\ntop of the atmosphere "
4796 "altitude around latitude ", latt,
" deg.")
4798 ppath.
pos(0, 0) = rte_pos[0];
4799 ppath.
pos(0, 1) = rte_pos[1];
4800 ppath.
r[0] = r_e + rte_pos[0];
4801 ppath.
los(0, 0) = rte_los[0];
4804 out1 <<
" ------- WARNING -------: path is totally outside "
4805 <<
"of the model atmosphere\n";
4808 ppath.
pos(0, 0) =
interp(itwt, z_field(lp, joker, 0), gp_latt);
4811 ppath.
pos(0, 1) =
interp(itwt, lat_grid, gp_latt);
4817 ppath.
gp_p[0].idx = lp - 1;
4818 ppath.
gp_p[0].fd[0] = 1;
4819 ppath.
gp_p[0].fd[1] = 0;
4825 if (cloudbox_on && cloudbox_limits[1] == lp) {
4827 if (fgp >= (Numeric)cloudbox_limits[2] &&
4828 fgp <= (Numeric)cloudbox_limits[3]) {
4840 const Index llat = lat_grid.nelem() - 1;
4841 const Index llon = lon_grid.nelem() - 1;
4844 Numeric lon2use = rte_pos[2];
4845 resolve_lon(lon2use, lon_grid[0], lon_grid[llon]);
4848 ppath.
end_pos[0] = rte_pos[0];
4849 ppath.
end_pos[1] = rte_pos[1];
4851 ppath.
end_los[0] = rte_los[0];
4852 ppath.
end_los[1] = rte_los[1];
4858 bool islatlonin =
false;
4860 Numeric z_toa = -99e99;
4862 if (rte_pos[1] >= lat_grid[0] && rte_pos[1] <= lat_grid[llat] &&
4863 lon2use >= lon_grid[0] && lon2use <= lon_grid[llon]) {
4865 gridpos(gp_lat, lat_grid, rte_pos[1]);
4866 gridpos(gp_lon, lon_grid, lon2use);
4868 z_toa =
interp(itw, z_field(lp, joker, joker), gp_lat, gp_lon);
4869 r_e =
refell2d(refellipsoid, lat_grid, gp_lat);
4871 r_e =
refell2r(refellipsoid, rte_pos[1]);
4875 if (islatlonin && rte_pos[0] < z_toa) {
4876 const Numeric z_s =
interp(itw, z_surface, gp_lat, gp_lon);
4880 "The ppath starting point is placed ", (z_s - rte_pos[0]) / 1e3,
4881 " km below the surface.")
4885 ppath.
r[0] = r_e + rte_pos[0];
4905 if (ppath.
pos(0, 0) <= z_s) {
4931 if (cloudbox_on && !ppath_inside_cloudbox_do) {
4935 if (fgp >= (Numeric)cloudbox_limits[0] &&
4936 fgp <= (Numeric)cloudbox_limits[1] &&
4937 fgl >= (Numeric)cloudbox_limits[2] &&
4938 fgl <= (Numeric)cloudbox_limits[3] &&
4939 fgo >= (Numeric)cloudbox_limits[4] &&
4940 fgo <= (Numeric)cloudbox_limits[5]) {
4951 fgp <= (Numeric)cloudbox_limits[1] &&
4952 fgl >= (Numeric)cloudbox_limits[2] &&
4953 fgl <= (Numeric)cloudbox_limits[3] &&
4954 fgo >= (Numeric)cloudbox_limits[4] &&
4955 fgo <= (Numeric)cloudbox_limits[5]);
4964 const Numeric r_p = r_e + rte_pos[0];
4968 Matrix r_toa(llat + 1, llon + 1);
4969 Numeric r_toa_min = 99e99, r_toa_max = -1;
4970 for (Index ilat = 0; ilat <= llat; ilat++) {
4971 const Numeric r_lat =
refell2r(refellipsoid, lat_grid[ilat]);
4972 for (Index ilon = 0; ilon <= llon; ilon++) {
4973 r_toa(ilat, ilon) = r_lat + z_field(lp, ilat, ilon);
4974 if (r_toa(ilat, ilon) < r_toa_min) {
4975 r_toa_min = r_toa(ilat, ilon);
4977 if (r_toa(ilat, ilon) > r_toa_max) {
4978 r_toa_max = r_toa(ilat, ilon);
4984 "r_p=",setprecision(18),r_p,
"\n"
4985 "r_toa_max=",r_toa_max,
"\n"
4986 "The sensor is horizontally outside of "
4987 "the model\natmosphere, but is at a radius smaller than "
4988 "the maximum value of\nthe top-of-the-atmosphere radii. "
4989 "This is not allowed. Make the\nmodel atmosphere larger "
4990 "to also cover the sensor position?")
4993 if (rte_los[0] <= 90) {
4994 ppath.
pos(0, 0) = rte_pos[0];
4995 ppath.
pos(0, 1) = rte_pos[1];
4996 ppath.
pos(0, 2) = rte_pos[2];
4997 ppath.
r[0] = r_e + rte_pos[0];
4998 ppath.
los(0, 0) = rte_los[0];
4999 ppath.
los(0, 1) = rte_los[1];
5002 out1 <<
" ------- WARNING -------: path is totally outside of "
5003 <<
"the model atmosphere\n";
5008 bool above =
false, ready =
false, failed =
false;
5009 Numeric rt = -1, latt, lont, lt, lt_old =
L_NOT_FOUND;
5018 if (islatlonin || ppath.
constant > r_toa_min) {
5026 Numeric x, y, z, dx, dy, dz;
5041 while (!ready && !failed) {
5068 if (latt < lat_grid[0] || latt > lat_grid[llat] ||
5069 lont < lon_grid[0] || lont > lon_grid[llon]) {
5076 if (abs(lt - lt_old) <
LACC) {
5082 gridpos(gp_latt, lat_grid, latt);
5083 gridpos(gp_lont, lon_grid, lont);
5085 rt =
interp(itwt, r_toa, gp_latt, gp_lont);
5091 "The path does not enter the model atmosphere. It\n"
5092 "reaches the top of the atmosphere altitude around:\n"
5093 " lat: ", latt,
" deg.\n lon: ", lont,
" deg.")
5095 ppath.
pos(0, 0) = rte_pos[0];
5096 ppath.
pos(0, 1) = rte_pos[1];
5097 ppath.
pos(0, 2) = rte_pos[2];
5099 ppath.
r[0] = r_e + rte_pos[0];
5100 ppath.
los(0, 0) = rte_los[0];
5101 ppath.
los(0, 1) = rte_los[1];
5104 out1 <<
" ------- WARNING -------: path is totally outside "
5105 <<
"of the model atmosphere\n";
5133 interp(itwt, z_field(lp, joker, joker), gp_latt, gp_lont);
5137 ppath.
gp_p[0].idx = lp - 1;
5138 ppath.
gp_p[0].fd[0] = 1;
5139 ppath.
gp_p[0].fd[1] = 0;
5146 if (cloudbox_on && cloudbox_limits[1] == lp) {
5149 if (fgp1 >= (Numeric)cloudbox_limits[2] &&
5150 fgp1 <= (Numeric)cloudbox_limits[3] &&
5151 fgp2 >= (Numeric)cloudbox_limits[4] &&
5152 fgp2 <= (Numeric)cloudbox_limits[5]) {
5164 const Agenda& ppath_step_agenda,
5165 const Index& atmosphere_dim,
5166 const Vector& p_grid,
5167 const Vector& lat_grid,
5168 const Vector& lon_grid,
5169 const Tensor3& z_field,
5170 const Vector& f_grid,
5171 const Vector& refellipsoid,
5172 const Matrix& z_surface,
5173 const Index& cloudbox_on,
5175 const Vector& rte_pos,
5176 const Vector& rte_los,
5177 const Numeric& ppath_lmax,
5178 const Numeric& ppath_lraytrace,
5179 const bool& ppath_inside_cloudbox_do,
5190 "The WSV *ppath_inside_cloudbox_do* can only be set "
5191 "to 1 if also *cloudbox_on* is 1.");
5215 ppath_inside_cloudbox_do,
5224 const Numeric end_lstep = ppath_step.
end_lstep;
5225 const Vector end_pos = ppath_step.
end_pos;
5226 const Vector end_los = ppath_step.
end_los;
5238 const Index imax_p = p_grid.nelem() - 1;
5239 const Index imax_lat = lat_grid.nelem() - 1;
5240 const Index imax_lon = lon_grid.nelem() - 1;
5260 const Index n = ppath_step.
np;
5266 "10 000 path points have been reached. Is this an infinite loop?");
5273 if (!ppath_inside_cloudbox_do) {
5278 (abs(ppath_step.
los(n - 1, 0)) < 90 &&
5284 if (atmosphere_dim == 2) {
5287 "The path exits the atmosphere through the lower "
5288 "latitude end face.\nThe exit point is at an altitude"
5289 " of ", ppath_step.
pos(n - 1, 0) / 1e3,
" km.")
5291 "The path exits the atmosphere through the upper "
5292 "latitude end face.\nThe exit point is at an altitude"
5293 " of ", ppath_step.
pos(n - 1, 0) / 1e3,
" km.")
5295 if (atmosphere_dim == 3) {
5299 "The path exits the atmosphere through the lower "
5300 "latitude end face.\nThe exit point is at an altitude"
5301 " of ", ppath_step.
pos(n - 1, 0) / 1e3,
" km.")
5304 "The path exits the atmosphere through the upper "
5305 "latitude end face.\nThe exit point is at an altitude"
5306 " of ", ppath_step.
pos(n - 1, 0) / 1e3,
" km.")
5312 ppath_step.
los(n - 1, 1) < 0 &&
5313 abs(ppath_step.
pos(n - 1, 1)) < 90) {
5315 if (lon_grid[imax_lon] - lon_grid[0] >= 360) {
5316 ppath_step.
pos(n - 1, 2) = ppath_step.
pos(n - 1, 2) + 360;
5318 ppath_step.
gp_lon[n - 1], lon_grid, ppath_step.
pos(n - 1, 2));
5321 "The path exits the atmosphere through the lower "
5322 "longitude end face.\nThe exit point is at an "
5323 "altitude of ", ppath_step.
pos(n - 1, 0) / 1e3,
" km.")
5326 ppath_step.
los(n - 1, 1) > 0 &&
5327 abs(ppath_step.
pos(n - 1, 1)) < 90) {
5329 if (lon_grid[imax_lon] - lon_grid[0] >= 360) {
5330 ppath_step.
pos(n - 1, 2) = ppath_step.
pos(n - 1, 2) - 360;
5332 ppath_step.
gp_lon[n - 1], lon_grid, ppath_step.
pos(n - 1, 2));
5335 "The path exits the atmosphere through the upper "
5336 "longitude end face.\nThe exit point is at an "
5337 "altitude of ", ppath_step.
pos(n - 1, 0) / 1e3,
" km.")
5346 if (ipos >= Numeric(cloudbox_limits[0]) &&
5347 ipos <= Numeric(cloudbox_limits[1])) {
5348 if (atmosphere_dim == 1) {
5353 if (ipos >= Numeric(cloudbox_limits[2]) &&
5354 ipos <= Numeric(cloudbox_limits[3])) {
5355 if (atmosphere_dim == 2) {
5360 if (ipos >= Numeric(cloudbox_limits[4]) &&
5361 ipos <= Numeric(cloudbox_limits[5])) {
5381 ARTS_ASSERT(ipos1 >= (Numeric)cloudbox_limits[0]);
5382 ARTS_ASSERT(ipos1 <= (Numeric)cloudbox_limits[1]);
5383 if (ipos1 <= (Numeric)cloudbox_limits[0] && ipos1 < ipos2) {
5387 else if (ipos1 >= Numeric(cloudbox_limits[1]) && ipos1 > ipos2) {
5391 else if (atmosphere_dim > 1) {
5395 ARTS_ASSERT(ipos1 >= (Numeric)cloudbox_limits[2]);
5396 ARTS_ASSERT(ipos1 <= (Numeric)cloudbox_limits[3]);
5397 if (ipos1 <= Numeric(cloudbox_limits[2]) && ipos1 < ipos2) {
5401 else if (ipos1 >= Numeric(cloudbox_limits[3]) && ipos1 > ipos2) {
5405 else if (atmosphere_dim > 2) {
5409 ARTS_ASSERT(ipos1 >= (Numeric)cloudbox_limits[4]);
5410 ARTS_ASSERT(ipos1 <= (Numeric)cloudbox_limits[5]);
5411 if (ipos1 <= Numeric(cloudbox_limits[4]) && ipos1 < ipos2) {
5415 else if (ipos1 >= Numeric(cloudbox_limits[5]) && ipos1 > ipos2) {
5433 ppath_array.push_back(ppath_step);
5441 const Index na = ppath_array.
nelem();
5466 for (Index i = 0; i < na; i++) {
5471 Index n = ppath_array[i].np;
5480 ppath.
r[Range(np, n - i1)] = ppath_array[i].r[Range(i1, n - i1)];
5481 ppath.
pos(Range(np, n - i1), joker) =
5482 ppath_array[i].pos(Range(i1, n - i1), joker);
5483 ppath.
los(Range(np, n - i1), joker) =
5484 ppath_array[i].los(Range(i1, n - i1), joker);
5485 ppath.
nreal[Range(np, n - i1)] = ppath_array[i].nreal[Range(i1, n - i1)];
5486 ppath.
ngroup[Range(np, n - i1)] =
5487 ppath_array[i].ngroup[Range(i1, n - i1)];
5488 ppath.
lstep[Range(np - i1, n - 1)] = ppath_array[i].lstep;
5491 for (Index j = i1; j < n; j++) {
5492 ppath.
gp_p[np + j - i1] = ppath_array[i].gp_p[j];
5494 if (atmosphere_dim >= 2) {
5495 for (Index j = i1; j < n; j++) {
5496 ppath.
gp_lat[np + j - i1] = ppath_array[i].gp_lat[j];
5499 if (atmosphere_dim == 3) {
5500 for (Index j = i1; j < n; j++) {
5501 ppath.
gp_lon[np + j - i1] = ppath_array[i].gp_lon[j];
Declarations for agendas.
This file contains the definition of Array.
base max(const Array< base > &x)
Max function.
base min(const Array< base > &x)
Min function.
Header file for helper functions for OpenMP.
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &f_grid, const Agenda &input_agenda)
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
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
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
void latlon_at_aa(Numeric &lat2, Numeric &lon2, const Numeric &lat1, const Numeric &lon1, const Numeric &aa, const Numeric &ddeg)
latlon_at_aa
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
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
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Index gridpos2gridrange(const GridPos &gp, const bool &upwards)
gridpos2gridrange
void gridpos_force_end_fd(GridPos &gp, const Index &n)
gridpos_force_end_fd
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gridpos_check_fd(GridPos &gp)
gridpos_check_fd
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
Numeric fractional_gp(const GridPos &gp)
fractional_gp
Numeric sign(const Numeric &x)
sign
Declarations having to do with the four output streams.
This file contains the definition of String, the ARTS string class.
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
int poly_root_solve(Matrix &roots, Vector &coeffs)
Contains the code to determine roots of polynomials.
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.
void ppath_end_1d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView za_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView z_field, ConstVectorView refellipsoid, const Index &ip, const Index &endface, const Numeric &ppc)
Internal help function for 1D path calculations.
void ppath_end_3d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView lon_v, ConstVectorView za_v, ConstVectorView aa_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, const Index &ip, const Index &ilat, const Index &ilon, const Index &endface, const Numeric &ppc)
Internal help function for 3D path calculations.
Numeric rsurf_at_lat(const Numeric &lat1, const Numeric &lat3, const Numeric &r1, const Numeric &r3, const Numeric &lat)
Determines the radius of a pressure level or the surface given the radius at the corners of a 2D grid...
bool is_los_downwards(const Numeric &za, const Numeric &tilt)
Determines if a line-of-sight is downwards compared to the angular tilt of the surface or a pressure ...
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.
void ppath_start_3d(Numeric &r_start, Numeric &lat_start, Numeric &lon_start, Numeric &za_start, Numeric &aa_start, Index &ip, Index &ilat, Index &ilon, Numeric &lat1, Numeric &lat3, Numeric &lon5, Numeric &lon6, Numeric &r15a, Numeric &r35a, Numeric &r36a, Numeric &r16a, Numeric &r15b, Numeric &r35b, Numeric &r36b, Numeric &r16b, Numeric &rsurface15, Numeric &rsurface35, Numeric &rsurface36, Numeric &rsurface16, Ppath &ppath, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface)
Internal help function for 3D path calculations.
const Numeric LON_NOT_FOUND
void zaaa2cart(Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &za, const Numeric &aa)
Converts zenith and azimuth angles to a cartesian unit vector.
void do_gridrange_1d(Vector &r_v, Vector &lat_v, Vector &za_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc, const Numeric &lmax, const Numeric &ra, const Numeric &rb, const Numeric &rsurface)
Calculates the geometrical path through a 1D grid range.
Numeric geompath_r_at_l(const Numeric &ppc, const Numeric &l)
Calculates the radius for a distance from the tangent point.
void raytrace_3d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &lon_array, Array< Numeric > &za_array, Array< Numeric > &aa_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView refellipsoid, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &rsurface15, const Numeric &rsurface35, const Numeric &rsurface36, const Numeric &rsurface16, const Numeric &r15a, const Numeric &r35a, const Numeric &r36a, const Numeric &r16a, const Numeric &r15b, const Numeric &r35b, const Numeric &r36b, const Numeric &r16b, Numeric r, Numeric lat, Numeric lon, Numeric za, Numeric aa)
Performs ray tracing for 3D with linear steps.
void ppath_start_1d(Numeric &r_start, Numeric &lat_start, Numeric &za_start, Index &ip, const Ppath &ppath)
Internal help function for 1D path calculations.
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
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.
const Numeric LAT_NOT_FOUND
void raytrace_1d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &za_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &rsurface, const Numeric &r1, const Numeric &r3, Numeric r, Numeric lat, Numeric za)
Performs ray tracing for 1D with linear steps.
void plevel_crossing_2d(Numeric &r, Numeric &lat, Numeric &l, const Numeric &r_start0, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc, const Numeric &lat1, const Numeric &lat3, const Numeric &r1, const Numeric &r3, const bool &above)
Handles the crossing with a geometric ppaths step and a atmospheric grid box level for 2D.
void do_gridcell_3d_byltest(Vector &r_v, Vector &lat_v, Vector &lon_v, Vector &za_v, Vector &aa_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start0, const Numeric &lon_start0, const Numeric &za_start, const Numeric &aa_start, const Numeric &l_start, const Index &icall, const Numeric &ppc, const Numeric &lmax, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15a, const Numeric &r35a, const Numeric &r36a, const Numeric &r16a, const Numeric &r15b, const Numeric &r35b, const Numeric &r36b, const Numeric &r16b, const Numeric &rsurface15, const Numeric &rsurface35, const Numeric &rsurface36, const Numeric &rsurface16)
See ATD for a description of the algorithm.
constexpr Numeric DEG2RAD
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Index first_pos_before_altitude(const Ppath &p, const Numeric &alt)
Determines ppath position just below an altitude.
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.
Numeric geompath_r_at_lat(const Numeric &ppc, const Numeric &lat0, const Numeric &za0, const Numeric &lat)
Calculates the radius for a given latitude.
void ppath_append(Ppath &ppath1, const Ppath &ppath2)
Combines two Ppath structures.
void rotationmat3D(Matrix &R, ConstVectorView vrot, const Numeric &a)
Creates a 3D rotation matrix.
Numeric geometrical_ppc(const Numeric &r, const Numeric &za)
Calculates the propagation path constant for pure geometrical calculations.
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.
void zaaa2enu(Numeric &de, Numeric &dn, Numeric &du, const Numeric &za, const Numeric &aa)
Converts zenith and azimuth angles to ENU unit vector.
void ppath_copy(Ppath &ppath1, const Ppath &ppath2, const Index &ncopy)
Copy the content in ppath2 to ppath1.
void enu2zaaa(Numeric &za, Numeric &aa, const Numeric &de, const Numeric &dn, const Numeric &du)
Converts ENU unit vector vector to zenith and azimuth.
void r_crossing_3d(Numeric &lat, Numeric &lon, Numeric &l, const Numeric &r_hit, const Numeric &r_start, const Numeric &lat_start, const Numeric &lon_start, const Numeric &za_start, const Numeric &ppc, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
Calculates where a 3D LOS crosses the specified radius.
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.
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.
void cart2zaaa(Numeric &za, Numeric &aa, const Numeric &dx, const Numeric &dy, const Numeric &dz)
Converts a cartesian directional vector to zenith and azimuth.
constexpr Numeric RAD2DEG
void ppath_end_2d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView za_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, const Index &ip, const Index &ilat, const Index &endface, const Numeric &ppc)
Internal help function for 2D path calculations.
Numeric rslope_crossing3d(const Numeric &rp, const Numeric &za, const Numeric &r0, Numeric c1, Numeric c2)
3D version of rslope_crossing2d.
void find_tanpoint(Index &it, const Ppath &ppath)
Identifies the tangent point of a propagation path.
void ppath_set_background(Ppath &ppath, const Index &case_nr)
Sets the background field of a Ppath structure.
void plevel_slope_3d(Numeric &c1, Numeric &c2, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15, const Numeric &r35, const Numeric &r36, const Numeric &r16, const Numeric &lat, const Numeric &lon, const Numeric &aa)
void r_crossing_2d(Numeric &lat, Numeric &l, const Numeric &r_hit, const Numeric &r_start, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc)
Calculates where a 2D LOS crosses the specified radius.
void resolve_lon(Numeric &lon, const Numeric &lon5, const Numeric &lon6)
Resolves which longitude angle that shall be used.
Numeric refraction_ppc(const Numeric &r, const Numeric &za, const Numeric &refr_index_air)
Calculates the propagation path constant for cases with refraction.
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.
void do_gridcell_2d_byltest(Vector &r_v, Vector &lat_v, Vector &za_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start0, const Numeric &za_start, const Numeric &l_start, const Index &icall, const Numeric &ppc, const Numeric &lmax, const Numeric &lat1, const Numeric &lat3, const Numeric &r1a, const Numeric &r3a, const Numeric &r3b, const Numeric &r1b, const Numeric &rsurface1, const Numeric &rsurface3)
Works as do_gridcell_3d_byltest, but downscaled to 2D.
Numeric geompath_r_at_za(const Numeric &ppc, const Numeric &za)
Calculates the zenith angle for a given radius along a geometrical propagation path.
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.
const Numeric R_NOT_FOUND
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.
Numeric rsurf_at_latlon(const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15, const Numeric &r35, const Numeric &r36, const Numeric &r16, const Numeric &lat, const Numeric &lon)
Determines the radius of a pressure level or the surface given the radius at the corners of a 3D grid...
void ppath_init_structure(Ppath &ppath, const Index &atmosphere_dim, const Index &np)
Initiates a Ppath structure to hold the given number of points.
Numeric plevel_angletilt(const Numeric &r, const Numeric &c1)
Calculates the angular tilt of the surface or a pressure level.
const Numeric L_NOT_FOUND
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.
void ppath_start_stepping(Ppath &ppath, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const bool &ppath_inside_cloudbox_do, ConstVectorView rte_pos, ConstVectorView rte_los, const Verbosity &verbosity)
Initiates a Ppath structure for calculation of a path with ppath_step.
void ppath_start_2d(Numeric &r_start, Numeric &lat_start, Numeric &za_start, Index &ip, Index &ilat, Numeric &lat1, Numeric &lat3, Numeric &r1a, Numeric &r3a, Numeric &r3b, Numeric &r1b, Numeric &rsurface1, Numeric &rsurface3, Ppath &ppath, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, ConstVectorView z_surface)
Internal help function for 2D path calculations.
Numeric geompath_l_at_r(const Numeric &ppc, const Numeric &r)
Calculates the length from the tangent point for the given radius.
void geompath_from_r1_to_r2(Vector &r, Vector &lat, Vector &za, Numeric &lstep, const Numeric &ppc, const Numeric &r1, const Numeric &lat1, const Numeric &za1, const Numeric &r2, const bool &tanpoint, const Numeric &lmax)
Determines radii, latitudes and zenith angles between two points of a propagation path.
void plevel_slope_2d(Numeric &c1, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstVectorView z_surf, const GridPos &gp, const Numeric &za)
Calculates the radial slope of the surface or a pressure level for 2D.
void raytrace_2d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &za_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &lat1, const Numeric &lat3, const Numeric &rsurface1, const Numeric &rsurface3, const Numeric &r1a, const Numeric &r3a, const Numeric &r3b, const Numeric &r1b, Numeric r, Numeric lat, Numeric za)
Performs ray tracing for 2D with linear steps.
Numeric rslope_crossing2d(const Numeric &rp, const Numeric &za, const Numeric &r0, Numeric c1)
Calculates the angular distance to a crossing with a level having a radial slope.
Propagation path structure and functions.
constexpr Numeric POLELAT
Size of north and south poles.
constexpr Numeric ANGTOL
Width of zenith and nadir directions.
void refr_gradients_3d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, Numeric &dndlat, Numeric &dndlon, 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)
refr_gradients_3d
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)
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
void refr_gradients_2d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, Numeric &dndlat, 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)
refr_gradients_2d
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
void refr_gradients_1d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, 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)
refr_gradients_1d
void adjust_los(VectorView los, const Index &atmosphere_dim)
Ensures that the zenith and azimuth angles of a line-of-sight vector are inside defined ranges.
Declaration of functions in rte.cc.
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.
The structure to describe a propagation path and releated quantities.
Matrix los
Line-of-sight at each ppath point.
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
String background
Radiative background.
Index np
Number of points describing the ppath.
Matrix pos
The distance between start pos and the last position in pos.
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Numeric end_lstep
The distance between end pos and the first position in pos.
Vector start_pos
Start position.
Vector lstep
The length between ppath points.
Numeric start_lstep
Length between sensor and atmospheric boundary.
Numeric constant
The propagation path constant (only used for 1D)
Vector r
Radius of each ppath point.
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Vector ngroup
The group index of refraction.
Index dim
Atmospheric dimensionality.
Vector end_pos
End position.
Vector start_los
Start line-of-sight.
Vector nreal
The real part of the refractive index at each path position.
Vector end_los
End line-of-sight.