geodist.c

Go to the documentation of this file.
00001 #include <math.h>
00002 #include <grass/gis.h>
00003 #include "pi.h"
00004 
00005 /* WARNING: this code is preliminary and may be changed,
00006  * including calling sequences to any of the functions
00007  * defined here
00008  */
00009 
00010 /* distance from point to point along a geodesic 
00011  * code from
00012  *   Paul D. Thomas
00013  *   "Spheroidal Geodesics, Reference Systems, and Local Geometry"
00014  *   U.S. Naval Oceanographic Office, p. 162
00015  *   Engineering Library 526.3 T36s
00016  */
00017 
00018 static double boa;
00019 static double f;
00020 static double ff64;
00021 static double al;
00022 
00023 /* must be called once to establish the ellipsoid */
00024 
00037 int G_begin_geodesic_distance(double a,double e2)
00038 {
00039     al = a;
00040     boa = sqrt (1 - e2);
00041     f = 1 - boa;
00042     ff64 = f*f/64;
00043 
00044     return 0;
00045 }
00046 
00047 static double t1,t2,t3,t4,t1r,t2r;
00048 
00049 /* must be called first */
00050 
00060 int G_set_geodesic_distance_lat1(double lat1)
00061 {
00062     t1r = atan(boa*tan(Radians(lat1)));
00063 
00064     return 0;
00065 }
00066 /* must be called second */
00067 
00077 int G_set_geodesic_distance_lat2( double lat2)
00078 {
00079     double stm,ctm,sdtm,cdtm;
00080     double tm, dtm;
00081 
00082     t2r = atan(boa*tan(Radians(lat2)));
00083 
00084     tm  = (t1r+t2r)/2;
00085     dtm = (t2r-t1r)/2;
00086 
00087     stm = sin(tm);
00088     ctm = cos(tm);
00089     sdtm = sin(dtm);
00090     cdtm = cos(dtm);
00091 
00092     t1 = stm*cdtm;
00093     t1 = t1 * t1 * 2;
00094 
00095     t2 = sdtm*ctm;
00096     t2 = t2 * t2 * 2;
00097 
00098     t3 = sdtm*sdtm;
00099     t4 = cdtm*cdtm-stm*stm;
00100 
00101     return 0;
00102 }
00103 
00104 
00119 double G_geodesic_distance_lon_to_lon (double lon1,double lon2)
00120 {
00121     double a, cd, d, e, /*dl,*/
00122            q, sd, sdlmr, 
00123            t, u, v, x, y;
00124 
00125 
00126     sdlmr = sin(Radians(lon2-lon1)/2);
00127 
00128     /* special case - shapiro */
00129     if (sdlmr == 0.0 && t1r == t2r) return 0.0;
00130 
00131     q = t3+sdlmr*sdlmr*t4;
00132 
00133     /* special case - shapiro */
00134     if (q == 1.0) return PI *al;
00135 
00136 /* Mod: shapiro
00137  * cd=1-2q is ill-conditioned if q is small O(10**-23)
00138  *   (for high lats? with lon1-lon2 < .25 degrees?)
00139  *   the computation of cd = 1-2*q will give cd==1.0.
00140  * However, note that t=dl/sd is dl/sin(dl) which approaches 1 as dl->0.
00141  * So the first step is to compute a good value for sd without using sin()
00142  *   and then check cd && q to see if we got cd==1.0 when we shouldn't.
00143  * Note that dl isn't used except to get t,
00144  *   but both cd and sd are used later
00145  */
00146 
00147 /* original code
00148     cd=1-2*q;
00149     dl=acos(cd);
00150     sd=sin(dl);
00151     t=dl/sd;
00152 */
00153 
00154     cd = 1-2*q;                 /* ill-conditioned subtraction for small q */
00155 /* mod starts here */
00156     sd = 2* sqrt (q - q * q);   /* sd^2 = 1 - cd^2 */
00157     if (q != 0.0 && cd == 1.0)  /* test for small q */
00158         t = 1.0;
00159     else if (sd == 0.0)
00160         t = 1.0;
00161     else
00162         t = acos(cd)/sd;          /* don't know how to fix acos(1-2*q) yet */
00163 /* mod ends here */
00164 
00165     u = t1/(1-q);
00166     v = t2/q;
00167     d = 4*t*t;
00168     x = u+v;
00169     e = -2*cd;
00170     y = u-v;
00171     a = -d*e;
00172 
00173     return ( al * sd *
00174                 ( t - f/4 * (t*x-y) +
00175                   ff64 *
00176                   (
00177                     x * ( a + (t - (a+e)/2) * x) +
00178                     y* (-2 * d + e*y)
00179                     + d*x*y
00180                   )
00181                 )
00182             );
00183 }
00184 
00185 
00202 double G_geodesic_distance (double lon1,double lat1,double lon2,double lat2)
00203 {
00204     G_set_geodesic_distance_lat1 (lat1);
00205     G_set_geodesic_distance_lat2 (lat2);
00206     return G_geodesic_distance_lon_to_lon (lon1, lon2);
00207 }

Generated on Fri Nov 21 11:02:18 2008 for GRASS by  doxygen 1.5.1