63 r = sqrt(x * x + z * z);
66 const Numeric absza = abs(za0);
67 if (absza < ANGTOL || absza > 180 -
ANGTOL) {
74 lat = lat - 360.0 * Numeric(round((lat - lat0) / 360.0));
105 const Numeric& za0) {
106 r = sqrt(x * x + z * z);
109 const Numeric absza = abs(za0);
110 if (absza < ANGTOL || absza > 180 -
ANGTOL) {
118 const Numeric latrad =
DEG2RAD * lat;
119 const Numeric coslat = cos(latrad);
120 const Numeric sinlat = sin(latrad);
121 const Numeric dr = coslat * dx + sinlat * dz;
133 if (std::isnan(za)) {
139 if (std::isnan(za)) {
169 const Numeric& lat2) {
172 Numeric x1, z1, x2, z2;
176 const Numeric dx = x2 - x1;
177 const Numeric dz = z2 - z1;
178 l = sqrt(dx * dx + dz * dz);
270 const Numeric
a = dx * dx + dz * dz;
271 const Numeric
b = 2 * (dx * (xl - xc) + dz * (zl - zc));
273 xc * xc + zc * zc + xl * xl + zl * zl - 2 * (xc * xl + zc * zl) - r * r;
275 Numeric
d =
b *
b - 4 *
a *
c;
278 const Numeric a2 = 2 *
a;
279 const Numeric b2 = -
b / a2;
280 const Numeric e = sqrt(
d) / a2;
282 const Numeric l1 = b2 + e;
283 const Numeric l2 = b2 - e;
314void pol2cart(Numeric& x, Numeric& z,
const Numeric& r,
const Numeric& lat) {
317 const Numeric latrad =
DEG2RAD * lat;
348 const Numeric latrad =
DEG2RAD * lat;
349 const Numeric zarad =
DEG2RAD * za;
351 const Numeric coslat = cos(latrad);
352 const Numeric sinlat = sin(latrad);
353 const Numeric cosza = cos(zarad);
354 const Numeric sinza = sin(zarad);
360 const Numeric dr = cosza;
361 const Numeric dlat = sinza;
363 dx = coslat * dr - sinlat * dlat;
364 dz = sinlat * dr + coslat * dlat;
422 const Numeric& aa0) {
424 r = sqrt(x * x + y * y + z * z);
427 if (za0 < ANGTOL || za0 > 180 -
ANGTOL) {
438 bool ns_case =
false;
439 bool lon_flip =
false;
447 if (abs(abs(lon - lon0) - 180) < 5) {
459 const Numeric latrad =
DEG2RAD * lat;
460 const Numeric lonrad =
DEG2RAD * lon;
461 const Numeric coslat = cos(latrad);
462 const Numeric sinlat = sin(latrad);
463 const Numeric coslon = cos(lonrad);
464 const Numeric sinlon = sin(lonrad);
472 if (std::isnan(za)) {
479 sqrt(pow(x - x0, 2.0) + pow(y - y0, 2.0) + pow(z - z0, 2.0));
480 const Numeric r0 = sqrt(x0 * x0 + y0 * y0 + z0 * z0);
507 const Numeric dlat = -sinlat * coslon / r * dx -
508 sinlat * sinlon / r * dy + coslat / r * dz;
509 const Numeric dlon = -sinlon / coslat / r * dx + coslon / coslat / r * dy;
513 if (std::isnan(aa)) {
519 }
else if (dlon < 0) {
539 const Numeric latrad =
DEG2RAD * lat;
540 const Numeric lonrad =
DEG2RAD * lon;
541 const Numeric coslat = cos(latrad);
542 const Numeric sinlat = sin(latrad);
543 const Numeric coslon = cos(lonrad);
544 const Numeric sinlon = sin(lonrad);
546 const Numeric dr = coslat*coslon*dx + sinlat*dz + coslat*sinlon*dy;
547 const Numeric dlat = -sinlat*coslon/r*dx + coslat/r*dz - sinlat*sinlon/r*dy;
548 const Numeric dlon = -sinlon/coslat/r*dx + coslon/coslat/r*dy;
551 aa =
RAD2DEG * acos( r * dlat / sin( za ) );
555 if (std::isnan(aa)) {
560 }
else if (dlon < 0) {
594 const Numeric& aa0) {
595 r = sqrt(x * x + y * y + z * z);
598 if (za0 < ANGTOL || za0 > 180 -
ANGTOL) {
612 if (abs(lon - lon0) < 1) {
630 r = sqrt(x * x + y * y + z * z);
658 const Numeric& lon2) {
659 Numeric x1, y1, z1, x2, y2, z2;
660 sph2cart(x1, y1, z1, r1, lat1, lon1);
661 sph2cart(x2, y2, z2, r2, lat2, lon2);
663 const Numeric dx = x2 - x1;
664 const Numeric dy = y2 - y1;
665 const Numeric dz = z2 - z1;
666 l = sqrt(dx * dx + dy * dy + dz * dz);
698 const Numeric& ppc) {
702 Numeric x, y, z, dx, dy, dz;
704 poslos2cart(x, y, z, dx, dy, dz, r, lat, lon, za, aa);
706 l_tan = sqrt(r * r - ppc * ppc);
858 const Vector& refellipsoid,
868 if (refellipsoid[1] < 1e-7) {
869 const Numeric p = x*dx + y*dy + z*dz;
870 const Numeric pp = p*p;
871 const Numeric q = x*x + y*y + z*z - refellipsoid[0]*refellipsoid[0];
875 const Numeric sq = sqrt(pp - q);
886 const Numeric
a = refellipsoid[0];
887 const Numeric
b = refellipsoid[0];
888 const Numeric
c = refellipsoid[0] * sqrt(1-refellipsoid[1]*refellipsoid[1]);
889 const Numeric a2 =
a*
a;
890 const Numeric b2 =
b*
b;
891 const Numeric c2 =
c*
c;
892 const Numeric x2 = x*x;
893 const Numeric y2 = y*y;
894 const Numeric z2 = z*z;
895 const Numeric dx2 = dx*dx;
896 const Numeric dy2 = dy*dy;
897 const Numeric dz2 = dz*dz;
898 const Numeric rad = a2*b2*dz2 + a2*c2*dy2 - a2*dy2*z2 + 2*a2*dy*dz*y*z -
899 a2*dz2*y2 + b2*c2*dx2 - b2*dx2*z2 + 2*b2*dx*dz*x*z -
900 b2*dz2*x2 - c2*dx2*y2 + 2*c2*dx*dy*x*y - c2*dy2*x2;
904 const Numeric val = -a2*b2*dz*z - a2*c2*dy*y - b2*c2*dx*x;
905 const Numeric mag = a2*b2*dz2 + a2*c2*dy2 + b2*c2*dx2;
906 const Numeric abc =
a*
b*
c*sqrt(rad);
908 l = (val - abc) / mag;
910 l = (val + abc) / mag;
949 const Numeric
a = dx * dx + dy * dy + dz * dz;
950 const Numeric
b = 2 * (dx * (xl - xc) + dy * (yl - yc) + dz * (zl - zc));
951 const Numeric
c = xc * xc + yc * yc + zc * zc + xl * xl + yl * yl + zl * zl -
952 2 * (xc * xl + yc * yl + zc * zl) - r * r;
954 Numeric
d =
b *
b - 4 *
a *
c;
957 const Numeric a2 = 2 *
a;
958 const Numeric b2 = -
b / a2;
959 const Numeric e = sqrt(
d) / a2;
961 const Numeric l1 = b2 + e;
962 const Numeric l2 = b2 - e;
1001 const Numeric& ddeg) {
1006 const Numeric dang =
DEG2RAD * ddeg;
1007 const Numeric cosdang = cos(dang);
1008 const Numeric sindang = sin(dang);
1009 const Numeric latrad =
DEG2RAD * lat1;
1010 const Numeric coslat = cos(latrad);
1011 const Numeric sinlat = sin(latrad);
1012 const Numeric aarad =
DEG2RAD * aa;
1014 lat2 = sinlat * cosdang + coslat * sindang * cos(aarad);
1015 lon2 = lon1 +
RAD2DEG * (atan2(sin(aarad) * sindang * coslat,
1016 cosdang - sinlat * lat2));
1045 const Numeric& lat1,
1046 const Numeric& lon1,
1052 const Numeric& z2) {
1053 Numeric dx = x2 - x1, dy = y2 - y1, dz = z2 - z1;
1054 const Numeric ldxyz = sqrt(dx * dx + dy * dy + dz * dz);
1060 const Numeric latrad =
DEG2RAD * lat1;
1061 const Numeric lonrad =
DEG2RAD * lon1;
1062 const Numeric coslat = cos(latrad);
1063 const Numeric sinlat = sin(latrad);
1064 const Numeric coslon = cos(lonrad);
1065 const Numeric sinlon = sin(lonrad);
1067 const Numeric dr = coslat * coslon * dx + coslat * sinlon * dy + sinlat * dz;
1068 const Numeric dlat =
1069 -sinlat * coslon / r1 * dx - sinlat * sinlon / r1 * dy + coslat / r1 * dz;
1070 const Numeric dlon = -sinlon / coslat / r1 * dx + coslon / coslat / r1 * dy;
1074 if (std::isnan(aa)) {
1080 }
else if (dlon < 0) {
1120 const Numeric& aa) {
1130 const Numeric s =
sign(lat);
1143 const Numeric latrad =
DEG2RAD * lat;
1144 const Numeric lonrad =
DEG2RAD * lon;
1145 const Numeric zarad =
DEG2RAD * za;
1146 const Numeric aarad =
DEG2RAD * aa;
1148 const Numeric coslat = cos(latrad);
1149 const Numeric sinlat = sin(latrad);
1150 const Numeric coslon = cos(lonrad);
1151 const Numeric sinlon = sin(lonrad);
1152 const Numeric cosza = cos(zarad);
1153 const Numeric sinza = sin(zarad);
1154 const Numeric cosaa = cos(aarad);
1155 const Numeric sinaa = sin(aarad);
1163 const Numeric dr = cosza;
1164 const Numeric dlat = sinza * cosaa;
1165 const Numeric dlon = sinza * sinaa / coslat;
1167 dx = coslat * coslon * dr - sinlat * coslon * dlat - coslat * sinlon * dlon;
1168 dz = sinlat * dr + coslat * dlat;
1169 dy = coslat * sinlon * dr - sinlat * sinlon * dlat + coslat * coslon * dlon;
1195 ConstVectorView refellipsoid,
1196 ConstVectorView lat_grid,
1197 ConstVectorView lon_grid,
1198 ConstVectorView rte_pos) {
1199 if (atmosphere_dim == 1) {
1200 return refellipsoid[0];
1206 if (rte_pos[1] < lat_grid[0] || rte_pos[1] >
last(lat_grid)) {
1208 }
else if (atmosphere_dim == 3) {
1210 if (rte_pos[2] < lon_grid[0] || rte_pos[2] >
last(lon_grid)) {
1217 gridpos(gp_lat, lat_grid, rte_pos[1]);
1218 return refell2d(refellipsoid, lat_grid, gp_lat);
1220 return refell2r(refellipsoid, rte_pos[1]);
1248Numeric
refell2r(ConstVectorView refellipsoid,
const Numeric& lat) {
1254 if (refellipsoid[1] < 1e-7)
1256 return refellipsoid[0];
1260 const Numeric
c = 1 - refellipsoid[1] * refellipsoid[1];
1261 const Numeric
b = refellipsoid[0] * sqrt(
c);
1263 const Numeric ct = cos(
v);
1264 const Numeric st = sin(
v);
1266 return b / sqrt(
c * ct * ct + st * st);
1288 ConstVectorView lat_grid,
1292 else if (gp.
fd[0] == 1)
1293 return refell2r(refellipsoid, lat_grid[gp.
idx + 1]);
1295 return gp.
fd[1] *
refell2r(refellipsoid, lat_grid[gp.
idx]) +
1319 const Numeric& lon1,
1320 const Numeric& lat2,
1321 const Numeric& lon2) {
1323 const Numeric slat = sin(
DEG2RAD * (lat2 - lat1) / 2.0);
1324 const Numeric slon = sin(
DEG2RAD * (lon2 - lon1) / 2.0);
1326 slat * slat + cos(
DEG2RAD * lat1) * cos(
DEG2RAD * lat2) * slon * slon;
1328 return RAD2DEG * 2 * atan2(sqrt(
a), sqrt(1 -
a));
1354 const Numeric& lon) {
1359 const Numeric latrad =
DEG2RAD * lat;
1360 const Numeric lonrad =
DEG2RAD * lon;
1362 x = r * cos(latrad);
1363 y = x * sin(lonrad);
1364 x = x * cos(lonrad);
1365 z = r * sin(latrad);
1392 ConstVectorView longrid_in,
1393 const Numeric lon) {
1394 longrid_out = longrid_in;
1395 if (longrid_in[longrid_in.nelem() - 1] >= lon + 360.)
1396 longrid_out += -360.;
1397 else if (longrid_in[0] <= lon - 360.)
1398 longrid_out += 360.;
1417 "Latitude values < -180.0 are not supported.");
1419 "Latitude values > 180.0 are not supported.");
1433 while (lon > 360.0) {
1463 const Vector& refellipsoid ) {
1465 if (refellipsoid[1] < 1e-7)
1468 h = r - refellipsoid[0];
1474 const Numeric sq = sqrt(x*x+y*y);
1475 Numeric B0 = atan2(z,sq);
1476 Numeric B = B0-1, N;
1477 const Numeric e2 = refellipsoid[1]*refellipsoid[1];
1479 while (abs(B-B0)>1e-10) {
1480 N = refellipsoid[0] / sqrt(1-e2*sin(B0)*sin(B0));
1481 h = sq / cos(B0) - N;
1483 B0 = atan((z/sq) * 1/(1-e2*N/(N+h)));
1510 const Vector& refellipsoid ) {
1512 if (refellipsoid[1] < 1e-7)
1513 {
sph2cart(x, y, z, h+refellipsoid[0], lat, lon); }
1516 const Numeric
a = refellipsoid[0];
1517 const Numeric e2 = refellipsoid[1]*refellipsoid[1];
1518 const Numeric sinlat = sin(
DEG2RAD*lat);
1519 const Numeric coslat = cos(
DEG2RAD*lat);
1520 const Numeric N =
a / sqrt(1 - e2*sinlat*sinlat);
1522 x = (N + h) * coslat * cos(
DEG2RAD*lon);
1523 y = (N + h) * coslat * sin(
DEG2RAD*lon);
1524 z = (N*(1 - e2) + h) * sinlat;
1560 const Vector& refellipsoid ) {
1569 const Numeric s =
sign(lat);
1573 z = s * (h + refellipsoid[0]*sqrt(1-refellipsoid[1]*refellipsoid[1]));
1582 const Numeric coslat = cos(
DEG2RAD * lat);
1583 const Numeric sinlat = sin(
DEG2RAD * lat);
1584 const Numeric coslon = cos(
DEG2RAD * lon);
1585 const Numeric sinlon = sin(
DEG2RAD * lon);
1594 dx = -sinlon*de - sinlat*coslon*dn + coslat*coslon*du;
1595 dy = coslon*de - sinlat*sinlon*dn + coslat*sinlon*du;
1596 dz = coslat* dn + sinlat* du;
1632 const Vector& refellipsoid ) {
1636 const Numeric latrad =
DEG2RAD * lat;
1637 const Numeric lonrad =
DEG2RAD * lon;
1638 const Numeric coslat = cos(latrad);
1639 const Numeric sinlat = sin(latrad);
1640 const Numeric coslon = cos(lonrad);
1641 const Numeric sinlon = sin(lonrad);
1645 const Numeric de = -sinlon*dx + coslon*dy;
1646 const Numeric dn = -sinlat*coslon*dx - sinlat*sinlon*dy + coslat*dz;
1647 const Numeric du = coslat*coslon*dx + coslat*sinlon*dy + sinlat*dz;
base min(const Array< base > &x)
Min function.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR_IF(condition,...)
void distance3D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &r2, const Numeric &lat2, const Numeric &lon2)
distance3D
void distance2D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
distance2D
void cart2geodeticposlos(Numeric &h, Numeric &lat, Numeric &lon, Numeric &za, Numeric &aa, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz, const Vector &refellipsoid)
cart2geodeticposlos
void sph2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &r, const Numeric &lat, const Numeric &lon)
sph2cart
void los2xyz(Numeric &za, Numeric &aa, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &x1, const Numeric &y1, const Numeric &z1, const Numeric &x2, const Numeric &y2, const Numeric &z2)
los2xyz
void cart2poslos_plain(Numeric &r, Numeric &lat, Numeric &lon, Numeric &za, Numeric &aa, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
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
constexpr Numeric DEG2RAD
void pol2cart(Numeric &x, Numeric &z, const Numeric &r, const Numeric &lat)
pol2cart
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
void geodeticposlos2cart(Numeric &x, Numeric &y, Numeric &z, Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &h, const Numeric &lat, const Numeric &lon, const Numeric &za, const Numeric &aa, const Vector &refellipsoid)
geodeticposlos2cart
void cart2geodetic(Numeric &h, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z, const Vector &refellipsoid)
Conversion from cartesian to geodetic coordinates.
void cycle_lat_lon(Numeric &lat, Numeric &lon)
Cyclic latitude longitude coordinates.
Numeric pos2refell_r(const Index &atmosphere_dim, ConstVectorView refellipsoid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView rte_pos)
pos2refell_r
void line_circle_intersect(Numeric &x, Numeric &z, const Numeric &xl, const Numeric &zl, const Numeric &dx, const Numeric &dz, const Numeric &xc, const Numeric &zc, const Numeric &r)
geomtanpoint2d
void line_refellipsoid_intersect(Numeric &l, const Vector &refellipsoid, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
geomtanpoint
constexpr Numeric RAD2DEG
Numeric sphdist(const Numeric &lat1, const Numeric &lon1, const Numeric &lat2, const Numeric &lon2)
sphdist
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
void lon_shiftgrid(Vector &longrid_out, ConstVectorView longrid_in, const Numeric lon)
lon_shiftgrid
void geodetic2cart(Numeric &x, Numeric &y, Numeric &z, const Numeric &h, const Numeric &lat, const Numeric &lon, const Vector &refellipsoid)
Conversion from geodetic to cartesian coordinates.
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
void cart2sph_plain(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z)
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
void geompath_tanpos_3d(Numeric &r_tan, Numeric &lat_tan, Numeric &lon_tan, Numeric &l_tan, const Numeric &r, const Numeric &lat, const Numeric &lon, const Numeric &za, const Numeric &aa, const Numeric &ppc)
geompath_tanpos_3d
void line_sphere_intersect(Numeric &x, Numeric &y, Numeric &z, const Numeric &xl, const Numeric &yl, const Numeric &zl, const Numeric &dx, const Numeric &dy, const Numeric &dz, const Numeric &xc, const Numeric &yc, const Numeric &zc, const Numeric &r)
line_sphere_intersect
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Numeric last(ConstVectorView x)
last
Numeric sign(const Numeric &x)
sign
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
void zaaa2enu(Numeric &de, Numeric &dn, Numeric &du, const Numeric &za, const Numeric &aa)
Converts zenith and azimuth angles to ENU unit vector.
void enu2zaaa(Numeric &za, Numeric &aa, const Numeric &de, const Numeric &dn, const Numeric &du)
Converts ENU unit vector vector to zenith and azimuth.
Numeric geompath_l_at_r(const Numeric &ppc, const Numeric &r)
Calculates the length from the tangent point for the given radius.
Propagation path structure and functions.
constexpr Numeric POLELAT
Size of north and south poles.
constexpr Numeric ANGTOL
Width of zenith and nadir directions.
Structure to store a grid position.
std::array< Numeric, 2 > fd