1 /*
2 * Copyright (C) 2018 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.navigation.geodesic;
17
18 /**
19 * A geodesic line
20 * GeodesicLine facilitates the determination of a series of points on a single geodesic. The
21 * starting point (<i>lat1</i>, <i>lon1</i>) and the azimuth <i>azi1</i> are specified in the
22 * constructor; alternatively, the {@link Geodesic#line} method can be used to create a GeodesicLine.
23 * {@link #position} returns the location of point 2 a distance <i>s12</i> along the geodesic.
24 * Alternatively {@link #arcPosition} gives the position of point 2 an arc length <i>a12</i> along
25 * the geodesic.
26 * You can register the position of a reference point 3 a distance (arc length), <i>s13</i>
27 * (<i>a13</i>) along the geodesic with the {@link #setDistance} ({@link #setArc}) functions. Points
28 * a fractional distance along the line can be found by providing, for example, 0.5 * {@link #getDistance}
29 * as an argument to {@link #position}. The {@link Geodesic#inverseLine} or {@link Geodesic#directLine}
30 * methods return GeodesicLine objects with point 3 set to the point 2 of the corresponding geodesic
31 * problem. GeodesicLine objects created with the public constructor or with {@link Geodesic#line}
32 * have <i>s13</i> and <i>a13</i> set to NaNs.
33 * The calculations are accurate to better than 15 nm (15 nanometers). See Sec. 9 of
34 * <a href="https://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> for details. The algorithms used
35 * by this class are based on series expansions using the flattening <i>f</i> as a small parameter.
36 * These are only accurate for |<i>f</i> < 0.02; however reasonably accurate results will be
37 * obtained for |<i>f</i> < 0.2.
38 * The algorithms are described in
39 * <ul>
40 * <li>
41 * C. F. F. Karney, <a href="https://doi.org/10.1007/s00190-012-0578-z">
42 * Algorithms for geodesics</a>,
43 * J. Geodesy <b>87</b>, 43–55 (2013)
44 * (<a href="https://geographiclib.sourceforge.io/geod-addenda.html">addenda</a>).
45 * </li>
46 * </ul>
47 * Here's an example of using this class
48 * <pre>
49 * {@code
50 * import net.sf.geographiclib.*;
51 * public class GeodesicLineTest {
52 * public static void main(String[] args) {
53 * // Print waypoints between JFK and SIN
54 * Geodesic geod = Geodesic.WGS84;
55 * double
56 * lat1 = 40.640, lon1 = -73.779, // JFK
57 * lat2 = 1.359, lon2 = 103.989; // SIN
58 * GeodesicLine line = geod.InverseLine(lat1, lon1, lat2, lon2,
59 * GeodesicMask.DISTANCE_IN |
60 * GeodesicMask.LATITUDE |
61 * GeodesicMask.LONGITUDE);
62 * double ds0 = 500e3; // Nominal distance between points = 500 km
63 * // The number of intervals
64 * int num = (int)(Math.ceil(line.Distance() / ds0));
65 * {
66 * // Use intervals of equal length
67 * double ds = line.Distance() / num;
68 * for (int i = 0; i <= num; ++i) {
69 * GeodesicData g = line.Position(i * ds,
70 * GeodesicMask.LATITUDE |
71 * GeodesicMask.LONGITUDE);
72 * System.out.println(i + " " + g.lat2 + " " + g.lon2);
73 * }
74 * }
75 * {
76 * // Slightly faster, use intervals of equal arc length
77 * double da = line.Arc() / num;
78 * for (int i = 0; i <= num; ++i) {
79 * GeodesicData g = line.ArcPosition(i * da,
80 * GeodesicMask.LATITUDE |
81 * GeodesicMask.LONGITUDE);
82 * System.out.println(i + " " + g.lat2 + " " + g.lon2);
83 * }
84 * }
85 * }
86 * }}</pre>
87 */
88 public class GeodesicLine {
89
90 private static final int NC1 = Geodesic.NC1;
91 private static final int NC1P = Geodesic.NC1P;
92 private static final int NC2 = Geodesic.NC2;
93 private static final int NC3 = Geodesic.NC3;
94 private static final int NC4 = Geodesic.NC4;
95
96 private double lat1;
97 private double lon1;
98 private double azi1;
99
100 private double a;
101 private double f;
102 private double b;
103 private double c2;
104 private double f1;
105 private double salp0;
106 private double calp0;
107 private double k2;
108 private double salp1;
109 private double calp1;
110 private double ssig1;
111 private double csig1;
112 private double dn1;
113 private double stau1;
114 private double ctau1;
115 private double somg1;
116 private double comg1;
117 private double a1m1;
118 private double a2m1;
119 private double a3c;
120 private double b11;
121 private double b21;
122 private double b31;
123 private double a4;
124 private double b41;
125
126 private double a13;
127 private double s13;
128
129 // index zero elements of mC1a, mC1pa, mC2a, mC3a are unused
130 // all the elements of mC4a are used.
131 private double[] c1a;
132 private double[] c1pa;
133 private double[] c2a;
134 private double[] c3a;
135 private double[] c4a;
136
137 private int caps;
138
139 /**
140 * Constructor for a geodesic line staring at latitude <i>lat1</i>, longitude <i>lon1</i>,
141 * and azimuth <i>azi1</i> (all in degrees).
142 * If the point is at a pole, the azimuth is defined by keeping <i>lon1</i> fixed, writing
143 * <i>lat1</i> = ±(90° − ε), and taking the limit ε
144 * → 0+.
145 *
146 * @param g a {@link Geodesic} object used to compute the necessary information about the
147 * GeodesicLine
148 * @param lat1 latitude of point 1 (degrees). <i>lat1</i> should be in the range [−90°, 90°].
149 * @param lon1 longitude of point 1 (degrees).
150 * @param azi1 azimuth at point 1 (degrees).
151 */
152 public GeodesicLine(final Geodesic g, final double lat1, final double lon1, final double azi1) {
153 this(g, lat1, lon1, azi1, GeodesicMask.ALL);
154 }
155
156 /**
157 * Constructor for a geodesic line starting at latitude <i>lat1</i>, longitude <i>lon1</i>,
158 * and azimuth <i>azi1</i> (all in degrees) with a subset of the capabilities included.
159 * The {@link GeodesicMask} values are:
160 * <ul>
161 * <li>
162 * <i>caps</i> |= {@link GeodesicMask#LATITUDE} for the latitude <i>lat2</i>; this is
163 * added automatically
164 * </li>
165 * <li>
166 * <i>caps</i> |= {@link GeodesicMask#LONGITUDE} for the longitude <i>lon2</i>
167 * </li>
168 * <li>
169 * <i>caps</i> |= {@link GeodesicMask#AZIMUTH} for the azimuth <i>azi2</i>; this is
170 * added automatically
171 * </li>
172 * <li>
173 * <i>caps</i> |= {@link GeodesicMask#DISTANCE} for the distance <i>s12</i>
174 * </li>
175 * <li>
176 * <i>caps</i> |= {@link GeodesicMask#REDUCED_LENGTH} for the reduced length <i>m12</i>
177 * </li>
178 * <li>
179 * <i>caps</i> |= {@link GeodesicMask#GEODESIC_SCALE} for the geodesic scales <i>M12</i>
180 * and <i>M21</i>
181 * </li>
182 * <li>
183 * <i>caps</i> |= {@link GeodesicMask#AREA} for the area <i>S12</i>
184 * </li>
185 * <li>
186 * <i>caps</i> |= {@link GeodesicMask#DISTANCE_IN} permits the length of the geodesic
187 * to be given in terms of <i>s12</i>; without this capability the length can only be
188 * specified in terms of arc length
189 * </li>
190 * <li>
191 * <i>caps</i> |= {@link GeodesicMask#ALL} for all of the above.
192 * </li>
193 * </ul>
194 *
195 * @param g a {@link Geodesic} object used to compute the necessary information about the
196 * GeodesicLine.
197 * @param lat1 latitude of point 1 (degrees).
198 * @param lon1 longitude of point 1 (degrees).
199 * @param azi1 azimuth at point 1 (degrees).
200 * @param caps bitor'ed combination of {@link GeodesicMask} values specifying the capabilities
201 * the GeodesicLine object should possess, i.e., which quantities can be returned in
202 * calls to {@link #position}.
203 */
204 public GeodesicLine(final Geodesic g, final double lat1, final double lon1, double azi1, final int caps) {
205 azi1 = GeoMath.angNormalize(azi1);
206
207 final var p = GeoMath.sincosd(GeoMath.angRound(azi1));
208 final var pSalp1 = p.getFirst();
209 final var pCalp1 = p.getSecond();
210
211 lineInit(g, lat1, lon1, azi1, pSalp1, pCalp1, caps);
212 }
213
214 protected GeodesicLine(
215 final Geodesic g, final double lat1, final double lon1, final double azi1,
216 final double salp1, final double calp1, final int caps, final boolean arcmode, final double s13A13) {
217 lineInit(g, lat1, lon1, azi1, salp1, calp1, caps);
218 genSetDistance(arcmode, s13A13);
219 }
220
221 /**
222 * Compute the position of point 2 which is a distance <i>s12</i> (meters) from point 1.
223 * The values of <i>lon2</i> and <i>azi2</i> returned are in the range [−180°, 180°].
224 * The GeodesicLine object <i>must</i> have been constructed with <i>caps</i>
225 * |= {@link GeodesicMask#DISTANCE_IN}; otherwise no parameters are set.
226 *
227 * @param s12 distance from point 1 to point 2 (meters); it can be negative.
228 * @return a {@link GeodesicData} object with the following fields: <i>lat1</i>, <i>lon1</i>,
229 * <i>azi1</i>, <i>lat2</i>, <i>lon2</i>, <i>azi2</i>, <i>s12</i>, <i>a12</i>. Some of these
230 * results may be missing if the GeodesicLine did not include the relevant capability.
231 */
232 public GeodesicData position(final double s12) {
233 return position(false, s12, GeodesicMask.STANDARD);
234 }
235
236 /**
237 * Compute the position of point 2 which is a distance <i>s12</i> (meters) from point 1 and
238 * with a subset of the geodesic results returned.
239 * The GeodesicLine object <i>must</i> have been constructed with <i>caps</i> |=
240 * {@link GeodesicMask#DISTANCE_IN}; otherwise no parameters are set.
241 * Requesting a value which the GeodesicLine object is not capable of computing is not an error
242 * (no parameters will be set). The value of <i>lon2</i> returned is normally in the range
243 * [−180°, 180°]; however if the <i>outmask</i> includes the
244 * {@link GeodesicMask#LONG_UNROLL} flag, the longitude is "unrolled" so that the quantity
245 * <i>lon2</i> − <i>lon1</i> indicates how many times and in what sense the geodesic
246 * encircles the ellipsoid.
247 *
248 * @param s12 distance from point 1 to point 2 (meters); it can be negative.
249 * @param outmask a bitor'ed combination of {@link GeodesicMask} values specifying which results
250 * should be returned.
251 * @return a {@link GeodesicData} object including the requested results.
252 */
253 public GeodesicData position(final double s12, final int outmask) {
254 return position(false, s12, outmask);
255 }
256
257 /**
258 * Compute the position of point 2 which is an arc length <i>a12</i> (degrees) from point 1.
259 * The values of <i>lon2</i> and <i>azi2</i> returned are in the range [−180°, 180°].
260 * The GeodesicLine object <i>must</i> have been constructed with <i>caps</i> |=
261 * {@link GeodesicMask#DISTANCE_IN}; otherwise no parameters are set.
262 *
263 * @param a12 arc length from point 1 to point 2 (degrees); it can be negative.
264 * @return a {@link GeodesicData} object with the following fields: <i>lat1</i>, <i>lon1</i>,
265 * <i>azi1</i>, <i>lat2</i>, <i>lon2</i>, <i>azi2</i>, <i>s12</i>, <i>a12</i>. Some of these
266 * results may be missing if the GeodesicLine did not include the relevant capability.
267 */
268 public GeodesicData arcPosition(final double a12) {
269 return position(true, a12, GeodesicMask.STANDARD);
270 }
271
272 /**
273 * Compute the position of point 2 which is an arc length <i>a12</i> (degrees) from point 1 and
274 * with a subset of th geodesic results returned.
275 * Requesting a value which the GeodesicLine object is not capable of computing is not an error
276 * (no parameters will be set). The value of <i>lon2</i> returned is in the range
277 * [−180°, 180°], unless the <i>outmask</i> includes the
278 * {@link GeodesicMask#LONG_UNROLL} flags.
279 *
280 * @param a12 arc length from point 1 to point 2 (degrees); it can be negative.
281 * @param outmask a bitor'ed combination of {@link GeodesicMask} values specifying which results
282 * should be returned.
283 * @return a {@link GeodesicData} object giving <i>lat1</i>, <i>lon2</i>, <i>azi2</i>, and
284 * <i>a12</i>.
285 */
286 public GeodesicData arcPosition(final double a12, final int outmask) {
287 return position(true, a12, outmask);
288 }
289
290 /**
291 * The general position function. {@link #position(double, int)} and {@link #arcPosition(double, int)}
292 * are defined in terms of this function.
293 * The {@link GeodesicMask} values possible for <i>outmask</i> are
294 * <ul>
295 * <li>
296 * <i>outmask</i> |= {@link GeodesicMask#LATITUDE} for the latitude <i>lat2</i>;
297 * </li>
298 * <li>
299 * <i>outmask</i> |= {@link GeodesicMask#LONGITUDE} for the latitude <i>lon2</i>;
300 * </li>
301 * <li>
302 * <i>outmask</i> |= {@link GeodesicMask#AZIMUTH} for the latitude <i>azi2</i>;
303 * </li>
304 * <li>
305 * <i>outmask</i> |= {@link GeodesicMask#DISTANCE} for the distance <i>s12</i>;
306 * </li>
307 * <li>
308 * <i>outmask</i> |= {@link GeodesicMask#REDUCED_LENGTH} for the reduced length
309 * <i>m12</i>;
310 * </li>
311 * <li>
312 * <i>outmask</i> |= {@link GeodesicMask#GEODESIC_SCALE} for the geodesic scales
313 * <i>M12</i> and <i>M21</i>;
314 * </li>
315 * <li>
316 * <i>outmask</i> |= {@link GeodesicMask#ALL} for all of the above;
317 * </li>
318 * <li>
319 * <i>outmask</i> |= {@link GeodesicMask#LONG_UNROLL} to unroll <i>lon2</i> (instead of
320 * reducing it to the range [−180°, 180°]).
321 * </li>
322 * </ul>
323 *
324 * @param arcmode boolean flag determining the meaning of the second parameter; if arcmode is false,
325 * then the GeodesicLine object must have been constructed with <i>caps</i> |=
326 * {@link GeodesicMask#DISTANCE_IN}.
327 * @param s12A12 if <i>arcmode</i> is false, this is the distance between point 1 and point 2 (meters);
328 * otherwise it is the arc length between point 1 and point 2 (degrees); it can be
329 * negative.
330 * @param outmask a bitor'ed combination of {@link GeodesicMask} values specifying which results
331 * should be returned.
332 * @return a {@link GeodesicData} object with the requested results. Requesting a value which the
333 * GeodesicLine object is not capable of computing is not an error; Double.NaN is returned instead.
334 */
335 public GeodesicData position(final boolean arcmode, final double s12A12, int outmask) {
336 outmask &= caps & GeodesicMask.OUT_MASK;
337 final var r = new GeodesicData();
338 if (!(init() && (arcmode || (caps & (GeodesicMask.OUT_MASK & GeodesicMask.DISTANCE_IN)) != 0))) {
339 // uninitialized or impossible distance calculation requested
340 return r;
341 }
342
343 r.setLat1(lat1);
344 r.setAzi1(azi1);
345 r.setLon1(((outmask & GeodesicMask.LONG_UNROLL) != 0) ? lon1 : GeoMath.angNormalize(lon1));
346
347 // avoid warning about uninitialized b12
348 double sig12;
349 double ssig12;
350 double csig12;
351 var b12 = 0.0;
352 var ab1 = 0.0;
353 if (arcmode) {
354 // interpret s12A12 as spherical arc length
355 r.setA12(s12A12);
356 sig12 = Math.toRadians(s12A12);
357
358 final var p = GeoMath.sincosd(s12A12);
359 ssig12 = p.getFirst();
360 csig12 = p.getSecond();
361 } else {
362 // interpret s12A12 as distance
363 r.setS12(s12A12);
364
365 final var tau12 = s12A12 / (b * (1 + a1m1));
366 final var s = Math.sin(tau12);
367 final var c = Math.cos(tau12);
368
369 // tau2 = tau1 + tau12
370 b12 = -Geodesic.sinCosSeries(true, stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, c1pa);
371 sig12 = tau12 - (b12 - b11);
372 ssig12 = Math.sin(sig12);
373 csig12 = Math.cos(sig12);
374
375 if (Math.abs(f) > 0.01) {
376 // reverted distance series is inaccurate for |f| > 1/100, so correct sig12 with 1
377 // Newton iteration. The following table shows the approximate maximum error for
378 // a = WGSa() and various f relative to GeodesicExact.
379 // erri = the error in the inverse solution (nm)
380 // errd = the error in the direct solution (series only) (nm)
381 // errda = the error in the direct solution (series + 1 Newton) (nm)
382
383 // f erri errd errda
384 // -1/5 12e6 1.2e9 69e6
385 // -1/10 123e3 12e6 765e3
386 // -1/20 1110 108e3 7155
387 // -1/50 18.63 200.9 27.12
388 // -1/100 18.63 23.78 23.37
389 // -1/150 18.63 21.05 20.26
390 // 1/150 22.35 24.73 25.83
391 // 1/100 22.35 25.03 25.31
392 // 1/50 29.80 231.9 30.44
393 // 1/20 5376 146e3 10e3
394 // 1/10 829e3 22e6 1.5e6
395 // 1/5 157e6 3.8e9 280e6
396
397 final var ssig2 = ssig1 * csig12 + csig1 * ssig12;
398 final var csig2 = csig1 * csig12 - ssig1 * ssig12;
399 b12 = Geodesic.sinCosSeries(true, ssig2, csig2, c1a);
400
401 final var serr = (1 + a1m1) * (sig12 + (b12 - b11)) - s12A12 / b;
402 sig12 = sig12 - serr / Math.sqrt(1 + k2 * GeoMath.sq(ssig2));
403 ssig12 = Math.sin(sig12);
404 csig12 = Math.cos(sig12);
405 // update b12 below
406 }
407 r.setA12(Math.toDegrees(sig12));
408 }
409
410 final var ssig2 = ssig1 * csig12 + csig1 * ssig12;
411 var csig2 = csig1 * csig12 - ssig1 * ssig12;
412 final double sbet2;
413 double cbet2;
414 final double salp2;
415 final double calp2;
416 // sig2 = sig1 + sig12
417
418
419 final var dn2 = Math.sqrt(1 + k2 * GeoMath.sq(ssig2));
420 if ((outmask & (GeodesicMask.DISTANCE | GeodesicMask.REDUCED_LENGTH | GeodesicMask.GEODESIC_SCALE)) != 0) {
421 if (arcmode || Math.abs(f) > 0.01) {
422 b12 = Geodesic.sinCosSeries(true, ssig2, csig2, c1a);
423 }
424 ab1 = (1 + a1m1) * (b12 - b11);
425 }
426
427 // sin(bet2) = cos(alp0) * sin(sig2)
428 sbet2 = calp0 * ssig2;
429 // alt: cbet2 = hypot(csig2, salp0 * ssig2)
430 cbet2 = GeoMath.hypot(salp0, calp0 * csig2);
431 if (cbet2 == 0) {
432 // i.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
433 cbet2 = csig2 = Geodesic.TINY;
434 }
435
436 // tan(alp0) = cos(sig2) * tan(alp2)
437 salp2 = salp0;
438 // no need to normalize
439 calp2 = calp0 * csig2;
440
441 if ((outmask & GeodesicMask.DISTANCE) != 0 && arcmode) {
442 r.setS12(b * ((1 + a1m1) * sig12 + ab1));
443 }
444
445 if ((outmask & GeodesicMask.LONGITUDE) != 0) {
446 // tan(omg2) = sin(alp0) * tan(sig2)
447
448 // no need to normalize east or west going?
449 //noinspection all
450 final var somg2 = salp0 * ssig2;
451 final var e = GeoMath.copysign(1, salp0);
452
453 // omg12 = omg2 - omg1
454 final var omg12 = ((outmask & GeodesicMask.LONG_UNROLL) != 0)
455 ? e * (sig12 - (Math.atan2(ssig2, csig2) - Math.atan2(ssig1, csig1))
456 + (Math.atan2(e * somg2, csig2) - Math.atan2(e * somg1, comg1)))
457 : Math.atan2(somg2 * comg1 - csig2 * somg1, csig2 * comg1 + somg2 * somg1);
458
459 final var lam12 = omg12 + a3c * (sig12 + (Geodesic.sinCosSeries(true, ssig2, csig2, c3a) - b31));
460 final var lon12 = Math.toDegrees(lam12);
461 r.setLon2(((outmask & GeodesicMask.LONG_UNROLL) != 0)
462 ? lon1 + lon12 : GeoMath.angNormalize(r.getLon1() + GeoMath.angNormalize(lon12)));
463 }
464
465 if ((outmask & GeodesicMask.LATITUDE) != 0) {
466 r.setLat2(GeoMath.atan2d(sbet2, f1 * cbet2));
467 }
468
469 if ((outmask & GeodesicMask.AZIMUTH) != 0) {
470 r.setAzi2(GeoMath.atan2d(salp2, calp2));
471 }
472
473 if ((outmask & (GeodesicMask.REDUCED_LENGTH | GeodesicMask.GEODESIC_SCALE)) != 0) {
474 final var b22 = Geodesic.sinCosSeries(true, ssig2, csig2, c2a);
475 final var ab2 = (1 + a2m1) * (b22 - b21);
476 final var j12 = (a1m1 - a2m1) * sig12 + (ab1 - ab2);
477 if ((outmask & GeodesicMask.REDUCED_LENGTH) != 0) {
478 // add parens around (mCsig1 * ssig2) and (mSsig1 * csig2) to ensure
479 // accurate cancellation in the case of coincident points
480 r.setM12(b * ((dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2)) - csig1 * csig2 * j12));
481 }
482 if ((outmask & GeodesicMask.GEODESIC_SCALE) != 0) {
483 final var t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2);
484 r.setScaleM12(csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1);
485 r.setScaleM21(csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2);
486 }
487 }
488
489 if ((outmask & GeodesicMask.AREA) != 0) {
490 final var b42 = Geodesic.sinCosSeries(false, ssig2, csig2, c4a);
491 final double salp12;
492 final double calp12;
493 if (calp0 == 0 || salp0 == 0) {
494 //alp12 = alp2 - alp1, used in atan2 so no need to normalize
495 salp12 = salp2 * calp1 - calp2 * salp1;
496 calp12 = calp2 * calp1 + salp2 * salp1;
497 } else {
498 // tan(alp) = tan(alp0) * sec(sig)
499 // tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
500 // = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
501 // if csig12 > 0, write
502 // csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
503 // else
504 // csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
505 // no need to normalize
506 salp12 = calp0 * salp0 * (csig12 <= 0
507 ? csig1 * (1 - csig12) + ssig12 * ssig1
508 : ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1));
509 calp12 = GeoMath.sq(salp0) + GeoMath.sq(calp0) * csig1 * csig2;
510 }
511 r.setAreaS12(c2 * Math.atan2(salp12, calp12) + a4 * (b42 - b41));
512 }
513
514 return r;
515 }
516
517 /**
518 * Specify position of point 3 in terms of distance.
519 * This is only useful if the GeodesicLine object has been constructed with <i>caps</i> |=
520 * {@link GeodesicMask#DISTANCE_IN}.
521 *
522 * @param s13 the distance from point 1 to point 3 (meters); it can be negative.
523 */
524 public void setDistance(final double s13) {
525 this.s13 = s13;
526 final var g = position(false, this.s13, 0);
527 a13 = g.getA12();
528 }
529
530 /**
531 * Specify position of point 3 in terms of either distance or arc length.
532 *
533 * @param arcmode boolean flag determining the meaning of the second parameter; if <i>arcmode</i>
534 * is false, then the GeodesicLine object must have been constructed with
535 * <i>caps</i> |= {@link GeodesicMask#DISTANCE_IN}.
536 * @param s13A13 if <i>arcmode</i> is false, this is the distance from point 1 to point 3 (meters);
537 * otherwise it is the arc length from point 1 to point 3 (degrees); it can be
538 * negative.
539 */
540 public void genSetDistance(final boolean arcmode, final double s13A13) {
541 if (arcmode) {
542 setArc(s13A13);
543 } else {
544 setDistance(s13A13);
545 }
546 }
547
548 /**
549 * Gets the latitude of point 1 (degrees).
550 *
551 * @return <i>lat1</i> the latitude of point 1 (degrees).
552 */
553 public double getLatitude() {
554 return init() ? lat1 : Double.NaN;
555 }
556
557 /**
558 * Gets the longitude of point 1 (degrees).
559 *
560 * @return <i>lon1</i> the longitude of point 1 (degrees).
561 */
562 public double getLongitude() {
563 return init() ? lon1 : Double.NaN;
564 }
565
566 /**
567 * Gets the azimuth (degrees) of the geodesic line at point 1.
568 *
569 * @return <i>azi1</i> the azimuth (degrees) of the geodesic line at point 1.
570 */
571 public double getAzimuth() {
572 return init() ? azi1 : Double.NaN;
573 }
574
575 /**
576 * Gets a pair of sine and cosine of <i>azi1</i> the azimuth (degrees) of the geodesic line at
577 * point 1.
578 *
579 * @return pair of sine and cosine of <i>azi1</i> the azimuth (degrees) of the geodesic line
580 * at point 1.
581 */
582 public Pair getAzimuthCosines() {
583 return new Pair(init() ? salp1 : Double.NaN, init() ? calp1 : Double.NaN);
584 }
585
586 /**
587 * Gets the azimuth (degrees) of the geodesic line as it crosses the equator in a northward
588 * direction.
589 *
590 * @return <i>azi0</i> the azimuth (degrees) of the geodesic line as it crosses the equator in a
591 * northward direction.
592 */
593 public double getEquatorialAzimuth() {
594 return init() ? GeoMath.atan2d(salp0, calp0) : Double.NaN;
595 }
596
597 /**
598 * Gets a pair of sine and cosine of <i>azi0</i> the azimuth of the godesic line as it crosses the
599 * equator in a northward direction.
600 *
601 * @return pair of sine and cosine of <i>azi0</i> the azimuth of the godesic line as it crosses
602 * the equator in a northward direction.
603 */
604 public Pair getEquatorialAzimuthCosines() {
605 return new Pair(init() ? salp0 : Double.NaN, init() ? calp0 : Double.NaN);
606 }
607
608 /**
609 * Gets the arc length (degrees) between the northward equatorial crossing and point 1.
610 *
611 * @return <i>a1</i> the arc length (degrees) between the northward equatorial crossing and
612 * point 1.
613 */
614 public double getEquatorialArc() {
615 return init() ? GeoMath.atan2d(ssig1, csig1) : Double.NaN;
616 }
617
618 /**
619 * Gets the equatorial radius of the ellipsoid (meters). This is the value inherited from the
620 * Geodesic object used in the constructor.
621 *
622 * @return <i>a</i> the equatorial radius of the ellipsoid (meters).
623 */
624 public double getMajorRadius() {
625 return init() ? a : Double.NaN;
626 }
627
628 /**
629 * Gets the flattening of the ellipsoid. This is the value inherited from the Geodesic object
630 * used in the constructor.
631 *
632 * @return <i>f</i> the flattening of the ellipsoid.
633 */
634 public double getFlattening() {
635 return init() ? f : Double.NaN;
636 }
637
638 /**
639 * Gets the computational capabilities that this object was constructed with. LATITUDE and AZIMUTH
640 * are always included.
641 *
642 * @return <i>caps</i> the computation capabilities that this object was constructed with.
643 */
644 public int getCapabilities() {
645 return caps;
646 }
647
648 /**
649 * Indicates whether this GeodesicLine object has all tested capabilities
650 *
651 * @param testcaps a set of bitor'ed {@link GeodesicMask} values.
652 * @return true if the GeodesicLine object has all these capabilities.
653 */
654 public boolean capabilities(int testcaps) {
655 testcaps &= GeodesicMask.OUT_ALL;
656 return (caps & testcaps) == testcaps;
657 }
658
659 /**
660 * The distance or arc length to point 3.
661 *
662 * @param arcmode boolean flag determining the meaning of returned value.
663 * @return <i>s13</i> if <i>arcmode</i> is false; <i>a13</i> if <i>arcmode</i> is true.
664 */
665 public double genDistance(final boolean arcmode) {
666 final var tmp = arcmode ? a13 : s13;
667 return init() ? tmp : Double.NaN;
668 }
669
670 /**
671 * Gets the distance to point 3 (meters).
672 *
673 * @return <i>s13</i> the disance to point 3 (meters).
674 */
675 public double getDistance() {
676 return genDistance(false);
677 }
678
679 /**
680 * Gets the arc length to point 3 (degrees).
681 *
682 * @return <i>a13</i> the arc length to point 3 (degrees).
683 */
684 public double getArc() {
685 return genDistance(true);
686 }
687
688 /**
689 * Specify position of point 3 in terms of arc length.
690 * The distance <i>s13</i> is only set if the GeodesicLine object has been constructed with
691 * <i>caps</i> |= {@link GeodesicMask#DISTANCE}.
692 *
693 * @param a13 the arc length from point 1 to point 3 (degrees); it can be negative.
694 */
695 void setArc(final double a13) {
696 this.a13 = a13;
697 final var g = position(true, this.a13, GeodesicMask.DISTANCE);
698 s13 = g.getS12();
699 }
700
701 /**
702 * @return true if the object has been initialized.
703 */
704 private boolean init() {
705 return caps != 0;
706 }
707
708 private void lineInit(
709 final Geodesic g, final double lat1, final double lon1, final double azi1, final double salp1,
710 final double calp1, final int caps) {
711 a = g.a;
712 f = g.f;
713 b = g.b;
714 c2 = g.c2;
715 f1 = g.f1;
716
717 // always allow latitude and azimuth and unrolling the longitude
718 this.caps = caps | GeodesicMask.LATITUDE | GeodesicMask.AZIMUTH | GeodesicMask.LONG_UNROLL;
719
720 this.lat1 = GeoMath.latFix(lat1);
721 this.lon1 = lon1;
722 this.azi1 = azi1;
723 this.salp1 = salp1;
724 this.calp1 = calp1;
725
726 var p = GeoMath.sincosd(GeoMath.angRound(this.lat1));
727 var sbet1 = f1 * p.getFirst();
728 var cbet1 = p.getSecond();
729
730 // ensure cbet1 = +epsilon at poles
731 p = GeoMath.norm(sbet1, cbet1);
732 sbet1 = p.getFirst();
733 cbet1 = Math.max(Geodesic.TINY, p.getSecond());
734
735 dn1 = Math.sqrt(1 + g.ep2 * GeoMath.sq(sbet1));
736
737 // evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
738
739 // alp0 in [0, pi/2 - |bet1|]
740 salp0 = this.salp1 * cbet1;
741
742 // alt: calp0 = hypot(sbet1, calp1 * cbet1). The following is slightly
743 // better (consider the case salp1 = 0).
744 calp0 = GeoMath.hypot(this.calp1, this.salp1 * sbet1);
745
746 // Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
747 // sig = 0 is nearest northward crossing of the equator.
748 // With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
749 // With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
750 // With bet1 = -pi/2, alp1 = 0, sig1 = -pi/2
751 // Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
752 // With alp0 in (0, pi/2], quadrants for sig and omg coincide.
753 // No atan2(0,0) ambiguity at poles since cbet1 = +epsilon
754 // With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
755 ssig1 = sbet1;
756 somg1 = salp0 * sbet1;
757 csig1 = comg1 = sbet1 != 0 || this.calp1 != 0 ? cbet1 * this.calp1 : 1;
758
759 p = GeoMath.norm(ssig1, csig1);
760 ssig1 = p.getFirst();
761 // sig 1 in (-pi, pi]
762 csig1 = p.getSecond();
763
764 // GeoMath.norm(mSomg1, mComg1); -- don't need to normalize!
765
766 k2 = GeoMath.sq(calp0) * g.ep2;
767 final var eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2);
768
769 if ((this.caps & GeodesicMask.CAP_C1) != 0) {
770 a1m1 = Geodesic.a1m1f(eps);
771 c1a = new double[NC1 + 1];
772 Geodesic.c1f(eps, c1a);
773 b11 = Geodesic.sinCosSeries(true, ssig1, csig1, c1a);
774 final var s = Math.sin(b11);
775 final var c = Math.cos(b11);
776 // tau1 = sig1 + b11
777 stau1 = ssig1 * c + csig1 * s;
778 ctau1 = csig1 * c - ssig1 * s;
779 // not necessary because c1pa rverts c1a
780 // mB11 = -sinCosSeries(true, mStau1, mCtau1, mC1pa, NC1P)
781 }
782
783 if ((this.caps & GeodesicMask.CAP_C1P) != 0) {
784 c1pa = new double[NC1P + 1];
785 Geodesic.c1pf(eps, c1pa);
786 }
787
788 if ((this.caps & GeodesicMask.CAP_C2) != 0) {
789 c2a = new double[NC2 + 1];
790 a2m1 = Geodesic.a2m1f(eps);
791 Geodesic.c2f(eps, c2a);
792 b21 = Geodesic.sinCosSeries(true, ssig1, csig1, c2a);
793 }
794
795 if ((this.caps & GeodesicMask.CAP_C3) != 0) {
796 c3a = new double[NC3];
797 g.c3f(eps, c3a);
798 a3c = -f * salp0 * g.a3f(eps);
799 b31 = Geodesic.sinCosSeries(true, ssig1, csig1, c3a);
800 }
801
802 if ((this.caps & GeodesicMask.CAP_C4) != 0) {
803 c4a = new double[NC4];
804 g.c4f(eps, c4a);
805 // multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
806 a4 = GeoMath.sq(a) * calp0 * salp0 * g.e2;
807 b41 = Geodesic.sinCosSeries(false, ssig1, csig1, c4a);
808 }
809 }
810 }