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