00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include <math.h>
00020 #include <grass/gis.h>
00021 #include "pi.h"
00022
00023 static double TAN_A, TAN1, TAN2, L;
00024 static int parallel;
00025 static int adjust_lat(double *);
00026 static int adjust_lon(double *);
00027
00028 int G_begin_rhumbline_equation (
00029 double lon1,double lat1,double lon2,double lat2)
00030 {
00031 adjust_lat (&lat1);
00032 adjust_lat (&lat2);
00033
00034 if (lon1 == lon2)
00035 {
00036 parallel = 1;
00037 L = lat1;
00038 return 0;
00039 }
00040 if (lat1 == lat2)
00041 {
00042 parallel = 1;
00043 L = lat1;
00044 return 1;
00045 }
00046 parallel = 0;
00047 lon1 = Radians(lon1);
00048 lon2 = Radians(lon2);
00049 lat1 = Radians(lat1);
00050 lat2 = Radians(lat2);
00051
00052 TAN1 = tan (PI/4 + lat1/2.0);
00053 TAN2 = tan (PI/4 + lat2/2.0);
00054 TAN_A = (lon2 - lon1) / (log(TAN2) - log(TAN1));
00055 L = lon1;
00056
00057 return 1;
00058 }
00059
00060
00061
00062 double G_rhumbline_lat_from_lon (double lon)
00063 {
00064 if (parallel) return L;
00065
00066 lon = Radians(lon);
00067
00068 return Degrees (2 * atan(exp((lon-L)/TAN_A) * TAN1) - PI/2.0) ;
00069 }
00070
00071 static int adjust_lon(double *lon)
00072 {
00073 while (*lon > 180.0)
00074 *lon -= 360.0;
00075 while (*lon < -180.0)
00076 *lon += 360.0;
00077
00078 return 0;
00079 }
00080
00081 static int adjust_lat(double *lat)
00082 {
00083 if (*lat > 90.0) *lat = 90.0;
00084 if (*lat < -90.0) *lat = -90.0;
00085
00086 return 0;
00087 }