GRASS Programmer's Manual  6.4.1(2011)
geodist.c
Go to the documentation of this file.
00001 
00026 #include <math.h>
00027 #include <grass/gis.h>
00028 #include "pi.h"
00029 
00030 
00031 static double boa;
00032 static double f;
00033 static double ff64;
00034 static double al;
00035 static double t1, t2, t3, t4, t1r, t2r;
00036 
00037 
00052 int G_begin_geodesic_distance(double a, double e2)
00053 {
00054     al = a;
00055     boa = sqrt(1 - e2);
00056     f = 1 - boa;
00057     ff64 = f * f / 64;
00058 
00059     return 0;
00060 }
00061 
00062 
00063 
00075 int G_set_geodesic_distance_lat1(double lat1)
00076 {
00077     t1r = atan(boa * tan(Radians(lat1)));
00078 
00079     return 0;
00080 }
00081 
00082 
00094 int G_set_geodesic_distance_lat2(double lat2)
00095 {
00096     double stm, ctm, sdtm, cdtm;
00097     double tm, dtm;
00098 
00099     t2r = atan(boa * tan(Radians(lat2)));
00100 
00101     tm = (t1r + t2r) / 2;
00102     dtm = (t2r - t1r) / 2;
00103 
00104     stm = sin(tm);
00105     ctm = cos(tm);
00106     sdtm = sin(dtm);
00107     cdtm = cos(dtm);
00108 
00109     t1 = stm * cdtm;
00110     t1 = t1 * t1 * 2;
00111 
00112     t2 = sdtm * ctm;
00113     t2 = t2 * t2 * 2;
00114 
00115     t3 = sdtm * sdtm;
00116     t4 = cdtm * cdtm - stm * stm;
00117 
00118     return 0;
00119 }
00120 
00121 
00135 double G_geodesic_distance_lon_to_lon(double lon1, double lon2)
00136 {
00137     double a, cd, d, e,         /*dl, */
00138       q, sd, sdlmr, t, u, v, x, y;
00139 
00140 
00141     sdlmr = sin(Radians(lon2 - lon1) / 2);
00142 
00143     /* special case - shapiro */
00144     if (sdlmr == 0.0 && t1r == t2r)
00145         return 0.0;
00146 
00147     q = t3 + sdlmr * sdlmr * t4;
00148 
00149     /* special case - shapiro */
00150     if (q == 1.0)
00151         return M_PI * al;
00152 
00153     /* Mod: shapiro
00154      * cd=1-2q is ill-conditioned if q is small O(10**-23)
00155      *   (for high lats? with lon1-lon2 < .25 degrees?)
00156      *   the computation of cd = 1-2*q will give cd==1.0.
00157      * However, note that t=dl/sd is dl/sin(dl) which approaches 1 as dl->0.
00158      * So the first step is to compute a good value for sd without using sin()
00159      *   and then check cd && q to see if we got cd==1.0 when we shouldn't.
00160      * Note that dl isn't used except to get t,
00161      *   but both cd and sd are used later
00162      */
00163 
00164     /* original code
00165        cd=1-2*q;
00166        dl=acos(cd);
00167        sd=sin(dl);
00168        t=dl/sd;
00169      */
00170 
00171     cd = 1 - 2 * q;             /* ill-conditioned subtraction for small q */
00172     /* mod starts here */
00173     sd = 2 * sqrt(q - q * q);   /* sd^2 = 1 - cd^2 */
00174     if (q != 0.0 && cd == 1.0)  /* test for small q */
00175         t = 1.0;
00176     else if (sd == 0.0)
00177         t = 1.0;
00178     else
00179         t = acos(cd) / sd;      /* don't know how to fix acos(1-2*q) yet */
00180     /* mod ends here */
00181 
00182     u = t1 / (1 - q);
00183     v = t2 / q;
00184     d = 4 * t * t;
00185     x = u + v;
00186     e = -2 * cd;
00187     y = u - v;
00188     a = -d * e;
00189 
00190     return (al * sd *
00191             (t - f / 4 * (t * x - y) +
00192              ff64 * (x * (a + (t - (a + e) / 2) * x) + y * (-2 * d + e * y)
00193                      + d * x * y)
00194             )
00195         );
00196 }
00197 
00198 
00212 double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
00213 {
00214     G_set_geodesic_distance_lat1(lat1);
00215     G_set_geodesic_distance_lat2(lat2);
00216     return G_geodesic_distance_lon_to_lon(lon1, lon2);
00217 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines