geodesic.for File Reference

Go to the source code of this file.

Functions/Subroutines

*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid 
at (\e lat1,\e *!lon1) and(\e lat2
*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid e lon2 is called the
geodesic Its length is *e s12
and the geodesic from point to
point has forward azimuths *e
azi1 and e azi2 at the two end
points **!Traditionally two
geodesic problems are e e and
e *!determine e e and e azi2
This is solved by the
*!subroutine 
direct ().*!-the inverse problem--given\e lat1
*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid e lon2 is called the
geodesic Its length is *e s12
and the geodesic from point to
point has forward azimuths *e
azi1 and e azi2 at the two end
points **!Traditionally two
geodesic problems are e e and
e *!determine e e and e azi2
This is solved by the
*!subroutine e e e *!determine
e e and e azi2 This is solved
by the *!subroutine 
invers ().*!*!The ellipsoid is specified by its equatorial radius\e a(typically *!in meters) and flattening\e f.The routines are accurate to round *!off with double precision arithmetic provided that|< i >f</i >|&lt
*for the WGS84 the errors are
less than * 
!nanometers (Reasonably accurate results are obtained for|< i >f</i >|*!&lt;1/5.) For a prolate ellipsoid
i it is the measured counter
*!of the geodesic
quadrilateral with 
corners (\e lat1,\e lon1)
i it is the measured counter
*!of the geodesic
quadrilateral with e e 
and (\e lat2,\e lon2).*!-\e m12
i it is the measured counter
*!of the geodesic
quadrilateral with e e the
reduced length of the geodesic
is defined such that if *!the
initial azimuth is perturbed
by e 
dazi1 (radians) then the *!second point is displaced by\e m12\e dazi1 in the direction *!perpendicular to the geodesic.On a curved surface the reduced *!length obeys a symmetry relation
i it is the measured counter
*!of the geodesic
quadrilateral with e e the
reduced length of the geodesic
is defined such that if *!the
initial azimuth is perturbed
by e e m12 e we have e then
*!they are separated by a
distance e MM12 e dt at point
e MM21 *!is defined 
similarly (with the geodesics being parallel to one *!another at point 2).On a flat surface
e MM23 e MM32 e *!m12 e m23
**!The shortest distance
returned by the solution of
the inverse problem * 
!is (obviously) uniquely defined.However
e lat2 (with neither point at a pole).If\e azi1
e SS12 (This occurs when the longitude difference is near *!&plusmn;180 &deg;for oblate ellipsoids.)*!-\e lon2
e SS12 (This occurs when\e lat2 is near *!&minus;\e lat1 for prolate ellipsoids.)*!-Points 1 and 2 at opposite poles.There are infinitely many *!geodesics which can be generated by setting[\e azi1
e for arbitrary e d (For *!spheres, this prescription applies when points 1 and 2 are *!antipodal.)*!-\e s12=0(coincident points).There are infinitely many *!geodesics which can be generated by setting[\e azi1
e for arbitrary e e for
arbitrary e d **!These
routines are a simple
transcription of the
corresponding C *!classes in
< a href="http:*! Because of
the limitations of Fortran
77, the classes have been*!
replaced by simple subroutines
with no attempt to save "state"
across*! subroutine calls.
Most of the internal comments
have been retained.*! However,
in the process of
transcription some
documentation has been*! lost
and the documentation for the
C++ classes,*!
GeographicLib::Geodesic,
GeographicLib::GeodesicLine,
and*!
GeographicLib::PolygonArea,
should be consulted. The C++
code*! remains the "reference
implementation". Think twice
about*! restructuring the
internals of the Fortran code
since this may make*! porting
fixes from the C++ code more
difficult.*!*! Copyright (c)
Charles Karney (2012-2013)
<charles@karney.com> and*!
licensed under the MIT/X11
License. For more information,
see*! http:*!*! This library
was distributed with*! <a href="../index.html">
GeographicLib</a>
1.31.*> Solve the direct
geodesic problem* * param [in]
a the equatorial 
radius (meters).*!@param[in] f the flattening of the ellipsoid.Setting\e f=0gives *!a sphere.Negative\e f gives a prolate ellipsoid.*!@param[in] lat1 latitude of point 1(degrees).*!@param[in] lon1 longitude of point 1(degrees).*!@param[in] azi1 azimuth at point 1(degrees).*!@param[in] s12a12 if\e arcmod is false
e for arbitrary e e for
arbitrary e d **!These
routines are a simple
transcription of the
corresponding C *!classes in
< a href="http:*! Because of
the limitations of Fortran
77, the classes have been*!
replaced by simple subroutines
with no attempt to save "state"
across*! subroutine calls.
Most of the internal comments
have been retained.*! However,
in the process of
transcription some
documentation has been*! lost
and the documentation for the
C++ classes,*!
GeographicLib::Geodesic,
GeographicLib::GeodesicLine,
and*!
GeographicLib::PolygonArea,
should be consulted. The C++
code*! remains the "reference
implementation". Think twice
about*! restructuring the
internals of the Fortran code
since this may make*! porting
fixes from the C++ code more
difficult.*!*! Copyright (c)
Charles Karney (2012-2013)
<charles@karney.com> and*!
licensed under the MIT/X11
License. For more information,
see*! http:*!*! This library
was distributed with*! <a href="../index.html">
GeographicLib</a>
1.31.*> Solve the direct
geodesic problem* * param [in]
a the equatorial this is the
distance between* ! point and 
point (meters)
otherwise it is the arc length
*!between point and 
point (degrees)

Variables

*The subroutines in this files
are documented at * 
http
*The subroutines in this files
are documented at *< ahref="http:*!Algorithmsforgeodesics
</a> * ! J Geodesy<b></b> * 
!DOI
!addenda
 Vincenty
are *accurate to round off for
< i > f</i > & 
lt
*the solution of the inverse
problem is always 
found
*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid e lon2 is called the
geodesic Its length is *e s12
and the geodesic from point to
point has forward azimuths *e
azi1 and e azi2 at the two end
points **!Traditionally two
geodesic problems are 
considered
*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid e lon2 is called the
geodesic Its length is *e s12
and the geodesic from point to
point has forward azimuths *e
azi1 and e azi2 at the two end
points **!Traditionally two
geodesic problems are e 
lon1
*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid e lon2 is called the
geodesic Its length is *e s12
and the geodesic from point to
point has forward azimuths *e
azi1 and e azi2 at the two end
points **!Traditionally two
geodesic problems are e e 
s12
*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid e lon2 is called the
geodesic Its length is *e s12
and the geodesic from point to
point has forward azimuths *e
azi1 and e azi2 at the two end
points **!Traditionally two
geodesic problems are e e and
e 
azi1 = 0&deg
*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid e lon2 is called the
geodesic Its length is *e s12
and the geodesic from point to
point has forward azimuths *e
azi1 and e azi2 at the two end
points **!Traditionally two
geodesic problems are e e and
e *!determine e 
lat2
*differential and integral
properties of geodesics are
computed **!The shortest path
between two points on the
ellipsoid e lon2 is called the
geodesic Its length is *e s12
and the geodesic from point to
point has forward azimuths *e
azi1 and e azi2 at the two end
points **!Traditionally two
geodesic problems are e e and
e *!determine e e 
lon2
*for the WGS84 ellipsoid
**!The routines also calculate
several other quantities of
interest *e SS12 is the area
between the geodesic from
point to point *!and the 
equator
e
i it is the area
i it is the measured counter clockwise
i it is the measured counter
*!of the geodesic
quadrilateral with e
!lon1
i it is the measured counter
*!of the geodesic
quadrilateral with e e the
reduced length of the geodesic
is defined such that if *!the
initial azimuth is perturbed
by e e m12 e 
m21
i it is the measured counter
*!of the geodesic
quadrilateral with e e the
reduced length of the geodesic
is defined such that if *!the
initial azimuth is perturbed
by e e m12 e we have e 
m12
i it is the measured counter
*!of the geodesic
quadrilateral with e e the
reduced length of the geodesic
is defined such that if *!the
initial azimuth is perturbed
by e e m12 e we have e then
*!they are separated by a
distance e MM12 e dt at point
e MM21 *!is defined we have e 
MM12
**!If points
**!If and lie on a single geodesic
**!If and lie on a single then
the following *!addition rules 
hold
minus
e MM12 e MM21 e *!m23 e m12 *e MM31 = \e MM32 \e MM21 &minus
e MM23 e MM32 e *!m12 e m23
**!The shortest distance
returned by the solution of
the inverse problem *in a few
special cases *!there are
multiple azimuths which yield
the same shortest distance
*!Here is a catalog of those 
cases
e the geodesic is unique
Otherwise there are two
geodesics *!and the second one
is obtained by 
setting [\e azi1,\e azi2]
e the geodesic is unique
Otherwise there are two
geodesics *!and the second one
is obtained by e 
SS12
deg
or & plusmn
the geodesic is unique
*!Otherwise there are two
geodesics and the second one
is obtained by * 
!setting [\e azi1,\e azi2] = [&minus
e azi2
e e!SS12 = &minus
e d
otherwise it is the distance
*!between point and whose
binary bits are interpreted
*!as follows *return e a12
*return e m12 *return e MM12
and e MM21 *return e SS12 **e
lat1 should be in the 
range [&minus;90 &deg;, 90 &deg;]
e lon1 and *e azi1 should be
in the range[&minus;540 &deg;,
540 &deg;).The *!values of\e
lon2 and\e azi2 returned are
in the range *![&minus;180
&deg;, 180 &deg;).*!*!If
either point is at a pole, the
azimuth is defined by keeping
the *!longitude fixed, writing\e
lat=\e lat=&plusmn;(90 &deg;&minus;*!&epsilon;),
and taking the limit &epsilon;&rarr;0+.An
arc length *!greater that
180 &deg;signifies a geodesic
which is not a shortest
*!path.(For a prolate
ellipsoid, an additional
condition is necessary *!for a
shortest path:the longitudinal
extent must not exceed of
*!180 &deg;.) subroutine
direct(a, f, lat1, lon1, azi1,
s12a12, arcmod,+lat2, lon2,
azi2, omask, a12s12, m12, MM12,
MM21, SS12)*input double
precision a, f, lat1, lon1,
azi1, s12a12 logical arcmod
integer omask *output double
precision lat2, lon2, azi2
*optional output double
precision a12s12, m12, MM12,
MM21, SS12 integer ord, nC1,
nC1p, nC2, nA3, nA3x, nC3,
nC3x, nC4, nC4x parameter(ord=6,
nC1=ord, nC1p=ord,+nC2=ord,
nA3=ord, nA3x=nA3,+nC3=ord,
nC3x=(nC3 *(nC3-1))/2,+nC4=ord,
nC4x=(nC4 *(nC4+1))/2) double
precision A3x(0:nA3x-1), C3x(0:nC3x-1),
C4x(0:nC4x-1),+C1a(nC1), C1pa(nC1p),
C2a(nC2), C3a(nC3-1), C4a(0:nC4-1)
double precision csmgt, atanhx,
hypotx,+AngNm, AngNm2, AngRnd,
TrgSum, A1m1f, A2m1f, A3f
logical arcp, redlp, scalp,
areap double precision e2, f1,
ep2, n, b, c2,+lon1x, azi1x,
phi, alp1, salp0, calp0, k2,
eps,+salp1, calp1, ssig1,
csig1, cbet1, sbet1, dn1,
somg1, comg1,+salp2, calp2,
ssig2, csig2, sbet2, cbet2,
dn2, somg2, comg2,+ssig12,
csig12, salp12, calp12, omg12,
lam12, lon12,+sig12, stau1,
ctau1, tau12, s12a, t, s, c,
serr,+A1m1, A2m1, A3c, A4, AB1,
AB2,+B11, B12, B21, B22, B31,
B41, B42, J12 double precision
dblmin, dbleps, pi, degree,
tiny,+tol0, tol1, tol2, tolb,
xthrsh integer digits, maxit1,
maxit2 logical init common/geocom/dblmin,
dbleps, pi, degree, tiny,+tol0,
tol1, tol2, tolb, xthrsh,
digits, maxit1, maxit2, init
if(.not.init) call geoini e2=f
*(2-f) ep2=e2/(1-e2) f1=1-f n=f/(2-f)
b=a *f1 c2=0 arcp=mod(omask/1,
2)==1 redlp=mod(omask/2,
2)==1 scalp=mod(omask/4,
2)==1 areap=mod(omask/8,
2)==1 if(areap) then if(e2.eq.0)
then c2=a **2 else if(e2.gt.0)
then c2=(a **2+b **2 *atanhx(sqrt(e2))/sqrt(e2))/2
else c2=(a **2+b **2 *atan(sqrt(abs(e2)))/sqrt(abs(e2)))/2
end if end if call A3cof(n,
A3x) call C3cof(n, C3x) if(areap)
call C4cof(n, C4x)*Guard
against underflow in salp0
azi1x=AngRnd(AngNm(azi1))
lon1x=AngNm(lon1)*alp1 is in 
!k2 [0, pi] alp1=azi1x *degree *Enforce sin(pi)==0 and cos(pi/2)==0.Better to face the ensuing *problems directly than to skirt them.salp1=csmgt(0d0, sin(alp1), azi1x.eq.-180) calp1=csmgt(0d0, cos(alp1), abs(azi1x).eq.90) phi=lat1 *degree *Ensure cbet1=+dbleps at poles sbet1=f1 *sin(phi) cbet1=csmgt(tiny, cos(phi), abs(lat1).eq.90) call Norm(sbet1, cbet1) dn1=sqrt(1+ep2 *sbet1 **2)*Evaluate alp0 from sin(alp1)*cos(bet1)=sin(alp0),*alp0 in[0, pi/2-|bet1|] salp0=salp1 *cbet1 *Alt:calp0=hypot(sbet1, calp1 *cbet1).The following *is slightly better(consider the case salp1=0).calp0=hypotx(calp1, salp1 *sbet1)*Evaluate sig with tan(bet1)=tan(sig1)*cos(alp1).*sig=0 is nearest northward crossing of equator.*With bet1=0, alp1=pi/2, we have sig1=0(equatorial line).*With bet1=pi/2, alp1=-pi, sig1=pi/2 *With bet1=-pi/2, alp1=0, sig1=-pi/2 *Evaluate omg1 with tan(omg1)=sin(alp0)*tan(sig1).*With alp0 in(0, pi/2], quadrants for sig and omg coincide.*No atan2(0, 0) ambiguity at poles since cbet1=+dbleps.*With alp0=0, omg1=0 for alp1=0, omg1=pi for alp1=pi.ssig1=sbet1 somg1=salp0 *sbet1 csig1=csmgt(cbet1 *calp1, 1d0, sbet1.ne.0.or.calp1.ne.0) comg1=csig1 *sig1 in(-pi, pi] call Norm(ssig1, csig1)*Geodesic don t need to normalize
e lon1 and *e azi1 should be
in the range[&minus;540 &deg;,
540 &deg;).The *!values of\e
lon2 and\e azi2 returned are
in the range *![&minus;180
&deg;, 180 &deg;).*!*!If
either point is at a pole, the
azimuth is defined by keeping
the *!longitude fixed, writing\e
lat=\e lat=&plusmn;(90 &deg;&minus;*!&epsilon;),
and taking the limit &epsilon;&rarr;0+.An
arc length *!greater that
180 &deg;signifies a geodesic
which is not a shortest
*!path.(For a prolate
ellipsoid, an additional
condition is necessary *!for a
shortest path:the longitudinal
extent must not exceed of
*!180 &deg;.) subroutine
direct(a, f, lat1, lon1, azi1,
s12a12, arcmod,+lat2, lon2,
azi2, omask, a12s12, m12, MM12,
MM21, SS12)*input double
precision a, f, lat1, lon1,
azi1, s12a12 logical arcmod
integer omask *output double
precision lat2, lon2, azi2
*optional output double
precision a12s12, m12, MM12,
MM21, SS12 integer ord, nC1,
nC1p, nC2, nA3, nA3x, nC3,
nC3x, nC4, nC4x parameter(ord=6,
nC1=ord, nC1p=ord,+nC2=ord,
nA3=ord, nA3x=nA3,+nC3=ord,
nC3x=(nC3 *(nC3-1))/2,+nC4=ord,
nC4x=(nC4 *(nC4+1))/2) double
precision A3x(0:nA3x-1), C3x(0:nC3x-1),
C4x(0:nC4x-1),+C1a(nC1), C1pa(nC1p),
C2a(nC2), C3a(nC3-1), C4a(0:nC4-1)
double precision csmgt, atanhx,
hypotx,+AngNm, AngNm2, AngRnd,
TrgSum, A1m1f, A2m1f, A3f
logical arcp, redlp, scalp,
areap double precision e2, f1,
ep2, n, b, c2,+lon1x, azi1x,
phi, alp1, salp0, calp0, k2,
eps,+salp1, calp1, ssig1,
csig1, cbet1, sbet1, dn1,
somg1, comg1,+salp2, calp2,
ssig2, csig2, sbet2, cbet2,
dn2, somg2, comg2,+ssig12,
csig12, salp12, calp12, omg12,
lam12, lon12,+sig12, stau1,
ctau1, tau12, s12a, t, s, c,
serr,+A1m1, A2m1, A3c, A4, AB1,
AB2,+B11, B12, B21, B22, B31,
B41, B42, J12 double precision
dblmin, dbleps, pi, degree,
tiny,+tol0, tol1, tol2, tolb,
xthrsh integer digits, maxit1,
maxit2 logical init common/geocom/dblmin,
dbleps, pi, degree, tiny,+tol0,
tol1, tol2, tolb, xthrsh,
digits, maxit1, maxit2, init
if(.not.init) call geoini e2=f
*(2-f) ep2=e2/(1-e2) f1=1-f n=f/(2-f)
b=a *f1 c2=0 arcp=mod(omask/1,
2)==1 redlp=mod(omask/2,
2)==1 scalp=mod(omask/4,
2)==1 areap=mod(omask/8,
2)==1 if(areap) then if(e2.eq.0)
then c2=a **2 else if(e2.gt.0)
then c2=(a **2+b **2 *atanhx(sqrt(e2))/sqrt(e2))/2
else c2=(a **2+b **2 *atan(sqrt(abs(e2)))/sqrt(abs(e2)))/2
end if end if call A3cof(n,
A3x) call C3cof(n, C3x) if(areap)
call C4cof(n, C4x)*Guard
against underflow in salp0
azi1x=AngRnd(AngNm(azi1))
lon1x=AngNm(lon1)*alp1 is in 
salp12 [0, pi] alp1=azi1x *degree *Enforce sin(pi)==0 and cos(pi/2)==0.Better to face the ensuing *problems directly than to skirt them.salp1=csmgt(0d0, sin(alp1), azi1x.eq.-180) calp1=csmgt(0d0, cos(alp1), abs(azi1x).eq.90) phi=lat1 *degree *Ensure cbet1=+dbleps at poles sbet1=f1 *sin(phi) cbet1=csmgt(tiny, cos(phi), abs(lat1).eq.90) call Norm(sbet1, cbet1) dn1=sqrt(1+ep2 *sbet1 **2)*Evaluate alp0 from sin(alp1)*cos(bet1)=sin(alp0),*alp0 in[0, pi/2-|bet1|] salp0=salp1 *cbet1 *Alt:calp0=hypot(sbet1, calp1 *cbet1).The following *is slightly better(consider the case salp1=0).calp0=hypotx(calp1, salp1 *sbet1)*Evaluate sig with tan(bet1)=tan(sig1)*cos(alp1).*sig=0 is nearest northward crossing of equator.*With bet1=0, alp1=pi/2, we have sig1=0(equatorial line).*With bet1=pi/2, alp1=-pi, sig1=pi/2 *With bet1=-pi/2, alp1=0, sig1=-pi/2 *Evaluate omg1 with tan(omg1)=sin(alp0)*tan(sig1).*With alp0 in(0, pi/2], quadrants for sig and omg coincide.*No atan2(0, 0) ambiguity at poles since cbet1=+dbleps.*With alp0=0, omg1=0 for alp1=0, omg1=pi for alp1=pi.ssig1=sbet1 somg1=salp0 *sbet1 csig1=csmgt(cbet1 *calp1, 1d0, sbet1.ne.0.or.calp1.ne.0) comg1=csig1 *sig1 in(-pi, pi] call Norm(ssig1, csig1)*Geodesic don t need to normalize converts to used in atan2 so no need to normalized

Function Documentation

* differential and integral properties of geodesics are computed* * ! The shortest path between two points on the ellipsoid at ( \e  lat1,
\e *!  lon1 
)
* differential and integral properties of geodesics are computed* * ! The shortest path between two points on the ellipsoid e lon2 is called the geodesic Its length is* e s12 and the geodesic from point to point has forward azimuths* e azi1 and e azi2 at the two end points* * ! Traditionally two geodesic problems are e e and e * ! determine e e and e azi2 This is solved by the* ! subroutine direct (  ) 
Type Constraints
* differential and integral properties of geodesics are computed* * ! The shortest path between two points on the ellipsoid e lon2 is called the geodesic Its length is* e s12 and the geodesic from point to point has forward azimuths* e azi1 and e azi2 at the two end points* * ! Traditionally two geodesic problems are e e and e * ! determine e e and e azi2 This is solved by the* ! subroutine e e e * ! determine e e and e azi2 This is solved by the* ! subroutine invers (  ) 
Type Constraints
* for the WGS84 the errors are less than* !nanometers ( Reasonably accurate results are obtained for|< i >f</i >|*!&lt;1/  5.  ) 
Type Constraints
i it is the measured counter * ! of the geodesic quadrilateral with corners ( \e  lat1,
\e  lon1 
)
Type Constraints
i it is the measured counter * ! of the geodesic quadrilateral with e e and ( \e  lat2,
\e  lon2 
)
Type Constraints
i it is the measured counter * ! of the geodesic quadrilateral with e e the reduced length of the geodesic is defined such that if* ! the initial azimuth is perturbed by e dazi1 ( radians   ) 
i it is the measured counter * ! of the geodesic quadrilateral with e e the reduced length of the geodesic is defined such that if* ! the initial azimuth is perturbed by e e m12 e we have e then* ! they are separated by a distance e MM12 e dt at point e MM21* ! is defined similarly ( with the geodesics being parallel to one *!another at point  2  ) 
Type Constraints
e MM23 e MM32 e* ! m12 e m23* * ! The shortest distance returned by the solution of the inverse problem* !is ( obviously   ) 
e lat2 ( with neither point at a  pole  ) 
e SS12 ( This occurs when the longitude difference is near *!&plusmn;180 &deg;for oblate  ellipsoids.  ) 
e SS12 ( This occurs when\e lat2 is near *!&minus;\e lat1 for prolate  ellipsoids.  ) 
e for arbitrary e d ( For *!  spheres,
this prescription applies when points 1 and 2 are *!  antipodal. 
) [pure virtual]
Type Constraints
e for arbitrary e e for arbitrary e d* * ! These routines are a simple transcription of the corresponding C* ! classes in<a href="http:*! Because of the limitations of Fortran 77, the classes have been*! replaced by simple subroutines with no attempt to save "state" across*! subroutine calls. Most of the internal comments have been retained.*! However, in the process of transcription some documentation has been*! lost and the documentation for the C++ classes,*! GeographicLib::Geodesic, GeographicLib::GeodesicLine, and*! GeographicLib::PolygonArea, should be consulted. The C++ code*! remains the "reference implementation". Think twice about*! restructuring the internals of the Fortran code since this may make*! porting fixes from the C++ code more difficult.*!*! Copyright (c) Charles Karney (2012-2013) <charles@karney.com> and*! licensed under the MIT/X11 License. For more information, see*! http:*!*! This library was distributed with*! <a href="../index.html">GeographicLib</a> 1.31.*> Solve the direct geodesic problem* * param [in] a the equatorial radius ( meters   )  [pure virtual]
Type Constraints
otherwise it is the distance *between point and point ( meters   ) 
otherwise it is the arc length* ! between point and point ( degrees   ) 

Variable Documentation

* The subroutines in this files are documented at* http

Definition at line 9 of file geodesic.for.

* The subroutines in this files are documented at *<ahref="http:*!Algorithmsforgeodesics</a> * ! J Geodesy<b></b> * !DOI

Definition at line 9 of file geodesic.for.

* !addenda

Definition at line 19 of file geodesic.for.

Definition at line 19 of file geodesic.for.

*for the WGS84 the errors are less than *specify e f & lt

Definition at line 19 of file geodesic.for.

* the solution of the inverse problem is always found

Definition at line 20 of file geodesic.for.

* differential and integral properties of geodesics are computed* * ! The shortest path between two points on the ellipsoid e lon2 is called the geodesic Its length is* e s12 and the geodesic from point to point has forward azimuths* e azi1 and e azi2 at the two end points* * ! Traditionally two geodesic problems are considered

Definition at line 30 of file geodesic.for.

* differential and integral properties of geodesics are computed* * ! The shortest path between two points on the ellipsoid e lon2 is called the geodesic Its length is* e s12 and the geodesic from point to point has forward azimuths* e azi1 and e azi2 at the two end points* * ! Traditionally two geodesic problems are e e and e * ! determine e e and e azi2 This is solved by the* ! subroutine e lon1

Definition at line 30 of file geodesic.for.

* differential and integral properties of geodesics are computed* * ! The shortest path between two points on the ellipsoid e lon2 is called the geodesic Its length is* e s12 and the geodesic from point to point has forward azimuths* e azi1 and e azi2 at the two end points* * ! Traditionally two geodesic problems are e e and e * ! determine e e and e azi2 This is solved by the* ! subroutine e e e * ! determine e s12

Definition at line 30 of file geodesic.for.

e azi1 = 0&deg

Definition at line 30 of file geodesic.for.

* differential and integral properties of geodesics are computed* * ! The shortest path between two points on the ellipsoid e lon2 is called the geodesic Its length is* e s12 and the geodesic from point to point has forward azimuths* e azi1 and e azi2 at the two end points* * ! Traditionally two geodesic problems are e e and e * ! determine e e and e azi2 This is solved by the* ! subroutine e e lat2

Definition at line 30 of file geodesic.for.

i it is the measured counter * ! of the geodesic quadrilateral with e e lon2

Definition at line 30 of file geodesic.for.

* for the WGS84 ellipsoid

Definition at line 40 of file geodesic.for.

* * ! The routines also calculate several other quantities of interest* e SS12 is the area between the geodesic from point to point* ! and the equator

Definition at line 46 of file geodesic.for.

i e

Definition at line 46 of file geodesic.for.

i it is the area

Definition at line 46 of file geodesic.for.

i it is the measured counter clockwise

Definition at line 46 of file geodesic.for.

i it is the measured counter * ! of the geodesic quadrilateral with e* !lon1

Definition at line 47 of file geodesic.for.

i it is the measured counter * ! of the geodesic quadrilateral with e e the reduced length of the geodesic is defined such that if* ! the initial azimuth is perturbed by e e m12 e m21
Initial value:
 0.  On a flat
*!   surface

Definition at line 53 of file geodesic.for.

i it is the measured counter * ! of the geodesic quadrilateral with e e the reduced length of the geodesic is defined such that if* ! the initial azimuth is perturbed by e e m12 e we have e m12
Initial value:
 \e s12.
*! - \e MM12 and \e MM21 are geodesic scales.  If two geodesics are
*!   parallel at point 1 and separated by a small distance \e dt

Definition at line 54 of file geodesic.for.

i it is the measured counter * ! of the geodesic quadrilateral with e e the reduced length of the geodesic is defined such that if* ! the initial azimuth is perturbed by e e m12 e we have e then* ! they are separated by a distance e MM12 e dt at point e MM21* ! is defined we have e MM12
Initial value:
 \e MM21
*!   = 1.
*! - \e a12 is the arc length on the auxiliary sphere.  This is a
*!   construct for converting the problem to one in spherical
*!   trigonometry.  \e a12 is measured in degrees.  The spherical arc
*!   length from one equator crossing to the next is always 180&deg

Definition at line 59 of file geodesic.for.

* * ! If points

Definition at line 66 of file geodesic.for.

* * ! If and lie on a single geodesic

Definition at line 66 of file geodesic.for.

* * ! If and lie on a single then the following* ! addition rules hold

Definition at line 66 of file geodesic.for.

e & minus

Definition at line 72 of file geodesic.for.

e MM12 e MM21 e* ! m23 e m12* e MM31 = \e MM32 \e MM21 &minus

Definition at line 74 of file geodesic.for.

e MM23 e MM32 e* ! m12 e m23* * ! The shortest distance returned by the solution of the inverse problem* in a few special cases* ! there are multiple azimuths which yield the same shortest distance* ! Here is a catalog of those cases

Definition at line 81 of file geodesic.for.

e the geodesic is unique Otherwise there are two geodesics* ! and the second one is obtained by setting[\e azi1,\e azi2]
Initial value:
 [\e
*!   azi2, \e azi1]

Definition at line 83 of file geodesic.for.

e the geodesic is unique Otherwise there are two geodesics* ! and the second one is obtained by e SS12
Initial value:
*!   &minus

Definition at line 84 of file geodesic.for.

& deg

Definition at line 87 of file geodesic.for.

or& plusmn

Definition at line 88 of file geodesic.for.

the geodesic is unique* ! Otherwise there are two geodesics and the second one is obtained by* !setting[\e azi1,\e azi2] = [&minus

Definition at line 90 of file geodesic.for.

e lon1 and* e azi1 should be in the range [&minus;540&deg;, 540&deg;). The*! values of \e lon2 and \e azi2 returned are in the range*! [&minus;180&deg;, 180&deg;).*!*! If either point is at a pole, the azimuth is defined by keeping the*! longitude fixed, writing \e lat = \e lat = &plusmn;(90&deg; &minus;*! &epsilon;), and taking the limit &epsilon; &rarr; 0+. An arc length*! greater that 180&deg; signifies a geodesic which is not a shortest*! path. (For a prolate ellipsoid, an additional condition is necessary*! for a shortest path: the longitudinal extent must not exceed of*! 180&deg;.) subroutine direct(a, f, lat1, lon1, azi1, s12a12, arcmod, + lat2, lon2, azi2, omask, a12s12, m12, MM12, MM21, SS12)* input double precision a, f, lat1, lon1, azi1, s12a12 logical arcmod integer omask* output double precision lat2, lon2, azi2* optional output double precision a12s12, m12, MM12, MM21, SS12 integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x parameter (ord = 6, nC1 = ord, nC1p = ord, + nC2 = ord, nA3 = ord, nA3x = nA3, + nC3 = ord, nC3x = (nC3 * (nC3 - 1)) / 2, + nC4 = ord, nC4x = (nC4 * (nC4 + 1)) / 2) double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1), + C1a(nC1), C1pa(nC1p), C2a(nC2), C3a(nC3-1), C4a(0:nC4-1) double precision csmgt, atanhx, hypotx, + AngNm, AngNm2, AngRnd, TrgSum, A1m1f, A2m1f, A3f logical arcp, redlp, scalp, areap double precision e2, f1, ep2, n, b, c2, + lon1x, azi1x, phi, alp1, salp0, calp0, k2, eps, + salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1, + salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2, + ssig12, csig12, salp12, calp12, omg12, lam12, lon12, + sig12, stau1, ctau1, tau12, s12a, t, s, c, serr, + A1m1, A2m1, A3c, A4, AB1, AB2, + B11, B12, B21, B22, B31, B41, B42, J12 double precision dblmin, dbleps, pi, degree, tiny, + tol0, tol1, tol2, tolb, xthrsh integer digits, maxit1, maxit2 logical init common /geocom/ dblmin, dbleps, pi, degree, tiny, + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init if (.not.init) call geoini e2 = f * (2 - f) ep2 = e2 / (1 - e2) f1 = 1 - f n = f / (2 - f) b = a * f1 c2 = 0 arcp = mod(omask/1, 2) == 1 redlp = mod(omask/2, 2) == 1 scalp = mod(omask/4, 2) == 1 areap = mod(omask/8, 2) == 1 if (areap) then if (e2 .eq. 0) then c2 = a**2 else if (e2 .gt. 0) then c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2 else c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2 end if end if call A3cof(n, A3x) call C3cof(n, C3x) if (areap) call C4cof(n, C4x)* Guard against underflow in salp0 azi1x = AngRnd(AngNm(azi1)) lon1x = AngNm(lon1)* alp1 is in azi2[0, pi] alp1 = azi1x * degree* Enforce sin(pi) == 0 and cos(pi/2) == 0. Better to face the ensuing* problems directly than to skirt them. salp1 = csmgt(0d0, sin(alp1), azi1x .eq. -180) calp1 = csmgt(0d0, cos(alp1), abs(azi1x) .eq. 90) phi = lat1 * degree* Ensure cbet1 = +dbleps at poles sbet1 = f1 * sin(phi) cbet1 = csmgt(tiny, cos(phi), abs(lat1) .eq. 90) call Norm(sbet1, cbet1) dn1 = sqrt(1 + ep2 * sbet1**2)* Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),* alp0 in [0, pi/2 - |bet1|] salp0 = salp1 * cbet1* Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following* is slightly better (consider the case salp1 = 0). calp0 = hypotx(calp1, salp1 * sbet1)* Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).* sig = 0 is nearest northward crossing of equator.* With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).* With bet1 = pi/2, alp1 = -pi, sig1 = pi/2* With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2* Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).* With alp0 in (0, pi/2], quadrants for sig and omg coincide.* No atan2(0,0) ambiguity at poles since cbet1 = +dbleps.* With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. ssig1 = sbet1 somg1 = salp0 * sbet1 csig1 = csmgt(cbet1 * calp1, 1d0, sbet1 .ne. 0 .or. calp1 .ne. 0) comg1 = csig1* sig1 in (-pi, pi] call Norm (ssig1, csig1) * Geodesic don t need to normalize converts to
Initial value:
*!   [\e azi1, \e azi2] + [\e d, &minus

Definition at line 90 of file geodesic.for.

e e* !SS12 = &minus

Definition at line 91 of file geodesic.for.

e d

Definition at line 95 of file geodesic.for.

e lon1 and *e azi1 should be in the range[&minus540 &deg, 540 &deg).The *!values of\e lon2 and\e azi2 returned are in the range *![&minus180 &deg, 180 &deg).*!*!If either point is at a pole, the azimuth is defined by keeping the *!longitude fixed, writing\e lat=\e lat=&plusmn(90 &deg &minus *!&epsilon) salp12, and taking the limit &epsilon &rarr 0+.An arc length *!greater that 180 &deg signifies a geodesic which is not a shortest *!path.(For a prolate ellipsoid, an additional condition is necessary *!for a shortest path:the longitudinal extent must not exceed of *!180 &deg.) subroutine direct(a, f, lat1, lon1, azi1, s12a12, arcmod,+lat2, lon2, azi2, omask, a12s12, m12, MM12, MM21, SS12)*input double precision a, f, lat1, lon1, azi1, s12a12 logical arcmod integer omask *output double precision lat2, lon2, azi2 *optional output double precision a12s12, m12, MM12, MM21, SS12 integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x parameter(ord=6, nC1=ord, nC1p=ord,+nC2=ord, nA3=ord, nA3x=nA3,+nC3=ord, nC3x=(nC3 *(nC3-1))/2,+nC4=ord, nC4x=(nC4 *(nC4+1))/2) double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1),+C1a(nC1), C1pa(nC1p), C2a(nC2), C3a(nC3-1), C4a(0:nC4-1) double precision csmgt, atanhx, hypotx,+AngNm, AngNm2, AngRnd, TrgSum, A1m1f, A2m1f, A3f logical arcp, redlp, scalp, areap double precision e2, f1, ep2, n, b, c2,+lon1x, azi1x, phi, alp1, salp0, calp0, k2, eps,+salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1,+salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2,+ssig12, csig12, salp12, calp12, omg12, lam12, lon12,+sig12, stau1, ctau1, tau12, s12a, t, s, c, serr,+A1m1, A2m1, A3c, A4, AB1, AB2,+B11, B12, B21, B22, B31, B41, B42, J12 double precision dblmin, dbleps, pi, degree, tiny,+tol0, tol1, tol2, tolb, xthrsh integer digits, maxit1, maxit2 logical init common/geocom/dblmin, dbleps, pi, degree, tiny,+tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init if(.not.init) call geoini e2=f *(2-f) ep2=e2/(1-e2) f1=1-f n=f/(2-f) b=a *f1 c2=0 arcp=mod(omask/1, 2)==1 redlp=mod(omask/2, 2)==1 scalp=mod(omask/4, 2)==1 areap=mod(omask/8, 2)==1 if(areap) then if(e2.eq.0) then c2=a **2 else if(e2.gt.0) then c2=(a **2+b **2 *atanhx(sqrt(e2))/sqrt(e2))/2 else c2=(a **2+b **2 *atan(sqrt(abs(e2)))/sqrt(abs(e2)))/2 end if end if call A3cof(n, A3x) call C3cof(n, C3x) if(areap) call C4cof(n, C4x)*Guard against underflow in salp0 azi1x=AngRnd(AngNm(azi1)) lon1x=AngNm(lon1)*alp1 is in[0, pi] alp1=azi1x *degree *Enforce sin(pi)==0 and cos(pi/2)==0.Better to face the ensuing *problems directly than to skirt them.salp1=csmgt(0d0, sin(alp1), azi1x.eq.-180) calp1=csmgt(0d0, cos(alp1), abs(azi1x).eq.90) phi=lat1 *degree *Ensure cbet1=+dbleps at poles sbet1=f1 *sin(phi) cbet1=csmgt(tiny, cos(phi), abs(lat1).eq.90) call Norm(sbet1, cbet1) dn1=sqrt(1+ep2 *sbet1 **2)*Evaluate alp0 from sin(alp1)*cos(bet1)=sin(alp0),*alp0 in[0, pi/2-|bet1|] salp0=salp1 *cbet1 *Alt:calp0=hypot(sbet1, calp1 *cbet1).The following *is slightly better(consider the case salp1=0).calp0=hypotx(calp1, salp1 *sbet1)*Evaluate sig with tan(bet1)=tan(sig1)*cos(alp1).*sig=0 is nearest northward crossing of equator.*With bet1=0, alp1=pi/2, we have sig1=0(equatorial line).*With bet1=pi/2, alp1=-pi, sig1=pi/2 *With bet1=-pi/2, alp1=0, sig1=-pi/2 *Evaluate omg1 with tan(omg1)=sin(alp0)*tan(sig1).*With alp0 in(0, pi/2], quadrants for sig and omg coincide.*No atan2(0, 0) ambiguity at poles since cbet1=+dbleps.*With alp0=0, omg1=0 for alp1=0, omg1=pi for alp1=pi.ssig1=sbet1 somg1=salp0 *sbet1 csig1=csmgt(cbet1 *calp1, 1d0, sbet1.ne.0.or.calp1.ne.0) comg1=csig1 *sig1 in(-pi, pi] call Norm(ssig1, csig1)*Geodesic don t need to normalize converts to used in atan2 so no need to normalized viz * (  ) 
Initial value:
 0 - atan2(-salp2, calp2) / degree

      if (redlp .or. scalp) then
        B22 = TrgSum(.true., ssig2, csig2, C2a, nC2)
        AB2 = (1 + A2m1) * (B22 - B21)
        J12 = (A1m1 - A2m1) * sig12 + (AB1 - AB2)
      end if
* Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
* accurate cancellation in the case of coincident points.
      if (redlp) m12 = b * ((dn2 * (csig1 * ssig2) -
     +    dn1 * (ssig1 * csig2)) - csig1 * csig2 * J12)
      if (scalp) then
        t = k2 * (ssig2 - ssig1) * (ssig2 + ssig1) / (dn1 + dn2)
        MM12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1
        MM21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2
      end if

      if (areap) then
        B42 = TrgSum(.false., ssig2, csig2, C4a, nC4)
        if (calp0 .eq. 0 .or. salp0 .eq. 0) then
* alp12 = alp2 - alp1

Definition at line 157 of file geodesic.for.

e lon1 and* e azi1 should be in the range [&minus;540&deg;, 540&deg;). The*! values of \e lon2 and \e azi2 returned are in the range*! [&minus;180&deg;, 180&deg;).*!*! If either point is at a pole, the azimuth is defined by keeping the*! longitude fixed, writing \e lat = \e lat = &plusmn;(90&deg; &minus;*! &epsilon;), and taking the limit &epsilon; &rarr; 0+. An arc length*! greater that 180&deg; signifies a geodesic which is not a shortest*! path. (For a prolate ellipsoid, an additional condition is necessary*! for a shortest path: the longitudinal extent must not exceed of*! 180&deg;.) subroutine direct(a, f, lat1, lon1, azi1, s12a12, arcmod, + lat2, lon2, azi2, omask, a12s12, m12, MM12, MM21, SS12)* input double precision a, f, lat1, lon1, azi1, s12a12 logical arcmod integer omask* output double precision lat2, lon2, azi2* optional output double precision a12s12, m12, MM12, MM21, SS12 integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x parameter (ord = 6, nC1 = ord, nC1p = ord, + nC2 = ord, nA3 = ord, nA3x = nA3, + nC3 = ord, nC3x = (nC3 * (nC3 - 1)) / 2, + nC4 = ord, nC4x = (nC4 * (nC4 + 1)) / 2) double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1), + C1a(nC1), C1pa(nC1p), C2a(nC2), C3a(nC3-1), C4a(0:nC4-1) double precision csmgt, atanhx, hypotx, + AngNm, AngNm2, AngRnd, TrgSum, A1m1f, A2m1f, A3f logical arcp, redlp, scalp, areap double precision e2, f1, ep2, n, b, c2, + lon1x, azi1x, phi, alp1, salp0, calp0, k2, eps, + salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1, + salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2, + ssig12, csig12, salp12, calp12, omg12, lam12, lon12, + sig12, stau1, ctau1, tau12, s12a, t, s, c, serr, + A1m1, A2m1, A3c, A4, AB1, AB2, + B11, B12, B21, B22, B31, B41, B42, J12 double precision dblmin, dbleps, pi, degree, tiny, + tol0, tol1, tol2, tolb, xthrsh integer digits, maxit1, maxit2 logical init common /geocom/ dblmin, dbleps, pi, degree, tiny, + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init if (.not.init) call geoini e2 = f * (2 - f) ep2 = e2 / (1 - e2) f1 = 1 - f n = f / (2 - f) b = a * f1 c2 = 0 arcp = mod(omask/1, 2) == 1 redlp = mod(omask/2, 2) == 1 scalp = mod(omask/4, 2) == 1 areap = mod(omask/8, 2) == 1 if (areap) then if (e2 .eq. 0) then c2 = a**2 else if (e2 .gt. 0) then c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2 else c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2 end if end if call A3cof(n, A3x) call C3cof(n, C3x) if (areap) call C4cof(n, C4x)* Guard against underflow in salp0 azi1x = AngRnd(AngNm(azi1)) lon1x = AngNm(lon1)* alp1 is in !k2[0, pi] alp1 = azi1x * degree* Enforce sin(pi) == 0 and cos(pi/2) == 0. Better to face the ensuing* problems directly than to skirt them. salp1 = csmgt(0d0, sin(alp1), azi1x .eq. -180) calp1 = csmgt(0d0, cos(alp1), abs(azi1x) .eq. 90) phi = lat1 * degree* Ensure cbet1 = +dbleps at poles sbet1 = f1 * sin(phi) cbet1 = csmgt(tiny, cos(phi), abs(lat1) .eq. 90) call Norm(sbet1, cbet1) dn1 = sqrt(1 + ep2 * sbet1**2)* Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),* alp0 in [0, pi/2 - |bet1|] salp0 = salp1 * cbet1* Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following* is slightly better (consider the case salp1 = 0). calp0 = hypotx(calp1, salp1 * sbet1)* Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).* sig = 0 is nearest northward crossing of equator.* With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).* With bet1 = pi/2, alp1 = -pi, sig1 = pi/2* With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2* Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).* With alp0 in (0, pi/2], quadrants for sig and omg coincide.* No atan2(0,0) ambiguity at poles since cbet1 = +dbleps.* With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. ssig1 = sbet1 somg1 = salp0 * sbet1 csig1 = csmgt(cbet1 * calp1, 1d0, sbet1 .ne. 0 .or. calp1 .ne. 0) comg1 = csig1* sig1 in (-pi, pi] call Norm (ssig1, csig1) * Geodesic don t need to normalize

Definition at line 277 of file geodesic.for.

e lon1 and* e azi1 should be in the range [&minus;540&deg;, 540&deg;). The*! values of \e lon2 and \e azi2 returned are in the range*! [&minus;180&deg;, 180&deg;).*!*! If either point is at a pole, the azimuth is defined by keeping the*! longitude fixed, writing \e lat = \e lat = &plusmn;(90&deg; &minus;*! &epsilon;), and taking the limit &epsilon; &rarr; 0+. An arc length*! greater that 180&deg; signifies a geodesic which is not a shortest*! path. (For a prolate ellipsoid, an additional condition is necessary*! for a shortest path: the longitudinal extent must not exceed of*! 180&deg;.) subroutine direct(a, f, lat1, lon1, azi1, s12a12, arcmod, + lat2, lon2, azi2, omask, a12s12, m12, MM12, MM21, SS12)* input double precision a, f, lat1, lon1, azi1, s12a12 logical arcmod integer omask* output double precision lat2, lon2, azi2* optional output double precision a12s12, m12, MM12, MM21, SS12 integer ord, nC1, nC1p, nC2, nA3, nA3x, nC3, nC3x, nC4, nC4x parameter (ord = 6, nC1 = ord, nC1p = ord, + nC2 = ord, nA3 = ord, nA3x = nA3, + nC3 = ord, nC3x = (nC3 * (nC3 - 1)) / 2, + nC4 = ord, nC4x = (nC4 * (nC4 + 1)) / 2) double precision A3x(0:nA3x-1), C3x(0:nC3x-1), C4x(0:nC4x-1), + C1a(nC1), C1pa(nC1p), C2a(nC2), C3a(nC3-1), C4a(0:nC4-1) double precision csmgt, atanhx, hypotx, + AngNm, AngNm2, AngRnd, TrgSum, A1m1f, A2m1f, A3f logical arcp, redlp, scalp, areap double precision e2, f1, ep2, n, b, c2, + lon1x, azi1x, phi, alp1, salp0, calp0, k2, eps, + salp1, calp1, ssig1, csig1, cbet1, sbet1, dn1, somg1, comg1, + salp2, calp2, ssig2, csig2, sbet2, cbet2, dn2, somg2, comg2, + ssig12, csig12, salp12, calp12, omg12, lam12, lon12, + sig12, stau1, ctau1, tau12, s12a, t, s, c, serr, + A1m1, A2m1, A3c, A4, AB1, AB2, + B11, B12, B21, B22, B31, B41, B42, J12 double precision dblmin, dbleps, pi, degree, tiny, + tol0, tol1, tol2, tolb, xthrsh integer digits, maxit1, maxit2 logical init common /geocom/ dblmin, dbleps, pi, degree, tiny, + tol0, tol1, tol2, tolb, xthrsh, digits, maxit1, maxit2, init if (.not.init) call geoini e2 = f * (2 - f) ep2 = e2 / (1 - e2) f1 = 1 - f n = f / (2 - f) b = a * f1 c2 = 0 arcp = mod(omask/1, 2) == 1 redlp = mod(omask/2, 2) == 1 scalp = mod(omask/4, 2) == 1 areap = mod(omask/8, 2) == 1 if (areap) then if (e2 .eq. 0) then c2 = a**2 else if (e2 .gt. 0) then c2 = (a**2 + b**2 * atanhx(sqrt(e2)) / sqrt(e2)) / 2 else c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2 end if end if call A3cof(n, A3x) call C3cof(n, C3x) if (areap) call C4cof(n, C4x)* Guard against underflow in salp0 azi1x = AngRnd(AngNm(azi1)) lon1x = AngNm(lon1)* alp1 is in salp12[0, pi] alp1 = azi1x * degree* Enforce sin(pi) == 0 and cos(pi/2) == 0. Better to face the ensuing* problems directly than to skirt them. salp1 = csmgt(0d0, sin(alp1), azi1x .eq. -180) calp1 = csmgt(0d0, cos(alp1), abs(azi1x) .eq. 90) phi = lat1 * degree* Ensure cbet1 = +dbleps at poles sbet1 = f1 * sin(phi) cbet1 = csmgt(tiny, cos(phi), abs(lat1) .eq. 90) call Norm(sbet1, cbet1) dn1 = sqrt(1 + ep2 * sbet1**2)* Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),* alp0 in [0, pi/2 - |bet1|] salp0 = salp1 * cbet1* Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following* is slightly better (consider the case salp1 = 0). calp0 = hypotx(calp1, salp1 * sbet1)* Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).* sig = 0 is nearest northward crossing of equator.* With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).* With bet1 = pi/2, alp1 = -pi, sig1 = pi/2* With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2* Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).* With alp0 in (0, pi/2], quadrants for sig and omg coincide.* No atan2(0,0) ambiguity at poles since cbet1 = +dbleps.* With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi. ssig1 = sbet1 somg1 = salp0 * sbet1 csig1 = csmgt(cbet1 * calp1, 1d0, sbet1 .ne. 0 .or. calp1 .ne. 0) comg1 = csig1* sig1 in (-pi, pi] call Norm (ssig1, csig1) * Geodesic don t need to normalize converts to used in atan2 so no need to normalized viz* ( pi   ) 
Initial value:
 salp2 * calp1 - calp2 * salp1
          calp12 = calp2 * calp1 + salp2 * salp1
* The right thing appears to happen if alp1 = +/-180 and alp2 = 0

Definition at line 431 of file geodesic.for.

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends

Generated on 6 Oct 2014 for Fortran library for Geodesics by  doxygen 1.6.1