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 >|< |
*for the WGS84 the errors are less than * | !nanometers (Reasonably accurate results are obtained for|< i >f</i >|*!<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 *!±180 °for oblate ellipsoids.)*!-\e lon2 |
e | SS12 (This occurs when\e lat2 is near *!−\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° |
*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 |
i | 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 [−90 °, 90 °] |
e lon1 and *e azi1 should be in the range[−540 °, 540 °).The *!values of\e lon2 and\e azi2 returned are in the range *![−180 °, 180 °).*!*!If either point is at a pole, the azimuth is defined by keeping the *!longitude fixed, writing\e lat=\e lat=±(90 °−*!ε), and taking the limit ε→0+.An arc length *!greater that 180 °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 °.) 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[−540 °, 540 °).The *!values of\e lon2 and\e azi2 returned are in the range *![−180 °, 180 °).*!*!If either point is at a pole, the azimuth is defined by keeping the *!longitude fixed, writing\e lat=\e lat=±(90 °−*!ε), and taking the limit ε→0+.An arc length *!greater that 180 °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 °.) 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 |
* 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 | ( | ) |
* 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 | ( | ) |
* for the WGS84 the errors are less than* !nanometers | ( | Reasonably accurate results are obtained for|< i >f</i >|*!<1/ | 5. | ) |
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 | ) |
e MM23 e MM32 e* ! m12 e m23* * ! The shortest distance returned by the solution of the inverse problem* !is | ( | obviously | ) |
e SS12 | ( | This occurs when the longitude difference is near *!±180 °for oblate | ellipsoids. | ) |
e for arbitrary e d | ( | For *! | spheres, | |
this prescription applies when points 1 and 2 are *! | antipodal. | |||
) | [pure virtual] |
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] |
otherwise it is the distance *between point and point | ( | meters | ) |
otherwise it is the arc length* ! between point and point | ( | degrees | ) |
* 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.
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.
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.
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.
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 |
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 then* ! they are separated by a distance e MM12 e dt at point e MM21* ! is defined we have e MM12 |
\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°
Definition at line 59 of file geodesic.for.
* * ! If points |
Definition at line 66 of file geodesic.for.
Definition at line 66 of file geodesic.for.
Definition at line 66 of file geodesic.for.
Definition at line 72 of file geodesic.for.
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 e SS12 |
*! &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 [−540°, 540°). The*! values of \e lon2 and \e azi2 returned are in the range*! [−180°, 180°).*!*! If either point is at a pole, the azimuth is defined by keeping the*! longitude fixed, writing \e lat = \e lat = ±(90° −*! ε), and taking the limit ε → 0+. An arc length*! greater that 180° 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°.) 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 |
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 °, 540 °).The *!values of\e lon2 and\e azi2 returned are in the range *![&minus180 °, 180 °).*!*!If either point is at a pole, the azimuth is defined by keeping the *!longitude fixed, writing\e lat=\e lat=±(90 ° &minus *!&epsilon) salp12, and taking the limit &epsilon &rarr 0+.An arc length *!greater that 180 ° 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 °.) 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 * | ( | ) |
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 [−540°, 540°). The*! values of \e lon2 and \e azi2 returned are in the range*! [−180°, 180°).*!*! If either point is at a pole, the azimuth is defined by keeping the*! longitude fixed, writing \e lat = \e lat = ±(90° −*! ε), and taking the limit ε → 0+. An arc length*! greater that 180° 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°.) 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 [−540°, 540°). The*! values of \e lon2 and \e azi2 returned are in the range*! [−180°, 180°).*!*! If either point is at a pole, the azimuth is defined by keeping the*! longitude fixed, writing \e lat = \e lat = ±(90° −*! ε), and taking the limit ε → 0+. An arc length*! greater that 180° 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°.) 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 | ) |
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.