GRASS Programmer's Manual
6.4.2(2012)
|
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 }