1+ /**
2+ * Implementation of the net.sf.geographiclib.GeoMath class
3+ *
4+ * Copyright (c) Charles Karney (2013-2020) <[email protected] > and licensed 5+ * under the MIT/X11 License. For more information, see
6+ * https://geographiclib.sourceforge.io/
7+ **********************************************************************/
8+ package org .locationtech .proj4j .geodesic ;
9+
10+ /**
11+ * Mathematical functions needed by GeographicLib.
12+ * <p>
13+ * Define mathematical functions and constants so that any version of Java
14+ * can be used.
15+ **********************************************************************/
16+ public class GeoMath {
17+ /**
18+ * The number of binary digits in the fraction of a double precision
19+ * number (equivalent to C++'s {@code numeric_limits<double>::digits}).
20+ **********************************************************************/
21+ public static final int digits = 53 ;
22+
23+ /**
24+ * Square a number.
25+ * <p>
26+ * @param x the argument.
27+ * @return <i>x</i><sup>2</sup>.
28+ **********************************************************************/
29+ public static double sq (double x ) { return x * x ; }
30+
31+ /**
32+ * The inverse hyperbolic tangent function. This is defined in terms of
33+ * Math.log1p(<i>x</i>) in order to maintain accuracy near <i>x</i> = 0.
34+ * In addition, the odd parity of the function is enforced.
35+ * <p>
36+ * @param x the argument.
37+ * @return atanh(<i>x</i>).
38+ **********************************************************************/
39+ public static double atanh (double x ) {
40+ double y = Math .abs (x ); // Enforce odd parity
41+ y = Math .log1p (2 * y /(1 - y ))/2 ;
42+ return x > 0 ? y : (x < 0 ? -y : x );
43+ }
44+
45+ /**
46+ * Normalize a sine cosine pair.
47+ * <p>
48+ * @param p return parameter for normalized quantities with sinx<sup>2</sup>
49+ * + cosx<sup>2</sup> = 1.
50+ * @param sinx the sine.
51+ * @param cosx the cosine.
52+ **********************************************************************/
53+ public static void norm (Pair p , double sinx , double cosx ) {
54+ double r = Math .hypot (sinx , cosx );
55+ p .first = sinx /r ; p .second = cosx /r ;
56+ }
57+
58+ /**
59+ * The error-free sum of two numbers.
60+ * <p>
61+ * @param u the first number in the sum.
62+ * @param v the second number in the sum.
63+ * @param p output Pair(<i>s</i>, <i>t</i>) with <i>s</i> = round(<i>u</i> +
64+ * <i>v</i>) and <i>t</i> = <i>u</i> + <i>v</i> - <i>s</i>.
65+ * <p>
66+ * See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B.
67+ **********************************************************************/
68+ public static void sum (Pair p , double u , double v ) {
69+ double s = u + v ;
70+ double up = s - v ;
71+ double vpp = s - up ;
72+ up -= u ;
73+ vpp -= v ;
74+ double t = s != 0 ? 0.0 - (up + vpp ) : s ;
75+ // u + v = s + t
76+ // = round(u + v) + t
77+ p .first = s ; p .second = t ;
78+ }
79+
80+ /**
81+ * Evaluate a polynomial.
82+ * <p>
83+ * @param N the order of the polynomial.
84+ * @param p the coefficient array (of size <i>N</i> + <i>s</i> + 1 or more).
85+ * @param s starting index for the array.
86+ * @param x the variable.
87+ * @return the value of the polynomial.
88+ * <p>
89+ * Evaluate <i>y</i> = ∑<sub><i>n</i>=0..<i>N</i></sub>
90+ * <i>p</i><sub><i>s</i>+<i>n</i></sub>
91+ * <i>x</i><sup><i>N</i>−<i>n</i></sup>. Return 0 if <i>N</i> < 0.
92+ * Return <i>p</i><sub><i>s</i></sub>, if <i>N</i> = 0 (even if <i>x</i> is
93+ * infinite or a nan). The evaluation uses Horner's method.
94+ **********************************************************************/
95+ public static double polyval (int N , double p [], int s , double x ) {
96+ double y = N < 0 ? 0 : p [s ++];
97+ while (--N >= 0 ) y = y * x + p [s ++];
98+ return y ;
99+ }
100+
101+ /**
102+ * Coarsen a value close to zero.
103+ * <p>
104+ * @param x the argument
105+ * @return the coarsened value.
106+ * <p>
107+ * This makes the smallest gap in <i>x</i> = 1/16 − nextafter(1/16, 0)
108+ * = 1/2<sup>57</sup> for reals = 0.7 pm on the earth if <i>x</i> is an angle
109+ * in degrees. (This is about 1000 times more resolution than we get with
110+ * angles around 90 degrees.) We use this to avoid having to deal with near
111+ * singular cases when <i>x</i> is non-zero but tiny (e.g.,
112+ * 10<sup>−200</sup>). This converts −0 to +0; however tiny
113+ * negative numbers get converted to −0.
114+ **********************************************************************/
115+ public static double AngRound (double x ) {
116+ final double z = 1 /16.0 ;
117+ double y = Math .abs (x );
118+ // The compiler mustn't "simplify" z - (z - y) to y
119+ y = y < z ? z - (z - y ) : y ;
120+ return Math .copySign (y , x );
121+ }
122+
123+ /**
124+ * Normalize an angle.
125+ * <p>
126+ * @param x the angle in degrees.
127+ * @return the angle reduced to the range [−180°, 180°).
128+ * <p>
129+ * The range of <i>x</i> is unrestricted.
130+ **********************************************************************/
131+ public static double AngNormalize (double x ) {
132+ double y = Math .IEEEremainder (x , 360.0 );
133+ return Math .abs (y ) == 180 ? Math .copySign (180.0 , x ) : y ;
134+ }
135+
136+ /**
137+ * Normalize a latitude.
138+ * <p>
139+ * @param x the angle in degrees.
140+ * @return x if it is in the range [−90°, 90°], otherwise
141+ * return NaN.
142+ **********************************************************************/
143+ public static double LatFix (double x ) {
144+ return Math .abs (x ) > 90 ? Double .NaN : x ;
145+ }
146+
147+ /**
148+ * The exact difference of two angles reduced to [−180°, 180°].
149+ * <p>
150+ * @param x the first angle in degrees.
151+ * @param y the second angle in degrees.
152+ * @param p output Pair(<i>d</i>, <i>e</i>) with <i>d</i> being the rounded
153+ * difference and <i>e</i> being the error.
154+ * <p>
155+ * This computes <i>z</i> = <i>y</i> − <i>x</i> exactly, reduced to
156+ * [−180°, 180°]; and then sets <i>z</i> = <i>d</i> + <i>e</i>
157+ * where <i>d</i> is the nearest representable number to <i>z</i> and
158+ * <i>e</i> is the truncation error. If <i>z</i> = ±0° or
159+ * ±180°, then the sign of <i>d</i> is given by the sign of
160+ * <i>y</i> − <i>x</i>. The maximum absolute value of <i>e</i> is
161+ * 2<sup>−26</sup> (for doubles).
162+ **********************************************************************/
163+ public static void AngDiff (Pair p , double x , double y ) {
164+ sum (p , Math .IEEEremainder (-x , 360.0 ), Math .IEEEremainder (y , 360.0 ));
165+ sum (p , Math .IEEEremainder (p .first , 360.0 ), p .second );
166+ if (p .first == 0 || Math .abs (p .first ) == 180 )
167+ // p = [d, e]...
168+ // If e == 0, take sign from y - x
169+ // else (e != 0, implies d = +/-180), d and e must have opposite signs
170+ p .first = Math .copySign (p .first , p .second == 0 ? y - x : -p .second );
171+ }
172+
173+ /**
174+ * Evaluate the sine and cosine function with the argument in degrees
175+ *
176+ * @param p return Pair(<i>s</i>, <i>t</i>) with <i>s</i> = sin(<i>x</i>) and
177+ * <i>c</i> = cos(<i>x</i>).
178+ * @param x in degrees.
179+ * <p>
180+ * The results obey exactly the elementary properties of the trigonometric
181+ * functions, e.g., sin 9° = cos 81° = − sin 123456789°.
182+ **********************************************************************/
183+ public static void sincosd (Pair p , double x ) {
184+ // In order to minimize round-off errors, this function exactly reduces
185+ // the argument to the range [-45, 45] before converting it to radians.
186+ double r ; int q ;
187+ r = x % 360.0 ;
188+ q = (int )Math .round (r / 90 ); // If r is NaN this returns 0
189+ r -= 90 * q ;
190+ // now abs(r) <= 45
191+ r = Math .toRadians (r );
192+ // Possibly could call the gnu extension sincos
193+ double s = Math .sin (r ), c = Math .cos (r );
194+ double sinx , cosx ;
195+ switch (q & 3 ) {
196+ case 0 : sinx = s ; cosx = c ; break ;
197+ case 1 : sinx = c ; cosx = -s ; break ;
198+ case 2 : sinx = -s ; cosx = -c ; break ;
199+ default : sinx = -c ; cosx = s ; break ; // case 3
200+ }
201+ if (sinx == 0 ) sinx = Math .copySign (sinx , x );
202+ p .first = sinx ; p .second = 0.0 + cosx ;
203+ }
204+
205+ /**
206+ * Evaluate the sine and cosine function with reduced argument plus correction
207+ *
208+ * @param p return Pair(<i>s</i>, <i>t</i>) with <i>s</i> =
209+ * sin(<i>x</i> + <i>t</i>) and <i>c</i> = cos(<i>x</i> + <i>t</i>).
210+ * @param x reduced angle in degrees.
211+ * @param t correction in degrees.
212+ * <p>
213+ * This is a variant of GeoMath.sincosd allowing a correction to the angle to
214+ * be supplied. <i>x</i> x must be in [−180°, 180°] and
215+ * <i>t</i> is assumed to be a <i>small</i> correction. GeoMath.AngRound is
216+ * applied to the reduced angle to prevent problems with <i>x</i> + <i>t</i>
217+ * being extremely close but not exactly equal to one of the four cardinal
218+ * directions.
219+ **********************************************************************/
220+ public static void sincosde (Pair p , double x , double t ) {
221+ // In order to minimize round-off errors, this function exactly reduces
222+ // the argument to the range [-45, 45] before converting it to radians.
223+ double r ; int q ;
224+ q = (int )Math .round (x / 90 ); // If r is NaN this returns 0
225+ r = x - 90 * q ;
226+ // now abs(r) <= 45
227+ r = Math .toRadians (GeoMath .AngRound (r + t ));
228+ // Possibly could call the gnu extension sincos
229+ double s = Math .sin (r ), c = Math .cos (r );
230+ double sinx , cosx ;
231+ switch (q & 3 ) {
232+ case 0 : sinx = s ; cosx = c ; break ;
233+ case 1 : sinx = c ; cosx = -s ; break ;
234+ case 2 : sinx = -s ; cosx = -c ; break ;
235+ default : sinx = -c ; cosx = s ; break ; // case 3
236+ }
237+ if (sinx == 0 ) sinx = Math .copySign (sinx , x );
238+ p .first = sinx ; p .second = 0.0 + cosx ;
239+ }
240+
241+ /**
242+ * Evaluate the atan2 function with the result in degrees
243+ *
244+ * @param y the sine of the angle
245+ * @param x the cosine of the angle
246+ * @return atan2(<i>y</i>, <i>x</i>) in degrees.
247+ * <p>
248+ * The result is in the range [−180° 180°]. N.B.,
249+ * atan2d(±0, −1) = ±180°.
250+ **********************************************************************/
251+ public static double atan2d (double y , double x ) {
252+ // In order to minimize round-off errors, this function rearranges the
253+ // arguments so that result of atan2 is in the range [-pi/4, pi/4] before
254+ // converting it to degrees and mapping the result to the correct
255+ // quadrant.
256+ int q = 0 ;
257+ if (Math .abs (y ) > Math .abs (x )) { double t ; t = x ; x = y ; y = t ; q = 2 ; }
258+ if (x < 0 ) { x = -x ; ++q ; }
259+ // here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4]
260+ double ang = Math .toDegrees (Math .atan2 (y , x ));
261+ switch (q ) {
262+ // Note that atan2d(-0.0, 1.0) will return -0. However, we expect that
263+ // atan2d will not be called with y = -0. If need be, include
264+ //
265+ // case 0: ang = 0 + ang; break;
266+ //
267+ // and handle mpfr as in AngRound.
268+ case 1 : ang = Math .copySign (180.0 , y ) - ang ; break ;
269+ case 2 : ang = 90 - ang ; break ;
270+ case 3 : ang = -90 + ang ; break ;
271+ default : break ;
272+ }
273+ return ang ;
274+ }
275+
276+ private GeoMath () {}
277+ }
0 commit comments