Fortran library for Geodesics  1.49
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-2017) <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.49.
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 *! Only simple polygons (which are not self-intersecting) are allowed.
965 *! There's no need to "close" the polygon by repeating the first vertex.
966 *! The area returned is signed with counter-clockwise traversal being
967 *! treated as positive.
968 
969  subroutine area(a, f, lats, lons, n, AA, PP)
970 * input
971  integer n
972  double precision a, f, lats(n), lons(n)
973 * output
974  double precision AA, PP
975 
976  integer i, omask, cross, trnsit
977  double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
978  + atanhx, aacc(2), pacc(2)
979 
980  double precision dblmin, dbleps, pi, degree, tiny,
981  + tol0, tol1, tol2, tolb, xthrsh
982  integer digits, maxit1, maxit2
983  logical init
984  common /geocom/ dblmin, dbleps, pi, degree, tiny,
985  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
986 
987  omask = 8
988  call accini(aacc)
989  call accini(pacc)
990  cross = 0
991  do 10 i = 0, n-1
992  call invers(a, f, lats(i+1), lons(i+1),
993  + lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
994  + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, ss12)
995  call accadd(pacc, s12)
996  call accadd(aacc, -ss12)
997  cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
998  10 continue
999  pp = pacc(1)
1000  b = a * (1 - f)
1001  e2 = f * (2 - f)
1002  if (e2 .eq. 0) then
1003  c2 = a**2
1004  else if (e2 .gt. 0) then
1005  c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
1006  else
1007  c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
1008  end if
1009  area0 = 4 * pi * c2
1010  if (mod(abs(cross), 2) .eq. 1) then
1011  if (aacc(1) .lt. 0) then
1012  call accadd(aacc, +area0/2)
1013  else
1014  call accadd(aacc, -area0/2)
1015  end if
1016  end if
1017  if (aacc(1) .gt. area0/2) then
1018  call accadd(aacc, -area0)
1019  else if (aacc(1) .le. -area0/2) then
1020  call accadd(aacc, +area0)
1021  end if
1022  aa = aacc(1)
1023 
1024  return
1025  end
1026 
1027 *> Return the version numbers for this package.
1028 *!
1029 *! @param[out] major the major version number.
1030 *! @param[out] minor the minor version number.
1031 *! @param[out] patch the patch number.
1032 *!
1033 *! This subroutine was added with version 1.44.
1034 
1035  subroutine geover(major, minor, patch)
1036 * output
1037  integer major, minor, patch
1038 
1039  major = 1
1040  minor = 49
1041  patch = 0
1042 
1043  return
1044  end
1045 
1046 *> @cond SKIP
1047 
1048  block data geodat
1049  double precision dblmin, dbleps, pi, degree, tiny,
1050  + tol0, tol1, tol2, tolb, xthrsh
1051  integer digits, maxit1, maxit2
1052  logical init
1053  data init /.false./
1054  common /geocom/ dblmin, dbleps, pi, degree, tiny,
1055  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1056  end
1057 
1058  subroutine geoini
1059  double precision dblmin, dbleps, pi, degree, tiny,
1060  + tol0, tol1, tol2, tolb, xthrsh
1061  integer digits, maxit1, maxit2
1062  logical init
1063  common /geocom/ dblmin, dbleps, pi, degree, tiny,
1064  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1065 
1066  digits = 53
1067  dblmin = 0.5d0**1022
1068  dbleps = 0.5d0**(digits-1)
1069 
1070  pi = atan2(0d0, -1d0)
1071  degree = pi/180
1072 * This is about cbrt(dblmin). With other implementations, sqrt(dblmin)
1073 * is used. The larger value is used here to avoid complaints about a
1074 * IEEE_UNDERFLOW_FLAG IEEE_DENORMAL signal. This is triggered when
1075 * invers is called with points at opposite poles.
1076  tiny = 0.5d0**((1022+1)/3)
1077  tol0 = dbleps
1078 * Increase multiplier in defn of tol1 from 100 to 200 to fix inverse
1079 * case 52.784459512564 0 -52.784459512563990912 179.634407464943777557
1080 * which otherwise failed for Visual Studio 10 (Release and Debug)
1081  tol1 = 200 * tol0
1082  tol2 = sqrt(tol0)
1083 * Check on bisection interval
1084  tolb = tol0 * tol2
1085  xthrsh = 1000 * tol2
1086  maxit1 = 20
1087  maxit2 = maxit1 + digits + 10
1088 
1089  init = .true.
1090 
1091  return
1092  end
1093 
1094  subroutine lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1095  + cbet1, cbet2, omask,
1096  + s12b, m12b, m0, MM12, MM21, ep2, Ca)
1097 * input
1098  double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1099  + cbet1, cbet2, ep2
1100  integer omask
1101 * optional output
1102  double precision s12b, m12b, m0, MM12, MM21
1103 * temporary storage
1104  double precision Ca(*)
1105 
1106  integer ord, nC1, nC2
1107  parameter (ord = 6, nc1 = ord, nc2 = ord)
1108 
1109  double precision A1m1f, A2m1f, TrgSum
1110  double precision m0x, J12, A1, A2, B1, B2, csig12, t, Cb(nC2)
1111  logical distp, redlp, scalp
1112  integer l
1113 
1114 * Return m12b = (reduced length)/b; also calculate s12b = distance/b,
1115 * and m0 = coefficient of secular term in expression for reduced length.
1116 
1117  distp = (mod(omask/16, 2) .eq. 1)
1118  redlp = (mod(omask/2, 2) .eq. 1)
1119  scalp = (mod(omask/4, 2) .eq. 1)
1120 
1121 * Suppress compiler warnings
1122  m0x = 0
1123  j12 = 0
1124  a1 = 0
1125  a2 = 0
1126  if (distp .or. redlp .or. scalp) then
1127  a1 = a1m1f(eps)
1128  call c1f(eps, ca)
1129  if (redlp .or. scalp) then
1130  a2 = a2m1f(eps)
1131  call c2f(eps, cb)
1132  m0x = a1 - a2
1133  a2 = 1 + a2
1134  end if
1135  a1 = 1 + a1
1136  end if
1137  if (distp) then
1138  b1 = trgsum(.true., ssig2, csig2, ca, nc1) -
1139  + trgsum(.true., ssig1, csig1, ca, nc1)
1140 * Missing a factor of b
1141  s12b = a1 * (sig12 + b1)
1142  if (redlp .or. scalp) then
1143  b2 = trgsum(.true., ssig2, csig2, cb, nc2) -
1144  + trgsum(.true., ssig1, csig1, cb, nc2)
1145  j12 = m0x * sig12 + (a1 * b1 - a2 * b2)
1146  end if
1147  else if (redlp .or. scalp) then
1148 * Assume here that nC1 >= nC2
1149  do 10 l = 1, nc2
1150  cb(l) = a1 * ca(l) - a2 * cb(l)
1151  10 continue
1152  j12 = m0x * sig12 + (trgsum(.true., ssig2, csig2, cb, nc2) -
1153  + trgsum(.true., ssig1, csig1, cb, nc2))
1154  end if
1155  if (redlp) then
1156  m0 = m0x
1157 * Missing a factor of b.
1158 * Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
1159 * accurate cancellation in the case of coincident points.
1160  m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
1161  + csig1 * csig2 * j12
1162  end if
1163  if (scalp) then
1164  csig12 = csig1 * csig2 + ssig1 * ssig2
1165  t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
1166  mm12 = csig12 + (t * ssig2 - csig2 * j12) * ssig1 / dn1
1167  mm21 = csig12 - (t * ssig1 - csig1 * j12) * ssig2 / dn2
1168  end if
1169 
1170  return
1171  end
1172 
1173  double precision function astrd(x, y)
1174 * Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
1175 * This solution is adapted from Geocentric::Reverse.
1176 * input
1177  double precision x, y
1178 
1179  double precision cbrt
1180  double precision k, p, q, r, S, r2, r3, disc, u,
1181  + t3, t, ang, v, uv, w
1182 
1183  p = x**2
1184  q = y**2
1185  r = (p + q - 1) / 6
1186  if ( .not. (q .eq. 0 .and. r .lt. 0) ) then
1187 * Avoid possible division by zero when r = 0 by multiplying equations
1188 * for s and t by r^3 and r, resp.
1189 * S = r^3 * s
1190  s = p * q / 4
1191  r2 = r**2
1192  r3 = r * r2
1193 * The discriminant of the quadratic equation for T3. This is zero on
1194 * the evolute curve p^(1/3)+q^(1/3) = 1
1195  disc = s * (s + 2 * r3)
1196  u = r
1197  if (disc .ge. 0) then
1198  t3 = s + r3
1199 * Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
1200 * of precision due to cancellation. The result is unchanged because
1201 * of the way the T is used in definition of u.
1202 * T3 = (r * t)^3
1203  if (t3 .lt. 0) then
1204  disc = -sqrt(disc)
1205  else
1206  disc = sqrt(disc)
1207  end if
1208  t3 = t3 + disc
1209 * N.B. cbrt always returns the real root. cbrt(-8) = -2.
1210 * T = r * t
1211  t = cbrt(t3)
1212 * T can be zero; but then r2 / T -> 0.
1213  if (t .ne. 0) u = u + t + r2 / t
1214  else
1215 * T is complex, but the way u is defined the result is real.
1216  ang = atan2(sqrt(-disc), -(s + r3))
1217 * There are three possible cube roots. We choose the root which
1218 * avoids cancellation. Note that disc < 0 implies that r < 0.
1219  u = u + 2 * r * cos(ang / 3)
1220  end if
1221 * guaranteed positive
1222  v = sqrt(u**2 + q)
1223 * Avoid loss of accuracy when u < 0.
1224 * u+v, guaranteed positive
1225  if (u .lt. 0) then
1226  uv = q / (v - u)
1227  else
1228  uv = u + v
1229  end if
1230 * positive?
1231  w = (uv - q) / (2 * v)
1232 * Rearrange expression for k to avoid loss of accuracy due to
1233 * subtraction. Division by 0 not possible because uv > 0, w >= 0.
1234 * guaranteed positive
1235  k = uv / (sqrt(uv + w**2) + w)
1236  else
1237 * q == 0 && r <= 0
1238 * y = 0 with |x| <= 1. Handle this case directly.
1239 * for y small, positive root is k = abs(y)/sqrt(1-x^2)
1240  k = 0
1241  end if
1242  astrd = k
1243 
1244  return
1245  end
1246 
1247  double precision function invsta(sbet1, cbet1, dn1,
1248  + sbet2, cbet2, dn2, lam12, slam12, clam12, f, A3x,
1249  + salp1, calp1, salp2, calp2, dnm,
1250  + Ca)
1251 * Return a starting point for Newton's method in salp1 and calp1
1252 * (function value is -1). If Newton's method doesn't need to be used,
1253 * return also salp2, calp2, and dnm and function value is sig12.
1254 * input
1255  double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1256  + lam12, slam12, clam12, f, A3x(*)
1257 * output
1258  double precision salp1, calp1, salp2, calp2, dnm
1259 * temporary
1260  double precision Ca(*)
1261 
1262  double precision hypotx, A3f, Astrd
1263  logical shortp
1264  double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
1265  + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
1266  + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
1267  + k, omg12a, sbetm2, lam12x
1268 
1269  double precision dblmin, dbleps, pi, degree, tiny,
1270  + tol0, tol1, tol2, tolb, xthrsh
1271  integer digits, maxit1, maxit2
1272  logical init
1273  common /geocom/ dblmin, dbleps, pi, degree, tiny,
1274  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1275 
1276  f1 = 1 - f
1277  e2 = f * (2 - f)
1278  ep2 = e2 / (1 - e2)
1279  n = f / (2 - f)
1280 * The sig12 threshold for "really short". Using the auxiliary sphere
1281 * solution with dnm computed at (bet1 + bet2) / 2, the relative error in
1282 * the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
1283 * (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a
1284 * given f and sig12, the max error occurs for lines near the pole. If
1285 * the old rule for computing dnm = (dn1 + dn2)/2 is used, then the error
1286 * increases by a factor of 2.) Setting this equal to epsilon gives
1287 * sig12 = etol2. Here 0.1 is a safety factor (error decreased by 100)
1288 * and max(0.001, abs(f)) stops etol2 getting too large in the nearly
1289 * spherical case.
1290  etol2 = 0.1d0 * tol2 /
1291  + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
1292 
1293 * Return value
1294  sig12 = -1
1295 * bet12 = bet2 - bet1 in [0, pi); bt12a = bet2 + bet1 in (-pi, 0]
1296  sbet12 = sbet2 * cbet1 - cbet2 * sbet1
1297  cbet12 = cbet2 * cbet1 + sbet2 * sbet1
1298  sbt12a = sbet2 * cbet1 + cbet2 * sbet1
1299 
1300  shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
1301  + cbet2 * lam12 .lt. 0.5d0
1302 
1303  if (shortp) then
1304  sbetm2 = (sbet1 + sbet2)**2
1305 * sin((bet1+bet2)/2)^2
1306 * = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
1307  sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
1308  dnm = sqrt(1 + ep2 * sbetm2)
1309  omg12 = lam12 / (f1 * dnm)
1310  somg12 = sin(omg12)
1311  comg12 = cos(omg12)
1312  else
1313  somg12 = slam12
1314  comg12 = clam12
1315  end if
1316 
1317  salp1 = cbet2 * somg12
1318  if (comg12 .ge. 0) then
1319  calp1 = sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12)
1320  else
1321  calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1322  end if
1323 
1324  ssig12 = hypotx(salp1, calp1)
1325  csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
1326 
1327  if (shortp .and. ssig12 .lt. etol2) then
1328 * really short lines
1329  salp2 = cbet1 * somg12
1330  if (comg12 .ge. 0) then
1331  calp2 = somg12**2 / (1 + comg12)
1332  else
1333  calp2 = 1 - comg12
1334  end if
1335  calp2 = sbet12 - cbet1 * sbet2 * calp2
1336  call norm2x(salp2, calp2)
1337 * Set return value
1338  sig12 = atan2(ssig12, csig12)
1339  else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
1340  + ssig12 .ge. 6 * abs(n) * pi * cbet1**2) then
1341 * Nothing to do, zeroth order spherical approximation is OK
1342  continue
1343  else
1344 * lam12 - pi
1345  lam12x = atan2(-slam12, -clam12)
1346 * Scale lam12 and bet2 to x, y coordinate system where antipodal point
1347 * is at origin and singular point is at y = 0, x = -1.
1348  if (f .ge. 0) then
1349 * x = dlong, y = dlat
1350  k2 = sbet1**2 * ep2
1351  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1352  lamscl = f * cbet1 * a3f(eps, a3x) * pi
1353  betscl = lamscl * cbet1
1354  x = lam12x / lamscl
1355  y = sbt12a / betscl
1356  else
1357 * f < 0: x = dlat, y = dlong
1358  cbt12a = cbet2 * cbet1 - sbet2 * sbet1
1359  bt12a = atan2(sbt12a, cbt12a)
1360 * In the case of lon12 = 180, this repeats a calculation made in
1361 * Inverse.
1362  call lengs(n, pi + bt12a,
1363  + sbet1, -cbet1, dn1, sbet2, cbet2, dn2, cbet1, cbet2, 2,
1364  + dummy, m12b, m0, dummy, dummy, ep2, ca)
1365  x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
1366  if (x .lt. -0.01d0) then
1367  betscl = sbt12a / x
1368  else
1369  betscl = -f * cbet1**2 * pi
1370  end if
1371  lamscl = betscl / cbet1
1372  y = lam12x / lamscl
1373  end if
1374 
1375  if (y .gt. -tol1 .and. x .gt. -1 - xthrsh) then
1376 * strip near cut
1377  if (f .ge. 0) then
1378  salp1 = min(1d0, -x)
1379  calp1 = - sqrt(1 - salp1**2)
1380  else
1381  if (x .gt. -tol1) then
1382  calp1 = 0
1383  else
1384  calp1 = 1
1385  end if
1386  calp1 = max(calp1, x)
1387  salp1 = sqrt(1 - calp1**2)
1388  end if
1389  else
1390 * Estimate alp1, by solving the astroid problem.
1391 *
1392 * Could estimate alpha1 = theta + pi/2, directly, i.e.,
1393 * calp1 = y/k; salp1 = -x/(1+k); for f >= 0
1394 * calp1 = x/(1+k); salp1 = -y/k; for f < 0 (need to check)
1395 *
1396 * However, it's better to estimate omg12 from astroid and use
1397 * spherical formula to compute alp1. This reduces the mean number of
1398 * Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
1399 * (min 0 max 5). The changes in the number of iterations are as
1400 * follows:
1401 *
1402 * change percent
1403 * 1 5
1404 * 0 78
1405 * -1 16
1406 * -2 0.6
1407 * -3 0.04
1408 * -4 0.002
1409 *
1410 * The histogram of iterations is (m = number of iterations estimating
1411 * alp1 directly, n = number of iterations estimating via omg12, total
1412 * number of trials = 148605):
1413 *
1414 * iter m n
1415 * 0 148 186
1416 * 1 13046 13845
1417 * 2 93315 102225
1418 * 3 36189 32341
1419 * 4 5396 7
1420 * 5 455 1
1421 * 6 56 0
1422 *
1423 * Because omg12 is near pi, estimate work with omg12a = pi - omg12
1424  k = astrd(x, y)
1425  if (f .ge. 0) then
1426  omg12a = -x * k/(1 + k)
1427  else
1428  omg12a = -y * (1 + k)/k
1429  end if
1430  omg12a = lamscl * omg12a
1431  somg12 = sin(omg12a)
1432  comg12 = -cos(omg12a)
1433 * Update spherical estimate of alp1 using omg12 instead of lam12
1434  salp1 = cbet2 * somg12
1435  calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
1436  end if
1437  end if
1438 * Sanity check on starting guess. Backwards check allows NaN through.
1439  if (.not. (salp1 .le. 0)) then
1440  call norm2x(salp1, calp1)
1441  else
1442  salp1 = 1
1443  calp1 = 0
1444  end if
1445  invsta = sig12
1446 
1447  return
1448  end
1449 
1450  double precision function lam12f(sbet1, cbet1, dn1,
1451  + sbet2, cbet2, dn2, salp1, calp1, slm120, clm120, f, A3x, C3x,
1452  + salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps,
1453  + domg12, diffp, dlam12, Ca)
1454 * input
1455  double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
1456  + salp1, calp1, slm120, clm120, f, A3x(*), C3x(*)
1457  logical diffp
1458 * output
1459  double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
1460  + eps, domg12
1461 * optional output
1462  double precision dlam12
1463 * temporary
1464  double precision Ca(*)
1465 
1466  integer ord, nC3
1467  parameter(ord = 6, nc3 = ord)
1468 
1469  double precision hypotx, A3f, TrgSum
1470 
1471  double precision f1, e2, ep2, salp0, calp0,
1472  + somg1, comg1, somg2, comg2, somg12, comg12,
1473  + lam12, eta, b312, k2, dummy
1474 
1475  double precision dblmin, dbleps, pi, degree, tiny,
1476  + tol0, tol1, tol2, tolb, xthrsh
1477  integer digits, maxit1, maxit2
1478  logical init
1479  common /geocom/ dblmin, dbleps, pi, degree, tiny,
1480  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
1481 
1482  f1 = 1 - f
1483  e2 = f * (2 - f)
1484  ep2 = e2 / (1 - e2)
1485 * Break degeneracy of equatorial line. This case has already been
1486 * handled.
1487  if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
1488 
1489 * sin(alp1) * cos(bet1) = sin(alp0)
1490  salp0 = salp1 * cbet1
1491 * calp0 > 0
1492  calp0 = hypotx(calp1, salp1 * sbet1)
1493 
1494 * tan(bet1) = tan(sig1) * cos(alp1)
1495 * tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
1496  ssig1 = sbet1
1497  somg1 = salp0 * sbet1
1498  csig1 = calp1 * cbet1
1499  comg1 = csig1
1500  call norm2x(ssig1, csig1)
1501 * norm2x(somg1, comg1); -- don't need to normalize!
1502 
1503 * Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful
1504 * about this case, since this can yield singularities in the Newton
1505 * iteration.
1506 * sin(alp2) * cos(bet2) = sin(alp0)
1507  if (cbet2 .ne. cbet1) then
1508  salp2 = salp0 / cbet2
1509  else
1510  salp2 = salp1
1511  end if
1512 * calp2 = sqrt(1 - sq(salp2))
1513 * = sqrt(sq(calp0) - sq(sbet2)) / cbet2
1514 * and subst for calp0 and rearrange to give (choose positive sqrt
1515 * to give alp2 in [0, pi/2]).
1516  if (cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1) then
1517  if (cbet1 .lt. -sbet1) then
1518  calp2 = (cbet2 - cbet1) * (cbet1 + cbet2)
1519  else
1520  calp2 = (sbet1 - sbet2) * (sbet1 + sbet2)
1521  end if
1522  calp2 = sqrt((calp1 * cbet1)**2 + calp2) / cbet2
1523  else
1524  calp2 = abs(calp1)
1525  end if
1526 * tan(bet2) = tan(sig2) * cos(alp2)
1527 * tan(omg2) = sin(alp0) * tan(sig2).
1528  ssig2 = sbet2
1529  somg2 = salp0 * sbet2
1530  csig2 = calp2 * cbet2
1531  comg2 = csig2
1532  call norm2x(ssig2, csig2)
1533 * norm2x(somg2, comg2); -- don't need to normalize!
1534 
1535 * sig12 = sig2 - sig1, limit to [0, pi]
1536  sig12 = atan2(0d0 + max(0d0, csig1 * ssig2 - ssig1 * csig2),
1537  + csig1 * csig2 + ssig1 * ssig2)
1538 
1539 * omg12 = omg2 - omg1, limit to [0, pi]
1540  somg12 = 0d0 + max(0d0, comg1 * somg2 - somg1 * comg2)
1541  comg12 = comg1 * comg2 + somg1 * somg2
1542 * eta = omg12 - lam120
1543  eta = atan2(somg12 * clm120 - comg12 * slm120,
1544  + comg12 * clm120 + somg12 * slm120);
1545  k2 = calp0**2 * ep2
1546  eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
1547  call c3f(eps, c3x, ca)
1548  b312 = (trgsum(.true., ssig2, csig2, ca, nc3-1) -
1549  + trgsum(.true., ssig1, csig1, ca, nc3-1))
1550  domg12 = -f * a3f(eps, a3x) * salp0 * (sig12 + b312)
1551  lam12 = eta + domg12
1552 
1553  if (diffp) then
1554  if (calp2 .eq. 0) then
1555  dlam12 = - 2 * f1 * dn1 / sbet1
1556  else
1557  call lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
1558  + cbet1, cbet2, 2,
1559  + dummy, dlam12, dummy, dummy, dummy, ep2, ca)
1560  dlam12 = dlam12 * f1 / (calp2 * cbet2)
1561  end if
1562  end if
1563  lam12f = lam12
1564 
1565  return
1566  end
1567 
1568  double precision function a3f(eps, A3x)
1569 * Evaluate A3
1570  integer ord, nA3, nA3x
1571  parameter(ord = 6, na3 = ord, na3x = na3)
1572 
1573 * input
1574  double precision eps
1575 * output
1576  double precision A3x(0: nA3x-1)
1577 
1578  double precision polval
1579  A3f = polval(na3 - 1, a3x, eps)
1580 
1581  return
1582  end
1583 
1584  subroutine c3f(eps, C3x, c)
1585 * Evaluate C3 coeffs
1586 * Elements c[1] thru c[nC3-1] are set
1587  integer ord, nC3, nC3x
1588  parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1589 
1590 * input
1591  double precision eps, C3x(0:nC3x-1)
1592 * output
1593  double precision c(nC3-1)
1594 
1595  integer o, m, l
1596  double precision mult, polval
1597 
1598  mult = 1
1599  o = 0
1600  do 10 l = 1, nc3 - 1
1601  m = nc3 - l - 1
1602  mult = mult * eps
1603  c(l) = mult * polval(m, c3x(o), eps)
1604  o = o + m + 1
1605  10 continue
1606 
1607  return
1608  end
1609 
1610  subroutine c4f(eps, C4x, c)
1611 * Evaluate C4
1612 * Elements c[0] thru c[nC4-1] are set
1613  integer ord, nC4, nC4x
1614  parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1615 
1616 * input
1617  double precision eps, C4x(0:nC4x-1)
1618 *output
1619  double precision c(0:nC4-1)
1620 
1621  integer o, m, l
1622  double precision mult, polval
1623 
1624  mult = 1
1625  o = 0
1626  do 10 l = 0, nc4 - 1
1627  m = nc4 - l - 1
1628  c(l) = mult * polval(m, c4x(o), eps)
1629  o = o + m + 1
1630  mult = mult * eps
1631  10 continue
1632 
1633  return
1634  end
1635 
1636  double precision function a1m1f(eps)
1637 * The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
1638 * input
1639  double precision eps
1640 
1641  double precision t
1642  integer ord, nA1, o, m
1643  parameter(ord = 6, na1 = ord)
1644  double precision polval, coeff(nA1/2 + 2)
1645  data coeff /1, 4, 64, 0, 256/
1646 
1647  o = 1
1648  m = na1/2
1649  t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1650  a1m1f = (t + eps) / (1 - eps)
1651 
1652  return
1653  end
1654 
1655  subroutine c1f(eps, c)
1656 * The coefficients C1[l] in the Fourier expansion of B1
1657  integer ord, nC1
1658  parameter(ord = 6, nc1 = ord)
1659 
1660 * input
1661  double precision eps
1662 * output
1663  double precision c(nC1)
1664 
1665  double precision eps2, d
1666  integer o, m, l
1667  double precision polval, coeff((nC1**2 + 7*nC1 - 2*(nC1/2))/4)
1668  data coeff /
1669  + -1, 6, -16, 32,
1670  + -9, 64, -128, 2048,
1671  + 9, -16, 768,
1672  + 3, -5, 512,
1673  + -7, 1280,
1674  + -7, 2048/
1675 
1676  eps2 = eps**2
1677  d = eps
1678  o = 1
1679  do 10 l = 1, nc1
1680  m = (nc1 - l) / 2
1681  c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1682  o = o + m + 2
1683  d = d * eps
1684  10 continue
1685 
1686  return
1687  end
1688 
1689  subroutine c1pf(eps, c)
1690 * The coefficients C1p[l] in the Fourier expansion of B1p
1691  integer ord, nC1p
1692  parameter(ord = 6, nc1p = ord)
1693 
1694 * input
1695  double precision eps
1696 * output
1697  double precision c(nC1p)
1698 
1699  double precision eps2, d
1700  integer o, m, l
1701  double precision polval, coeff((nC1p**2 + 7*nC1p - 2*(nC1p/2))/4)
1702  data coeff /
1703  + 205, -432, 768, 1536,
1704  + 4005, -4736, 3840, 12288,
1705  + -225, 116, 384,
1706  + -7173, 2695, 7680,
1707  + 3467, 7680,
1708  + 38081, 61440/
1709 
1710  eps2 = eps**2
1711  d = eps
1712  o = 1
1713  do 10 l = 1, nc1p
1714  m = (nc1p - l) / 2
1715  c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1716  o = o + m + 2
1717  d = d * eps
1718  10 continue
1719 
1720  return
1721  end
1722 
1723 * The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
1724  double precision function a2m1f(eps)
1725 * input
1726  double precision eps
1727 
1728  double precision t
1729  integer ord, nA2, o, m
1730  parameter(ord = 6, na2 = ord)
1731  double precision polval, coeff(nA2/2 + 2)
1732  data coeff /-11, -28, -192, 0, 256/
1733 
1734  o = 1
1735  m = na2/2
1736  t = polval(m, coeff(o), eps**2) / coeff(o + m + 1)
1737  a2m1f = (t - eps) / (1 + eps)
1738 
1739  return
1740  end
1741 
1742  subroutine c2f(eps, c)
1743 * The coefficients C2[l] in the Fourier expansion of B2
1744  integer ord, nC2
1745  parameter(ord = 6, nc2 = ord)
1746 
1747 * input
1748  double precision eps
1749 * output
1750  double precision c(nC2)
1751 
1752  double precision eps2, d
1753  integer o, m, l
1754  double precision polval, coeff((nC2**2 + 7*nC2 - 2*(nC2/2))/4)
1755  data coeff /
1756  + 1, 2, 16, 32,
1757  + 35, 64, 384, 2048,
1758  + 15, 80, 768,
1759  + 7, 35, 512,
1760  + 63, 1280,
1761  + 77, 2048/
1762 
1763  eps2 = eps**2
1764  d = eps
1765  o = 1
1766  do 10 l = 1, nc2
1767  m = (nc2 - l) / 2
1768  c(l) = d * polval(m, coeff(o), eps2) / coeff(o + m + 1)
1769  o = o + m + 2
1770  d = d * eps
1771  10 continue
1772 
1773  return
1774  end
1775 
1776  subroutine a3cof(n, A3x)
1777 * The scale factor A3 = mean value of (d/dsigma)I3
1778  integer ord, nA3, nA3x
1779  parameter(ord = 6, na3 = ord, na3x = na3)
1780 
1781 * input
1782  double precision n
1783 * output
1784  double precision A3x(0:nA3x-1)
1785 
1786  integer o, m, k, j
1787  double precision polval, coeff((nA3**2 + 7*nA3 - 2*(nA3/2))/4)
1788  data coeff /
1789  + -3, 128,
1790  + -2, -3, 64,
1791  + -1, -3, -1, 16,
1792  + 3, -1, -2, 8,
1793  + 1, -1, 2,
1794  + 1, 1/
1795 
1796  o = 1
1797  k = 0
1798  do 10 j = na3 - 1, 0, -1
1799  m = min(na3 - j - 1, j)
1800  a3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1801  k = k + 1
1802  o = o + m + 2
1803  10 continue
1804 
1805  return
1806  end
1807 
1808  subroutine c3cof(n, C3x)
1809 * The coefficients C3[l] in the Fourier expansion of B3
1810  integer ord, nC3, nC3x
1811  parameter(ord = 6, nc3 = ord, nc3x = (nc3 * (nc3 - 1)) / 2)
1812 
1813 * input
1814  double precision n
1815 * output
1816  double precision C3x(0:nC3x-1)
1817 
1818  integer o, m, l, j, k
1819  double precision polval,
1820  + coeff(((nc3-1)*(nc3**2 + 7*nc3 - 2*(nc3/2)))/8)
1821  data coeff /
1822  + 3, 128,
1823  + 2, 5, 128,
1824  + -1, 3, 3, 64,
1825  + -1, 0, 1, 8,
1826  + -1, 1, 4,
1827  + 5, 256,
1828  + 1, 3, 128,
1829  + -3, -2, 3, 64,
1830  + 1, -3, 2, 32,
1831  + 7, 512,
1832  + -10, 9, 384,
1833  + 5, -9, 5, 192,
1834  + 7, 512,
1835  + -14, 7, 512,
1836  + 21, 2560/
1837 
1838  o = 1
1839  k = 0
1840  do 20 l = 1, nc3 - 1
1841  do 10 j = nc3 - 1, l, -1
1842  m = min(nc3 - j - 1, j)
1843  c3x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1844  k = k + 1
1845  o = o + m + 2
1846  10 continue
1847  20 continue
1848 
1849  return
1850  end
1851 
1852  subroutine c4cof(n, C4x)
1853 * The coefficients C4[l] in the Fourier expansion of I4
1854  integer ord, nC4, nC4x
1855  parameter(ord = 6, nc4 = ord, nc4x = (nc4 * (nc4 + 1)) / 2)
1856 
1857 * input
1858  double precision n
1859 * output
1860  double precision C4x(0:nC4x-1)
1861 
1862  integer o, m, l, j, k
1863  double precision polval, coeff((nC4 * (nC4 + 1) * (nC4 + 5)) / 6)
1864  data coeff /
1865  + 97, 15015, 1088, 156, 45045, -224, -4784, 1573, 45045,
1866  + -10656, 14144, -4576, -858, 45045,
1867  + 64, 624, -4576, 6864, -3003, 15015,
1868  + 100, 208, 572, 3432, -12012, 30030, 45045,
1869  + 1, 9009, -2944, 468, 135135, 5792, 1040, -1287, 135135,
1870  + 5952, -11648, 9152, -2574, 135135,
1871  + -64, -624, 4576, -6864, 3003, 135135,
1872  + 8, 10725, 1856, -936, 225225, -8448, 4992, -1144, 225225,
1873  + -1440, 4160, -4576, 1716, 225225,
1874  + -136, 63063, 1024, -208, 105105,
1875  + 3584, -3328, 1144, 315315,
1876  + -128, 135135, -2560, 832, 405405, 128, 99099/
1877 
1878  o = 1
1879  k = 0
1880  do 20 l = 0, nc4 - 1
1881  do 10 j = nc4 - 1, l, -1
1882  m = nc4 - j - 1
1883  c4x(k) = polval(m, coeff(o), n) / coeff(o + m + 1)
1884  k = k + 1
1885  o = o + m + 2
1886  10 continue
1887  20 continue
1888 
1889  return
1890  end
1891 
1892  double precision function sumx(u, v, t)
1893 * input
1894  double precision u, v
1895 * output
1896  double precision t
1897 
1898  double precision up, vpp
1899  sumx = u + v
1900  up = sumx - v
1901  vpp = sumx - up
1902  up = up - u
1903  vpp = vpp - v
1904  t = -(up + vpp)
1905 
1906  return
1907  end
1908 
1909  double precision function angnm(x)
1910 * input
1911  double precision x
1912 
1913  angnm = mod(x, 360d0)
1914  if (angnm .le. -180) then
1915  angnm = angnm + 360
1916  else if (angnm .gt. 180) then
1917  angnm = angnm - 360
1918  end if
1919 
1920  return
1921  end
1922 
1923  double precision function latfix(x)
1924 * input
1925  double precision x
1926 
1927  latfix = x
1928  if (.not. (abs(x) .gt. 90)) return
1929 * concoct a NaN
1930  latfix = sqrt(90 - abs(x))
1931 
1932  return
1933  end
1934 
1935  double precision function angdif(x, y, e)
1936 * Compute y - x. x and y must both lie in [-180, 180]. The result is
1937 * equivalent to computing the difference exactly, reducing it to (-180,
1938 * 180] and rounding the result. Note that this prescription allows -180
1939 * to be returned (e.g., if x is tiny and negative and y = 180). The
1940 * error in the difference is returned in e
1941 * input
1942  double precision x, y
1943 * output
1944  double precision e
1945 
1946  double precision d, t, sumx, AngNm
1947  d = angnm(sumx(angnm(-x), angnm(y), t))
1948  if (d .eq. 180 .and. t .gt. 0) then
1949  d = -180
1950  end if
1951  angdif = sumx(d, t, e)
1952 
1953  return
1954  end
1955 
1956  double precision function angrnd(x)
1957 * The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
1958 * for reals = 0.7 pm on the earth if x is an angle in degrees. (This is
1959 * about 1000 times more resolution than we get with angles around 90
1960 * degrees.) We use this to avoid having to deal with near singular
1961 * cases when x is non-zero but tiny (e.g., 1.0e-200). This also
1962 * converts -0 to +0.
1963 * input
1964  double precision x
1965 
1966  double precision y, z
1967  z = 1/16d0
1968  y = abs(x)
1969 * The compiler mustn't "simplify" z - (z - y) to y
1970  if (y .lt. z) y = z - (z - y)
1971  angrnd = sign(y, x)
1972  if (x .eq. 0) angrnd = 0
1973 
1974  return
1975  end
1976 
1977  subroutine swap(x, y)
1978 * input/output
1979  double precision x, y
1980 
1981  double precision z
1982  z = x
1983  x = y
1984  y = z
1985 
1986  return
1987  end
1988 
1989  double precision function hypotx(x, y)
1990 * input
1991  double precision x, y
1992 
1993 * With Fortran 2008, this becomes: hypotx = hypot(x, y)
1994  hypotx = sqrt(x**2 + y**2)
1995 
1996  return
1997  end
1998 
1999  subroutine norm2x(x, y)
2000 * input/output
2001  double precision x, y
2002 
2003  double precision hypotx, r
2004  r = hypotx(x, y)
2005  x = x/r
2006  y = y/r
2007 
2008  return
2009  end
2010 
2011  double precision function log1px(x)
2012 * input
2013  double precision x
2014 
2015  double precision y, z
2016  y = 1 + x
2017  z = y - 1
2018  if (z .eq. 0) then
2019  log1px = x
2020  else
2021  log1px = x * log(y) / z
2022  end if
2023 
2024  return
2025  end
2026 
2027  double precision function atanhx(x)
2028 * input
2029  double precision x
2030 
2031 * With Fortran 2008, this becomes: atanhx = atanh(x)
2032  double precision log1px, y
2033  y = abs(x)
2034  y = log1px(2 * y/(1 - y))/2
2035  atanhx = sign(y, x)
2036 
2037  return
2038  end
2039 
2040  double precision function cbrt(x)
2041 * input
2042  double precision x
2043 
2044  cbrt = sign(abs(x)**(1/3d0), x)
2045 
2046  return
2047  end
2048 
2049  double precision function trgsum(sinp, sinx, cosx, c, n)
2050 * Evaluate
2051 * y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
2052 * sum(c[i] * cos((2*i-1) * x), i, 1, n)
2053 * using Clenshaw summation.
2054 * Approx operation count = (n + 5) mult and (2 * n + 2) add
2055 * input
2056  logical sinp
2057  integer n
2058  double precision sinx, cosx, c(n)
2059 
2060  double precision ar, y0, y1
2061  integer n2, k
2062 
2063 * 2 * cos(2 * x)
2064  ar = 2 * (cosx - sinx) * (cosx + sinx)
2065 * accumulators for sum
2066  if (mod(n, 2) .eq. 1) then
2067  y0 = c(n)
2068  n2 = n - 1
2069  else
2070  y0 = 0
2071  n2 = n
2072  end if
2073  y1 = 0
2074 * Now n2 is even
2075  do 10 k = n2, 1, -2
2076 * Unroll loop x 2, so accumulators return to their original role
2077  y1 = ar * y0 - y1 + c(k)
2078  y0 = ar * y1 - y0 + c(k-1)
2079  10 continue
2080  if (sinp) then
2081 * sin(2 * x) * y0
2082  trgsum = 2 * sinx * cosx * y0
2083  else
2084 * cos(x) * (y0 - y1)
2085  trgsum = cosx * (y0 - y1)
2086  end if
2087 
2088  return
2089  end
2090 
2091  integer function trnsit(lon1, lon2)
2092 * input
2093  double precision lon1, lon2
2094 
2095  double precision lon1x, lon2x, lon12, AngNm, AngDif, e
2096  lon1x = angnm(lon1)
2097  lon2x = angnm(lon2)
2098  lon12 = angdif(lon1x, lon2x, e)
2099  trnsit = 0
2100  if (lon1x .le. 0 .and. lon2x .gt. 0 .and. lon12 .gt. 0) then
2101  trnsit = 1
2102  else if (lon2x .le. 0 .and. lon1x .gt. 0 .and. lon12 .lt. 0) then
2103  trnsit = -1
2104  end if
2105 
2106  return
2107  end
2108 
2109  subroutine accini(s)
2110 * Initialize an accumulator; this is an array with two elements.
2111 * input/output
2112  double precision s(2)
2113 
2114  s(1) = 0
2115  s(2) = 0
2116 
2117  return
2118  end
2119 
2120  subroutine accadd(s, y)
2121 * Add y to an accumulator.
2122 * input
2123  double precision y
2124 * input/output
2125  double precision s(2)
2126 
2127  double precision z, u, sumx
2128  z = sumx(y, s(2), u)
2129  s(1) = sumx(z, s(1), s(2))
2130  if (s(1) .eq. 0) then
2131  s(1) = u
2132  else
2133  s(2) = s(2) + u
2134  end if
2135 
2136  return
2137  end
2138 
2139  subroutine sncsdx(x, sinx, cosx)
2140 * Compute sin(x) and cos(x) with x in degrees
2141 * input
2142  double precision x
2143 * input/output
2144  double precision sinx, cosx
2145 
2146  double precision dblmin, dbleps, pi, degree, tiny,
2147  + tol0, tol1, tol2, tolb, xthrsh
2148  integer digits, maxit1, maxit2
2149  logical init
2150  common /geocom/ dblmin, dbleps, pi, degree, tiny,
2151  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2152 
2153  double precision r, s, c
2154  integer q
2155  r = mod(x, 360d0)
2156  q = nint(r / 90)
2157  r = (r - 90 * q) * degree
2158  s = sin(r)
2159 * sin(-0) isn't reliably -0, so...
2160  if (x .eq. 0) s = x
2161  c = cos(r)
2162  q = mod(q + 4, 4)
2163  if (q .eq. 0) then
2164  sinx = s
2165  cosx = c
2166  else if (q .eq. 1) then
2167  sinx = c
2168  cosx = -s
2169  else if (q .eq. 2) then
2170  sinx = -s
2171  cosx = -c
2172  else
2173 * q.eq.3
2174  sinx = -c
2175  cosx = s
2176  end if
2177 
2178  if (x .ne. 0) then
2179  sinx = 0d0 + sinx
2180  cosx = 0d0 + cosx
2181  end if
2182 
2183  return
2184  end
2185 
2186  double precision function atn2dx(y, x)
2187 * input
2188  double precision x, y
2189 
2190  double precision dblmin, dbleps, pi, degree, tiny,
2191  + tol0, tol1, tol2, tolb, xthrsh
2192  integer digits, maxit1, maxit2
2193  logical init
2194  common /geocom/ dblmin, dbleps, pi, degree, tiny,
2195  + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
2196 
2197  double precision xx, yy
2198  integer q
2199  if (abs(y) .gt. abs(x)) then
2200  xx = y
2201  yy = x
2202  q = 2
2203  else
2204  xx = x
2205  yy = y
2206  q = 0
2207  end if
2208  if (xx .lt. 0) then
2209  xx = -xx
2210  q = q + 1
2211  end if
2212  atn2dx = atan2(yy, xx) / degree
2213  if (q .eq. 1) then
2214  if (yy .ge. 0) then
2215  atn2dx = 180 - atn2dx
2216  else
2217  atn2dx = -180 - atn2dx
2218  end if
2219  else if (q .eq. 2) then
2220  atn2dx = 90 - atn2dx
2221  else if (q .eq. 3) then
2222  atn2dx = -90 + atn2dx
2223  end if
2224 
2225  return
2226  end
2227 
2228  double precision function polval(N, p, x)
2229 * input
2230  integer N
2231  double precision p(0:N), x
2232 
2233  integer i
2234  if (N .lt. 0) then
2235  polval = 0
2236  else
2237  polval = p(0)
2238  end if
2239  do 10 i = 1, n
2240  polval = polval * x + p(i)
2241  10 continue
2242 
2243  return
2244  end
2245 
2246 * Table of name abbreviations to conform to the 6-char limit and
2247 * potential name conflicts.
2248 * A3coeff A3cof
2249 * C3coeff C3cof
2250 * C4coeff C4cof
2251 * AngNormalize AngNm
2252 * AngDiff AngDif
2253 * AngRound AngRnd
2254 * arcmode arcmod
2255 * Astroid Astrd
2256 * betscale betscl
2257 * lamscale lamscl
2258 * cbet12a cbt12a
2259 * sbet12a sbt12a
2260 * epsilon dbleps
2261 * realmin dblmin
2262 * geodesic geod
2263 * inverse invers
2264 * InverseStart InvSta
2265 * Lambda12 Lam12f
2266 * latsign latsgn
2267 * lonsign lonsgn
2268 * Lengths Lengs
2269 * meridian merid
2270 * outmask omask
2271 * shortline shortp
2272 * norm norm2x
2273 * SinCosSeries TrgSum
2274 * xthresh xthrsh
2275 * transit trnsit
2276 * polyval polval
2277 * LONG_UNROLL unroll
2278 * sincosd sncsdx
2279 * atan2d atn2dx
2280 *> @endcond SKIP