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,
00138 q, sd, sdlmr, t, u, v, x, y;
00139
00140
00141 sdlmr = sin(Radians(lon2 - lon1) / 2);
00142
00143
00144 if (sdlmr == 0.0 && t1r == t2r)
00145 return 0.0;
00146
00147 q = t3 + sdlmr * sdlmr * t4;
00148
00149
00150 if (q == 1.0)
00151 return M_PI * al;
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 cd = 1 - 2 * q;
00172
00173 sd = 2 * sqrt(q - q * q);
00174 if (q != 0.0 && cd == 1.0)
00175 t = 1.0;
00176 else if (sd == 0.0)
00177 t = 1.0;
00178 else
00179 t = acos(cd) / sd;
00180
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 }