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