Fortran library for Geodesics  1.50.1
geodesic.for
Go to the documentation of this file.
1 * The subroutines in this files are documented at
2 * https://geographiclib.sourceforge.io/html/Fortran/
3 *
4 *> @file geodesic.for
5 *! @brief Implementation of geodesic routines in Fortran
6 *!
7 *! This is a Fortran implementation of the geodesic algorithms described
8 *! in
9 *! - C. F. F. Karney,
10 *! <a href="https://doi.org/10.1007/s00190-012-0578-z">
11 *! Algorithms for geodesics</a>,
12 *! J. Geodesy <b>87</b>, 43--55 (2013);
13 *! DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
14 *! 10.1007/s00190-012-0578-z</a>;
15 *! addenda: <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
16 *! geod-addenda.html</a>.
17 *! .
18 *! The principal advantages of these algorithms over previous ones
19 *! (e.g., Vincenty, 1975) are
20 *! - accurate to round off for |<i>f</i>| &lt; 1/50;
21 *! - the solution of the inverse problem is always found;
22 *! - differential and integral properties of geodesics are computed.
23 *!
24 *! The shortest path between two points on the ellipsoid at (\e lat1, \e
25 *! lon1) and (\e lat2, \e lon2) is called the geodesic. Its length is
26 *! \e s12 and the geodesic from point 1 to point 2 has forward azimuths
27 *! \e azi1 and \e azi2 at the two end points.
28 *!
29 *! Traditionally two geodesic problems are considered:
30 *! - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1,
31 *! determine \e lat2, \e lon2, and \e azi2. This is solved by the
32 *! subroutine direct().
33 *! - the inverse problem -- given \e lat1, \e lon1, \e lat2, \e lon2,
34 *! determine \e s12, \e azi1, and \e azi2. This is solved by the
35 *! subroutine invers().
36 *!
37 *! The ellipsoid is specified by its equatorial radius \e a (typically
38 *! in meters) and flattening \e f. The routines are accurate to round
39 *! off with double precision arithmetic provided that |<i>f</i>| &lt;
40 *! 1/50; for the WGS84 ellipsoid, the errors are less than 15
41 *! nanometers. (Reasonably accurate results are obtained for |<i>f</i>|
42 *! &lt; 1/5.) For a prolate ellipsoid, specify \e f &lt; 0.
43 *!
44 *! The routines also calculate several other quantities of interest
45 *! - \e SS12 is the area between the geodesic from point 1 to point 2
46 *! and the equator; i.e., it is the area, measured counter-clockwise,
47 *! of the geodesic quadrilateral with corners (\e lat1,\e lon1), (0,\e
48 *! lon1), (0,\e lon2), and (\e lat2,\e lon2).
49 *! - \e m12, the reduced length of the geodesic is defined such that if
50 *! the initial azimuth is perturbed by \e dazi1 (radians) then the
51 *! second point is displaced by \e m12 \e dazi1 in the direction
52 *! perpendicular to the geodesic. On a curved surface the reduced
53 *! length obeys a symmetry relation, \e m12 + \e m21 = 0. On a flat
54 *! surface, we have \e m12 = \e s12.
55 *! - \e MM12 and \e MM21 are geodesic scales. If two geodesics are
56 *! parallel at point 1 and separated by a small distance \e dt, then
57 *! they are separated by a distance \e MM12 \e dt at point 2. \e MM21
58 *! is defined similarly (with the geodesics being parallel to one
59 *! another at point 2). On a flat surface, we have \e MM12 = \e MM21
60 *! = 1.
61 *! - \e a12 is the arc length on the auxiliary sphere. This is a
62 *! construct for converting the problem to one in spherical
63 *! trigonometry. \e a12 is measured in degrees. The spherical arc
64 *! length from one equator crossing to the next is always 180&deg;.
65 *!
66 *! If points 1, 2, and 3 lie on a single geodesic, then the following
67 *! addition rules hold:
68 *! - \e s13 = \e s12 + \e s23
69 *! - \e a13 = \e a12 + \e a23
70 *! - \e SS13 = \e SS12 + \e SS23
71 *! - \e m13 = \e m12 \e MM23 + \e m23 \e MM21
72 *! - \e MM13 = \e MM12 \e MM23 &minus; (1 &minus; \e MM12 \e MM21) \e
73 *! m23 / \e m12
74 *! - \e MM31 = \e MM32 \e MM21 &minus; (1 &minus; \e MM23 \e MM32) \e
75 *! m12 / \e m23
76 *!
77 *! The shortest distance returned by the solution of the inverse problem
78 *! is (obviously) uniquely defined. However, in a few special cases
79 *! there are multiple azimuths which yield the same shortest distance.
80 *! Here is a catalog of those cases:
81 *! - \e lat1 = &minus;\e lat2 (with neither point at a pole). If \e
82 *! azi1 = \e azi2, the geodesic is unique. Otherwise there are two
83 *! geodesics and the second one is obtained by setting [\e azi1, \e
84 *! azi2] &rarr; [\e azi2, \e azi1], [\e MM12, \e MM21] &rarr; [\e
85 *! MM21, \e MM12], \e SS12 &rarr; &minus;\e SS12. (This occurs when
86 *! the longitude difference is near &plusmn;180&deg; for oblate
87 *! ellipsoids.)
88 *! - \e lon2 = \e lon1 &plusmn; 180&deg; (with neither point at a pole).
89 *! If \e azi1 = 0&deg; or &plusmn;180&deg;, the geodesic is unique.
90 *! Otherwise there are two geodesics and the second one is obtained by
91 *! setting [\e azi1, \e azi2] &rarr; [&minus;\e azi1, &minus;\e azi2],
92 *! \e SS12 &rarr; &minus;\e SS12. (This occurs when \e lat2 is near
93 *! &minus;\e lat1 for prolate ellipsoids.)
94 *! - Points 1 and 2 at opposite poles. There are infinitely many
95 *! geodesics which can be generated by setting [\e azi1, \e azi2]
96 *! &rarr; [\e azi1, \e azi2] + [\e d, &minus;\e d], for arbitrary \e
97 *! d. (For spheres, this prescription applies when points 1 and 2 are
98 *! antipodal.)
99 *! - \e s12 = 0 (coincident points). There are infinitely many
100 *! geodesics which can be generated by setting [\e azi1, \e azi2]
101 *! &rarr; [\e azi1, \e azi2] + [\e d, \e d], for arbitrary \e d.
102 *!
103 *! These routines are a simple transcription of the corresponding C++
104 *! classes in <a href="https://geographiclib.sourceforge.io">
105 *! GeographicLib</a>. Because of the limitations of Fortran 77, the
106 *! classes have been replaced by simple subroutines with no attempt to
107 *! save "state" across subroutine calls. Most of the internal comments
108 *! have been retained. However, in the process of transcription some
109 *! documentation has been lost and the documentation for the C++
110 *! classes, GeographicLib::Geodesic, GeographicLib::GeodesicLine, and
111 *! GeographicLib::PolygonAreaT, should be consulted. The C++ code
112 *! remains the "reference implementation". Think twice about
113 *! restructuring the internals of the Fortran code since this may make
114 *! porting fixes from the C++ code more difficult.
115 *!
116 *! Copyright (c) Charles Karney (2012-2019) <charles@karney.com> and
117 *! licensed under the MIT/X11 License. For more information, see
118 *! https://geographiclib.sourceforge.io/
119 *!
120 *! This library was distributed with
121 *! <a href="../index.html">GeographicLib</a> 1.50.
122 
123 *> Solve the direct geodesic problem
124 *!
125 *! @param[in] a the equatorial radius (meters).
126 *! @param[in] f the flattening of the ellipsoid. Setting \e f = 0 gives
127 *! a sphere. Negative \e f gives a prolate ellipsoid.
128 *! @param[in] lat1 latitude of point 1 (degrees).
129 *! @param[in] lon1 longitude of point 1 (degrees).
130 *! @param[in] azi1 azimuth at point 1 (degrees).
131 *! @param[in] s12a12 if \e arcmode is not set, this is the distance
132 *! from point 1 to point 2 (meters); otherwise it is the arc
133 *! length from point 1 to point 2 (degrees); it can be negative.
134 *! @param[in] flags a bitor'ed combination of the \e arcmode and \e
135 *! unroll flags.
136 *! @param[out] lat2 latitude of point 2 (degrees).
137 *! @param[out] lon2 longitude of point 2 (degrees).
138 *! @param[out] azi2 (forward) azimuth at point 2 (degrees).
139 *! @param[in] omask a bitor'ed combination of mask values
140 *! specifying which of the following parameters should be set.
141 *! @param[out] a12s12 if \e arcmode is not set, this is the arc length
142 *! from point 1 to point 2 (degrees); otherwise it is the distance
143 *! from point 1 to point 2 (meters).
144 *! @param[out] m12 reduced length of geodesic (meters).
145 *! @param[out] MM12 geodesic scale of point 2 relative to point 1
146 *! (dimensionless).
147 *! @param[out] MM21 geodesic scale of point 1 relative to point 2
148 *! (dimensionless).
149 *! @param[out] SS12 area under the geodesic (meters<sup>2</sup>).
150 *!
151 *! \e flags is an integer in [0, 4) whose binary bits are interpreted
152 *! as follows
153 *! - 1 the \e arcmode flag
154 *! - 2 the \e unroll flag
155 *! .
156 *! If \e arcmode is not set, \e s12a12 is \e s12 and \e a12s12 is \e
157 *! a12; otherwise, \e s12a12 is \e a12 and \e a12s12 is \e s12. It \e
158 *! unroll is not set, the value \e lon2 returned is in the range
159 *! [&minus;180&deg;, 180&deg;]; if unroll is set, the longitude variable
160 *! is "unrolled" so that \e lon2 &minus; \e lon1 indicates how many
161 *! times and in what sense the geodesic encircles the ellipsoid.
162 *!
163 *! \e omask is an integer in [0, 16) whose binary bits are interpreted
164 *! as follows
165 *! - 1 return \e a12
166 *! - 2 return \e m12
167 *! - 4 return \e MM12 and \e MM21
168 *! - 8 return \e SS12
169 *!
170 *! \e lat1 should be in the range [&minus;90&deg;, 90&deg;]. The value
171 *! \e azi2 returned is in the range [&minus;180&deg;, 180&deg;].
172 *!
173 *! If either point is at a pole, the azimuth is defined by keeping the
174 *! longitude fixed, writing \e lat = \e lat = &plusmn;(90&deg; &minus;
175 *! &epsilon;), and taking the limit &epsilon; &rarr; 0+. An arc length
176 *! greater that 180&deg; signifies a geodesic which is not a shortest
177 *! path. (For a prolate ellipsoid, an additional condition is necessary
178 *! for a shortest path: the longitudinal extent must not exceed of
179 *! 180&deg;.)
180 *!
181 *! Example of use:
182 *! \include geoddirect.for
183 
184  subroutine direct(a, f, lat1, lon1, azi1, s12a12, flags,
185  + lat2, lon2, azi2, omask, a12s12, m12, MM12, MM21, SS12)
186 * input
187  double precision a, f, lat1, lon1, azi1, s12a12
188  integer flags, omask
189 * output
190  double precision lat2, lon2, azi2
191 * optional output
192  double precision a12s12, m12, MM12, MM21, SS12
193 
194  integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
195  parameter(ord = 6, nc1 = ord, nc1p = ord,
196  + nc2 = ord, na3 = ord, na3x = na3,
197  + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
198  + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
199  double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1),
200  + c1a(nc1), c1pa(nc1p), c2a(nc2), c3a(nc3-1), c4a(0:nc4-1)
201 
202  double precision atanhx, hypotx,
203  + angnm, angrnd, trgsum, a1m1f, a2m1f, a3f, atn2dx, latfix
204  logical arcmod, unroll, arcp, redlp, scalp, areap
205  double precision e2, f1, ep2, n, b, c2,
206  + salp0, calp0, k2, eps,
207  + salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1,
208  + salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2,
209  + ssig12, csig12, salp12, calp12, omg12, lam12, lon12,
210  + sig12, stau1, ctau1, tau12, t, s, c, serr, e,
211  + a1m1, a2m1, a3c, a4, ab1, ab2,
212  + b11, b12, b21, b22, b31, b41, b42, j12
213 
214  double precision dblmin, dbleps, pi, degree, tiny,
215  + tol0, tol1, tol2, tolb, xthrsh
216  integer digits, maxit1, maxit2
217  logical init
218  common /geocom/ dblmin, dbleps, pi, degree, tiny,
219  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
220 
221  if (.not.init) call geoini
222 
223  e2 = f * (2 - f)
224  ep2 = e2 / (1 - e2)
225  f1 = 1 - f
226  n = f / (2 - f)
227  b = a * f1
228  c2 = 0
229 
230  arcmod = mod(flags/1, 2) .eq. 1
231  unroll = mod(flags/2, 2) .eq. 1
232 
233  arcp = mod(omask/1, 2) .eq. 1
234  redlp = mod(omask/2, 2) .eq. 1
235  scalp = mod(omask/4, 2) .eq. 1
236  areap = mod(omask/8, 2) .eq. 1
237 
238  if (areap) then
239  if (e2 .eq. 0) then
240  c2 = a**2
241  else if (e2 .gt. 0) then
242  c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
243  else
244  c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
245  end if
246  end if
247 
248  call a3cof(n, a3x)
249  call c3cof(n, c3x)
250  if (areap) call c4cof(n, c4x)
251 
252 * Guard against underflow in salp0
253  call sncsdx(angrnd(azi1), salp1, calp1)
254 
255  call sncsdx(angrnd(latfix(lat1)), sbet1, cbet1)
256  sbet1 = f1 * sbet1
257  call norm2x(sbet1, cbet1)
258 * Ensure cbet1 = +dbleps at poles
259  cbet1 = max(tiny, cbet1)
260  dn1 = sqrt(1 + ep2 * sbet1**2)
261 
262 * Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
263 * alp0 in [0, pi/2 - |bet1|]
264  salp0 = salp1 * cbet1
265 * Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following
266 * is slightly better (consider the case salp1 = 0).
267  calp0 = hypotx(calp1, salp1 * sbet1)
268 * Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
269 * sig = 0 is nearest northward crossing of equator.
270 * With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
271 * With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
272 * With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2
273 * Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
274 * With alp0 in (0, pi/2], quadrants for sig and omg coincide.
275 * No atan2(0,0) ambiguity at poles since cbet1 = +dbleps.
276 * With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
277  ssig1 = sbet1
278  somg1 = salp0 * sbet1
279  if (sbet1 .ne. 0 .or. calp1 .ne. 0) then
280  csig1 = cbet1 * calp1
281  else
282  csig1 = 1
283  end if
284  comg1 = csig1
285 * sig1 in (-pi, pi]
286  call norm2x(ssig1, csig1)
287 * norm2x(somg1, comg1); -- don't need to normalize!
288 
289  k2 = calp0**2 * ep2
290  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
291 
292  a1m1 = a1m1f(eps)
293  call c1f(eps, c1a)
294  b11 = trgsum(.true., ssig1, csig1, c1a, nc1)
295  s = sin(b11)
296  c = cos(b11)
297 * tau1 = sig1 + B11
298  stau1 = ssig1 * c + csig1 * s
299  ctau1 = csig1 * c - ssig1 * s
300 * Not necessary because C1pa reverts C1a
301 * B11 = -TrgSum(true, stau1, ctau1, C1pa, nC1p)
302 
303  if (.not. arcmod) call c1pf(eps, c1pa)
304 
305  if (redlp .or. scalp) then
306  a2m1 = a2m1f(eps)
307  call c2f(eps, c2a)
308  b21 = trgsum(.true., ssig1, csig1, c2a, nc2)
309  else
310 * Suppress bogus warnings about unitialized variables
311  a2m1 = 0
312  b21 = 0
313  end if
314 
315  call c3f(eps, c3x, c3a)
316  a3c = -f * salp0 * a3f(eps, a3x)
317  b31 = trgsum(.true., ssig1, csig1, c3a, nc3-1)
318 
319  if (areap) then
320  call c4f(eps, c4x, c4a)
321 * Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
322  a4 = a**2 * calp0 * salp0 * e2
323  b41 = trgsum(.false., ssig1, csig1, c4a, nc4)
324  else
325 * Suppress bogus warnings about unitialized variables
326  a4 = 0
327  b41 = 0
328  end if
329 
330  if (arcmod) then
331 * Interpret s12a12 as spherical arc length
332  sig12 = s12a12 * degree
333  call sncsdx(s12a12, ssig12, csig12)
334 * Suppress bogus warnings about unitialized variables
335  b12 = 0
336  else
337 * Interpret s12a12 as distance
338  tau12 = s12a12 / (b * (1 + a1m1))
339  s = sin(tau12)
340  c = cos(tau12)
341 * tau2 = tau1 + tau12
342  b12 = - trgsum(.true.,
343  + stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, c1pa, nc1p)
344  sig12 = tau12 - (b12 - b11)
345  ssig12 = sin(sig12)
346  csig12 = cos(sig12)
347  if (abs(f) .gt. 0.01d0) then
348 * Reverted distance series is inaccurate for |f| > 1/100, so correct
349 * sig12 with 1 Newton iteration. The following table shows the
350 * approximate maximum error for a = WGS_a() and various f relative to
351 * GeodesicExact.
352 * erri = the error in the inverse solution (nm)
353 * errd = the error in the direct solution (series only) (nm)
354 * errda = the error in the direct solution (series + 1 Newton) (nm)
355 *
356 * f erri errd errda
357 * -1/5 12e6 1.2e9 69e6
358 * -1/10 123e3 12e6 765e3
359 * -1/20 1110 108e3 7155
360 * -1/50 18.63 200.9 27.12
361 * -1/100 18.63 23.78 23.37
362 * -1/150 18.63 21.05 20.26
363 * 1/150 22.35 24.73 25.83
364 * 1/100 22.35 25.03 25.31
365 * 1/50 29.80 231.9 30.44
366 * 1/20 5376 146e3 10e3
367 * 1/10 829e3 22e6 1.5e6
368 * 1/5 157e6 3.8e9 280e6
369  ssig2 = ssig1 * csig12 + csig1 * ssig12
370  csig2 = csig1 * csig12 - ssig1 * ssig12
371  b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
372  serr = (1 + a1m1) * (sig12 + (b12 - b11)) - s12a12 / b
373  sig12 = sig12 - serr / sqrt(1 + k2 * ssig2**2)
374  ssig12 = sin(sig12)
375  csig12 = cos(sig12)
376 * Update B12 below
377  end if
378  end if
379 
380 * sig2 = sig1 + sig12
381  ssig2 = ssig1 * csig12 + csig1 * ssig12
382  csig2 = csig1 * csig12 - ssig1 * ssig12
383  dn2 = sqrt(1 + k2 * ssig2**2)
384  if (arcmod .or. abs(f) .gt. 0.01d0)
385  + b12 = trgsum(.true., ssig2, csig2, c1a, nc1)
386  ab1 = (1 + a1m1) * (b12 - b11)
387 
388 * sin(bet2) = cos(alp0) * sin(sig2)
389  sbet2 = calp0 * ssig2
390 * Alt: cbet2 = hypot(csig2, salp0 * ssig2)
391  cbet2 = hypotx(salp0, calp0 * csig2)
392  if (cbet2 .eq. 0) then
393 * I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
394  cbet2 = tiny
395  csig2 = cbet2
396  end if
397 * tan(omg2) = sin(alp0) * tan(sig2)
398 * No need to normalize
399  somg2 = salp0 * ssig2
400  comg2 = csig2
401 * tan(alp0) = cos(sig2)*tan(alp2)
402 * No need to normalize
403  salp2 = salp0
404  calp2 = calp0 * csig2
405 * East or west going?
406  e = sign(1d0, salp0)
407 * omg12 = omg2 - omg1
408  if (unroll) then
409  omg12 = e * (sig12
410  + - (atan2( ssig2, csig2) - atan2( ssig1, csig1))
411  + + (atan2(e * somg2, comg2) - atan2(e * somg1, comg1)))
412  else
413  omg12 = atan2(somg2 * comg1 - comg2 * somg1,
414  + comg2 * comg1 + somg2 * somg1)
415  end if
416 
417  lam12 = omg12 + a3c *
418  + ( sig12 + (trgsum(.true., ssig2, csig2, c3a, nc3-1)
419  + - b31))
420  lon12 = lam12 / degree
421  if (unroll) then
422  lon2 = lon1 + lon12
423  else
424  lon2 = angnm(angnm(lon1) + angnm(lon12))
425  end if
426  lat2 = atn2dx(sbet2, f1 * cbet2)
427  azi2 = atn2dx(salp2, calp2)
428 
429  if (redlp .or. scalp) then
430  b22 = trgsum(.true., ssig2, csig2, c2a, nc2)
431  ab2 = (1 + a2m1) * (b22 - b21)
432  j12 = (a1m1 - a2m1) * sig12 + (ab1 - ab2)
433  end if
434 * Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
435 * accurate cancellation in the case of coincident points.
436  if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
437  + dn1 * (ssig1 * csig2)) - csig1 * csig2 * j12)
438  if (scalp) then
439  t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
440  mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
441  mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
442  end if
443 
444  if (areap) then
445  b42 = trgsum(.false., ssig2, csig2, c4a, nc4)
446  if (calp0 .eq. 0 .or. salp0 .eq. 0) then
447 * alp12 = alp2 - alp1, used in atan2 so no need to normalize
448  salp12 = salp2 * calp1 - calp2 * salp1
449  calp12 = calp2 * calp1 + salp2 * salp1
450  else
451 * tan(alp) = tan(alp0) * sec(sig)
452 * tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
453 * = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
454 * If csig12 > 0, write
455 * csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
456 * else
457 * csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
458 * No need to normalize
459  if (csig12 .le. 0) then
460  salp12 = csig1 * (1 - csig12) + ssig12 * ssig1
461  else
462  salp12 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
463  end if
464  salp12 = calp0 * salp0 * salp12
465  calp12 = salp0**2 + calp0**2 * csig1 * csig2
466  end if
467  ss12 = c2 * atan2(salp12, calp12) + a4 * (b42 - b41)
468  end if
469 
470  if (arcp) then
471  if (arcmod) then
472  a12s12 = b * ((1 + a1m1) * sig12 + ab1)
473  else
474  a12s12 = sig12 / degree
475  end if
476  end if
477 
478  return
479  end
480 
481 *> Solve the inverse geodesic problem.
482 *!
483 *! @param[in] a the equatorial radius (meters).
484 *! @param[in] f the flattening of the ellipsoid. Setting \e f = 0 gives
485 *! a sphere. Negative \e f gives a prolate ellipsoid.
486 *! @param[in] lat1 latitude of point 1 (degrees).
487 *! @param[in] lon1 longitude of point 1 (degrees).
488 *! @param[in] lat2 latitude of point 2 (degrees).
489 *! @param[in] lon2 longitude of point 2 (degrees).
490 *! @param[out] s12 distance from point 1 to point 2 (meters).
491 *! @param[out] azi1 azimuth at point 1 (degrees).
492 *! @param[out] azi2 (forward) azimuth at point 2 (degrees).
493 *! @param[in] omask a bitor'ed combination of mask values
494 *! specifying which of the following parameters should be set.
495 *! @param[out] a12 arc length from point 1 to point 2 (degrees).
496 *! @param[out] m12 reduced length of geodesic (meters).
497 *! @param[out] MM12 geodesic scale of point 2 relative to point 1
498 *! (dimensionless).
499 *! @param[out] MM21 geodesic scale of point 1 relative to point 2
500 *! (dimensionless).
501 *! @param[out] SS12 area under the geodesic (meters<sup>2</sup>).
502 *!
503 *! \e omask is an integer in [0, 16) whose binary bits are interpreted
504 *! as follows
505 *! - 1 return \e a12
506 *! - 2 return \e m12
507 *! - 4 return \e MM12 and \e MM21
508 *! - 8 return \e SS12
509 *!
510 *! \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;].
511 *! The values of \e azi1 and \e azi2 returned are in the range
512 *! [&minus;180&deg;, 180&deg;].
513 *!
514 *! If either point is at a pole, the azimuth is defined by keeping the
515 *! longitude fixed, writing \e lat = &plusmn;(90&deg; &minus;
516 *! &epsilon;), and taking the limit &epsilon; &rarr; 0+.
517 *!
518 *! The solution to the inverse problem is found using Newton's method.
519 *! If this fails to converge (this is very unlikely in geodetic
520 *! applications but does occur for very eccentric ellipsoids), then the
521 *! bisection method is used to refine the solution.
522 *!
523 *! Example of use:
524 *! \include geodinverse.for
525 
526  subroutine invers(a, f, lat1, lon1, lat2, lon2,
527  + s12, azi1, azi2, omask, a12, m12, MM12, MM21, SS12)
528 * input
529  double precision a, f, lat1, lon1, lat2, lon2
530  integer omask
531 * output
532  double precision s12, azi1, azi2
533 * optional output
534  double precision a12, m12, MM12, MM21, SS12
535 
536  integer ord, nA3, nA3x, nC3, nC3x, nC4, nC4x, nC
537  parameter (ord = 6, na3 = ord, na3x = na3,
538  + nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2,
539  + nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2,
540  + nc = ord)
541  double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1),
542  + Ca(nC)
543 
544  double precision atanhx, hypotx,
545  + AngDif, AngRnd, TrgSum, Lam12f, InvSta, atn2dx, LatFix
546  integer latsgn, lonsgn, swapp, numit
547  logical arcp, redlp, scalp, areap, merid, tripn, tripb
548 
549  double precision e2, f1, ep2, n, b, c2,
550  + lat1x, lat2x, salp0, calp0, k2, eps,
551  + salp1, calp1, ssig1, csig1, cbet1, sbet1, dbet1, dn1,
552  + salp2, calp2, ssig2, csig2, sbet2, cbet2, dbet2, dn2,
553  + slam12, clam12, salp12, calp12, omg12, lam12, lon12, lon12s,
554  + salp1a, calp1a, salp1b, calp1b,
555  + dalp1, sdalp1, cdalp1, nsalp1, alp12, somg12, comg12, domg12,
556  + sig12, v, dv, dnm, dummy,
557  + a4, b41, b42, s12x, m12x, a12x, sdomg12, cdomg12
558 
559  double precision dblmin, dbleps, pi, degree, tiny,
560  + tol0, tol1, tol2, tolb, xthrsh
561  integer digits, maxit1, maxit2, lmask
562  logical init
563  common /geocom/ dblmin, dbleps, pi, degree, tiny,
564  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
565 
566  if (.not.init) call geoini
567 
568  f1 = 1 - f
569  e2 = f * (2 - f)
570  ep2 = e2 / f1**2
571  n = f / ( 2 - f)
572  b = a * f1
573  c2 = 0
574 
575  arcp = mod(omask/1, 2) .eq. 1
576  redlp = mod(omask/2, 2) .eq. 1
577  scalp = mod(omask/4, 2) .eq. 1
578  areap = mod(omask/8, 2) .eq. 1
579  if (scalp) then
580  lmask = 16 + 2 + 4
581  else
582  lmask = 16 + 2
583  end if
584 
585  if (areap) then
586  if (e2 .eq. 0) then
587  c2 = a**2
588  else if (e2 .gt. 0) then
589  c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
590  else
591  c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
592  end if
593  end if
594 
595  call a3cof(n, a3x)
596  call c3cof(n, c3x)
597  if (areap) call c4cof(n, c4x)
598 
599 * Compute longitude difference (AngDiff does this carefully). Result is
600 * in [-180, 180] but -180 is only for west-going geodesics. 180 is for
601 * east-going and meridional geodesics.
602 * If very close to being on the same half-meridian, then make it so.
603  lon12 = angdif(lon1, lon2, lon12s)
604 * Make longitude difference positive.
605  if (lon12 .ge. 0) then
606  lonsgn = 1
607  else
608  lonsgn = -1
609  end if
610  lon12 = lonsgn * angrnd(lon12)
611  lon12s = angrnd((180 - lon12) - lonsgn * lon12s)
612  lam12 = lon12 * degree
613  if (lon12 .gt. 90) then
614  call sncsdx(lon12s, slam12, clam12)
615  clam12 = -clam12
616  else
617  call sncsdx(lon12, slam12, clam12)
618  end if
619 
620 * If really close to the equator, treat as on equator.
621  lat1x = angrnd(latfix(lat1))
622  lat2x = angrnd(latfix(lat2))
623 * Swap points so that point with higher (abs) latitude is point 1
624 * If one latitude is a nan, then it becomes lat1.
625  if (abs(lat1x) .lt. abs(lat2x)) then
626  swapp = -1
627  else
628  swapp = 1
629  end if
630  if (swapp .lt. 0) then
631  lonsgn = -lonsgn
632  call swap(lat1x, lat2x)
633  end if
634 * Make lat1 <= 0
635  if (lat1x .lt. 0) then
636  latsgn = 1
637  else
638  latsgn = -1
639  end if
640  lat1x = lat1x * latsgn
641  lat2x = lat2x * latsgn
642 * Now we have
643 *
644 * 0 <= lon12 <= 180
645 * -90 <= lat1 <= 0
646 * lat1 <= lat2 <= -lat1
647 *
648 * longsign, swapp, latsgn register the transformation to bring the
649 * coordinates to this canonical form. In all cases, 1 means no change
650 * was made. We make these transformations so that there are few cases
651 * to check, e.g., on verifying quadrants in atan2. In addition, this
652 * enforces some symmetries in the results returned.
653 
654  call sncsdx(lat1x, sbet1, cbet1)
655  sbet1 = f1 * sbet1
656  call norm2x(sbet1, cbet1)
657 * Ensure cbet1 = +dbleps at poles
658  cbet1 = max(tiny, cbet1)
659 
660  call sncsdx(lat2x, sbet2, cbet2)
661  sbet2 = f1 * sbet2
662  call norm2x(sbet2, cbet2)
663 * Ensure cbet2 = +dbleps at poles
664  cbet2 = max(tiny, cbet2)
665 
666 * If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
667 * |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1
668 * is a better measure. This logic is used in assigning calp2 in
669 * Lambda12. Sometimes these quantities vanish and in that case we force
670 * bet2 = +/- bet1 exactly. An example where is is necessary is the
671 * inverse problem 48.522876735459 0 -48.52287673545898293
672 * 179.599720456223079643 which failed with Visual Studio 10 (Release and
673 * Debug)
674 
675  if (cbet1 .lt. -sbet1) then
676  if (cbet2 .eq. cbet1) sbet2 = sign(sbet1, sbet2)
677  else
678  if (abs(sbet2) .eq. -sbet1) cbet2 = cbet1
679  end if
680 
681  dn1 = sqrt(1 + ep2 * sbet1**2)
682  dn2 = sqrt(1 + ep2 * sbet2**2)
683 
684 * Suppress bogus warnings about unitialized variables
685  a12x = 0
686  merid = lat1x .eq. -90 .or. slam12 .eq. 0
687 
688  if (merid) then
689 
690 * Endpoints are on a single full meridian, so the geodesic might lie on
691 * a meridian.
692 
693 * Head to the target longitude
694  calp1 = clam12
695  salp1 = slam12
696 * At the target we're heading north
697  calp2 = 1
698  salp2 = 0
699 
700 * tan(bet) = tan(sig) * cos(alp)
701  ssig1 = sbet1
702  csig1 = calp1 * cbet1
703  ssig2 = sbet2
704  csig2 = calp2 * cbet2
705 
706 * sig12 = sig2 - sig1
707  sig12 = atan2(0d0 + max(0d0, csig1 * ssig2 - ssig1 * csig2),
708  + csig1 * csig2 + ssig1 * ssig2)
709  call lengs(n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
710  + cbet1, cbet2, lmask,
711  + s12x, m12x, dummy, mm12, mm21, ep2, ca)
712 
713 * Add the check for sig12 since zero length geodesics might yield m12 <
714 * 0. Test case was
715 *
716 * echo 20.001 0 20.001 0 | GeodSolve -i
717 *
718 * In fact, we will have sig12 > pi/2 for meridional geodesic which is
719 * not a shortest path.
720  if (sig12 .lt. 1 .or. m12x .ge. 0) then
721  if (sig12 .lt. 3 * tiny) then
722  sig12 = 0
723  m12x = 0
724  s12x = 0
725  end if
726  m12x = m12x * b
727  s12x = s12x * b
728  a12x = sig12 / degree
729  else
730 * m12 < 0, i.e., prolate and too close to anti-podal
731  merid = .false.
732  end if
733  end if
734 
735  omg12 = 0
736 * somg12 > 1 marks that it needs to be calculated
737  somg12 = 2
738  comg12 = 0
739  if (.not. merid .and. sbet1 .eq. 0 .and.
740  + (f .le. 0 .or. lon12s .ge. f * 180)) then
741 
742 * Geodesic runs along equator
743  calp1 = 0
744  calp2 = 0
745  salp1 = 1
746  salp2 = 1
747  s12x = a * lam12
748  sig12 = lam12 / f1
749  omg12 = sig12
750  m12x = b * sin(sig12)
751  if (scalp) then
752  mm12 = cos(sig12)
753  mm21 = mm12
754  end if
755  a12x = lon12 / f1
756  else if (.not. merid) then
757 * Now point1 and point2 belong within a hemisphere bounded by a
758 * meridian and geodesic is neither meridional or equatorial.
759 
760 * Figure a starting point for Newton's method
761  sig12 = invsta(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
762  + slam12, clam12, f, a3x, salp1, calp1, salp2, calp2, dnm, ca)
763 
764  if (sig12 .ge. 0) then
765 * Short lines (InvSta sets salp2, calp2, dnm)
766  s12x = sig12 * b * dnm
767  m12x = dnm**2 * b * sin(sig12 / dnm)
768  if (scalp) then
769  mm12 = cos(sig12 / dnm)
770  mm21 = mm12
771  end if
772  a12x = sig12 / degree
773  omg12 = lam12 / (f1 * dnm)
774  else
775 
776 * Newton's method. This is a straightforward solution of f(alp1) =
777 * lambda12(alp1) - lam12 = 0 with one wrinkle. f(alp) has exactly one
778 * root in the interval (0, pi) and its derivative is positive at the
779 * root. Thus f(alp) is positive for alp > alp1 and negative for alp <
780 * alp1. During the course of the iteration, a range (alp1a, alp1b) is
781 * maintained which brackets the root and with each evaluation of
782 * f(alp) the range is shrunk, if possible. Newton's method is
783 * restarted whenever the derivative of f is negative (because the new
784 * value of alp1 is then further from the solution) or if the new
785 * estimate of alp1 lies outside (0,pi); in this case, the new starting
786 * guess is taken to be (alp1a + alp1b) / 2.
787 
788 * Bracketing range
789  salp1a = tiny
790  calp1a = 1
791  salp1b = tiny
792  calp1b = -1
793  tripn = .false.
794  tripb = .false.
795  do 10 numit = 0, maxit2-1
796 * the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
797 * WGS84 and random input: mean = 2.85, sd = 0.60
798  v = lam12f(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
799  + salp1, calp1, slam12, clam12, f, a3x, c3x, salp2, calp2,
800  + sig12, ssig1, csig1, ssig2, csig2,
801  + eps, domg12, numit .lt. maxit1, dv, ca)
802 * Reversed test to allow escape with NaNs
803  if (tripn) then
804  dummy = 8
805  else
806  dummy = 1
807  end if
808  if (tripb .or. .not. (abs(v) .ge. dummy * tol0))
809  + go to 20
810 * Update bracketing values
811  if (v .gt. 0 .and. (numit .gt. maxit1 .or.
812  + calp1/salp1 .gt. calp1b/salp1b)) then
813  salp1b = salp1
814  calp1b = calp1
815  else if (v .lt. 0 .and. (numit .gt. maxit1 .or.
816  + calp1/salp1 .lt. calp1a/salp1a)) then
817  salp1a = salp1
818  calp1a = calp1
819  end if
820  if (numit .lt. maxit1 .and. dv .gt. 0) then
821  dalp1 = -v/dv
822  sdalp1 = sin(dalp1)
823  cdalp1 = cos(dalp1)
824  nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
825  if (nsalp1 .gt. 0 .and. abs(dalp1) .lt. pi) then
826  calp1 = calp1 * cdalp1 - salp1 * sdalp1
827  salp1 = nsalp1
828  call norm2x(salp1, calp1)
829 * In some regimes we don't get quadratic convergence because
830 * slope -> 0. So use convergence conditions based on dbleps
831 * instead of sqrt(dbleps).
832  tripn = abs(v) .le. 16 * tol0
833  go to 10
834  end if
835  end if
836 * Either dv was not positive or updated value was outside legal
837 * range. Use the midpoint of the bracket as the next estimate.
838 * This mechanism is not needed for the WGS84 ellipsoid, but it does
839 * catch problems with more eccentric ellipsoids. Its efficacy is
840 * such for the WGS84 test set with the starting guess set to alp1 =
841 * 90deg:
842 * the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
843 * WGS84 and random input: mean = 4.74, sd = 0.99
844  salp1 = (salp1a + salp1b)/2
845  calp1 = (calp1a + calp1b)/2
846  call norm2x(salp1, calp1)
847  tripn = .false.
848  tripb = abs(salp1a - salp1) + (calp1a - calp1) .lt. tolb
849  + .or. abs(salp1 - salp1b) + (calp1 - calp1b) .lt. tolb
850  10 continue
851  20 continue
852  call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
853  + cbet1, cbet2, lmask,
854  + s12x, m12x, dummy, mm12, mm21, ep2, ca)
855  m12x = m12x * b
856  s12x = s12x * b
857  a12x = sig12 / degree
858  if (areap) then
859  sdomg12 = sin(domg12)
860  cdomg12 = cos(domg12)
861  somg12 = slam12 * cdomg12 - clam12 * sdomg12
862  comg12 = clam12 * cdomg12 + slam12 * sdomg12
863  end if
864  end if
865  end if
866 
867 * Convert -0 to 0
868  s12 = 0 + s12x
869  if (redlp) m12 = 0 + m12x
870 
871  if (areap) then
872 * From Lambda12: sin(alp1) * cos(bet1) = sin(alp0)
873  salp0 = salp1 * cbet1
874  calp0 = hypotx(calp1, salp1 * sbet1)
875  if (calp0 .ne. 0 .and. salp0 .ne. 0) then
876 * From Lambda12: tan(bet) = tan(sig) * cos(alp)
877  ssig1 = sbet1
878  csig1 = calp1 * cbet1
879  ssig2 = sbet2
880  csig2 = calp2 * cbet2
881  k2 = calp0**2 * ep2
882  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
883 * Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0).
884  a4 = a**2 * calp0 * salp0 * e2
885  call norm2x(ssig1, csig1)
886  call norm2x(ssig2, csig2)
887  call c4f(eps, c4x, ca)
888  b41 = trgsum(.false., ssig1, csig1, ca, nc4)
889  b42 = trgsum(.false., ssig2, csig2, ca, nc4)
890  ss12 = a4 * (b42 - b41)
891  else
892 * Avoid problems with indeterminate sig1, sig2 on equator
893  ss12 = 0
894  end if
895 
896  if (.not. merid .and. somg12 .gt. 1) then
897  somg12 = sin(omg12)
898  comg12 = cos(omg12)
899  end if
900 
901  if (.not. merid .and. comg12 .ge. 0.7071d0
902  + .and. sbet2 - sbet1 .lt. 1.75d0) then
903 * Use tan(Gamma/2) = tan(omg12/2)
904 * * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
905 * with tan(x/2) = sin(x)/(1+cos(x))
906  domg12 = 1 + comg12
907  dbet1 = 1 + cbet1
908  dbet2 = 1 + cbet2
909  alp12 = 2 * atan2(somg12 * (sbet1 * dbet2 + sbet2 * dbet1),
910  + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
911  else
912 * alp12 = alp2 - alp1, used in atan2 so no need to normalize
913  salp12 = salp2 * calp1 - calp2 * salp1
914  calp12 = calp2 * calp1 + salp2 * salp1
915 * The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
916 * salp12 = -0 and alp12 = -180. However this depends on the sign
917 * being attached to 0 correctly. The following ensures the correct
918 * behavior.
919  if (salp12 .eq. 0 .and. calp12 .lt. 0) then
920  salp12 = tiny * calp1
921  calp12 = -1
922  end if
923  alp12 = atan2(salp12, calp12)
924  end if
925  ss12 = ss12 + c2 * alp12
926  ss12 = ss12 * swapp * lonsgn * latsgn
927 * Convert -0 to 0
928  ss12 = 0 + ss12
929  end if
930 
931 * Convert calp, salp to azimuth accounting for lonsgn, swapp, latsgn.
932  if (swapp .lt. 0) then
933  call swap(salp1, salp2)
934  call swap(calp1, calp2)
935  if (scalp) call swap(mm12, mm21)
936  end if
937 
938  salp1 = salp1 * swapp * lonsgn
939  calp1 = calp1 * swapp * latsgn
940  salp2 = salp2 * swapp * lonsgn
941  calp2 = calp2 * swapp * latsgn
942 
943  azi1 = atn2dx(salp1, calp1)
944  azi2 = atn2dx(salp2, calp2)
945 
946  if (arcp) a12 = a12x
947 
948  return
949  end
950 
951 *> Determine the area of a geodesic polygon
952 *!
953 *! @param[in] a the equatorial radius (meters).
954 *! @param[in] f the flattening of the ellipsoid. Setting \e f = 0 gives
955 *! a sphere. Negative \e f gives a prolate ellipsoid.
956 *! @param[in] lats an array of the latitudes of the vertices (degrees).
957 *! @param[in] lons an array of the longitudes of the vertices (degrees).
958 *! @param[in] n the number of vertices.
959 *! @param[out] AA the (signed) area of the polygon (meters<sup>2</sup>).
960 *! @param[out] PP the perimeter of the polygon.
961 *!
962 *! \e lats should be in the range [&minus;90&deg;, 90&deg;].
963 *!
964 *! Arbitrarily complex polygons are allowed. In the case of
965 *! self-intersecting polygons the area is accumulated "algebraically",
966 *! e.g., the areas of the 2 loops in a figure-8 polygon will partially
967 *! cancel. There's no need to "close" the polygon by repeating the
968 *! first vertex. The area returned is signed with counter-clockwise
969 *! traversal being treated as positive.
970 
971  subroutine area(a, f, lats, lons, n, AA, PP)
972 * input
973  integer n
974  double precision a, f, lats(n), lons(n)
975 * output
976  double precision AA, PP
977 
978  integer i, omask, cross, trnsit
979  double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
980  + atanhx, aacc(2), pacc(2)
981 
982  double precision dblmin, dbleps, pi, degree, tiny,
983  + tol0, tol1, tol2, tolb, xthrsh
984  integer digits, maxit1, maxit2
985  logical init
986  common /geocom/ dblmin, dbleps, pi, degree, tiny,
987  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
988 
989  omask = 8
990  call accini(aacc)
991  call accini(pacc)
992  cross = 0
993  do 10 i = 0, n-1
994  call invers(a, f, lats(i+1), lons(i+1),
995  + lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
996  + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, ss12)
997  call accadd(pacc, s12)
998  call accadd(aacc, -ss12)
999  cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
1000  10 continue
1001  pp = pacc(1)
1002  b = a * (1 - f)
1003  e2 = f * (2 - f)
1004  if (e2 .eq. 0) then
1005  c2 = a**2
1006  else if (e2 .gt. 0) then
1007  c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
1008  else
1009  c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
1010  end if
1011  area0 = 4 * pi * c2
1012  call accrem(aacc, area0)
1013  if (mod(abs(cross), 2) .eq. 1) then
1014  if (aacc(1) .lt. 0) then
1015  call accadd(aacc, +area0/2)
1016  else
1017  call accadd(aacc, -area0/2)
1018  end if
1019  end if
1020  if (aacc(1) .gt. area0/2) then
1021  call accadd(aacc, -area0)
1022  else if (aacc(1) .le. -area0/2) then
1023  call accadd(aacc, +area0)
1024  end if
1025  aa = aacc(1)
1026 
1027  return
1028  end
1029 
1030 *> Return the version numbers for this package.
1031 *!
1032 *! @param[out] major the major version number.
1033 *! @param[out] minor the minor version number.
1034 *! @param[out] patch the patch number.
1035 *!
1036 *! This subroutine was added with version 1.44.
1037 
1038  subroutine geover(major, minor, patch)
1039 * output
1040  integer major, minor, patch
1041 
1042  major = 1
1043  minor = 50
1044  patch = 0
1045 
1046  return
1047  end
1048 
1049 *> @cond SKIP
1050 
1051  block data geodat
1052  double precision dblmin, dbleps, pi, degree, tiny,
1053  + tol0, tol1, tol2, tolb, xthrsh
1054  integer digits, maxit1, maxit2
1055  logical init
1056  data init /.false./
1057  common /geocom/ dblmin, dbleps, pi, degree, tiny,
1058  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1059  end
1060 
1061  subroutine geoini
1062  double precision dblmin, dbleps, pi, degree, tiny,
1063  + tol0, tol1, tol2, tolb, xthrsh
1064  integer digits, maxit1, maxit2
1065  logical init
1066  common /geocom/ dblmin, dbleps, pi, degree, tiny,
1067  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1068 
1069  digits = 53
1070  dblmin = 0.5d0**1022
1071  dbleps = 0.5d0**(digits-1)
1072 
1073  pi = atan2(0d0, -1d0)
1074  degree = pi/180
1075 * This is about cbrt(dblmin). With other implementations, sqrt(dblmin)
1076 * is used. The larger value is used here to avoid complaints about a
1077 * IEEE_UNDERFLOW_FLAG IEEE_DENORMAL signal. This is triggered when
1078 * invers is called with points at opposite poles.
1079  tiny = 0.5d0**((1022+1)/3)
1080  tol0 = dbleps
1081 * Increase multiplier in defn of tol1 from 100 to 200 to fix inverse
1082 * case 52.784459512564 0 -52.784459512563990912 179.634407464943777557
1083 * which otherwise failed for Visual Studio 10 (Release and Debug)
1084  tol1 = 200 * tol0
1085  tol2 = sqrt(tol0)
1086 * Check on bisection interval
1087  tolb = tol0 * tol2
1088  xthrsh = 1000 * tol2
1089  maxit1 = 20
1090  maxit2 = maxit1 + digits + 10
1091 
1092  init = .true.
1093 
1094  return
1095  end
1096 
1097  subroutine lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1098  + cbet1, cbet2, omask,
1099  + s12b, m12b, m0, MM12, MM21, ep2, Ca)
1100 * input
1101  double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1102  + cbet1, cbet2, ep2
1103  integer omask
1104 * optional output
1105  double precision s12b, m12b, m0, MM12, MM21
1106 * temporary storage
1107  double precision Ca(*)
1108 
1109  integer ord, nC1, nC2
1110  parameter (ord = 6, nc1 = ord, nc2 = ord)
1111 
1112  double precision A1m1f, A2m1f, TrgSum
1113  double precision m0x, J12, A1, A2, B1, B2, csig12, t, Cb(nC2)
1114  logical distp, redlp, scalp
1115  integer l
1116 
1117 * Return m12b = (reduced length)/b; also calculate s12b = distance/b,
1118 * and m0 = coefficient of secular term in expression for reduced length.
1119 
1120  distp = (mod(omask/16, 2) .eq. 1)
1121  redlp = (mod(omask/2, 2) .eq. 1)
1122  scalp = (mod(omask/4, 2) .eq. 1)
1123 
1124 * Suppress compiler warnings
1125  m0x = 0
1126  j12 = 0
1127  a1 = 0
1128  a2 = 0
1129  if (distp .or. redlp .or. scalp) then
1130  a1 = a1m1f(eps)
1131  call c1f(eps, ca)
1132  if (redlp .or. scalp) then
1133  a2 = a2m1f(eps)
1134  call c2f(eps, cb)
1135  m0x = a1 - a2
1136  a2 = 1 + a2
1137  end if
1138  a1 = 1 + a1
1139  end if
1140  if (distp) then
1141  b1 = trgsum(.true., ssig2, csig2, ca, nc1) -
1142  + trgsum(.true., ssig1, csig1, ca, nc1)
1143 * Missing a factor of b
1144  s12b = a1 * (sig12 + b1)
1145  if (redlp .or. scalp) then
1146  b2 = trgsum(.true., ssig2, csig2, cb, nc2) -
1147  + trgsum(.true., ssig1, csig1, cb, nc2)
1148  j12 = m0x * sig12 + (a1 * b1 - a2 * b2)
1149  end if
1150  else if (redlp .or. scalp) then
1151 * Assume here that nC1 >= nC2
1152  do 10 l = 1, nc2
1153  cb(l) = a1 * ca(l) - a2 * cb(l)
1154  10 continue
1155  j12 = m0x * sig12 + (trgsum(.true., ssig2, csig2, cb, nc2) -
1156  + trgsum(.true., ssig1, csig1, cb, nc2))
1157  end if
1158  if (redlp) then
1159  m0 = m0x
1160 * Missing a factor of b.
1161 * Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
1162 * accurate cancellation in the case of coincident points.
1163  m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1164  + csig1 * csig2 * j12
1165  end if
1166  if (scalp) then
1167  csig12 = csig1 * csig2 + ssig1 * ssig2
1168  t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
1169  mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
1170  mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
1171  end if
1172 
1173  return
1174  end
1175 
1176  double precision function astrd(x, y)
1177 * Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
1178 * This solution is adapted from Geocentric::Reverse.
1179 * input
1180  double precision x, y
1181 
1182  double precision cbrt
1183  double precision k, p, q, r, S, r2, r3, disc, u,
1184  + t3, t, ang, v, uv, w
1185 
1186  p = x**2
1187  q = y**2
1188  r = (p + q - 1) / 6
1189  if ( .not. (q .eq. 0 .and. r .lt. 0) ) then
1190 * Avoid possible division by zero when r = 0 by multiplying equations
1191 * for s and t by r^3 and r, resp.
1192 * S = r^3 * s
1193  s = p * q / 4
1194  r2 = r**2
1195  r3 = r * r2
1196 * The discriminant of the quadratic equation for T3. This is zero on
1197 * the evolute curve p^(1/3)+q^(1/3) = 1
1198  disc = s * (s + 2 * r3)
1199  u = r
1200  if (disc .ge. 0) then
1201  t3 = s + r3
1202 * Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
1203 * of precision due to cancellation. The result is unchanged because
1204 * of the way the T is used in definition of u.
1205 * T3 = (r * t)^3
1206  if (t3 .lt. 0) then
1207  disc = -sqrt(disc)
1208  else
1209  disc = sqrt(disc)
1210  end if
1211  t3 = t3 + disc
1212 * N.B. cbrt always returns the real root. cbrt(-8) = -2.
1213 * T = r * t
1214  t = cbrt(t3)
1215 * T can be zero; but then r2 / T -> 0.
1216  if (t .ne. 0) u = u + t + r2 / t
1217  else
1218 * T is complex, but the way u is defined the result is real.
1219  ang = atan2(sqrt(-disc), -(s + r3))
1220 * There are three possible cube roots. We choose the root which
1221 * avoids cancellation. Note that disc < 0 implies that r < 0.
1222  u = u + 2 * r * cos(ang / 3)
1223  end if
1224 * guaranteed positive
1225  v = sqrt(u**2 + q)
1226 * Avoid loss of accuracy when u < 0.
1227 * u+v, guaranteed positive
1228  if (u .lt. 0) then
1229  uv = q / (v - u)
1230  else
1231  uv = u + v
1232  end if
1233 * positive?
1234  w = (uv - q) / (2 * v)
1235 * Rearrange expression for k to avoid loss of accuracy due to
1236 * subtraction. Division by 0 not possible because uv > 0, w >= 0.
1237 * guaranteed positive
1238  k = uv / (sqrt(uv + w**2) + w)
1239  else
1240 * q == 0 && r <= 0
1241 * y = 0 with |x| <= 1. Handle this case directly.
1242 * for y small, positive root is k = abs(y)/sqrt(1-x^2)
1243  k = 0
1244  end if
1245  astrd = k
1246 
1247  return
1248  end
1249 
1250  double precision function invsta(sbet1, cbet1, dn1,
1251  + sbet2, cbet2, dn2, lam12, slam12, clam12, f, A3x,
1252  + salp1, calp1, salp2, calp2, dnm,
1253  + Ca)
1254 * Return a starting point for Newton's method in salp1 and calp1
1255 * (function value is -1). If Newton's method doesn't need to be used,
1256 * return also salp2, calp2, and dnm and function value is sig12.
1257 * input
1258  double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1259  + lam12, slam12, clam12, f, A3x(*)
1260 * output
1261  double precision salp1, calp1, salp2, calp2, dnm
1262 * temporary
1263  double precision Ca(*)
1264 
1265  double precision hypotx, A3f, Astrd
1266  logical shortp
1267  double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
1268  + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
1269  + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
1270  + k, omg12a, sbetm2, lam12x
1271 
1272  double precision dblmin, dbleps, pi, degree, tiny,
1273  + tol0, tol1, tol2, tolb, xthrsh
1274  integer digits, maxit1, maxit2
1275  logical init
1276  common /geocom/ dblmin, dbleps, pi, degree, tiny,
1277  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1278 
1279  f1 = 1 - f
1280  e2 = f * (2 - f)
1281  ep2 = e2 / (1 - e2)
1282  n = f / (2 - f)
1283 * The sig12 threshold for "really short". Using the auxiliary sphere
1284 * solution with dnm computed at (bet1 + bet2) / 2, the relative error in
1285 * the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
1286 * (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a
1287 * given f and sig12, the max error occurs for lines near the pole. If
1288 * the old rule for computing dnm = (dn1 + dn2)/2 is used, then the error
1289 * increases by a factor of 2.) Setting this equal to epsilon gives
1290 * sig12 = etol2. Here 0.1 is a safety factor (error decreased by 100)
1291 * and max(0.001, abs(f)) stops etol2 getting too large in the nearly
1292 * spherical case.
1293  etol2 = 0.1d0 * tol2 /
1294  + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
1295 
1296 * Return value
1297  sig12 = -1
1298 * bet12 = bet2 - bet1 in [0, pi); bt12a = bet2 + bet1 in (-pi, 0]
1299  sbet12 = sbet2 * cbet1 - cbet2 * sbet1
1300  cbet12 = cbet2 * cbet1 + sbet2 * sbet1
1301  sbt12a = sbet2 * cbet1 + cbet2 * sbet1
1302 
1303  shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
1304  + cbet2 * lam12 .lt. 0.5d0
1305 
1306  if (shortp) then
1307  sbetm2 = (sbet1 + sbet2)**2
1308 * sin((bet1+bet2)/2)^2
1309 * = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
1310  sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
1311  dnm = sqrt(1 + ep2 * sbetm2)
1312  omg12 = lam12 / (f1 * dnm)
1313  somg12 = sin(omg12)
1314  comg12 = cos(omg12)
1315  else
1316  somg12 = slam12
1317  comg12 = clam12
1318  end if
1319 
1320  salp1 = cbet2 * somg12
1321  if (comg12 .ge. 0) then
1322  calp1 = sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12)
1323  else
1324  calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1325  end if
1326 
1327  ssig12 = hypotx(salp1, calp1)
1328  csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
1329 
1330  if (shortp .and. ssig12 .lt. etol2) then
1331 * really short lines
1332  salp2 = cbet1 * somg12
1333  if (comg12 .ge. 0) then
1334  calp2 = somg12**2 / (1 + comg12)
1335  else
1336  calp2 = 1 - comg12
1337  end if
1338  calp2 = sbet12 - cbet1 * sbet2 * calp2
1339  call norm2x(salp2, calp2)
1340 * Set return value
1341  sig12 = atan2(ssig12, csig12)
1342  else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
1343  + ssig12 .ge. 6 * abs(n) * pi * cbet1**2) then
1344 * Nothing to do, zeroth order spherical approximation is OK
1345  continue
1346  else
1347 * lam12 - pi
1348  lam12x = atan2(-slam12, -clam12)
1349 * Scale lam12 and bet2 to x, y coordinate system where antipodal point
1350 * is at origin and singular point is at y = 0, x = -1.
1351  if (f .ge. 0) then
1352 * x = dlong, y = dlat
1353  k2 = sbet1**2 * ep2
1354  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1355  lamscl = f * cbet1 * a3f(eps, a3x) * pi
1356  betscl = lamscl * cbet1
1357  x = lam12x / lamscl
1358  y = sbt12a / betscl
1359  else
1360 * f < 0: x = dlat, y = dlong
1361  cbt12a = cbet2 * cbet1 - sbet2 * sbet1
1362  bt12a = atan2(sbt12a, cbt12a)
1363 * In the case of lon12 = 180, this repeats a calculation made in
1364 * Inverse.
1365  call lengs(n, pi + bt12a,
1366  + sbet1, -cbet1, dn1, sbet2, cbet2, dn2, cbet1, cbet2, 2,
1367  + dummy, m12b, m0, dummy, dummy, ep2, ca)
1368  x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
1369  if (x .lt. -0.01d0) then
1370  betscl = sbt12a / x
1371  else
1372  betscl = -f * cbet1**2 * pi
1373  end if
1374  lamscl = betscl / cbet1
1375  y = lam12x / lamscl
1376  end if
1377 
1378  if (y .gt. -tol1 .and. x .gt. -1 - xthrsh) then
1379 * strip near cut
1380  if (f .ge. 0) then
1381  salp1 = min(1d0, -x)
1382  calp1 = - sqrt(1 - salp1**2)
1383  else
1384  if (x .gt. -tol1) then
1385  calp1 = 0
1386  else
1387  calp1 = 1
1388  end if
1389  calp1 = max(calp1, x)
1390  salp1 = sqrt(1 - calp1**2)
1391  end if
1392  else
1393 * Estimate alp1, by solving the astroid problem.
1394 *
1395 * Could estimate alpha1 = theta + pi/2, directly, i.e.,
1396 * calp1 = y/k; salp1 = -x/(1+k); for f >= 0
1397 * calp1 = x/(1+k); salp1 = -y/k; for f < 0 (need to check)
1398 *
1399 * However, it's better to estimate omg12 from astroid and use
1400 * spherical formula to compute alp1. This reduces the mean number of
1401 * Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
1402 * (min 0 max 5). The changes in the number of iterations are as
1403 * follows:
1404 *
1405 * change percent
1406 * 1 5
1407 * 0 78
1408 * -1 16
1409 * -2 0.6
1410 * -3 0.04
1411 * -4 0.002
1412 *
1413 * The histogram of iterations is (m = number of iterations estimating
1414 * alp1 directly, n = number of iterations estimating via omg12, total
1415 * number of trials = 148605):
1416 *
1417 * iter m n
1418 * 0 148 186
1419 * 1 13046 13845
1420 * 2 93315 102225
1421 * 3 36189 32341
1422 * 4 5396 7
1423 * 5 455 1
1424 * 6 56 0
1425 *
1426 * Because omg12 is near pi, estimate work with omg12a = pi - omg12
1427  k = astrd(x, y)
1428  if (f .ge. 0) then
1429  omg12a = -x * k/(1 + k)
1430  else
1431  omg12a = -y * (1 + k)/k
1432  end if
1433  omg12a = lamscl * omg12a
1434  somg12 = sin(omg12a)
1435  comg12 = -cos(omg12a)
1436 * Update spherical estimate of alp1 using omg12 instead of lam12
1437  salp1 = cbet2 * somg12
1438  calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1439  end if
1440  end if
1441 * Sanity check on starting guess. Backwards check allows NaN through.
1442  if (.not. (salp1 .le. 0)) then
1443  call norm2x(salp1, calp1)
1444  else
1445  salp1 = 1
1446  calp1 = 0
1447  end if
1448  invsta = sig12
1449 
1450  return
1451  end
1452 
1453  double precision function lam12f(sbet1, cbet1, dn1,
1454  + sbet2, cbet2, dn2, salp1, calp1, slm120, clm120, f, A3x, C3x,
1455  + salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps,
1456  + domg12, diffp, dlam12, Ca)
1457 * input
1458  double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1459  + salp1, calp1, slm120, clm120, f, A3x(*), C3x(*)
1460  logical diffp
1461 * output
1462  double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1463  + eps, domg12
1464 * optional output
1465  double precision dlam12
1466 * temporary
1467  double precision Ca(*)
1468 
1469  integer ord, nC3
1470  parameter(ord = 6, nc3 = ord)
1471 
1472  double precision hypotx, A3f, TrgSum
1473 
1474  double precision f1, e2, ep2, salp0, calp0,
1475  + somg1, comg1, somg2, comg2, somg12, comg12,
1476  + lam12, eta, b312, k2, dummy
1477 
1478  double precision dblmin, dbleps, pi, degree, tiny,
1479  + tol0, tol1, tol2, tolb, xthrsh
1480  integer digits, maxit1, maxit2
1481  logical init
1482  common /geocom/ dblmin, dbleps, pi, degree, tiny,
1483  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1484 
1485  f1 = 1 - f
1486  e2 = f * (2 - f)
1487  ep2 = e2 / (1 - e2)
1488 * Break degeneracy of equatorial line. This case has already been
1489 * handled.
1490  if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
1491 
1492 * sin(alp1) * cos(bet1) = sin(alp0)
1493  salp0 = salp1 * cbet1
1494 * calp0 > 0
1495  calp0 = hypotx(calp1, salp1 * sbet1)
1496 
1497 * tan(bet1) = tan(sig1) * cos(alp1)
1498 * tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
1499  ssig1 = sbet1
1500  somg1 = salp0 * sbet1
1501  csig1 = calp1 * cbet1
1502  comg1 = csig1
1503  call norm2x(ssig1, csig1)
1504 * norm2x(somg1, comg1); -- don't need to normalize!
1505 
1506 * Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful
1507 * about this case, since this can yield singularities in the Newton
1508 * iteration.
1509 * sin(alp2) * cos(bet2) = sin(alp0)
1510  if (cbet2 .ne. cbet1) then
1511  salp2 = salp0 / cbet2
1512  else
1513  salp2 = salp1
1514  end if
1515 * calp2 = sqrt(1 - sq(salp2))
1516 * = sqrt(sq(calp0) - sq(sbet2)) / cbet2
1517 * and subst for calp0 and rearrange to give (choose positive sqrt
1518 * to give alp2 in [0, pi/2]).
1519  if (cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1) then
1520  if (cbet1 .lt. -sbet1) then
1521  calp2 = (cbet2 - cbet1) * (cbet1 + cbet2)
1522  else
1523  calp2 = (sbet1 - sbet2) * (sbet1 + sbet2)
1524  end if
1525  calp2 = sqrt((calp1 * cbet1)**2 + calp2) / cbet2
1526  else
1527  calp2 = abs(calp1)
1528  end if
1529 * tan(bet2) = tan(sig2) * cos(alp2)
1530 * tan(omg2) = sin(alp0) * tan(sig2).
1531  ssig2 = sbet2
1532  somg2 = salp0 * sbet2
1533  csig2 = calp2 * cbet2
1534  comg2 = csig2
1535  call norm2x(ssig2, csig2)
1536 * norm2x(somg2, comg2); -- don't need to normalize!
1537 
1538 * sig12 = sig2 - sig1, limit to [0, pi]
1539  sig12 = atan2(0d0 + max(0d0, csig1 * ssig2 - ssig1 * csig2),
1540  + csig1 * csig2 + ssig1 * ssig2)
1541 
1542 * omg12 = omg2 - omg1, limit to [0, pi]
1543  somg12 = 0d0 + max(0d0, comg1 * somg2 - somg1 * comg2)
1544  comg12 = comg1 * comg2 + somg1 * somg2
1545 * eta = omg12 - lam120
1546  eta = atan2(somg12 * clm120 - comg12 * slm120,
1547  + comg12 * clm120 + somg12 * slm120)
1548  k2 = calp0**2 * ep2
1549  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1550  call c3f(eps, c3x, ca)
1551  b312 = (trgsum(.true., ssig2, csig2, ca, nc3-1) -
1552  + trgsum(.true., ssig1, csig1, ca, nc3-1))
1553  domg12 = -f * a3f(eps, a3x) * salp0 * (sig12 + b312)
1554  lam12 = eta + domg12
1555 
1556  if (diffp) then
1557  if (calp2 .eq. 0) then
1558  dlam12 = - 2 * f1 * dn1 / sbet1
1559  else
1560  call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1561  + cbet1, cbet2, 2,
1562  + dummy, dlam12, dummy, dummy, dummy, ep2, ca)
1563  dlam12 = dlam12 * f1 / (calp2 * cbet2)
1564  end if
1565  end if
1566  lam12f = lam12
1567 
1568  return
1569  end
1570 
1571  double precision function a3f(eps, A3x)
1572 * Evaluate A3
1573  integer ord, nA3, nA3x
1574  parameter(ord = 6, na3 = ord, na3x = na3)
1575 
1576 * input
1577  double precision eps
1578 * output
1579  double precision A3x(0: nA3x-1)
1580 
1581  double precision polval
1582  A3f = polval(na3 - 1, a3x, eps)
1583 
1584  return
1585  end
1586 
1587  subroutine c3f(eps, C3x, c)
1588 * Evaluate C3 coeffs
1589 * Elements c[1] thru c[nC3-1] are set
1590  integer ord, nC3, nC3x
1591  parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1592 
1593 * input
1594  double precision eps, C3x(0:nC3x-1)
1595 * output
1596  double precision c(nC3-1)
1597 
1598  integer o, m, l
1599  double precision mult, polval
1600 
1601  mult = 1
1602  o = 0
1603  do 10 l = 1, nc3 - 1
1604  m = nc3 - l - 1
1605  mult = mult * eps
1606  c(l) = mult * polval(m, c3x(o), eps)
1607  o = o + m + 1
1608  10 continue
1609 
1610  return
1611  end
1612 
1613  subroutine c4f(eps, C4x, c)
1614 * Evaluate C4
1615 * Elements c[0] thru c[nC4-1] are set
1616  integer ord, nC4, nC4x
1617  parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1618 
1619 * input
1620  double precision eps, C4x(0:nC4x-1)
1621 *output
1622  double precision c(0:nC4-1)
1623 
1624  integer o, m, l
1625  double precision mult, polval
1626 
1627  mult = 1
1628  o = 0
1629  do 10 l = 0, nc4 - 1
1630  m = nc4 - l - 1
1631  c(l) = mult * polval(m, c4x(o), eps)
1632  o = o + m + 1
1633  mult = mult * eps
1634  10 continue
1635 
1636  return
1637  end
1638 
1639  double precision function a1m1f(eps)
1640 * The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
1641 * input
1642  double precision eps
1643 
1644  double precision t
1645  integer ord, nA1, o, m
1646  parameter(ord = 6, na1 = ord)
1647  double precision polval, coeff(nA1/2 + 2)
1648  data coeff /1, 4, 64, 0, 256/
1649 
1650  o = 1
1651  m = na1/2
1652  t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1653  a1m1f = (t + eps) / (1 - eps)
1654 
1655  return
1656  end
1657 
1658  subroutine c1f(eps, c)
1659 * The coefficients C1[l] in the Fourier expansion of B1
1660  integer ord, nC1
1661  parameter(ord = 6, nc1 = ord)
1662 
1663 * input
1664  double precision eps
1665 * output
1666  double precision c(nC1)
1667 
1668  double precision eps2, d
1669  integer o, m, l
1670  double precision polval, coeff((nC1**2 + 7*nC1 - 2*(nC1/2))/4)
1671  data coeff /
1672  + -1, 6, -16, 32,
1673  + -9, 64, -128, 2048,
1674  + 9, -16, 768,
1675  + 3, -5, 512,
1676  + -7, 1280,
1677  + -7, 2048/
1678 
1679  eps2 = eps**2
1680  d = eps
1681  o = 1
1682  do 10 l = 1, nc1
1683  m = (nc1 - l) / 2
1684  c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1685  o = o + m + 2
1686  d = d * eps
1687  10 continue
1688 
1689  return
1690  end
1691 
1692  subroutine c1pf(eps, c)
1693 * The coefficients C1p[l] in the Fourier expansion of B1p
1694  integer ord, nC1p
1695  parameter(ord = 6, nc1p = ord)
1696 
1697 * input
1698  double precision eps
1699 * output
1700  double precision c(nC1p)
1701 
1702  double precision eps2, d
1703  integer o, m, l
1704  double precision polval, coeff((nC1p**2 + 7*nC1p - 2*(nC1p/2))/4)
1705  data coeff /
1706  + 205, -432, 768, 1536,
1707  + 4005, -4736, 3840, 12288,
1708  + -225, 116, 384,
1709  + -7173, 2695, 7680,
1710  + 3467, 7680,
1711  + 38081, 61440/
1712 
1713  eps2 = eps**2
1714  d = eps
1715  o = 1
1716  do 10 l = 1, nc1p
1717  m = (nc1p - l) / 2
1718  c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1719  o = o + m + 2
1720  d = d * eps
1721  10 continue
1722 
1723  return
1724  end
1725 
1726 * The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
1727  double precision function a2m1f(eps)
1728 * input
1729  double precision eps
1730 
1731  double precision t
1732  integer ord, nA2, o, m
1733  parameter(ord = 6, na2 = ord)
1734  double precision polval, coeff(nA2/2 + 2)
1735  data coeff /-11, -28, -192, 0, 256/
1736 
1737  o = 1
1738  m = na2/2
1739  t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1740  a2m1f = (t - eps) / (1 + eps)
1741 
1742  return
1743  end
1744 
1745  subroutine c2f(eps, c)
1746 * The coefficients C2[l] in the Fourier expansion of B2
1747  integer ord, nC2
1748  parameter(ord = 6, nc2 = ord)
1749 
1750 * input
1751  double precision eps
1752 * output
1753  double precision c(nC2)
1754 
1755  double precision eps2, d
1756  integer o, m, l
1757  double precision polval, coeff((nC2**2 + 7*nC2 - 2*(nC2/2))/4)
1758  data coeff /
1759  + 1, 2, 16, 32,
1760  + 35, 64, 384, 2048,
1761  + 15, 80, 768,
1762  + 7, 35, 512,
1763  + 63, 1280,
1764  + 77, 2048/
1765 
1766  eps2 = eps**2
1767  d = eps
1768  o = 1
1769  do 10 l = 1, nc2
1770  m = (nc2 - l) / 2
1771  c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1772  o = o + m + 2
1773  d = d * eps
1774  10 continue
1775 
1776  return
1777  end
1778 
1779  subroutine a3cof(n, A3x)
1780 * The scale factor A3 = mean value of (d/dsigma)I3
1781  integer ord, nA3, nA3x
1782  parameter(ord = 6, na3 = ord, na3x = na3)
1783 
1784 * input
1785  double precision n
1786 * output
1787  double precision A3x(0:nA3x-1)
1788 
1789  integer o, m, k, j
1790  double precision polval, coeff((nA3**2 + 7*nA3 - 2*(nA3/2))/4)
1791  data coeff /
1792  + -3, 128,
1793  + -2, -3, 64,
1794  + -1, -3, -1, 16,
1795  + 3, -1, -2, 8,
1796  + 1, -1, 2,
1797  + 1, 1/
1798 
1799  o = 1
1800  k = 0
1801  do 10 j = na3 - 1, 0, -1
1802  m = min(na3 - j - 1, j)
1803  a3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1804  k = k + 1
1805  o = o + m + 2
1806  10 continue
1807 
1808  return
1809  end
1810 
1811  subroutine c3cof(n, C3x)
1812 * The coefficients C3[l] in the Fourier expansion of B3
1813  integer ord, nC3, nC3x
1814  parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1815 
1816 * input
1817  double precision n
1818 * output
1819  double precision C3x(0:nC3x-1)
1820 
1821  integer o, m, l, j, k
1822  double precision polval,
1823  + coeff(((nc3-1)*(nc3**2 + 7*nc3 - 2*(nc3/2)))/8)
1824  data coeff /
1825  + 3, 128,
1826  + 2, 5, 128,
1827  + -1, 3, 3, 64,
1828  + -1, 0, 1, 8,
1829  + -1, 1, 4,
1830  + 5, 256,
1831  + 1, 3, 128,
1832  + -3, -2, 3, 64,
1833  + 1, -3, 2, 32,
1834  + 7, 512,
1835  + -10, 9, 384,
1836  + 5, -9, 5, 192,
1837  + 7, 512,
1838  + -14, 7, 512,
1839  + 21, 2560/
1840 
1841  o = 1
1842  k = 0
1843  do 20 l = 1, nc3 - 1
1844  do 10 j = nc3 - 1, l, -1
1845  m = min(nc3 - j - 1, j)
1846  c3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1847  k = k + 1
1848  o = o + m + 2
1849  10 continue
1850  20 continue
1851 
1852  return
1853  end
1854 
1855  subroutine c4cof(n, C4x)
1856 * The coefficients C4[l] in the Fourier expansion of I4
1857  integer ord, nC4, nC4x
1858  parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1859 
1860 * input
1861  double precision n
1862 * output
1863  double precision C4x(0:nC4x-1)
1864 
1865  integer o, m, l, j, k
1866  double precision polval, coeff((nC4 * (nC4 + 1) * (nC4 + 5)) / 6)
1867  data coeff /
1868  + 97, 15015, 1088, 156, 45045, -224, -4784, 1573, 45045,
1869  + -10656, 14144, -4576, -858, 45045,
1870  + 64, 624, -4576, 6864, -3003, 15015,
1871  + 100, 208, 572, 3432, -12012, 30030, 45045,
1872  + 1, 9009, -2944, 468, 135135, 5792, 1040, -1287, 135135,
1873  + 5952, -11648, 9152, -2574, 135135,
1874  + -64, -624, 4576, -6864, 3003, 135135,
1875  + 8, 10725, 1856, -936, 225225, -8448, 4992, -1144, 225225,
1876  + -1440, 4160, -4576, 1716, 225225,
1877  + -136, 63063, 1024, -208, 105105,
1878  + 3584, -3328, 1144, 315315,
1879  + -128, 135135, -2560, 832, 405405, 128, 99099/
1880 
1881  o = 1
1882  k = 0
1883  do 20 l = 0, nc4 - 1
1884  do 10 j = nc4 - 1, l, -1
1885  m = nc4 - j - 1
1886  c4x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1887  k = k + 1
1888  o = o + m + 2
1889  10 continue
1890  20 continue
1891 
1892  return
1893  end
1894 
1895  double precision function sumx(u, v, t)
1896 * input
1897  double precision u, v
1898 * output
1899  double precision t
1900 
1901  double precision up, vpp
1902  sumx = u + v
1903  up = sumx - v
1904  vpp = sumx - up
1905  up = up - u
1906  vpp = vpp - v
1907  t = -(up + vpp)
1908 
1909  return
1910  end
1911 
1912  double precision function remx(x, y)
1913 * the remainder function but not worrying how ties are handled
1914 * y must be positive
1915 * input
1916  double precision x, y
1917 
1918  remx = mod(x, y)
1919  if (remx .lt. -y/2) then
1920  remx = remx + y
1921  else if (remx .gt. +y/2) then
1922  remx = remx - y
1923  end if
1924 
1925  return
1926  end
1927 
1928  double precision function angnm(x)
1929 * input
1930  double precision x
1931 
1932  double precision remx
1933  angnm = remx(x, 360d0)
1934  if (angnm .eq. -180) then
1935  angnm = 180
1936  end if
1937 
1938  return
1939  end
1940 
1941  double precision function latfix(x)
1942 * input
1943  double precision x
1944 
1945  latfix = x
1946  if (.not. (abs(x) .gt. 90)) return
1947 * concoct a NaN
1948  latfix = sqrt(90 - abs(x))
1949 
1950  return
1951  end
1952 
1953  double precision function angdif(x, y, e)
1954 * Compute y - x. x and y must both lie in [-180, 180]. The result is
1955 * equivalent to computing the difference exactly, reducing it to (-180,
1956 * 180] and rounding the result. Note that this prescription allows -180
1957 * to be returned (e.g., if x is tiny and negative and y = 180). The
1958 * error in the difference is returned in e
1959 * input
1960  double precision x, y
1961 * output
1962  double precision e
1963 
1964  double precision d, t, sumx, AngNm
1965  d = angnm(sumx(angnm(-x), angnm(y), t))
1966  if (d .eq. 180 .and. t .gt. 0) then
1967  d = -180
1968  end if
1969  angdif = sumx(d, t, e)
1970 
1971  return
1972  end
1973 
1974  double precision function angrnd(x)
1975 * The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
1976 * for reals = 0.7 pm on the earth if x is an angle in degrees. (This is
1977 * about 1000 times more resolution than we get with angles around 90
1978 * degrees.) We use this to avoid having to deal with near singular
1979 * cases when x is non-zero but tiny (e.g., 1.0e-200). This also
1980 * converts -0 to +0.
1981 * input
1982  double precision x
1983 
1984  double precision y, z
1985  z = 1/16d0
1986  y = abs(x)
1987 * The compiler mustn't "simplify" z - (z - y) to y
1988  if (y .lt. z) y = z - (z - y)
1989  angrnd = sign(y, x)
1990  if (x .eq. 0) angrnd = 0
1991 
1992  return
1993  end
1994 
1995  subroutine swap(x, y)
1996 * input/output
1997  double precision x, y
1998 
1999  double precision z
2000  z = x
2001  x = y
2002  y = z
2003 
2004  return
2005  end
2006 
2007  double precision function hypotx(x, y)
2008 * input
2009  double precision x, y
2010 
2011 * With Fortran 2008, this becomes: hypotx = hypot(x, y)
2012  hypotx = sqrt(x**2 + y**2)
2013 
2014  return
2015  end
2016 
2017  subroutine norm2x(x, y)
2018 * input/output
2019  double precision x, y
2020 
2021  double precision hypotx, r
2022  r = hypotx(x, y)
2023  x = x/r
2024  y = y/r
2025 
2026  return
2027  end
2028 
2029  double precision function log1px(x)
2030 * input
2031  double precision x
2032 
2033  double precision y, z
2034  y = 1 + x
2035  z = y - 1
2036  if (z .eq. 0) then
2037  log1px = x
2038  else
2039  log1px = x * log(y) / z
2040  end if
2041 
2042  return
2043  end
2044 
2045  double precision function atanhx(x)
2046 * input
2047  double precision x
2048 
2049 * With Fortran 2008, this becomes: atanhx = atanh(x)
2050  double precision log1px, y
2051  y = abs(x)
2052  y = log1px(2 * y/(1 - y))/2
2053  atanhx = sign(y, x)
2054 
2055  return
2056  end
2057 
2058  double precision function cbrt(x)
2059 * input
2060  double precision x
2061 
2062  cbrt = sign(abs(x)**(1/3d0), x)
2063 
2064  return
2065  end
2066 
2067  double precision function trgsum(sinp, sinx, cosx, c, n)
2068 * Evaluate
2069 * y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
2070 * sum(c[i] * cos((2*i-1) * x), i, 1, n)
2071 * using Clenshaw summation.
2072 * Approx operation count = (n + 5) mult and (2 * n + 2) add
2073 * input
2074  logical sinp
2075  integer n
2076  double precision sinx, cosx, c(n)
2077 
2078  double precision ar, y0, y1
2079  integer n2, k
2080 
2081 * 2 * cos(2 * x)
2082  ar = 2 * (cosx - sinx) * (cosx + sinx)
2083 * accumulators for sum
2084  if (mod(n, 2) .eq. 1) then
2085  y0 = c(n)
2086  n2 = n - 1
2087  else
2088  y0 = 0
2089  n2 = n
2090  end if
2091  y1 = 0
2092 * Now n2 is even
2093  do 10 k = n2, 2, -2
2094 * Unroll loop x 2, so accumulators return to their original role
2095  y1 = ar * y0 - y1 + c(k)
2096  y0 = ar * y1 - y0 + c(k-1)
2097  10 continue
2098  if (sinp) then
2099 * sin(2 * x) * y0
2100  trgsum = 2 * sinx * cosx * y0
2101  else
2102 * cos(x) * (y0 - y1)
2103  trgsum = cosx * (y0 - y1)
2104  end if
2105 
2106  return
2107  end
2108 
2109  integer function trnsit(lon1, lon2)
2110 * input
2111  double precision lon1, lon2
2112 
2113  double precision lon1x, lon2x, lon12, AngNm, AngDif, e
2114  lon1x = angnm(lon1)
2115  lon2x = angnm(lon2)
2116  lon12 = angdif(lon1x, lon2x, e)
2117  trnsit = 0
2118  if (lon1x .le. 0 .and. lon2x .gt. 0 .and. lon12 .gt. 0) then
2119  trnsit = 1
2120  else if (lon2x .le. 0 .and. lon1x .gt. 0 .and. lon12 .lt. 0) then
2121  trnsit = -1
2122  end if
2123 
2124  return
2125  end
2126 
2127  subroutine accini(s)
2128 * Initialize an accumulator; this is an array with two elements.
2129 * input/output
2130  double precision s(2)
2131 
2132  s(1) = 0
2133  s(2) = 0
2134 
2135  return
2136  end
2137 
2138  subroutine accadd(s, y)
2139 * Add y to an accumulator.
2140 * input
2141  double precision y
2142 * input/output
2143  double precision s(2)
2144 
2145  double precision z, u, sumx
2146  z = sumx(y, s(2), u)
2147  s(1) = sumx(z, s(1), s(2))
2148  if (s(1) .eq. 0) then
2149  s(1) = u
2150  else
2151  s(2) = s(2) + u
2152  end if
2153 
2154  return
2155  end
2156 
2157  subroutine accrem(s, y)
2158 * Reduce s to [-y/2, y/2].
2159 * input
2160  double precision y
2161 * input/output
2162  double precision s(2)
2163 
2164  double precision remx
2165  s(1) = remx(s(1), y)
2166  call accadd(s, 0d0)
2167 
2168  return
2169  end
2170 
2171  subroutine sncsdx(x, sinx, cosx)
2172 * Compute sin(x) and cos(x) with x in degrees
2173 * input
2174  double precision x
2175 * input/output
2176  double precision sinx, cosx
2177 
2178  double precision dblmin, dbleps, pi, degree, tiny,
2179  + tol0, tol1, tol2, tolb, xthrsh
2180  integer digits, maxit1, maxit2
2181  logical init
2182  common /geocom/ dblmin, dbleps, pi, degree, tiny,
2183  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2184 
2185  double precision r, s, c
2186  integer q
2187  r = mod(x, 360d0)
2188  q = nint(r / 90)
2189  r = (r - 90 * q) * degree
2190  s = sin(r)
2191 * sin(-0) isn't reliably -0, so...
2192  if (x .eq. 0) s = x
2193  c = cos(r)
2194  q = mod(q + 4, 4)
2195  if (q .eq. 0) then
2196  sinx = s
2197  cosx = c
2198  else if (q .eq. 1) then
2199  sinx = c
2200  cosx = -s
2201  else if (q .eq. 2) then
2202  sinx = -s
2203  cosx = -c
2204  else
2205 * q.eq.3
2206  sinx = -c
2207  cosx = s
2208  end if
2209 
2210  if (x .ne. 0) then
2211  sinx = 0d0 + sinx
2212  cosx = 0d0 + cosx
2213  end if
2214 
2215  return
2216  end
2217 
2218  double precision function atn2dx(y, x)
2219 * input
2220  double precision x, y
2221 
2222  double precision dblmin, dbleps, pi, degree, tiny,
2223  + tol0, tol1, tol2, tolb, xthrsh
2224  integer digits, maxit1, maxit2
2225  logical init
2226  common /geocom/ dblmin, dbleps, pi, degree, tiny,
2227  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2228 
2229  double precision xx, yy
2230  integer q
2231  if (abs(y) .gt. abs(x)) then
2232  xx = y
2233  yy = x
2234  q = 2
2235  else
2236  xx = x
2237  yy = y
2238  q = 0
2239  end if
2240  if (xx .lt. 0) then
2241  xx = -xx
2242  q = q + 1
2243  end if
2244  atn2dx = atan2(yy, xx) / degree
2245  if (q .eq. 1) then
2246  if (yy .ge. 0) then
2247  atn2dx = 180 - atn2dx
2248  else
2249  atn2dx = -180 - atn2dx
2250  end if
2251  else if (q .eq. 2) then
2252  atn2dx = 90 - atn2dx
2253  else if (q .eq. 3) then
2254  atn2dx = -90 + atn2dx
2255  end if
2256 
2257  return
2258  end
2259 
2260  double precision function polval(N, p, x)
2261 * input
2262  integer N
2263  double precision p(0:N), x
2264 
2265  integer i
2266  if (N .lt. 0) then
2267  polval = 0
2268  else
2269  polval = p(0)
2270  end if
2271  do 10 i = 1, n
2272  polval = polval * x + p(i)
2273  10 continue
2274 
2275  return
2276  end
2277 
2278 * Table of name abbreviations to conform to the 6-char limit and
2279 * potential name conflicts.
2280 * A3coeff A3cof
2281 * C3coeff C3cof
2282 * C4coeff C4cof
2283 * AngNormalize AngNm
2284 * AngDiff AngDif
2285 * AngRound AngRnd
2286 * arcmode arcmod
2287 * Astroid Astrd
2288 * betscale betscl
2289 * lamscale lamscl
2290 * cbet12a cbt12a
2291 * sbet12a sbt12a
2292 * epsilon dbleps
2293 * realmin dblmin
2294 * geodesic geod
2295 * inverse invers
2296 * InverseStart InvSta
2297 * Lambda12 Lam12f
2298 * latsign latsgn
2299 * lonsign lonsgn
2300 * Lengths Lengs
2301 * meridian merid
2302 * outmask omask
2303 * shortline shortp
2304 * norm norm2x
2305 * SinCosSeries TrgSum
2306 * xthresh xthrsh
2307 * transit trnsit
2308 * polyval polval
2309 * LONG_UNROLL unroll
2310 * sincosd sncsdx
2311 * atan2d atn2dx
2312 *> @endcond SKIP