00001 #include "gis.h"
00002 #include "pi.h"
00003
00004
00005
00006 double
00007 G_ellipsoid_grid_dist (a, e2, lon1, lat1, lon2, lat2)
00008 double a, e2, lon1, lat1, lon2, lat2;
00009 {
00010 double cos1, cos2, sin1, sin2;
00011 double M, m;
00012 double M2, m2;
00013 double M4, m4;
00014 double M6, m6;
00015 double A,B,C,D;
00016 double Ap, Bp, Cp, Dp;
00017 double distance;
00018 double e4, e6;
00019 double sin(), cos(), sqrt();
00020 double sc();
00021
00022
00023 if (lon1 > lon2)
00024 while ((lon1-lon2) > 180)
00025 lon2 += 360;
00026 else if (lon2 > lon1)
00027 while ((lon2-lon1) > 180)
00028 lon1 += 360;
00029
00030 lon1 = Radians(lon1);
00031 lon2 = Radians(lon2);
00032 lat1 = Radians(lat1);
00033 lat2 = Radians(lat2);
00034
00035 sin1 = sin(lat1);
00036 cos1 = cos(lon1);
00037 sin2 = sin(lat2);
00038 cos2 = cos(lon2);
00039
00040 e4 = e2 * e2;
00041 e6 = e2 * e4;
00042
00043
00044 if (lat1 == lat2)
00045 {
00046 distance = a * cos1 * (lon2 - lon1) / sqrt(1 - e2 * sin1 * sin1);
00047 }
00048
00049
00050 else if (lon1 == lon2)
00051 {
00052 A = 1 + e2 * (3.0/4.0) + e4 * (45.0/64.0) + e6 * (175.0/256.0);
00053 B = - e2 * (3.0/4.0) - e4 * (45.0/64.0) - e6 * (175.0/256.0);
00054 C = - e4 * (15.0/32.0) - e6 * (175.0/384.0);
00055 D = - e6 * (35.0/96.0);
00056
00057 distance = a * (1-e2) * (sc(A,B,C,D,lat2, 0., 0.)-sc(A,B,C,D,lat1, 0., 0.));
00058 }
00059
00060
00061 else
00062 {
00063 M = 1- e2;
00064 M2 = M*M;
00065 M4 = M2 * M2;
00066 M6 = M2 * M4;
00067
00068 m = (lat2-lat1)/(lon2-lon1);
00069 m2 = m*m;
00070 m4 = m2 * m2;
00071 m6 = m2 * m4;
00072
00073 A = sqrt(M2 + m2);
00074 B = (3.0*M2 + m2)/(2.0*A);
00075 C = (15.0*M4 + 22.0*M2*m2 + 3.0*m4)/(8.0*A*A*A);
00076 D = (35.0*M6 + 87.0*M4*m2 + 65.0*M2*m4 + 5.0*m6)/(16.0*A*A*A*A*A);
00077
00078 Ap = A + B*e2/2.0 + 3.0*C*e4/8.0 + 5.0*D*e6/16.0;
00079 Bp = - B*e2/2.0 - 3.0*C*e4/8.0 - 5.0*D*e6/16.0;
00080 Cp = - C*e4/4.0 - 5.0*D*e6/24.0;
00081 Dp = - D*e6/6.0;
00082
00083 distance = a * a *
00084 (sc(Ap,Bp,Cp,Dp,lat2,sin2,cos2)-sc(Ap,Bp,Cp,Dp,lat1,sin1,cos1));
00085 }
00086
00087 if (distance < 0.0)
00088 distance = -distance;
00089 return distance;
00090 }
00091
00092 static double sc(A,B,C,D,lat,s,c)
00093 double A,B,C,D,lat,s,c;
00094 {
00095 double s2;
00096
00097 s2 = s*s;
00098
00099 return A * lat + c*s*(B + s2 * (C + s2*D));
00100 }