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