ARTS 2.5.9 (git: 825fa5f2)
m_geodetic.cc
Go to the documentation of this file.
1/* Copyright (C) 2012
2 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3
4 This program is free software; you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by the
6 Free Software Foundation; either version 2, or (at your option) any
7 later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 USA. */
18
19/*===========================================================================
20 === File description
21 ===========================================================================*/
22
34/*===========================================================================
35 === External declarations
36 ===========================================================================*/
37
38#include <cmath>
39#include <stdexcept>
40#include "arts.h"
41#include "arts_conversions.h"
42#include "auto_md.h"
43#include "check_input.h"
44#include "geodetic.h"
45#include "matpackI.h"
46#include "messages.h"
47
49
50/*===========================================================================
51 === The functions (in alphabetical order)
52 ===========================================================================*/
53
54/* Workspace method: Doxygen documentation will be auto-generated */
55void refellipsoidForAzimuth(Vector& refellipsoid,
56 const Numeric& latitude,
57 const Numeric& azimuth,
58 const Verbosity&) {
59 ARTS_USER_ERROR_IF (refellipsoid.nelem() != 2,
60 "Input *refellispoid must be a vector of length 2*.");
61
62 if (refellipsoid[1] > 0) {
63 const Numeric e2 = refellipsoid[1] * refellipsoid[1];
64 const Numeric a = 1 - e2 * pow(sin(DEG2RAD * latitude), 2.0);
65
66 const Numeric rn = 1 / sqrt(a);
67 const Numeric rm = (1 - e2) * (rn / a);
68
69 const Numeric v = DEG2RAD * azimuth;
70
71 refellipsoid[0] =
72 refellipsoid[0] / (pow(cos(v), 2.0) / rm + pow(sin(v), 2.0) / rn);
73 refellipsoid[1] = 0;
74 }
75}
76
77
78
79/* Workspace method: Doxygen documentation will be auto-generated */
80void refellipsoidOrbitPlane(Vector& refellipsoid,
81 const Numeric& orbitinc,
82 const Verbosity&) {
83 ARTS_USER_ERROR_IF (refellipsoid.nelem() != 2,
84 "Input *refellispoid must be a vector of length 2*.");
85 chk_if_in_range("orbitinc", orbitinc, 0, 180);
86
87 // Radius at maximum latitude
88 const Numeric rp = refell2r(refellipsoid, orbitinc);
89
90 // New eccentricity
91 refellipsoid[1] = sqrt(1 - pow(rp / refellipsoid[0], 2.0));
92}
93
94
95
96/* Workspace method: Doxygen documentation will be auto-generated */
97void refellipsoidSet(Vector& refellipsoid,
98 const Numeric& re,
99 const Numeric& e,
100 const Verbosity&) {
101 refellipsoid.resize(2);
102
103 refellipsoid[0] = re;
104 refellipsoid[1] = e;
105}
106
107
108
109/* Workspace method: Doxygen documentation will be auto-generated */
111 Vector& rte_los,
112 const Matrix& sensor_pos_ecef,
113 const Matrix& sensor_los_ecef,
114 const Vector& refellipsoid,
115 const Verbosity& v) {
116 ARTS_USER_ERROR_IF (sensor_pos_ecef.nrows() != 1,
117 "For this WSM, *sensor_pos_ecef* can only have one row.");
118
119 Matrix sensor_pos, sensor_los;
120 sensor_poslosFromECEF(sensor_pos,
121 sensor_los,
122 sensor_pos_ecef,
123 sensor_los_ecef,
124 refellipsoid,
125 v );
126 rte_pos = sensor_pos(0,joker);
127 if(sensor_los.empty())
128 rte_los.resize(0);
129 else
130 rte_los = sensor_los(0,joker);
131}
132
133
134
135/* Workspace method: Doxygen documentation will be auto-generated */
137 Vector& rte_los,
138 const Matrix& sensor_pos_geodetic,
139 const Matrix& sensor_los_geodetic,
140 const Vector& refellipsoid,
141 const Verbosity& v) {
142 ARTS_USER_ERROR_IF (sensor_pos_geodetic.nrows() != 1,
143 "For this WSM, *sensor_pos_geodetic* can only have one row.");
144
145 Matrix sensor_pos, sensor_los;
146 sensor_poslosFromGeodetic(sensor_pos,
147 sensor_los,
148 sensor_pos_geodetic,
149 sensor_los_geodetic,
150 refellipsoid,
151 v );
152 rte_pos = sensor_pos(0,joker);
153 if(sensor_los.empty())
154 rte_los.resize(0);
155 else
156 rte_los = sensor_los(0,joker);
157}
158
159
160
161/* Workspace method: Doxygen documentation will be auto-generated */
163 Matrix& sensor_los,
164 const Matrix& sensor_pos_ecef,
165 const Matrix& sensor_los_ecef,
166 const Vector& refellipsoid,
167 const Verbosity&) {
168 const Index ncols = sensor_pos_ecef.ncols();
169 const Index nrows = sensor_pos_ecef.nrows();
170 ARTS_USER_ERROR_IF (ncols != 3,
171 "*sensor_pos_geodetic* must have three columns.");
172
173 sensor_pos.resize( nrows, ncols );
174
175 if (sensor_los_ecef.empty()) {
176 sensor_los = sensor_los_ecef;
177
178 for (Index i=0; i<nrows; i++) {
179 cart2sph_plain(sensor_pos(i,0),
180 sensor_pos(i,1),
181 sensor_pos(i,2),
182 sensor_pos_ecef(i,0),
183 sensor_pos_ecef(i,1),
184 sensor_pos_ecef(i,2));
185 sensor_pos(i,0) -= refell2r(refellipsoid,sensor_pos(i,1));
186 }
187 }
188 else {
189 const Index ncols2 = sensor_los_ecef.ncols();
190 ARTS_USER_ERROR_IF (ncols2 != 3,
191 "*sensor_los_ecef* must be empty or have "
192 "three columns.");
193 ARTS_USER_ERROR_IF (nrows != sensor_los_ecef.nrows(),
194 "*sensor_los_ecef* must be empty or have the "
195 "same number of rows as *sensor_pos_ecef*");
196
197 sensor_pos.resize( nrows, ncols );
198 sensor_los.resize( nrows, 2 );
199
200 for (Index i=0; i<nrows; i++) {
201 cart2poslos_plain(sensor_pos(i,0),
202 sensor_pos(i,1),
203 sensor_pos(i,2),
204 sensor_los(i,0),
205 sensor_los(i,1),
206 sensor_pos_ecef(i,0),
207 sensor_pos_ecef(i,1),
208 sensor_pos_ecef(i,2),
209 sensor_los_ecef(i,0),
210 sensor_los_ecef(i,1),
211 sensor_los_ecef(i,2));
212 sensor_pos(i,0) -= refell2r(refellipsoid,sensor_pos(i,1));
213 }
214 }
215}
216
217
218
219/* Workspace method: Doxygen documentation will be auto-generated */
221 Matrix& sensor_los,
222 const Matrix& sensor_pos_geodetic,
223 const Matrix& sensor_los_geodetic,
224 const Vector& refellipsoid,
225 const Verbosity&) {
226 const Index ncols = sensor_pos_geodetic.ncols();
227 const Index nrows = sensor_pos_geodetic.nrows();
228 ARTS_USER_ERROR_IF (ncols != 3,
229 "*sensor_pos_geodetic* must have three columns.");
230
231 // No conversion to do if geoid is spherical
232 if (refellipsoid[1] < 1e-7) {
233 sensor_pos = sensor_pos_geodetic;
234 sensor_los = sensor_los_geodetic;
235 }
236 else {
237 sensor_pos.resize( nrows, ncols );
238
239 if (sensor_los_geodetic.empty()) {
240 sensor_los = sensor_los_geodetic;
241 Numeric x, y, z;
242
243 for (Index i=0; i<nrows; i++) {
244 geodetic2cart(x, y, z,
245 sensor_pos_geodetic(i,0),
246 sensor_pos_geodetic(i,1),
247 sensor_pos_geodetic(i,2),
248 refellipsoid );
249 cart2sph_plain(sensor_pos(i,0),
250 sensor_pos(i,1),
251 sensor_pos(i,2),
252 x, y, z);
253 sensor_pos(i,0) -= refell2r(refellipsoid,sensor_pos(i,1));
254 }
255 } else {
256 const Index ncols2 = sensor_los_geodetic.ncols();
257 ARTS_USER_ERROR_IF (ncols2 != 2,
258 "*sensor_los_geodetic* must be empty or have "
259 "two columns.");
260 ARTS_USER_ERROR_IF (nrows != sensor_los_geodetic.nrows(),
261 "*sensor_los_geodetic* must be empty or have the "
262 "same number of rows as *sensor_pos_geodetic*");
263
264 sensor_pos.resize( nrows, ncols );
265 sensor_los.resize( nrows, ncols2 );
266 Numeric x, y, z, dx, dy, dz;
267
268 for (Index i=0; i<nrows; i++) {
269 geodeticposlos2cart(x, y, z, dx,dy, dz,
270 sensor_pos_geodetic(i,0),
271 sensor_pos_geodetic(i,1),
272 sensor_pos_geodetic(i,2),
273 sensor_los_geodetic(i,0),
274 sensor_los_geodetic(i,1),
275 refellipsoid );
276 cart2poslos_plain(sensor_pos(i,0),
277 sensor_pos(i,1),
278 sensor_pos(i,2),
279 sensor_los(i,0),
280 sensor_los(i,1),
281 x, y, z, dx, dy, dz);
282 sensor_pos(i,0) -= refell2r(refellipsoid,sensor_pos(i,1));
283 }
284 }
285 }
286}
287
288
289
290/* Workspace method: Doxygen documentation will be auto-generated */
292 const Vector& refellipsoid,
293 const Verbosity&) {
294 if (y_geo.empty() || refellipsoid[1] < 1e-7 || std::isnan(y_geo(0,0))) {
295 // Do nothing
296 } else {
297 for (Index i=0; i<y_geo.nrows(); i++) {
298 Numeric x, y, z, dx, dy, dz;
299 Numeric r = y_geo(i,0) + refell2r(refellipsoid,y_geo(i,1));
300 poslos2cart(x, y, z, dx, dy, dz, r,
301 y_geo(i,1), y_geo(i,2), y_geo(i,3), y_geo(i,4));
302 Numeric h, lat, lon, za, aa;
303 cart2geodeticposlos(h, lat, lon, za, aa, x, y, z, dx, dy, dz,refellipsoid);
304 y_geo(i,0) = h;
305 y_geo(i,1) = lat;
306 y_geo(i,2) = lon;
307 y_geo(i,3) = za;
308 y_geo(i,4) = aa;
309 }
310 }
311}
312
313
314/* Workspace method: Doxygen documentation will be auto-generated */
316 Matrix& los,
317 const Matrix& sensor_pos,
318 const Matrix& sensor_los,
319 const Vector& refellipsoid,
320 const Vector& lat_grid,
321 const Vector& lon_grid,
322 const Numeric& altitude,
323 const Verbosity&) {
324 //
325
326 ARTS_USER_ERROR_IF (sensor_pos.ncols() != 3,
327 "*sensor_pos* must have three columns.");
328 ARTS_USER_ERROR_IF (sensor_pos.nrows() == 0,
329 "*sensor_pos* must have at least one row.");
330 for (Index i=0; i<sensor_pos.nrows(); i++) {
331 ARTS_USER_ERROR_IF (sensor_pos(i,1) < -90 || sensor_pos(i,1) > 90,
332 "Unvalid latitude in *sensor_pos*.\n"
333 "Latitudes must be inside [-90,90],\n"
334 "but sensor_pos(",i,",1) is ", sensor_pos(i,1));
335 ARTS_USER_ERROR_IF (sensor_pos(i,2) < -180 || sensor_pos(i,1) > 360,
336 "Unvalid longitude in *sensor_pos*.\n"
337 "Longitudes must be inside [-180,360],\n"
338 "but sensor_pos(",i,",2) is ", sensor_pos(i,2));
339 }
340 ARTS_USER_ERROR_IF (sensor_los.ncols() != 2,
341 "*sensor_los* must have two columns.");
342 ARTS_USER_ERROR_IF (sensor_los.nrows() == 0,
343 "*sensor_los* must have at least one row.");
344 for (Index i=0; i<sensor_los.nrows(); i++) {
345 ARTS_USER_ERROR_IF (sensor_los(i,0) < 0 || sensor_los(i,0) > 180,
346 "Unvalid zenith angle in *sensor_los*.\n"
347 "Latitudes must be inside [0,180],\n"
348 "but sensor_los(",i,",0) is ", sensor_los(i,0));
349 ARTS_USER_ERROR_IF (sensor_los(i,1) < -180 || sensor_los(i,1) > 180,
350 "Unvalid azimuth angle in *sensor_los*.\n"
351 "Latitudes must be inside [-180,180],\n"
352 "but sensor_los(",i,",1) is ", sensor_los(i,1));
353 }
354 ARTS_USER_ERROR_IF (sensor_los.nrows() != sensor_pos.nrows(),
355 "*sensor_los* and *sensor_pos* must have the same number of rows.");
356
357 const Index n = sensor_pos.nrows();
358 pos.resize(n, 3);
359 los.resize(n, 2);
360
361 Numeric r1 = refellipsoid[0];
362 Numeric r2 = r1 * sqrt(1 - refellipsoid[1]*refellipsoid[1]);
363 r1 += altitude;
364 r2 += altitude;
365 Vector refell_target(2);
366 refell_target[0] = r1;
367 refell_target[1] = sqrt(1-pow(r2/r1,2.0));
368
369
370 const Index atmosphere_dim = 3;
371 for (Index i=0; i<n; i++) {
372 Numeric l, x, y, z, dx, dy, dz, r;
373 r = sensor_pos(i,0) + pos2refell_r(atmosphere_dim,
374 refellipsoid,
375 lat_grid,
376 lon_grid,
377 sensor_pos(i,joker));
378 poslos2cart(x, y, z, dx, dy, dz,
379 r, sensor_pos(i,1), sensor_pos(i,2),
380 sensor_los(i,0), sensor_los(i,1));
381 line_refellipsoid_intersect(l, refell_target, x, y, z, dx, dy, dz);
382
383 if (l<0) {
384 pos(i,joker) = NAN;
385 los(i,joker) = NAN;
386 } else {
387 cart2poslos_plain(pos(i,0),
388 pos(i,1),
389 pos(i,2),
390 los(i,0),
391 los(i,1),
392 x+l*dx,
393 y+l*dy,
394 z+l*dz,
395 dx,
396 dy,
397 dz);
398 pos(i,0) -= pos2refell_r(atmosphere_dim,
399 refellipsoid,
400 lat_grid,
401 lon_grid,
402 pos(i,joker));
403 // Check that pos matches expected altitude
404 //ARTS_ASSERT(abs(pos(i,0)-altitude) < 0.1);
405 }
406 }
407}
The global header file for ARTS.
Common ARTS conversions.
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:86
bool empty() const noexcept
Definition: matpackI.h:1086
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
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
Definition: geodetic.cc:1639
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)
Definition: geodetic.cc:543
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1266
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
Definition: geodetic.cc:1567
Numeric pos2refell_r(const Index &atmosphere_dim, ConstVectorView refellipsoid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView rte_pos)
pos2refell_r
Definition: geodetic.cc:1212
void line_refellipsoid_intersect(Numeric &l, const Vector &refellipsoid, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
geomtanpoint
Definition: geodetic.cc:875
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:356
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.
Definition: geodetic.cc:1522
void cart2sph_plain(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z)
Definition: geodetic.cc:642
void refellipsoidOrbitPlane(Vector &refellipsoid, const Numeric &orbitinc, const Verbosity &)
WORKSPACE METHOD: refellipsoidOrbitPlane.
Definition: m_geodetic.cc:80
void sensor_poslosFromGeodetic(Matrix &sensor_pos, Matrix &sensor_los, const Matrix &sensor_pos_geodetic, const Matrix &sensor_los_geodetic, const Vector &refellipsoid, const Verbosity &)
WORKSPACE METHOD: sensor_poslosFromGeodetic.
Definition: m_geodetic.cc:220
constexpr Numeric DEG2RAD
Definition: m_geodetic.cc:48
void y_geoToGeodetic(Matrix &y_geo, const Vector &refellipsoid, const Verbosity &)
WORKSPACE METHOD: y_geoToGeodetic.
Definition: m_geodetic.cc:291
void rte_poslosFromGeodetic(Vector &rte_pos, Vector &rte_los, const Matrix &sensor_pos_geodetic, const Matrix &sensor_los_geodetic, const Vector &refellipsoid, const Verbosity &v)
WORKSPACE METHOD: rte_poslosFromGeodetic.
Definition: m_geodetic.cc:136
void refellipsoidSet(Vector &refellipsoid, const Numeric &re, const Numeric &e, const Verbosity &)
WORKSPACE METHOD: refellipsoidSet.
Definition: m_geodetic.cc:97
void sensor_poslosFromECEF(Matrix &sensor_pos, Matrix &sensor_los, const Matrix &sensor_pos_ecef, const Matrix &sensor_los_ecef, const Vector &refellipsoid, const Verbosity &)
WORKSPACE METHOD: sensor_poslosFromECEF.
Definition: m_geodetic.cc:162
void rte_poslosFromECEF(Vector &rte_pos, Vector &rte_los, const Matrix &sensor_pos_ecef, const Matrix &sensor_los_ecef, const Vector &refellipsoid, const Verbosity &v)
WORKSPACE METHOD: rte_poslosFromECEF.
Definition: m_geodetic.cc:110
void refellipsoidForAzimuth(Vector &refellipsoid, const Numeric &latitude, const Numeric &azimuth, const Verbosity &)
WORKSPACE METHOD: refellipsoidForAzimuth.
Definition: m_geodetic.cc:55
void IntersectionGeometricalWithAltitude(Matrix &pos, Matrix &los, const Matrix &sensor_pos, const Matrix &sensor_los, const Vector &refellipsoid, const Vector &lat_grid, const Vector &lon_grid, const Numeric &altitude, const Verbosity &)
WORKSPACE METHOD: IntersectionGeometricalWithAltitude.
Definition: m_geodetic.cc:315
Implementation of Matrix, Vector, and such stuff.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
Declarations having to do with the four output streams.
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:705
#define v
#define a