grid_dist1.c

Go to the documentation of this file.
00001 #include "gis.h"
00002 #include "pi.h"
00003 
00004 /* I don't think the formula used here is correct */
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 /* adjust longitudes so that they are as close as can be */
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 /* along a parallel */
00044     if (lat1 == lat2)
00045     {
00046         distance = a * cos1 * (lon2 - lon1) / sqrt(1 - e2 * sin1 * sin1);
00047     }
00048 
00049 /* along a meridian */
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 /* other lines */
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 }

Generated on Sat Jul 22 22:06:14 2006 for GRASS by  doxygen 1.4.7