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