00001 #include <math.h>
00002 #include <grass/gis.h>
00003 #include "pi.h"
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 static double boa;
00019 static double f;
00020 static double ff64;
00021 static double al;
00022
00023
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
00050
00060 int G_set_geodesic_distance_lat1(double lat1)
00061 {
00062 t1r = atan(boa*tan(Radians(lat1)));
00063
00064 return 0;
00065 }
00066
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,
00122 q, sd, sdlmr,
00123 t, u, v, x, y;
00124
00125
00126 sdlmr = sin(Radians(lon2-lon1)/2);
00127
00128
00129 if (sdlmr == 0.0 && t1r == t2r) return 0.0;
00130
00131 q = t3+sdlmr*sdlmr*t4;
00132
00133
00134 if (q == 1.0) return PI *al;
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 cd = 1-2*q;
00155
00156 sd = 2* sqrt (q - q * q);
00157 if (q != 0.0 && cd == 1.0)
00158 t = 1.0;
00159 else if (sd == 0.0)
00160 t = 1.0;
00161 else
00162 t = acos(cd)/sd;
00163
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 }