function [dis az1 az2 ittime itnum ]= GeodesicInversec(lat1,lon1,lat2,lon2,Geod).
GeodesicLib1.rar is the zip file of function
[dis az1 az2] = GeodesicInversec(lat1,lon1,lat2,lon2,Geod) computes the geodesic distance and azimuths assuming that the points lie on the reference ellipsoid defined by the input Geod.
The input latitudes and longitudes, lat1, lon1, lat2, lon2, can be scalars or arrays of equal size and must be expressed in degrees. The Geoid vector is of the form Geod=[semimajor axis, eccentricity]. The output dis (distance) is expressed in the same distance units as the semimajor axis of the ellipsoid vector. The output az1 and az2 are the initial and final azimuths at the two endpoints in degrees. Additionally , the outputs ittime and itnum compute the time (second) and number of iteration of Newton's Method. If ellipsoid is omitted, the WGS84 Geod=[6378137 0.0818191908426215] is used.
This function duplicates some of the functionality of the DISTANCE function in the MATLAB mapping toolbox and GEODDISTANCE function in Karney's geographiclib.
The differences and equvalent abilities are:
(1) The series expansion of the rate of change of the longitude difference with respect to Clairaut constant provided here has the same coefficient matrices as provided in the inte-gral of longitude difference and is derived without the series expansion of integral of the distance. The simpler series expansion may reduce the computational cost in the iteration of Newton’s method for the inverse geodesic sailing.
(2) The accuracy of computation with sufficient terms in the series expansions is able to comply with the standard precision of 15 significant digits of double-precision floating-point format specified by IEEE 754 standard which has the almost same ability to Karney (2013).
(3) Any geodesic between two points can be determined after iterating Newton’s method with specified initial value to converge for all geodesic between two points after several steps by the algorithm given here.
(4) The differential quantity of longitude difference with respect to Clairaut constant or initial azimuth can be directly computed with the series expansion of function of longitude difference in terms of latitude and without the extra series expansion derived from the reduced length m12 and the final azimuth (Karney, 2013) which is the kernel in the solution of the inverse problem by Newton’s method.
(5) In cases of two points near antipodal regions, the rate of longitude difference is sensitively. The inadequate initial value will lead the process of Newton’s method to diverge or very slow converge (Karney, 2013). The adaptation (Karney, 2013) of the prescription given by Helmert (1880) can overcome this problem and accelerate the speed of conver-gence which can obtain the approximated initial azimuth by astroid plane. In our derivation following this method, the neighborhood of the antipodal point is to be treated as a plane. The derivation is a little straightforward and concise which gives the similar quartic equation to Karney’s derivation (2013).
GeodesicLib1.rar is the zip file of function
email:wilco0917@yahoo.com.tw