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
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #define SWAP(a,b) temp=a;a=b;b=temp
00032
00033 static int adjust_lat(double *);
00034 static int adjust_lon(double *);
00035
00036 static double A, B;
00037
00038
00039 int G_begin_geodesic_equation(double lon1, double lat1, double lon2,
00040 double lat2)
00041 {
00042 double sin21, tan1, tan2;
00043
00044 adjust_lon(&lon1);
00045 adjust_lon(&lon2);
00046 adjust_lat(&lat1);
00047 adjust_lat(&lat2);
00048 if (lon1 > lon2) {
00049 register double temp;
00050
00051 SWAP(lon1, lon2);
00052 SWAP(lat1, lat2);
00053 }
00054 if (lon1 == lon2) {
00055 A = B = 0.0;
00056 return 0;
00057 }
00058 lon1 = Radians(lon1);
00059 lon2 = Radians(lon2);
00060 lat1 = Radians(lat1);
00061 lat2 = Radians(lat2);
00062
00063 sin21 = sin(lon2 - lon1);
00064 tan1 = tan(lat1);
00065 tan2 = tan(lat2);
00066
00067 A = (tan2 * cos(lon1) - tan1 * cos(lon2)) / sin21;
00068 B = (tan2 * sin(lon1) - tan1 * sin(lon2)) / sin21;
00069
00070 return 1;
00071 }
00072
00073
00074
00075 double G_geodesic_lat_from_lon(double lon)
00076 {
00077 adjust_lon(&lon);
00078 lon = Radians(lon);
00079
00080 return Degrees(atan(A * sin(lon) - B * cos(lon)));
00081 }
00082
00083 static int adjust_lon(double *lon)
00084 {
00085 while (*lon > 180.0)
00086 *lon -= 360.0;
00087 while (*lon < -180.0)
00088 *lon += 360.0;
00089
00090 return 0;
00091 }
00092
00093 static int adjust_lat(double *lat)
00094 {
00095 if (*lat > 90.0)
00096 *lat = 90.0;
00097 if (*lat < -90.0)
00098 *lat = -90.0;
00099
00100 return 0;
00101 }