00001 * The subroutines in this files are documented at
00002 * http:
00003 *
00004 *> @file geodesic.for
00005 *! @brief Implementation of geodesic routines in Fortran
00006 *!
00007 *! This is a Fortran implementation of the geodesic algorithms described
00008 *! in
00009 *! - C. F. F. Karney,
00010 *! <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
00011 *! Algorithms for geodesics</a>,
00012 *! J. Geodesy <b>87</b>, 43--55 (2013);
00013 *! DOI: <a href="http://dx.doi.org/10.1007/s00190-012-0578-z">
00014 *! 10.1007/s00190-012-0578-z</a>;
00015 *! addenda: <a href="http://geographiclib.sf.net/geod-addenda.html">
00016 *! geod-addenda.html</a>.
00017 *! .
00018 *! The principal advantages of these algorithms over previous ones
00019 *! (e.g., Vincenty, 1975) are
00020 *! - accurate to round off for |<i>f</i>| < 1/50;
00021 *! - the solution of the inverse problem is always found;
00022 *! - differential and integral properties of geodesics are computed.
00023 *!
00024 *! The shortest path between two points on the ellipsoid at (\e lat1, \e
00025 *! lon1) and (\e lat2, \e lon2) is called the geodesic. Its length is
00026 *! \e s12 and the geodesic from point 1 to point 2 has forward azimuths
00027 *! \e azi1 and \e azi2 at the two end points.
00028 *!
00029 *! Traditionally two geodesic problems are considered:
00030 *! - the direct problem -- given \e lat1, \e lon1, \e s12, and \e azi1,
00031 *! determine \e lat2, \e lon2, and \e azi2. This is solved by the
00032 *! subroutine direct().
00033 *! - the inverse problem -- given \e lat1, \e lon1, \e lat2, \e lon2,
00034 *! determine \e s12, \e azi1, and \e azi2. This is solved by the
00035 *! subroutine invers().
00036 *!
00037 *! The ellipsoid is specified by its equatorial radius \e a (typically
00038 *! in meters) and flattening \e f. The routines are accurate to round
00039 *! off with double precision arithmetic provided that |<i>f</i>| <
00040 *! 1/50; for the WGS84 ellipsoid, the errors are less than 15
00041 *! nanometers. (Reasonably accurate results are obtained for |<i>f</i>|
00042 *! < 1/5.) For a prolate ellipsoid, specify \e f < 0.
00043 *!
00044 *! The routines also calculate several other quantities of interest
00045 *! - \e SS12 is the area between the geodesic from point 1 to point 2
00046 *! and the equator; i.e., it is the area, measured counter-clockwise,
00047 *! of the geodesic quadrilateral with corners (\e lat1,\e lon1), (0,\e
00048 *! lon1), (0,\e lon2), and (\e lat2,\e lon2).
00049 *! - \e m12, the reduced length of the geodesic is defined such that if
00050 *! the initial azimuth is perturbed by \e dazi1 (radians) then the
00051 *! second point is displaced by \e m12 \e dazi1 in the direction
00052 *! perpendicular to the geodesic. On a curved surface the reduced
00053 *! length obeys a symmetry relation, \e m12 + \e m21 = 0. On a flat
00054 *! surface, we have \e m12 = \e s12.
00055 *! - \e MM12 and \e MM21 are geodesic scales. If two geodesics are
00056 *! parallel at point 1 and separated by a small distance \e dt, then
00057 *! they are separated by a distance \e MM12 \e dt at point 2. \e MM21
00058 *! is defined similarly (with the geodesics being parallel to one
00059 *! another at point 2). On a flat surface, we have \e MM12 = \e MM21
00060 *! = 1.
00061 *! - \e a12 is the arc length on the auxiliary sphere. This is a
00062 *! construct for converting the problem to one in spherical
00063 *! trigonometry. \e a12 is measured in degrees. The spherical arc
00064 *! length from one equator crossing to the next is always 180°.
00065 *!
00066 *! If points 1, 2, and 3 lie on a single geodesic, then the following
00067 *! addition rules hold:
00068 *! - \e s13 = \e s12 + \e s23
00069 *! - \e a13 = \e a12 + \e a23
00070 *! - \e SS13 = \e SS12 + \e SS23
00071 *! - \e m13 = \e m12 \e MM23 + \e m23 \e MM21
00072 *! - \e MM13 = \e MM12 \e MM23 − (1 − \e MM12 \e MM21) \e
00073 *! m23 / \e m12
00074 *! - \e MM31 = \e MM32 \e MM21 − (1 − \e MM23 \e MM32) \e
00075 *! m12 / \e m23
00076 *!
00077 *! The shortest distance returned by the solution of the inverse problem
00078 *! is (obviously) uniquely defined. However, in a few special cases
00079 *! there are multiple azimuths which yield the same shortest distance.
00080 *! Here is a catalog of those cases:
00081 *! - \e lat1 = −\e lat2 (with neither point at a pole). If \e azi1 = \e
00082 *! azi2, the geodesic is unique. Otherwise there are two geodesics
00083 *! and the second one is obtained by setting [\e azi1, \e azi2] = [\e
00084 *! azi2, \e azi1], [\e MM12, \e MM21] = [\e MM21, \e MM12], \e SS12 =
00085 *! −\e SS12. (This occurs when the longitude difference is near
00086 *! ±180° for oblate ellipsoids.)
00087 *! - \e lon2 = \e lon1 ± 180° (with neither point at a pole). If
00088 *! \e azi1 = 0° or ±180°, the geodesic is unique.
00089 *! Otherwise there are two geodesics and the second one is obtained by
00090 *! setting [\e azi1, \e azi2] = [−\e azi1, −\e azi2], \e
00091 *! SS12 = −\e SS12. (This occurs when \e lat2 is near
00092 *! −\e lat1 for prolate ellipsoids.)
00093 *! - Points 1 and 2 at opposite poles. There are infinitely many
00094 *! geodesics which can be generated by setting [\e azi1, \e azi2] =
00095 *! [\e azi1, \e azi2] + [\e d, −\e d], for arbitrary \e d. (For
00096 *! spheres, this prescription applies when points 1 and 2 are
00097 *! antipodal.)
00098 *! - \e s12 = 0 (coincident points). There are infinitely many
00099 *! geodesics which can be generated by setting [\e azi1, \e azi2] =
00100 *! [\e azi1, \e azi2] + [\e d, \e d], for arbitrary \e d.
00101 *!
00102 *! These routines are a simple transcription of the corresponding C++
00103 *! classes in <a href="http:
00104 *! Because of the limitations of Fortran 77, the classes have been
00105 *! replaced by simple subroutines with no attempt to save "state" across
00106 *! subroutine calls. Most of the internal comments have been retained.
00107 *! However, in the process of transcription some documentation has been
00108 *! lost and the documentation for the C++ classes,
00109 *! GeographicLib::Geodesic, GeographicLib::GeodesicLine, and
00110 *! GeographicLib::PolygonArea, should be consulted. The C++ code
00111 *! remains the "reference implementation". Think twice about
00112 *! restructuring the internals of the Fortran code since this may make
00113 *! porting fixes from the C++ code more difficult.
00114 *!
00115 *! Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and
00116 *! licensed under the MIT/X11 License. For more information, see
00117 *! http:
00118 *!
00119 *! This library was distributed with
00120 *! <a href="../index.html">GeographicLib</a> 1.31.
00121
00122 *> Solve the direct geodesic problem
00123 *!
00124 *! @param[in] a the equatorial radius (meters).
00125 *! @param[in] f the flattening of the ellipsoid. Setting \e f = 0 gives
00126 *! a sphere. Negative \e f gives a prolate ellipsoid.
00127 *! @param[in] lat1 latitude of point 1 (degrees).
00128 *! @param[in] lon1 longitude of point 1 (degrees).
00129 *! @param[in] azi1 azimuth at point 1 (degrees).
00130 *! @param[in] s12a12 if \e arcmod is false, this is the distance between
00131 *! point 1 and point 2 (meters); otherwise it is the arc length
00132 *! between point 1 and point 2 (degrees); it can be negative.
00133 *! @param[in] arcmod logical flag determining the meaning of the \e
00134 *! s12a12.
00135 *! @param[out] lat2 latitude of point 2 (degrees).
00136 *! @param[out] lon2 longitude of point 2 (degrees).
00137 *! @param[out] azi2 (forward) azimuth at point 2 (degrees).
00138 *! @param[in] omask a bitor'ed combination of mask values
00139 *! specifying which of the following parameters should be set.
00140 *! @param[out] a12s12 if \e arcmod is false, this is the arc length
00141 *! between point 1 and point 2 (degrees); otherwise it is the distance
00142 *! between point 1 and point 2 (meters).
00143 *! @param[out] m12 reduced length of geodesic (meters).
00144 *! @param[out] MM12 geodesic scale of point 2 relative to point 1
00145 *! (dimensionless).
00146 *! @param[out] MM21 geodesic scale of point 1 relative to point 2
00147 *! (dimensionless).
00148 *! @param[out] SS12 area under the geodesic (meters<sup>2</sup>).
00149 *!
00150 *! \e omask is an integer in [0, 16) whose binary bits are interpreted
00151 *! as follows
00152 *! - 1 return \e a12
00153 *! - 2 return \e m12
00154 *! - 4 return \e MM12 and \e MM21
00155 *! - 8 return \e SS12
00156 *!
00157 *! \e lat1 should be in the range [−90°, 90°]; \e lon1 and
00158 *! \e azi1 should be in the range [−540°, 540°). The
00159 *! values of \e lon2 and \e azi2 returned are in the range
00160 *! [−180°, 180°).
00161 *!
00162 *! If either point is at a pole, the azimuth is defined by keeping the
00163 *! longitude fixed, writing \e lat = \e lat = ±(90° −
00164 *! ε), and taking the limit ε → 0+. An arc length
00165 *! greater that 180° signifies a geodesic which is not a shortest
00166 *! path. (For a prolate ellipsoid, an additional condition is necessary
00167 *! for a shortest path: the longitudinal extent must not exceed of
00168 *! 180°.)
00169
00170 subroutine direct(a, f, lat1, lon1, azi1, s12a12, arcmod,
00171 + lat2, lon2, azi2, omask, a12s12, m12, MM12, MM21, SS12)
00172 * input
00173 double precision a, f, lat1, lon1, azi1, s12a12
00174 logical arcmod
00175 integer omask
00176 * output
00177 double precision lat2, lon2, azi2
00178 * optional output
00179 double precision a12s12, m12, MM12, MM21, SS12
00180
00181 integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
00182 parameter (ord = 6, nC1 = ord, nC1p = ord,
00183 + nC2 = ord, nA3 = ord, nA3x = nA3,
00184 + nC3 = ord, nC3x = (nC3 * (nC3 - 1)) / 2,
00185 + nC4 = ord, nC4x = (nC4 * (nC4 + 1)) / 2)
00186 double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1),
00187 + C1a(nC1), C1pa(nC1p), C2a(nC2), C3a(nC3-1), C4a(0:nC4-1)
00188
00189 double precision csmgt, atanhx, hypotx,
00190 + AngNm, AngNm2, AngRnd, TrgSum, A1m1f, A2m1f, A3f
00191 logical arcp, redlp, scalp, areap
00192 double precision e2, f1, ep2, n, b, c2,
00193 + lon1x, azi1x, phi, alp1, salp0, calp0, k2, eps,
00194 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1,
00195 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2,
00196 + ssig12, csig12, salp12, calp12, omg12, lam12, lon12,
00197 + sig12, stau1, ctau1, tau12, s12a, t, s, c, serr,
00198 + A1m1, A2m1, A3c, A4, AB1, AB2,
00199 + B11, B12, B21, B22, B31, B41, B42, J12
00200
00201 double precision dblmin, dbleps, pi, degree, tiny,
00202 + tol0, tol1, tol2, tolb, xthrsh
00203 integer digits, maxit1, maxit2
00204 logical init
00205 common /geocom/ dblmin, dbleps, pi, degree, tiny,
00206 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
00207
00208 if (.not.init) call geoini
00209
00210 e2 = f * (2 - f)
00211 ep2 = e2 / (1 - e2)
00212 f1 = 1 - f
00213 n = f / (2 - f)
00214 b = a * f1
00215 c2 = 0
00216
00217 arcp = mod(omask/1, 2) == 1
00218 redlp = mod(omask/2, 2) == 1
00219 scalp = mod(omask/4, 2) == 1
00220 areap = mod(omask/8, 2) == 1
00221
00222 if (areap) then
00223 if (e2 .eq. 0) then
00224 c2 = a**2
00225 else if (e2 .gt. 0) then
00226 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
00227 else
00228 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
00229 end if
00230 end if
00231
00232 call A3cof(n, A3x)
00233 call C3cof(n, C3x)
00234 if (areap) call C4cof(n, C4x)
00235
00236 * Guard against underflow in salp0
00237 azi1x = AngRnd(AngNm(azi1))
00238 lon1x = AngNm(lon1)
00239
00240 * alp1 is in [0, pi]
00241 alp1 = azi1x * degree
00242 * Enforce sin(pi) == 0 and cos(pi/2) == 0. Better to face the ensuing
00243 * problems directly than to skirt them.
00244 salp1 = csmgt(0d0, sin(alp1), azi1x .eq. -180)
00245 calp1 = csmgt(0d0, cos(alp1), abs(azi1x) .eq. 90)
00246
00247 phi = lat1 * degree
00248 * Ensure cbet1 = +dbleps at poles
00249 sbet1 = f1 * sin(phi)
00250 cbet1 = csmgt(tiny, cos(phi), abs(lat1) .eq. 90)
00251 call Norm(sbet1, cbet1)
00252 dn1 = sqrt(1 + ep2 * sbet1**2)
00253
00254 * Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
00255 * alp0 in [0, pi/2 - |bet1|]
00256 salp0 = salp1 * cbet1
00257 * Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following
00258 * is slightly better (consider the case salp1 = 0).
00259 calp0 = hypotx(calp1, salp1 * sbet1)
00260 * Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
00261 * sig = 0 is nearest northward crossing of equator.
00262 * With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
00263 * With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
00264 * With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2
00265 * Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
00266 * With alp0 in (0, pi/2], quadrants for sig and omg coincide.
00267 * No atan2(0,0) ambiguity at poles since cbet1 = +dbleps.
00268 * With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
00269 ssig1 = sbet1
00270 somg1 = salp0 * sbet1
00271 csig1 = csmgt(cbet1 * calp1, 1d0, sbet1 .ne. 0 .or. calp1 .ne. 0)
00272 comg1 = csig1
00273 * sig1 in (-pi, pi]
00274 call Norm(ssig1, csig1)
00275 * Geodesic::Norm(somg1, comg1); -- don't need to normalize!
00276
00277 k2 = calp0**2 * ep2
00278 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
00279
00280 A1m1 = A1m1f(eps)
00281 call C1f(eps, C1a)
00282 B11 = TrgSum(.true., ssig1, csig1, C1a, nC1)
00283 s = sin(B11)
00284 c = cos(B11)
00285 * tau1 = sig1 + B11
00286 stau1 = ssig1 * c + csig1 * s
00287 ctau1 = csig1 * c - ssig1 * s
00288 * Not necessary because C1pa reverts C1a
00289 * B11 = -TrgSum(true, stau1, ctau1, C1pa, nC1p)
00290
00291 if (.not. arcmod) call C1pf(eps, C1pa)
00292
00293 if (redlp .or. scalp) then
00294 A2m1 = A2m1f(eps)
00295 call C2f(eps, C2a)
00296 B21 = TrgSum(.true., ssig1, csig1, C2a, nC2)
00297 else
00298 * Suppress bogus warnings about unitialized variables
00299 A2m1 = 0
00300 B21 = 0
00301 end if
00302
00303 call C3f(eps, C3x, C3a)
00304 A3c = -f * salp0 * A3f(eps, A3x)
00305 B31 = TrgSum(.true., ssig1, csig1, C3a, nC3-1)
00306
00307 if (areap) then
00308 call C4f(eps, C4x, C4a)
00309 * Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
00310 A4 = a**2 * calp0 * salp0 * e2
00311 B41 = TrgSum(.false., ssig1, csig1, C4a, nC4)
00312 else
00313 * Suppress bogus warnings about unitialized variables
00314 A4 = 0
00315 B41 = 0
00316 end if
00317
00318 if (arcmod) then
00319 * Interpret s12a12 as spherical arc length
00320 sig12 = s12a12 * degree
00321 s12a = abs(s12a12)
00322 s12a = s12a - 180 * aint(s12a / 180)
00323 ssig12 = csmgt(0d0, sin(sig12), s12a .eq. 0)
00324 csig12 = csmgt(0d0, cos(sig12), s12a .eq. 90)
00325 * Suppress bogus warnings about unitialized variables
00326 B12 = 0
00327 else
00328 * Interpret s12a12 as distance
00329 tau12 = s12a12 / (b * (1 + A1m1))
00330 s = sin(tau12)
00331 c = cos(tau12)
00332 * tau2 = tau1 + tau12
00333 B12 = - TrgSum(.true.,
00334 + stau1 * c + ctau1 * s, ctau1 * c - stau1 * s, C1pa, nC1p)
00335 sig12 = tau12 - (B12 - B11)
00336 ssig12 = sin(sig12)
00337 csig12 = cos(sig12)
00338 if (abs(f) .gt. 0.01d0) then
00339 * Reverted distance series is inaccurate for |f| > 1/100, so correct
00340 * sig12 with 1 Newton iteration. The following table shows the
00341 * approximate maximum error for a = WGS_a() and various f relative to
00342 * GeodesicExact.
00343 * erri = the error in the inverse solution (nm)
00344 * errd = the error in the direct solution (series only) (nm)
00345 * errda = the error in the direct solution (series + 1 Newton) (nm)
00346 *
00347 * f erri errd errda
00348 * -1/5 12e6 1.2e9 69e6
00349 * -1/10 123e3 12e6 765e3
00350 * -1/20 1110 108e3 7155
00351 * -1/50 18.63 200.9 27.12
00352 * -1/100 18.63 23.78 23.37
00353 * -1/150 18.63 21.05 20.26
00354 * 1/150 22.35 24.73 25.83
00355 * 1/100 22.35 25.03 25.31
00356 * 1/50 29.80 231.9 30.44
00357 * 1/20 5376 146e3 10e3
00358 * 1/10 829e3 22e6 1.5e6
00359 * 1/5 157e6 3.8e9 280e6
00360 ssig2 = ssig1 * csig12 + csig1 * ssig12
00361 csig2 = csig1 * csig12 - ssig1 * ssig12
00362 B12 = TrgSum(.true., ssig2, csig2, C1a, nC1)
00363 serr = (1 + A1m1) * (sig12 + (B12 - B11)) - s12a12 / b
00364 sig12 = sig12 - serr / sqrt(1 + k2 * ssig2**2)
00365 ssig12 = sin(sig12)
00366 csig12 = cos(sig12)
00367 * Update B12 below
00368 end if
00369 end if
00370
00371 * sig2 = sig1 + sig12
00372 ssig2 = ssig1 * csig12 + csig1 * ssig12
00373 csig2 = csig1 * csig12 - ssig1 * ssig12
00374 dn2 = sqrt(1 + k2 * ssig2**2)
00375 if (arcmod .or. abs(f) .gt. 0.01d0)
00376 + B12 = TrgSum(.true., ssig2, csig2, C1a, nC1)
00377 AB1 = (1 + A1m1) * (B12 - B11)
00378
00379 * sin(bet2) = cos(alp0) * sin(sig2)
00380 sbet2 = calp0 * ssig2
00381 * Alt: cbet2 = hypot(csig2, salp0 * ssig2)
00382 cbet2 = hypotx(salp0, calp0 * csig2)
00383 if (cbet2 .eq. 0) then
00384 * I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
00385 cbet2 = tiny
00386 csig2 = cbet2
00387 end if
00388 * tan(omg2) = sin(alp0) * tan(sig2)
00389 * No need to normalize
00390 somg2 = salp0 * ssig2
00391 comg2 = csig2
00392 * tan(alp0) = cos(sig2)*tan(alp2)
00393 * No need to normalize
00394 salp2 = salp0
00395 calp2 = calp0 * csig2
00396 * omg12 = omg2 - omg1
00397 omg12 = atan2(somg2 * comg1 - comg2 * somg1,
00398 + comg2 * comg1 + somg2 * somg1)
00399
00400 lam12 = omg12 + A3c *
00401 + ( sig12 + (TrgSum(.true., ssig2, csig2, C3a, nC3-1)
00402 + - B31))
00403 lon12 = lam12 / degree
00404 * Use Math::AngNm2 because longitude might have wrapped multiple
00405 * times.
00406 lon12 = AngNm2(lon12)
00407 lon2 = AngNm(lon1x + lon12)
00408 lat2 = atan2(sbet2, f1 * cbet2) / degree
00409 * minus signs give range [-180, 180). 0- converts -0 to +0.
00410 azi2 = 0 - atan2(-salp2, calp2) / degree
00411
00412 if (redlp .or. scalp) then
00413 B22 = TrgSum(.true., ssig2, csig2, C2a, nC2)
00414 AB2 = (1 + A2m1) * (B22 - B21)
00415 J12 = (A1m1 - A2m1) * sig12 + (AB1 - AB2)
00416 end if
00417 * Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
00418 * accurate cancellation in the case of coincident points.
00419 if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
00420 + dn1 * (ssig1 * csig2)) - csig1 * csig2 * J12)
00421 if (scalp) then
00422 t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
00423 MM12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1
00424 MM21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2
00425 end if
00426
00427 if (areap) then
00428 B42 = TrgSum(.false., ssig2, csig2, C4a, nC4)
00429 if (calp0 .eq. 0 .or. salp0 .eq. 0) then
00430 * alp12 = alp2 - alp1, used in atan2 so no need to normalized
00431 salp12 = salp2 * calp1 - calp2 * salp1
00432 calp12 = calp2 * calp1 + salp2 * salp1
00433 * The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
00434 * salp12 = -0 and alp12 = -180. However this depends on the sign being
00435 * attached to 0 correctly. The following ensures the correct behavior.
00436 if (salp12 .eq. 0 .and. calp12 .lt. 0) then
00437 salp12 = tiny * calp1
00438 calp12 = -1
00439 end if
00440 else
00441 * tan(alp) = tan(alp0) * sec(sig)
00442 * tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
00443 * = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
00444 * If csig12 > 0, write
00445 * csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
00446 * else
00447 * csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
00448 * No need to normalize
00449 salp12 = calp0 * salp0 *
00450 + csmgt(csig1 * (1 - csig12) + ssig12 * ssig1,
00451 + ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1),
00452 + csig12 .le. 0)
00453 calp12 = salp0**2 + calp0**2 * csig1 * csig2
00454 end if
00455 SS12 = c2 * atan2(salp12, calp12) + A4 * (B42 - B41)
00456 end if
00457
00458 if (arcp) a12s12 = csmgt(b * ((1 + A1m1) * sig12 + AB1),
00459 + sig12 / degree, arcmod)
00460
00461 return
00462 end
00463
00464 *> Solve the inverse geodesic problem.
00465 *!
00466 *! @param[in] a the equatorial radius (meters).
00467 *! @param[in] f the flattening of the ellipsoid. Setting \e f = 0 gives
00468 *! a sphere. Negative \e f gives a prolate ellipsoid.
00469 *! @param[in] lat1 latitude of point 1 (degrees).
00470 *! @param[in] lon1 longitude of point 1 (degrees).
00471 *! @param[in] lat2 latitude of point 2 (degrees).
00472 *! @param[in] lon2 longitude of point 2 (degrees).
00473 *! @param[out] s12 distance between point 1 and point 2 (meters).
00474 *! @param[out] azi1 azimuth at point 1 (degrees).
00475 *! @param[out] azi2 (forward) azimuth at point 2 (degrees).
00476 *! @param[in] omask a bitor'ed combination of mask values
00477 *! specifying which of the following parameters should be set.
00478 *! @param[out] a12 arc length of between point 1 and point 2 (degrees).
00479 *! @param[out] m12 reduced length of geodesic (meters).
00480 *! @param[out] MM12 geodesic scale of point 2 relative to point 1
00481 *! (dimensionless).
00482 *! @param[out] MM21 geodesic scale of point 1 relative to point 2
00483 *! (dimensionless).
00484 *! @param[out] SS12 area under the geodesic (meters<sup>2</sup>).
00485 *!
00486 *! \e omask is an integer in [0, 16) whose binary bits are interpreted
00487 *! as follows
00488 *! - 1 return \e a12
00489 *! - 2 return \e m12
00490 *! - 4 return \e MM12 and \e MM21
00491 *! - 8 return \e SS12
00492 *!
00493 *! \e lat1 and \e lat2 should be in the range [−90°, 90°];
00494 *! \e lon1 and \e lon2 should be in the range [−540°,
00495 *! 540°). The values of \e azi1 and \e azi2 returned are in the
00496 *! range [−180°, 180°).
00497 *!
00498 *! If either point is at a pole, the azimuth is defined by keeping the
00499 *! longitude fixed, writing \e lat = ±(90° −
00500 *! ε), and taking the limit ε → 0+.
00501 *!
00502 *! The solution to the inverse problem is found using Newton's method.
00503 *! If this fails to converge (this is very unlikely in geodetic
00504 *! applications but does occur for very eccentric ellipsoids), then the
00505 *! bisection method is used to refine the solution.
00506
00507 subroutine invers(a, f, lat1, lon1, lat2, lon2,
00508 + s12, azi1, azi2, omask, a12, m12, MM12, MM21, SS12)
00509 * input
00510 double precision a, f, lat1, lon1, lat2, lon2
00511 integer omask
00512 * output
00513 double precision s12, azi1, azi2
00514 * optional output
00515 double precision a12, m12, MM12, MM21, SS12
00516
00517 integer ord, nC1, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x
00518 parameter (ord = 6, nC1 = ord, nC2 = ord, nA3 = ord, nA3x = nA3,
00519 + nC3 = ord, nC3x = (nC3 * (nC3 - 1)) / 2,
00520 + nC4 = ord, nC4x = (nC4 * (nC4 + 1)) / 2)
00521 double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1),
00522 + C1a(nC1), C2a(nC2), C3a(nC3-1), C4a(0:nC4-1)
00523
00524 double precision csmgt, atanhx, hypotx,
00525 + AngNm, AngDif, AngRnd, TrgSum, Lam12f, InvSta
00526 integer latsgn, lonsgn, swapp, numit
00527 logical arcp, redlp, scalp, areap, merid, tripn, tripb
00528
00529 double precision e2, f1, ep2, n, b, c2,
00530 + lat1x, lat2x, phi, salp0, calp0, k2, eps,
00531 + salp1, calp1, ssig1, csig1, cbet1, sbet1, dbet1, dn1,
00532 + salp2, calp2, ssig2, csig2, sbet2, cbet2, dbet2, dn2,
00533 + slam12, clam12, salp12, calp12, omg12, lam12, lon12,
00534 + salp1a, calp1a, salp1b, calp1b,
00535 + dalp1, sdalp1, cdalp1, nsalp1, alp12, somg12, domg12,
00536 + sig12, v, dv, dnm, dummy,
00537 + A4, B41, B42, s12x, m12x, a12x
00538
00539 double precision dblmin, dbleps, pi, degree, tiny,
00540 + tol0, tol1, tol2, tolb, xthrsh
00541 integer digits, maxit1, maxit2
00542 logical init
00543 common /geocom/ dblmin, dbleps, pi, degree, tiny,
00544 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
00545
00546 if (.not.init) call geoini
00547
00548 f1 = 1 - f
00549 e2 = f * (2 - f)
00550 ep2 = e2 / f1**2
00551 n = f / ( 2 - f)
00552 b = a * f1
00553 c2 = 0
00554
00555 arcp = mod(omask/1, 2) == 1
00556 redlp = mod(omask/2, 2) == 1
00557 scalp = mod(omask/4, 2) == 1
00558 areap = mod(omask/8, 2) == 1
00559
00560 if (areap) then
00561 if (e2 .eq. 0) then
00562 c2 = a**2
00563 else if (e2 .gt. 0) then
00564 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
00565 else
00566 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
00567 end if
00568 end if
00569
00570 call A3cof(n, A3x)
00571 call C3cof(n, C3x)
00572 if (areap) call C4cof(n, C4x)
00573
00574 * Compute longitude difference (AngDiff does this carefully). Result is
00575 * in [-180, 180] but -180 is only for west-going geodesics. 180 is for
00576 * east-going and meridional geodesics.
00577 lon12 = AngDif(AngNm(lon1), AngNm(lon2))
00578 * If very close to being on the same half-meridian, then make it so.
00579 lon12 = AngRnd(lon12)
00580 * Make longitude difference positive.
00581 if (lon12 .ge. 0) then
00582 lonsgn = 1
00583 else
00584 lonsgn = -1
00585 end if
00586 lon12 = lon12 * lonsgn
00587 * If really close to the equator, treat as on equator.
00588 lat1x = AngRnd(lat1)
00589 lat2x = AngRnd(lat2)
00590 * Swap points so that point with higher (abs) latitude is point 1
00591 if (abs(lat1x) .ge. abs(lat2x)) then
00592 swapp = 1
00593 else
00594 swapp = -1
00595 end if
00596 if (swapp .lt. 0) then
00597 lonsgn = -lonsgn
00598 call swap(lat1x, lat2x)
00599 end if
00600 * Make lat1 <= 0
00601 if (lat1x .lt. 0) then
00602 latsgn = 1
00603 else
00604 latsgn = -1
00605 end if
00606 lat1x = lat1x * latsgn
00607 lat2x = lat2x * latsgn
00608 * Now we have
00609 *
00610 * 0 <= lon12 <= 180
00611 * -90 <= lat1 <= 0
00612 * lat1 <= lat2 <= -lat1
00613 *
00614 * longsign, swapp, latsgn register the transformation to bring the
00615 * coordinates to this canonical form. In all cases, 1 means no change
00616 * was made. We make these transformations so that there are few cases
00617 * to check, e.g., on verifying quadrants in atan2. In addition, this
00618 * enforces some symmetries in the results returned.
00619
00620 phi = lat1x * degree
00621 * Ensure cbet1 = +dbleps at poles
00622 sbet1 = f1 * sin(phi)
00623 cbet1 = csmgt(tiny, cos(phi), lat1x .eq. -90)
00624 call Norm(sbet1, cbet1)
00625
00626 phi = lat2x * degree
00627 * Ensure cbet2 = +dbleps at poles
00628 sbet2 = f1 * sin(phi)
00629 cbet2 = csmgt(tiny, cos(phi), abs(lat2x) .eq. 90)
00630 call Norm(sbet2, cbet2)
00631
00632 * If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
00633 * |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1
00634 * is a better measure. This logic is used in assigning calp2 in
00635 * Lambda12. Sometimes these quantities vanish and in that case we force
00636 * bet2 = +/- bet1 exactly. An example where is is necessary is the
00637 * inverse problem 48.522876735459 0 -48.52287673545898293
00638 * 179.599720456223079643 which failed with Visual Studio 10 (Release and
00639 * Debug)
00640
00641 if (cbet1 .lt. -sbet1) then
00642 if (cbet2 .eq. cbet1) sbet2 = sign(sbet1, sbet2)
00643 else
00644 if (abs(sbet2) .eq. -sbet1) cbet2 = cbet1
00645 end if
00646
00647 dn1 = sqrt(1 + ep2 * sbet1**2)
00648 dn2 = sqrt(1 + ep2 * sbet2**2)
00649
00650 lam12 = lon12 * degree
00651 slam12 = sin(lam12)
00652 if (lon12 .eq. 180) slam12 = 0
00653 * lon12 == 90 isn't interesting
00654 clam12 = cos(lam12)
00655
00656 * Suppress bogus warnings about unitialized variables
00657 a12x = 0
00658 merid = lat1x .eq. -90 .or. slam12 == 0
00659
00660 if (merid) then
00661
00662 * Endpoints are on a single full meridian, so the geodesic might lie on
00663 * a meridian.
00664
00665 * Head to the target longitude
00666 calp1 = clam12
00667 salp1 = slam12
00668 * At the target we're heading north
00669 calp2 = 1
00670 salp2 = 0
00671
00672 * tan(bet) = tan(sig) * cos(alp)
00673 ssig1 = sbet1
00674 csig1 = calp1 * cbet1
00675 ssig2 = sbet2
00676 csig2 = calp2 * cbet2
00677
00678 * sig12 = sig2 - sig1
00679 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
00680 + csig1 * csig2 + ssig1 * ssig2)
00681 call Lengs(n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
00682 + cbet1, cbet2, s12x, m12x, dummy,
00683 + scalp, MM12, MM21, ep2, C1a, C2a)
00684
00685 * Add the check for sig12 since zero length geodesics might yield m12 <
00686 * 0. Test case was
00687 *
00688 * echo 20.001 0 20.001 0 | GeodSolve -i
00689 *
00690 * In fact, we will have sig12 > pi/2 for meridional geodesic which is
00691 * not a shortest path.
00692 if (sig12 .lt. 1 .or. m12x .ge. 0) then
00693 m12x = m12x * b
00694 s12x = s12x * b
00695 a12x = sig12 / degree
00696 else
00697 * m12 < 0, i.e., prolate and too close to anti-podal
00698 merid = .false.
00699 end if
00700 end if
00701
00702 * Mimic the way Lambda12 works with calp1 = 0
00703 if (.not. merid .and. sbet1 .eq. 0 .and.
00704 + (f .le. 0 .or. lam12 .le. pi - f * pi)) then
00705
00706 * Geodesic runs along equator
00707 calp1 = 0
00708 calp2 = 0
00709 salp1 = 1
00710 salp2 = 1
00711 s12x = a * lam12
00712 sig12 = lam12 / f1
00713 omg12 = sig12
00714 m12x = b * sin(sig12)
00715 if (scalp) then
00716 MM12 = cos(sig12)
00717 MM21 = MM12
00718 end if
00719 a12x = lon12 / f1
00720 else if (.not. merid) then
00721 * Now point1 and point2 belong within a hemisphere bounded by a
00722 * meridian and geodesic is neither meridional or equatorial.
00723
00724 * Figure a starting point for Newton's method
00725 sig12 = InvSta(sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
00726 + f, A3x, salp1, calp1, salp2, calp2, dnm, C1a, C2a)
00727
00728 if (sig12 .ge. 0) then
00729 * Short lines (InvSta sets salp2, calp2, dnm)
00730 s12x = sig12 * b * dnm
00731 m12x = dnm**2 * b * sin(sig12 / dnm)
00732 if (scalp) then
00733 MM12 = cos(sig12 / dnm)
00734 MM21 = MM12
00735 end if
00736 a12x = sig12 / degree
00737 omg12 = lam12 / (f1 * dnm)
00738 else
00739
00740 * Newton's method. This is a straightforward solution of f(alp1) =
00741 * lambda12(alp1) - lam12 = 0 with one wrinkle. f(alp) has exactly one
00742 * root in the interval (0, pi) and its derivative is positive at the
00743 * root. Thus f(alp) is positive for alp > alp1 and negative for alp <
00744 * alp1. During the course of the iteration, a range (alp1a, alp1b) is
00745 * maintained which brackets the root and with each evaluation of
00746 * f(alp) the range is shrunk, if possible. Newton's method is
00747 * restarted whenever the derivative of f is negative (because the new
00748 * value of alp1 is then further from the solution) or if the new
00749 * estimate of alp1 lies outside (0,pi); in this case, the new starting
00750 * guess is taken to be (alp1a + alp1b) / 2.
00751
00752 * Bracketing range
00753 salp1a = tiny
00754 calp1a = 1
00755 salp1b = tiny
00756 calp1b = -1
00757 tripn = .false.
00758 tripb = .false.
00759 do 10 numit = 0, maxit2-1
00760 * the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
00761 * WGS84 and random input: mean = 2.85, sd = 0.60
00762 v = Lam12f(sbet1, cbet1, dn1, sbet2, cbet2, dn2,
00763 + salp1, calp1, f, A3x, C3x, salp2, calp2, sig12,
00764 + ssig1, csig1, ssig2, csig2,
00765 + eps, omg12, numit .lt. maxit1, dv,
00766 + C1a, C2a, C3a) - lam12
00767 * 2 * tol0 is approximately 1 ulp for a number in [0, pi].
00768 * Reversed test to allow escape with NaNs
00769 if (tripb .or.
00770 + .not. (abs(v) .ge. csmgt(8d0, 2d0, tripn) * tol0))
00771 + go to 20
00772 * Update bracketing values
00773 if (v .gt. 0 .and. (numit .gt. maxit1 .or.
00774 + calp1/salp1 .gt. calp1b/salp1b)) then
00775 salp1b = salp1
00776 calp1b = calp1
00777 else if (v .lt. 0 .and. (numit .gt. maxit1 .or.
00778 + calp1/salp1 .lt. calp1a/salp1a)) then
00779 salp1a = salp1
00780 calp1a = calp1
00781 end if
00782 if (numit .lt. maxit1 .and. dv .gt. 0) then
00783 dalp1 = -v/dv
00784 sdalp1 = sin(dalp1)
00785 cdalp1 = cos(dalp1)
00786 nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
00787 if (nsalp1 .gt. 0 .and. abs(dalp1) .lt. pi) then
00788 calp1 = calp1 * cdalp1 - salp1 * sdalp1
00789 salp1 = nsalp1
00790 call Norm(salp1, calp1)
00791 * In some regimes we don't get quadratic convergence because
00792 * slope -> 0. So use convergence conditions based on dbleps
00793 * instead of sqrt(dbleps).
00794 tripn = abs(v) .le. 16 * tol0
00795 go to 10
00796 end if
00797 end if
00798 * Either dv was not postive or updated value was outside legal
00799 * range. Use the midpoint of the bracket as the next estimate.
00800 * This mechanism is not needed for the WGS84 ellipsoid, but it does
00801 * catch problems with more eccentric ellipsoids. Its efficacy is
00802 * such for the WGS84 test set with the starting guess set to alp1 =
00803 * 90deg:
00804 * the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
00805 * WGS84 and random input: mean = 4.74, sd = 0.99
00806 salp1 = (salp1a + salp1b)/2
00807 calp1 = (calp1a + calp1b)/2
00808 call Norm(salp1, calp1)
00809 tripn = .false.
00810 tripb = abs(salp1a - salp1) + (calp1a - calp1) .lt. tolb
00811 + .or. abs(salp1 - salp1b) + (calp1 - calp1b) .lt. tolb
00812 10 continue
00813 20 continue
00814 call Lengs(eps, sig12, ssig1, csig1, dn1,
00815 + ssig2, csig2, dn2, cbet1, cbet2, s12x, m12x, dummy,
00816 + scalp, MM12, MM21, ep2, C1a, C2a)
00817 m12x = m12x * b
00818 s12x = s12x * b
00819 a12x = sig12 / degree
00820 omg12 = lam12 - omg12
00821 end if
00822 end if
00823
00824 * Convert -0 to 0
00825 s12 = 0 + s12x
00826 if (redlp) m12 = 0 + m12x
00827
00828 if (areap) then
00829 * From Lambda12: sin(alp1) * cos(bet1) = sin(alp0)
00830 salp0 = salp1 * cbet1
00831 calp0 = hypotx(calp1, salp1 * sbet1)
00832 if (calp0 .ne. 0 .and. salp0 .ne. 0) then
00833 * From Lambda12: tan(bet) = tan(sig) * cos(alp)
00834 ssig1 = sbet1
00835 csig1 = calp1 * cbet1
00836 ssig2 = sbet2
00837 csig2 = calp2 * cbet2
00838 k2 = calp0**2 * ep2
00839 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
00840 * Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0).
00841 A4 = a**2 * calp0 * salp0 * e2
00842 call Norm(ssig1, csig1)
00843 call Norm(ssig2, csig2)
00844 call C4f(eps, C4x, C4a)
00845 B41 = TrgSum(.false., ssig1, csig1, C4a, nC4)
00846 B42 = TrgSum(.false., ssig2, csig2, C4a, nC4)
00847 SS12 = A4 * (B42 - B41)
00848 else
00849 * Avoid problems with indeterminate sig1, sig2 on equator
00850 SS12 = 0
00851 end if
00852
00853 if (.not. merid .and. omg12 .lt. 0.75d0 * pi
00854 + .and. sbet2 - sbet1 .lt. 1.75d0) then
00855 * Use tan(Gamma/2) = tan(omg12/2)
00856 * * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
00857 * with tan(x/2) = sin(x)/(1+cos(x))
00858 somg12 = sin(omg12)
00859 domg12 = 1 + cos(omg12)
00860 dbet1 = 1 + cbet1
00861 dbet2 = 1 + cbet2
00862 alp12 = 2 * atan2(somg12 * (sbet1 * dbet2 + sbet2 * dbet1),
00863 + domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
00864 else
00865 * alp12 = alp2 - alp1, used in atan2 so no need to normalize
00866 salp12 = salp2 * calp1 - calp2 * salp1
00867 calp12 = calp2 * calp1 + salp2 * salp1
00868 * The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
00869 * salp12 = -0 and alp12 = -180. However this depends on the sign
00870 * being attached to 0 correctly. The following ensures the correct
00871 * behavior.
00872 if (salp12 .eq. 0 .and. calp12 .lt. 0) then
00873 salp12 = tiny * calp1
00874 calp12 = -1
00875 end if
00876 alp12 = atan2(salp12, calp12)
00877 end if
00878 SS12 = SS12 + c2 * alp12
00879 SS12 = SS12 * swapp * lonsgn * latsgn
00880 * Convert -0 to 0
00881 SS12 = 0 + SS12
00882 end if
00883
00884 * Convert calp, salp to azimuth accounting for lonsgn, swapp, latsgn.
00885 if (swapp .lt. 0) then
00886 call swap(salp1, salp2)
00887 call swap(calp1, calp2)
00888 if (scalp) call swap(MM12, MM21)
00889 end if
00890
00891 salp1 = salp1 * swapp * lonsgn
00892 calp1 = calp1 * swapp * latsgn
00893 salp2 = salp2 * swapp * lonsgn
00894 calp2 = calp2 * swapp * latsgn
00895
00896 * minus signs give range [-180, 180). 0- converts -0 to +0.
00897 azi1 = 0 - atan2(-salp1, calp1) / degree
00898 azi2 = 0 - atan2(-salp2, calp2) / degree
00899
00900 if (arcp) a12 = a12x
00901
00902 return
00903 end
00904
00905 *> Determine the area of a geodesic polygon
00906 *!
00907 *! @param[in] a the equatorial radius (meters).
00908 *! @param[in] f the flattening of the ellipsoid. Setting \e f = 0 gives
00909 *! a sphere. Negative \e f gives a prolate ellipsoid.
00910 *! @param[in] lats an array of the latitudes of the vertices (degrees).
00911 *! @param[in] lons an array of the longitudes of the vertices (degrees).
00912 *! @param[in] n the number of vertices
00913 *! @param[out] AA the (signed) area of the polygon (meters<sup>2</sup>)
00914 *! @param[out] PP the perimeter of the polygon
00915 *!
00916 *! \e lats should be in the range [−90°, 90°]; \e lons
00917 *! should be in the range [−540°, 540°).
00918 *!
00919 *! Only simple polygons (which are not self-intersecting) are allowed.
00920 *! There's no need to "close" the polygon by repeating the first vertex.
00921 *! The area returned is signed with counter-clockwise traversal being
00922 *! treated as positive.
00923
00924 subroutine area(a, f, lats, lons, n, AA, PP)
00925 * input
00926 integer n
00927 double precision a, f, lats(n), lons(n)
00928 * output
00929 double precision AA, PP
00930
00931 integer i, omask, cross, trnsit
00932 double precision s12, azi1, azi2, dummy, SS12, b, e2, c2, area0,
00933 + atanhx, Aacc(2), Pacc(2)
00934
00935 double precision dblmin, dbleps, pi, degree, tiny,
00936 + tol0, tol1, tol2, tolb, xthrsh
00937 integer digits, maxit1, maxit2
00938 logical init
00939 common /geocom/ dblmin, dbleps, pi, degree, tiny,
00940 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
00941
00942 omask = 8
00943 call accini(Aacc)
00944 call accini(Pacc)
00945 cross = 0
00946 do 10 i = 0, n-1
00947 call invers(a, f, lats(i+1), lons(i+1),
00948 + lats(mod(i+1,n)+1), lons(mod(i+1,n)+1),
00949 + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, SS12)
00950 call accadd(Pacc, s12)
00951 call accadd(Aacc, -SS12)
00952 cross = cross + trnsit(lons(i+1), lons(mod(i+1,n)+1))
00953 10 continue
00954 PP = Pacc(1)
00955 b = a * (1 - f)
00956 e2 = f * (2 - f)
00957 if (e2 .eq. 0) then
00958 c2 = a**2
00959 else if (e2 .gt. 0) then
00960 c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2
00961 else
00962 c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2
00963 end if
00964 area0 = 4 * pi * c2
00965 if (mod(abs(cross), 2) .eq. 1) then
00966 if (Aacc(1) .lt. 0) then
00967 call accadd(Aacc, +area0/2)
00968 else
00969 call accadd(Aacc, -area0/2)
00970 end if
00971 end if
00972 if (Aacc(1) .gt. area0/2) then
00973 call accadd(Aacc, -area0)
00974 else if (Aacc(1) .le. -area0/2) then
00975 call accadd(Aacc, +area0)
00976 end if
00977 AA = Aacc(1)
00978
00979 return
00980 end
00981
00982 *> @cond SKIP
00983
00984 block data geodat
00985 double precision dblmin, dbleps, pi, degree, tiny,
00986 + tol0, tol1, tol2, tolb, xthrsh
00987 integer digits, maxit1, maxit2
00988 logical init
00989 data init /.false./
00990 common /geocom/ dblmin, dbleps, pi, degree, tiny,
00991 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
00992 end
00993
00994 subroutine geoini
00995 double precision dblmin, dbleps, pi, degree, tiny,
00996 + tol0, tol1, tol2, tolb, xthrsh
00997 integer digits, maxit1, maxit2
00998 logical init
00999 common /geocom/ dblmin, dbleps, pi, degree, tiny,
01000 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
01001
01002 digits = 53
01003 dblmin = 0.5d0**1022
01004 dbleps = 0.5d0**(digits-1)
01005
01006 pi = atan2(0.0d0, -1.0d0)
01007 degree = pi/180
01008 tiny = sqrt(dblmin)
01009 tol0 = dbleps
01010 * Increase multiplier in defn of tol1 from 100 to 200 to fix inverse
01011 * case 52.784459512564 0 -52.784459512563990912 179.634407464943777557
01012 * which otherwise failed for Visual Studio 10 (Release and Debug)
01013 tol1 = 200 * tol0
01014 tol2 = sqrt(tol0)
01015 * Check on bisection interval
01016 tolb = tol0 * tol2
01017 xthrsh = 1000 * tol2
01018 maxit1 = 20
01019 maxit2 = maxit1 + digits + 10
01020
01021 init = .true.
01022
01023 return
01024 end
01025
01026 subroutine Lengs(eps, sig12,
01027 + ssig1, csig1, dn1, ssig2, csig2, dn2,
01028 + cbet1, cbet2, s12b, m12b, m0,
01029 + scalp, MM12, MM21, ep2, C1a, C2a)
01030 * input
01031 double precision eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
01032 + cbet1, cbet2, ep2
01033 logical scalp
01034 * output
01035 double precision s12b, m12b, m0
01036 * optional output
01037 double precision MM12, MM21
01038 * temporary storage
01039 double precision C1a(*), C2a(*)
01040
01041 integer ord, nC1, nC2
01042 parameter (ord = 6, nC1 = ord, nC2 = ord)
01043
01044 double precision A1m1f, A2m1f, TrgSum
01045 double precision A1m1, AB1, A2m1, AB2, J12, csig12, t
01046
01047 * Return m12b = (reduced length)/b; also calculate s12b = distance/b,
01048 * and m0 = coefficient of secular term in expression for reduced length.
01049 call C1f(eps, C1a)
01050 call C2f(eps, C2a)
01051
01052 A1m1 = A1m1f(eps)
01053 AB1 = (1 + A1m1) * (TrgSum(.true., ssig2, csig2, C1a, nC1) -
01054 + TrgSum(.true., ssig1, csig1, C1a, nC1))
01055 A2m1 = A2m1f(eps)
01056 AB2 = (1 + A2m1) * (TrgSum(.true., ssig2, csig2, C2a, nC2) -
01057 + TrgSum(.true., ssig1, csig1, C2a, nC2))
01058 m0 = A1m1 - A2m1
01059 J12 = m0 * sig12 + (AB1 - AB2)
01060 * Missing a factor of b.
01061 * Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
01062 * accurate cancellation in the case of coincident points.
01063 m12b = dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
01064 + csig1 * csig2 * J12
01065 * Missing a factor of b
01066 s12b = (1 + A1m1) * sig12 + AB1
01067 if (scalp) then
01068 csig12 = csig1 * csig2 + ssig1 * ssig2
01069 t = ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
01070 MM12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1
01071 MM21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2
01072 end if
01073
01074 return
01075 end
01076
01077 double precision function Astrd(x, y)
01078 * Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
01079 * This solution is adapted from Geocentric::Reverse.
01080 * input
01081 double precision x, y
01082
01083 double precision cbrt, csmgt
01084 double precision k, p, q, r, S, r2, r3, disc, u,
01085 + T3, T, ang, v, uv, w
01086
01087 p = x**2
01088 q = y**2
01089 r = (p + q - 1) / 6
01090 if ( .not. (q .eq. 0 .and. r .lt. 0) ) then
01091 * Avoid possible division by zero when r = 0 by multiplying equations
01092 * for s and t by r^3 and r, resp.
01093 * S = r^3 * s
01094 S = p * q / 4
01095 r2 = r**2
01096 r3 = r * r2
01097 * The discrimant of the quadratic equation for T3. This is zero on
01098 * the evolute curve p^(1/3)+q^(1/3) = 1
01099 disc = S * (S + 2 * r3)
01100 u = r
01101 if (disc .ge. 0) then
01102 T3 = S + r3
01103 * Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
01104 * of precision due to cancellation. The result is unchanged because
01105 * of the way the T is used in definition of u.
01106 * T3 = (r * t)^3
01107 T3 = T3 + csmgt(-sqrt(disc), sqrt(disc), T3 .lt. 0)
01108 * N.B. cbrt always returns the real root. cbrt(-8) = -2.
01109 * T = r * t
01110 T = cbrt(T3)
01111 * T can be zero; but then r2 / T -> 0.
01112 if (T .ne. 0) u = u + T + r2 / T
01113 else
01114 * T is complex, but the way u is defined the result is real.
01115 ang = atan2(sqrt(-disc), -(S + r3))
01116 * There are three possible cube roots. We choose the root which
01117 * avoids cancellation. Note that disc < 0 implies that r < 0.
01118 u = u + 2 * r * cos(ang / 3)
01119 end if
01120 * guaranteed positive
01121 v = sqrt(u**2 + q)
01122 * Avoid loss of accuracy when u < 0.
01123 * u+v, guaranteed positive
01124 uv = csmgt(q / (v - u), u + v, u .lt. 0)
01125 * positive?
01126 w = (uv - q) / (2 * v)
01127 * Rearrange expression for k to avoid loss of accuracy due to
01128 * subtraction. Division by 0 not possible because uv > 0, w >= 0.
01129 * guaranteed positive
01130 k = uv / (sqrt(uv + w**2) + w)
01131 else
01132 * q == 0 && r <= 0
01133 * y = 0 with |x| <= 1. Handle this case directly.
01134 * for y small, positive root is k = abs(y)/sqrt(1-x^2)
01135 k = 0
01136 end if
01137 Astrd = k
01138
01139 return
01140 end
01141
01142 double precision function InvSta(sbet1, cbet1, dn1,
01143 + sbet2, cbet2, dn2, lam12, f, A3x,
01144 + salp1, calp1, salp2, calp2, dnm,
01145 + C1a, C2a)
01146 * Return a starting point for Newton's method in salp1 and calp1
01147 * (function value is -1). If Newton's method doesn't need to be used,
01148 * return also salp2, calp2, and dnm and function value is sig12.
01149 * input
01150 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
01151 + f, A3x(*)
01152 * output
01153 double precision salp1, calp1, salp2, calp2, dnm
01154 * temporary
01155 double precision C1a(*), C2a(*)
01156
01157 double precision csmgt, hypotx, A3f, Astrd
01158 logical shortp
01159 double precision f1, e2, ep2, n, etol2, k2, eps, sig12,
01160 + sbet12, cbet12, sbt12a, omg12, somg12, comg12, ssig12, csig12,
01161 + x, y, lamscl, betscl, cbt12a, bt12a, m12b, m0, dummy,
01162 + k, omg12a, sbetm2
01163
01164 double precision dblmin, dbleps, pi, degree, tiny,
01165 + tol0, tol1, tol2, tolb, xthrsh
01166 integer digits, maxit1, maxit2
01167 logical init
01168 common /geocom/ dblmin, dbleps, pi, degree, tiny,
01169 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
01170
01171 f1 = 1 - f
01172 e2 = f * (2 - f)
01173 ep2 = e2 / (1 - e2)
01174 n = f / (2 - f)
01175 * The sig12 threshold for "really short". Using the auxiliary sphere
01176 * solution with dnm computed at (bet1 + bet2) / 2, the relative error in
01177 * the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
01178 * (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a
01179 * given f and sig12, the max error occurs for lines near the pole. If
01180 * the old rule for computing dnm = (dn1 + dn2)/2 is used, then the error
01181 * increases by a factor of 2.) Setting this equal to epsilon gives
01182 * sig12 = etol2. Here 0.1 is a safety factor (error decreased by 100)
01183 * and max(0.001, abs(f)) stops etol2 getting too large in the nearly
01184 * spherical case.
01185 etol2 = 0.1d0 * tol2 /
01186 + sqrt( max(0.001d0, abs(f)) * min(1d0, 1 - f/2) / 2 )
01187
01188 * Return value
01189 sig12 = -1
01190 * bet12 = bet2 - bet1 in [0, pi); bt12a = bet2 + bet1 in (-pi, 0]
01191 sbet12 = sbet2 * cbet1 - cbet2 * sbet1
01192 cbet12 = cbet2 * cbet1 + sbet2 * sbet1
01193 sbt12a = sbet2 * cbet1 + cbet2 * sbet1
01194
01195 shortp = cbet12 .ge. 0 .and. sbet12 .lt. 0.5d0 .and.
01196 + cbet2 * lam12 .lt. 0.5d0
01197
01198 omg12 = lam12
01199 if (shortp) then
01200 sbetm2 = (sbet1 + sbet2)**2
01201 * sin((bet1+bet2)/2)^2
01202 * = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
01203 sbetm2 = sbetm2 / (sbetm2 + (cbet1 + cbet2)**2)
01204 dnm = sqrt(1 + ep2 * sbetm2)
01205 omg12 = omg12 / (f1 * dnm)
01206 end if
01207 somg12 = sin(omg12)
01208 comg12 = cos(omg12)
01209
01210 salp1 = cbet2 * somg12
01211 calp1 = csmgt(sbet12 + cbet2 * sbet1 * somg12**2 / (1 + comg12),
01212 + sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12),
01213 + comg12 .ge. 0)
01214
01215 ssig12 = hypotx(salp1, calp1)
01216 csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
01217
01218 if (shortp .and. ssig12 .lt. etol2) then
01219 * really short lines
01220 salp2 = cbet1 * somg12
01221 calp2 = sbet12 - cbet1 * sbet2 *
01222 + csmgt(somg12**2 / (1 + comg12), 1 - comg12, comg12 .ge. 0)
01223 call Norm(salp2, calp2)
01224 * Set return value
01225 sig12 = atan2(ssig12, csig12)
01226 else if (abs(n) .gt. 0.1d0 .or. csig12 .ge. 0 .or.
01227 + ssig12 .ge. 6 * abs(n) * pi * cbet1**2) then
01228 * Nothing to do, zeroth order spherical approximation is OK
01229 continue
01230 else
01231 * Scale lam12 and bet2 to x, y coordinate system where antipodal point
01232 * is at origin and singular point is at y = 0, x = -1.
01233 if (f .ge. 0) then
01234 * x = dlong, y = dlat
01235 k2 = sbet1**2 * ep2
01236 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
01237 lamscl = f * cbet1 * A3f(eps, A3x) * pi
01238 betscl = lamscl * cbet1
01239 x = (lam12 - pi) / lamscl
01240 y = sbt12a / betscl
01241 else
01242 * f < 0: x = dlat, y = dlong
01243 cbt12a = cbet2 * cbet1 - sbet2 * sbet1
01244 bt12a = atan2(sbt12a, cbt12a)
01245 * In the case of lon12 = 180, this repeats a calculation made in
01246 * Inverse.
01247 call Lengs(n, pi + bt12a,
01248 + sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
01249 + cbet1, cbet2, dummy, m12b, m0, .false.,
01250 + dummy, dummy, ep2, C1a, C2a)
01251 x = -1 + m12b / (cbet1 * cbet2 * m0 * pi)
01252 betscl = csmgt(sbt12a / x, -f * cbet1**2 * pi,
01253 + x .lt. -0.01d0)
01254 lamscl = betscl / cbet1
01255 y = (lam12 - pi) / lamscl
01256 end if
01257
01258 if (y .gt. -tol1 .and. x .gt. -1 - xthrsh) then
01259 * strip near cut
01260 if (f .ge. 0) then
01261 salp1 = min(1d0, -x)
01262 calp1 = - sqrt(1 - salp1**2)
01263 else
01264 calp1 = max(csmgt(0d0, 1d0, x .gt. -tol1), x)
01265 salp1 = sqrt(1 - calp1**2)
01266 end if
01267 else
01268 * Estimate alp1, by solving the astroid problem.
01269 *
01270 * Could estimate alpha1 = theta + pi/2, directly, i.e.,
01271 * calp1 = y/k; salp1 = -x/(1+k); for f >= 0
01272 * calp1 = x/(1+k); salp1 = -y/k; for f < 0 (need to check)
01273 *
01274 * However, it's better to estimate omg12 from astroid and use
01275 * spherical formula to compute alp1. This reduces the mean number of
01276 * Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
01277 * (min 0 max 5). The changes in the number of iterations are as
01278 * follows:
01279 *
01280 * change percent
01281 * 1 5
01282 * 0 78
01283 * -1 16
01284 * -2 0.6
01285 * -3 0.04
01286 * -4 0.002
01287 *
01288 * The histogram of iterations is (m = number of iterations estimating
01289 * alp1 directly, n = number of iterations estimating via omg12, total
01290 * number of trials = 148605):
01291 *
01292 * iter m n
01293 * 0 148 186
01294 * 1 13046 13845
01295 * 2 93315 102225
01296 * 3 36189 32341
01297 * 4 5396 7
01298 * 5 455 1
01299 * 6 56 0
01300 *
01301 * Because omg12 is near pi, estimate work with omg12a = pi - omg12
01302 k = Astrd(x, y)
01303 omg12a = lamscl *
01304 + csmgt(-x * k/(1 + k), -y * (1 + k)/k, f .ge. 0)
01305 somg12 = sin(omg12a)
01306 comg12 = -cos(omg12a)
01307 * Update spherical estimate of alp1 using omg12 instead of lam12
01308 salp1 = cbet2 * somg12
01309 calp1 = sbt12a - cbet2 * sbet1 * somg12**2 / (1 - comg12)
01310 end if
01311 end if
01312 * Sanity check on starting guess
01313 if (salp1 .gt. 0) then
01314 call Norm(salp1, calp1)
01315 else
01316 salp1 = 1
01317 calp1 = 0
01318 end if
01319 InvSta = sig12
01320
01321 return
01322 end
01323
01324 double precision function Lam12f(sbet1, cbet1, dn1,
01325 + sbet2, cbet2, dn2, salp1, calp1, f, A3x, C3x, salp2, calp2,
01326 + sig12, ssig1, csig1, ssig2, csig2, eps, domg12, diffp, dlam12,
01327 + C1a, C2a, C3a)
01328 * input
01329 double precision sbet1, cbet1, dn1, sbet2, cbet2, dn2,
01330 + salp1, calp1, f, A3x(*), C3x(*)
01331 logical diffp
01332 * output
01333 double precision salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
01334 + eps, domg12
01335 * optional output
01336 double precision dlam12
01337 * temporary
01338 double precision C1a(*), C2a(*), C3a(*)
01339
01340 integer ord, nC3
01341 parameter (ord = 6, nC3 = ord)
01342
01343 double precision csmgt, hypotx, A3f, TrgSum
01344
01345 double precision f1, e2, ep2, salp0, calp0,
01346 + somg1, comg1, somg2, comg2, omg12, lam12, B312, h0, k2, dummy
01347
01348 double precision dblmin, dbleps, pi, degree, tiny,
01349 + tol0, tol1, tol2, tolb, xthrsh
01350 integer digits, maxit1, maxit2
01351 logical init
01352 common /geocom/ dblmin, dbleps, pi, degree, tiny,
01353 + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init
01354
01355 f1 = 1 - f
01356 e2 = f * (2 - f)
01357 ep2 = e2 / (1 - e2)
01358 * Break degeneracy of equatorial line. This case has already been
01359 * handled.
01360 if (sbet1 .eq. 0 .and. calp1 .eq. 0) calp1 = -tiny
01361
01362 * sin(alp1) * cos(bet1) = sin(alp0)
01363 salp0 = salp1 * cbet1
01364 * calp0 > 0
01365 calp0 = hypotx(calp1, salp1 * sbet1)
01366
01367 * tan(bet1) = tan(sig1) * cos(alp1)
01368 * tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
01369 ssig1 = sbet1
01370 somg1 = salp0 * sbet1
01371 csig1 = calp1 * cbet1
01372 comg1 = csig1
01373 call Norm(ssig1, csig1)
01374 * Norm(somg1, comg1); -- don't need to normalize!
01375
01376 * Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful
01377 * about this case, since this can yield singularities in the Newton
01378 * iteration.
01379 * sin(alp2) * cos(bet2) = sin(alp0)
01380 salp2 = csmgt(salp0 / cbet2, salp1, cbet2 .ne. cbet1)
01381 * calp2 = sqrt(1 - sq(salp2))
01382 * = sqrt(sq(calp0) - sq(sbet2)) / cbet2
01383 * and subst for calp0 and rearrange to give (choose positive sqrt
01384 * to give alp2 in [0, pi/2]).
01385 calp2 = csmgt(sqrt((calp1 * cbet1)**2 +
01386 + csmgt((cbet2 - cbet1) * (cbet1 + cbet2),
01387 + (sbet1 - sbet2) * (sbet1 + sbet2),
01388 + cbet1 .lt. -sbet1)) / cbet2,
01389 + abs(calp1), cbet2 .ne. cbet1 .or. abs(sbet2) .ne. -sbet1)
01390 * tan(bet2) = tan(sig2) * cos(alp2)
01391 * tan(omg2) = sin(alp0) * tan(sig2).
01392 ssig2 = sbet2
01393 somg2 = salp0 * sbet2
01394 csig2 = calp2 * cbet2
01395 comg2 = csig2
01396 call Norm(ssig2, csig2)
01397 * Norm(somg2, comg2); -- don't need to normalize!
01398
01399 * sig12 = sig2 - sig1, limit to [0, pi]
01400 sig12 = atan2(max(csig1 * ssig2 - ssig1 * csig2, 0d0),
01401 + csig1 * csig2 + ssig1 * ssig2)
01402
01403 * omg12 = omg2 - omg1, limit to [0, pi]
01404 omg12 = atan2(max(comg1 * somg2 - somg1 * comg2, 0d0),
01405 + comg1 * comg2 + somg1 * somg2)
01406 k2 = calp0**2 * ep2
01407 eps = k2 / (2 * (1 + sqrt(1 + k2)) + k2)
01408 call C3f(eps, C3x, C3a)
01409 B312 = (TrgSum(.true., ssig2, csig2, C3a, nC3-1) -
01410 + TrgSum(.true., ssig1, csig1, C3a, nC3-1))
01411 h0 = -f * A3f(eps, A3x)
01412 domg12 = salp0 * h0 * (sig12 + B312)
01413 lam12 = omg12 + domg12
01414
01415 if (diffp) then
01416 if (calp2 .eq. 0) then
01417 dlam12 = - 2 * f1 * dn1 / sbet1
01418 else
01419 call Lengs(eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
01420 + cbet1, cbet2, dummy, dlam12, dummy,
01421 + .false., dummy, dummy, ep2, C1a, C2a)
01422 dlam12 = dlam12 * f1 / (calp2 * cbet2)
01423 end if
01424 end if
01425 Lam12f = lam12
01426
01427 return
01428 end
01429
01430 double precision function A3f(eps, A3x)
01431 * Evaluate sum(A3x[k] * eps^k, k, 0, nA3x-1) by Horner's method
01432 integer ord, nA3, nA3x
01433 parameter (ord = 6, nA3 = ord, nA3x = nA3)
01434
01435 * input
01436 double precision eps
01437 * output
01438 double precision A3x(0: nA3x-1)
01439
01440 integer i
01441 A3f = 0
01442 do 10 i = nA3x-1, 0, -1
01443 A3f = eps * A3f + A3x(i)
01444 10 continue
01445
01446 return
01447 end
01448
01449 subroutine C3f(eps, C3x, c)
01450 * Evaluate C3 coeffs by Horner's method
01451 * Elements c[1] thru c[nC3-1] are set
01452 integer ord, nC3, nC3x
01453 parameter (ord = 6, nC3 = ord, nC3x = (nC3 * (nC3 - 1)) / 2)
01454
01455 * input
01456 double precision eps, C3x(0:nC3x-1)
01457 * output
01458 double precision c(nC3-1)
01459
01460 integer i, j, k
01461 double precision t, mult
01462
01463 j = nC3x
01464 do 20 k = nC3-1, 1 , -1
01465 t = 0
01466 do 10 i = nC3 - k, 1, -1
01467 j = j - 1
01468 t = eps * t + C3x(j)
01469 10 continue
01470 c(k) = t
01471 20 continue
01472
01473 mult = 1
01474 do 30 k = 1, nC3-1
01475 mult = mult * eps
01476 c(k) = c(k) * mult
01477 30 continue
01478
01479 return
01480 end
01481
01482 subroutine C4f(eps, C4x, c)
01483 * Evaluate C4 coeffs by Horner's method
01484 * Elements c[0] thru c[nC4-1] are set
01485 integer ord, nC4, nC4x
01486 parameter (ord = 6, nC4 = ord, nC4x = (nC4 * (nC4 + 1)) / 2)
01487
01488 * input
01489 double precision eps, C4x(0:nC4x-1)
01490 *output
01491 double precision c(0:nC4-1)
01492
01493 integer i, j, k
01494 double precision t, mult
01495
01496 j = nC4x
01497 do 20 k = nC4-1, 0, -1
01498 t = 0
01499 do 10 i = nC4 - k, 1, -1
01500 j = j - 1
01501 t = eps * t + C4x(j)
01502 10 continue
01503 c(k) = t
01504 20 continue
01505
01506 mult = 1
01507 do 30 k = 1, nC4-1
01508 mult = mult * eps
01509 c(k) = c(k) * mult
01510 30 continue
01511
01512 return
01513 end
01514
01515 * Generated by Maxima on 2010-09-04 10:26:17-04:00
01516
01517 double precision function A1m1f(eps)
01518 * The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
01519 * input
01520 double precision eps
01521
01522 double precision eps2, t
01523
01524 eps2 = eps**2
01525 t = eps2*(eps2*(eps2+4)+64)/256
01526 A1m1f = (t + eps) / (1 - eps)
01527
01528 return
01529 end
01530
01531 subroutine C1f(eps, c)
01532 * The coefficients C1[l] in the Fourier expansion of B1
01533 integer ord, nC1
01534 parameter (ord = 6, nC1 = ord)
01535
01536 * input
01537 double precision eps
01538 * output
01539 double precision c(nC1)
01540
01541 double precision eps2, d
01542
01543 eps2 = eps**2
01544 d = eps
01545 c(1) = d*((6-eps2)*eps2-16)/32
01546 d = d * eps
01547 c(2) = d*((64-9*eps2)*eps2-128)/2048
01548 d = d * eps
01549 c(3) = d*(9*eps2-16)/768
01550 d = d * eps
01551 c(4) = d*(3*eps2-5)/512
01552 d = d * eps
01553 c(5) = -7*d/1280
01554 d = d * eps
01555 c(6) = -7*d/2048
01556
01557 return
01558 end
01559
01560 subroutine C1pf(eps, c)
01561 * The coefficients C1p[l] in the Fourier expansion of B1p
01562 integer ord, nC1p
01563 parameter (ord = 6, nC1p = ord)
01564
01565 * input
01566 double precision eps
01567 * output
01568 double precision c(nC1p)
01569
01570 double precision eps2, d
01571
01572 eps2 = eps**2
01573 d = eps
01574 c(1) = d*(eps2*(205*eps2-432)+768)/1536
01575 d = d * eps
01576 c(2) = d*(eps2*(4005*eps2-4736)+3840)/12288
01577 d = d * eps
01578 c(3) = d*(116-225*eps2)/384
01579 d = d * eps
01580 c(4) = d*(2695-7173*eps2)/7680
01581 d = d * eps
01582 c(5) = 3467*d/7680
01583 d = d * eps
01584 c(6) = 38081*d/61440
01585
01586 return
01587 end
01588
01589 * The scale factor A2-1 = mean value of (d/dsigma)I2 - 1
01590 double precision function A2m1f(eps)
01591 * input
01592 double precision eps
01593
01594 double precision eps2, t
01595
01596 eps2 = eps**2
01597 t = eps2*(eps2*(25*eps2+36)+64)/256
01598 A2m1f = t * (1 - eps) - eps
01599
01600 return
01601 end
01602
01603 subroutine C2f(eps, c)
01604 * The coefficients C2[l] in the Fourier expansion of B2
01605 integer ord, nC2
01606 parameter (ord = 6, nC2 = ord)
01607
01608 * input
01609 double precision eps
01610 * output
01611 double precision c(nC2)
01612
01613 double precision eps2, d
01614
01615 eps2 = eps**2
01616 d = eps
01617 c(1) = d*(eps2*(eps2+2)+16)/32
01618 d = d * eps
01619 c(2) = d*(eps2*(35*eps2+64)+384)/2048
01620 d = d * eps
01621 c(3) = d*(15*eps2+80)/768
01622 d = d * eps
01623 c(4) = d*(7*eps2+35)/512
01624 d = d * eps
01625 c(5) = 63*d/1280
01626 d = d * eps
01627 c(6) = 77*d/2048
01628
01629 return
01630 end
01631
01632 subroutine A3cof(n, A3x)
01633 * The scale factor A3 = mean value of (d/dsigma)I3
01634 integer ord, nA3, nA3x
01635 parameter (ord = 6, nA3 = ord, nA3x = nA3)
01636
01637 * input
01638 double precision n
01639 * output
01640 double precision A3x(0:nA3x-1)
01641
01642 A3x(0) = 1
01643 A3x(1) = (n-1)/2
01644 A3x(2) = (n*(3*n-1)-2)/8
01645 A3x(3) = ((-n-3)*n-1)/16
01646 A3x(4) = (-2*n-3)/64
01647 A3x(5) = -3/128d0
01648
01649 return
01650 end
01651
01652 subroutine C3cof(n, C3x)
01653 * The coefficients C3[l] in the Fourier expansion of B3
01654 integer ord, nC3, nC3x
01655 parameter (ord = 6, nC3 = ord, nC3x = (nC3 * (nC3 - 1)) / 2)
01656
01657 * input
01658 double precision n
01659 * output
01660 double precision C3x(0:nC3x-1)
01661
01662 C3x(0) = (1-n)/4
01663 C3x(1) = (1-n*n)/8
01664 C3x(2) = ((3-n)*n+3)/64
01665 C3x(3) = (2*n+5)/128
01666 C3x(4) = 3/128d0
01667 C3x(5) = ((n-3)*n+2)/32
01668 C3x(6) = ((-3*n-2)*n+3)/64
01669 C3x(7) = (n+3)/128
01670 C3x(8) = 5/256d0
01671 C3x(9) = (n*(5*n-9)+5)/192
01672 C3x(10) = (9-10*n)/384
01673 C3x(11) = 7/512d0
01674 C3x(12) = (7-14*n)/512
01675 C3x(13) = 7/512d0
01676 C3x(14) = 21/2560d0
01677
01678 return
01679 end
01680
01681 * Generated by Maxima on 2012-10-19 08:02:34-04:00
01682
01683 subroutine C4cof(n, C4x)
01684 * The coefficients C4[l] in the Fourier expansion of I4
01685 integer ord, nC4, nC4x
01686 parameter (ord = 6, nC4 = ord, nC4x = (nC4 * (nC4 + 1)) / 2)
01687
01688 * input
01689 double precision n
01690 * output
01691 double precision C4x(0:nC4x-1)
01692
01693 C4x(0) = (n*(n*(n*(n*(100*n+208)+572)+3432)-12012)+30030)/45045
01694 C4x(1) = (n*(n*(n*(64*n+624)-4576)+6864)-3003)/15015
01695 C4x(2) = (n*((14144-10656*n)*n-4576)-858)/45045
01696 C4x(3) = ((-224*n-4784)*n+1573)/45045
01697 C4x(4) = (1088*n+156)/45045
01698 C4x(5) = 97/15015d0
01699 C4x(6) = (n*(n*((-64*n-624)*n+4576)-6864)+3003)/135135
01700 C4x(7) = (n*(n*(5952*n-11648)+9152)-2574)/135135
01701 C4x(8) = (n*(5792*n+1040)-1287)/135135
01702 C4x(9) = (468-2944*n)/135135
01703 C4x(10) = 1/9009d0
01704 C4x(11) = (n*((4160-1440*n)*n-4576)+1716)/225225
01705 C4x(12) = ((4992-8448*n)*n-1144)/225225
01706 C4x(13) = (1856*n-936)/225225
01707 C4x(14) = 8/10725d0
01708 C4x(15) = (n*(3584*n-3328)+1144)/315315
01709 C4x(16) = (1024*n-208)/105105
01710 C4x(17) = -136/63063d0
01711 C4x(18) = (832-2560*n)/405405
01712 C4x(19) = -128/135135d0
01713 C4x(20) = 128/99099d0
01714
01715 return
01716 end
01717
01718 double precision function sumx(u, v, t)
01719 * input
01720 double precision u, v
01721 * output
01722 double precision t
01723
01724 double precision up, vpp
01725 sumx = u + v
01726 up = sumx - v
01727 vpp = sumx - up
01728 up = up - u
01729 vpp = vpp - v
01730 t = -(up + vpp)
01731
01732 return
01733 end
01734
01735 double precision function AngNm(x)
01736 * input
01737 double precision x
01738
01739 if (x .ge. 180) then
01740 AngNm = x - 360
01741 else if (x .lt. -180) then
01742 AngNm = x + 360
01743 else
01744 AngNm = x
01745 end if
01746
01747 return
01748 end
01749
01750 double precision function AngNm2(x)
01751 * input
01752 double precision x
01753
01754 double precision AngNm
01755 AngNm2 = mod(x, 360d0)
01756 AngNm2 = AngNm(AngNm2)
01757
01758 return
01759 end
01760
01761 double precision function AngDif(x, y)
01762 * Compute y - x. x and y must both lie in [-180, 180]. The result is
01763 * equivalent to computing the difference exactly, reducing it to (-180,
01764 * 180] and rounding the result. Note that this prescription allows -180
01765 * to be returned (e.g., if x is tiny and negative and y = 180).
01766 * input
01767 double precision x, y
01768
01769 double precision d, t, sumx
01770 d = sumx(-x, y, t)
01771 if ((d - 180d0) + t .gt. 0d0) then
01772 d = d - 360d0
01773 else if ((d + 180d0) + t .le. 0d0) then
01774 d = d + 360d0
01775 end if
01776 AngDif = d + t
01777
01778 return
01779 end
01780
01781 double precision function AngRnd(x)
01782 * The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
01783 * for reals = 0.7 pm on the earth if x is an angle in degrees. (This
01784 * is about 1000 times more resolution than we get with angles around 90
01785 * degrees.) We use this to avoid having to deal with near singular
01786 * cases when x is non-zero but tiny (e.g., 1.0e-200).
01787 * input
01788 double precision x
01789
01790 double precision y, z
01791 z = 1/16d0
01792 y = abs(x)
01793 * The compiler mustn't "simplify" z - (z - y) to y
01794 if (y .lt. z) y = z - (z - y)
01795 AngRnd = sign(y, x)
01796
01797 return
01798 end
01799
01800 subroutine swap(x, y)
01801 * input/output
01802 double precision x, y
01803
01804 double precision z
01805 z = x
01806 x = y
01807 y = z
01808
01809 return
01810 end
01811
01812 double precision function hypotx(x, y)
01813 * input
01814 double precision x, y
01815
01816 hypotx = sqrt(x**2 + y**2)
01817
01818 return
01819 end
01820
01821 subroutine Norm(sinx, cosx)
01822 * input/output
01823 double precision sinx, cosx
01824
01825 double precision hypotx, r
01826 r = hypotx(sinx, cosx)
01827 sinx = sinx/r
01828 cosx = cosx/r
01829
01830 return
01831 end
01832
01833 double precision function log1px(x)
01834 * input
01835 double precision x
01836
01837 double precision csmgt, y, z
01838 y = 1 + x
01839 z = y - 1
01840 log1px = csmgt(x, x * log(y) / z, z .eq. 0)
01841
01842 return
01843 end
01844
01845 double precision function atanhx(x)
01846 * input
01847 double precision x
01848
01849 double precision log1px, y
01850 y = abs(x)
01851 y = log1px(2 * y/(1 - y))/2
01852 atanhx = sign(y, x)
01853
01854 return
01855 end
01856
01857 double precision function cbrt(x)
01858 * input
01859 double precision x
01860
01861 cbrt = sign(abs(x)**(1/3d0), x)
01862
01863 return
01864 end
01865
01866 double precision function csmgt(x, y, p)
01867 * input
01868 double precision x, y
01869 logical p
01870
01871 if (p) then
01872 csmgt = x
01873 else
01874 csmgt = y
01875 end if
01876
01877 return
01878 end
01879
01880 double precision function TrgSum(sinp, sinx, cosx, c, n)
01881 * Evaluate
01882 * y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
01883 * sum(c[i] * cos((2*i-1) * x), i, 1, n)
01884 * using Clenshaw summation.
01885 * Approx operation count = (n + 5) mult and (2 * n + 2) add
01886 * input
01887 logical sinp
01888 integer n
01889 double precision sinx, cosx, c(n)
01890
01891 double precision ar, y0, y1
01892 integer n2, k
01893
01894 * 2 * cos(2 * x)
01895 ar = 2 * (cosx - sinx) * (cosx + sinx)
01896 * accumulators for sum
01897 if (mod(n, 2) .eq. 1) then
01898 y0 = c(n)
01899 n2 = n - 1
01900 else
01901 y0 = 0
01902 n2 = n
01903 end if
01904 y1 = 0
01905 * Now n2 is even
01906 do 10 k = n2, 1, -2
01907 * Unroll loop x 2, so accumulators return to their original role
01908 y1 = ar * y0 - y1 + c(k)
01909 y0 = ar * y1 - y0 + c(k-1)
01910 10 continue
01911 if (sinp) then
01912 * sin(2 * x) * y0
01913 TrgSum = 2 * sinx * cosx * y0
01914 else
01915 * cos(x) * (y0 - y1)
01916 TrgSum = cosx * (y0 - y1)
01917 end if
01918
01919 return
01920 end
01921
01922 integer function trnsit(lon1, lon2)
01923 * input
01924 double precision lon1, lon2
01925
01926 double precision lon1x, lon2x, lon12, AngNm, AngDif
01927 lon1x = AngNm(lon1)
01928 lon2x = AngNm(lon2)
01929 lon12 = AngDif(lon1x, lon2x)
01930 trnsit = 0
01931 if (lon1x .lt. 0 .and. lon2x .ge. 0 .and. lon12 .gt. 0) then
01932 trnsit = 1
01933 else if (lon2x .lt. 0 .and. lon1x .ge. 0 .and. lon12 .lt. 0) then
01934 trnsit = -1
01935 end if
01936
01937 return
01938 end
01939
01940 subroutine accini(s)
01941 * Initialize an accumulator; this is an array with two elements.
01942 * input/output
01943 double precision s(2)
01944
01945 s(1) = 0
01946 s(2) = 0
01947
01948 return
01949 end
01950
01951 subroutine accadd(s, y)
01952 * Add y to an accumulator.
01953 * input
01954 double precision y
01955 * input/output
01956 double precision s(2)
01957
01958 double precision z, u, sumx
01959 z = sumx(y, s(2), u)
01960 s(1) = sumx(z, s(1), s(2))
01961 if (s(1) .eq. 0) then
01962 s(1) = u
01963 else
01964 s(2) = s(2) + u
01965 end if
01966
01967 return
01968 end
01969
01970 * Table of name abbreviations to conform to the 6-char limit
01971 * A3coeff A3cof
01972 * C3coeff C3cof
01973 * C4coeff C4cof
01974 * AngNormalize AngNm
01975 * AngNormalize2 AngNm2
01976 * AngDiff AngDif
01977 * AngRound AngRnd
01978 * arcmode arcmod
01979 * Astroid Astrd
01980 * betscale betscl
01981 * lamscale lamscl
01982 * cbet12a cbt12a
01983 * sbet12a sbt12a
01984 * epsilon dbleps
01985 * realmin dblmin
01986 * geodesic geod
01987 * inverse invers
01988 * InverseStart InvSta
01989 * Lambda12 Lam12f
01990 * latsign latsgn
01991 * lonsign lonsgn
01992 * Lengths Lengs
01993 * meridian merid
01994 * outmask omask
01995 * shortline shortp
01996 * SinCosNorm Norm
01997 * SinCosSeries TrgSum
01998 * xthresh xthrsh
01999 * transit trnsit
02000 *> @endcond SKIP