Peano
Loading...
Searching...
No Matches
ProjectionAEA.cpph
Go to the documentation of this file.
1// This file is part of the ExaHyPE2 project. For conditions of distribution and
2// use, please see the copyright notice at www.peano-framework.org
3
4#include <cmath>
5#include <utility>
6
7template <auto EllipsoidA, auto EllipsoidE, auto StdParallel1Deg, auto StdParallel2Deg, auto Lon0Deg, auto Lat0Deg>
8static inline auto projection::AEA::ellipsoid::forward(const auto lonDeg, const auto latDeg) {
9 constexpr auto a = EllipsoidA;
10 constexpr auto e = EllipsoidE;
11 constexpr auto stdParallel1 = StdParallel1Deg;
12 constexpr auto stdParallel2 = StdParallel2Deg;
13 constexpr auto lon0 = Lon0Deg;
14 constexpr auto lat0 = Lat0Deg; // Latitude of Origin (Grad)
15
16 constexpr auto toRad = M_PI / 180.0;
17
18 constexpr auto stdParallel1Rad = stdParallel1 * toRad;
19 constexpr auto stdParallel2Rad = stdParallel2 * toRad;
20 constexpr auto lat0Rad = lat0 * toRad;
21 constexpr auto lon0Rad = lon0 * toRad;
22
23 const auto lonRad = lonDeg * toRad;
24 const auto latRad = latDeg * toRad;
25
26 constexpr auto e2 = e * e;
27 constexpr auto oneMinusE2 = 1.0 - (e * e);
28
29 auto m = [](const auto phi) { return std::cos(phi) / std::sqrt(1.0 - (e2 * std::sin(phi) * std::sin(phi))); };
30
31 static const auto m1 = m(stdParallel1Rad);
32 static const auto m2 = m(stdParallel2Rad);
33
34 auto q = [](const auto phi) {
35 const auto sinPhi = std::sin(phi);
36 const auto sinPhi2 = sinPhi * sinPhi;
37 return oneMinusE2
38 * ((std::sin(phi) / (1.0 - (e2 * sinPhi2))) - ((1.0 / (2.0 * e)) * std::log((1.0 - (e * std::sin(phi))) / (1.0 + (e * std::sin(phi))))));
39 };
40
41 static const auto q0 = q(lat0Rad);
42 static const auto q1 = q(stdParallel1Rad);
43 static const auto q2 = q(stdParallel2Rad);
44
45 const auto qLat = q(latRad);
46
47 static const auto n = ((m1 * m1) - (m2 * m2)) / (q2 - q1);
48 static const auto C = (m1 * m1) + (n * q1);
49 static const auto p0 = a * std::sqrt(C - (n * q0)) / n;
50
51 const auto p = a * std::sqrt(C - (n * qLat)) / n;
52 const auto theta = n * (lonRad - lon0Rad);
53
54 const auto x = p * std::sin(theta);
55 const auto y = p0 - (p * std::cos(theta));
56
57 return std::pair{x, y};
58};
59
60template <
61 auto EllipsoidA,
62 auto EllipsoidE,
63 auto StdParallel1Deg,
64 auto StdParallel2Deg,
65 auto Lon0Deg,
66 auto Lat0Deg,
67 auto MaxIter,
68 auto Tolerance>
69static inline auto projection::AEA::ellipsoid::inverse(const auto x, const auto y) {
70 constexpr auto tol = Tolerance;
71
72 constexpr auto a = EllipsoidA; // semi-major axis
73 constexpr auto e = EllipsoidE; // eccentricity
74 constexpr auto stdParallel1 = StdParallel1Deg; // standard parallel 1 (deg)
75 constexpr auto stdParallel2 = StdParallel2Deg; // standard parallel 2 (deg)
76 constexpr auto lon0 = Lon0Deg; // central meridian (deg)
77 constexpr auto lat0 = Lat0Deg; // latitude of origin (deg) // Latitude of Origin (Grad)
78
79 constexpr auto toRad = M_PI / 180.0;
80
81 constexpr auto stdParallel1Rad = stdParallel1 * toRad;
82 constexpr auto stdParallel2Rad = stdParallel2 * toRad;
83 constexpr auto lat0Rad = lat0 * toRad;
84 constexpr auto lon0Rad = lon0 * toRad;
85
86 constexpr auto e2 = e * e;
87 constexpr auto oneMinusE2 = 1.0 - (e * e);
88
89 auto m = [](const auto phi) { return std::cos(phi) / std::sqrt(1.0 - (e2 * std::sin(phi) * std::sin(phi))); };
90
91 static const auto m1 = m(stdParallel1Rad);
92 static const auto m2 = m(stdParallel2Rad);
93
94 auto q = [](const auto phi) {
95 const auto sinPhi = std::sin(phi);
96 const auto sinPhi2 = sinPhi * sinPhi;
97 return oneMinusE2
98 * ((std::sin(phi) / (1.0 - (e2 * sinPhi2))) - ((1.0 / (2.0 * e)) * std::log((1.0 - (e * std::sin(phi))) / (1.0 + (e * std::sin(phi))))));
99 };
100
101 static const auto q0 = q(lat0Rad);
102 static const auto q1 = q(stdParallel1Rad);
103 static const auto q2 = q(stdParallel2Rad);
104
105 static const auto n = ((m1 * m1) - (m2 * m2)) / (q2 - q1);
106 static const auto C = (m1 * m1) + (n * q1);
107 static const auto p0 = a * std::sqrt(C - (n * q0)) / n;
108
109 const auto p = std::sqrt((x * x) + ((p0 - y) * (p0 - y)));
110 const auto qnew = (C - (p * p * n * n / (a * a))) / n;
111
112 auto phi = std::asin(qnew / 2);
113
114 for (int i = 0; i < static_cast<int>(MaxIter); ++i) {
115 const auto sinPhi = std::sin(phi);
116 const auto cosPhi = std::cos(phi);
117 const auto sinPhi2 = sinPhi * sinPhi;
118 const auto denom = 1.0 - (e2 * sinPhi2);
119
120 const auto num1 = qnew / oneMinusE2;
121 const auto num2 = sinPhi / denom;
122 const auto num3 = (1.0 / (2.0 * e)) * std::log((1.0 - e * sinPhi) / (1.0 + e * sinPhi));
123
124 const auto deltaPhi = (1.0 - (e2 * sinPhi2)) * (1.0 - (e2 * sinPhi2)) / (2.0 * cosPhi) * (num1 - num2 + num3);
125 const auto phiNew = phi + deltaPhi;
126
127 if (std::abs(phiNew - phi) < tol) {
128 break;
129 }
130 phi = phiNew;
131 }
132
133 const auto theta = std::atan2(x, p0 - y);
134 const auto lon = lon0Rad + (theta / n);
135
136 return std::pair{
137 lon / toRad,
138 phi / toRad,
139 };
140};
141
142template <auto Lon0Deg, auto Lat0Deg, auto StdParallel1Deg, auto StdParallel2Deg, auto EarthRadius>
143static inline auto projection::AEA::sphere::forward(const auto lon, const auto lat) {
144 constexpr auto lon0 = Lon0Deg;
145 constexpr auto lat0 = Lat0Deg;
146 constexpr auto stdParallel1 = StdParallel1Deg;
147 constexpr auto stdParallel2 = StdParallel2Deg;
148 constexpr auto R = EarthRadius;
149
150 constexpr auto toRad = M_PI / 180.0;
151
152 constexpr auto stdParallel1Rad = stdParallel1 * toRad;
153 constexpr auto stdParallel2Rad = stdParallel2 * toRad;
154 constexpr auto lat0Rad = lat0 * toRad;
155 constexpr auto lon0Rad = lon0 * toRad;
156
157 const auto lonRad = lon * toRad;
158 const auto latRad = lat * toRad;
159
160 static const auto n = (std::sin(stdParallel1Rad) + std::sin(stdParallel2Rad)) / 2;
161 static const auto C = (std::cos(stdParallel1Rad) * std::cos(stdParallel1Rad)) + (2 * n * std::sin(stdParallel1Rad));
162 static const auto p0 = R * std::sqrt(C - (2 * n * std::sin(lat0Rad))) / n;
163
164 const auto p = R * std::sqrt(C - (2 * n * std::sin(latRad))) / n;
165 const auto theta = n * (lonRad - lon0Rad);
166
167 const auto x = p * std::sin(theta);
168 const auto y = p0 - (p * std::cos(theta));
169
170 return std::pair{x, y};
171};
172
173template <auto Lon0Deg, auto Lat0Deg, auto StdParallel1Deg, auto StdParallel2Deg, auto EarthRadius>
174static inline auto projection::AEA::sphere::inverse(const auto x, const auto y) {
175 constexpr auto lon0 = Lon0Deg;
176 constexpr auto lat0 = Lat0Deg;
177 constexpr auto stdParallel1 = StdParallel1Deg;
178 constexpr auto stdParallel2 = StdParallel2Deg;
179 constexpr auto R = EarthRadius;
180
181 constexpr auto toRad = M_PI / 180.0;
182
183 constexpr auto stdParallel1Rad = stdParallel1 * toRad;
184 constexpr auto stdParallel2Rad = stdParallel2 * toRad;
185 constexpr auto lat0Rad = lat0 * toRad;
186 constexpr auto lon0Rad = lon0 * toRad;
187
188 static const auto n = (std::sin(stdParallel1Rad) + std::sin(stdParallel2Rad)) / 2;
189 static const auto C = (std::cos(stdParallel1Rad) * std::cos(stdParallel1Rad)) + (2 * n * std::sin(stdParallel1Rad));
190 static const auto p0 = R * std::sqrt(C - (2 * n * std::sin(lat0Rad))) / n;
191
192 const auto p = std::sqrt((x * x) + ((p0 - y) * (p0 - y)));
193 const auto theta = std::atan2(x, p0 - y);
194
195 const auto lat = std::asin((C - ((p * n / R) * (p * n / R))) / (2 * n));
196 const auto lon = lon0Rad + (theta / n);
197
198 return std::pair{lon / toRad, lat / toRad};
199};
csv_to_plot.y
list y
Definition csv_to_plot.py:98
projection::AEA::ellipsoid::forward
static auto forward(const auto lonDeg, const auto latDeg)
projection::AEA::ellipsoid::inverse
static auto inverse(const auto x, const auto y)
Inverts the Albers Equal-Area Conic projection for WGS_1984 and returns the geographic coordinates (l...
projection::AEA::sphere::forward
static auto forward(const auto lon, const auto lat)
Projects geographic coordinates (longitude, latitude) onto a planar surface using the Albers Equal-Ar...
projection::AEA::sphere::inverse
static auto inverse(const auto x, const auto y)
Inverts the spherical Albers Equal-Area Conic projection and returns the geographic coordinates (long...