00001 #include <grass/gis.h>
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "pi.h"
00027
00028 extern double sin(), cos(), tan(), atan();
00029
00030 static double A, B;
00031 #define swap(a,b) temp=a;a=b;b=temp
00032 static int adjust_lat(double *);
00033 static int adjust_lon(double *);
00034
00035 int G_begin_geodesic_equation(double lon1,double lat1,double lon2,double lat2)
00036 {
00037 double sin21, tan1,tan2;
00038
00039 adjust_lon (&lon1);
00040 adjust_lon (&lon2);
00041 adjust_lat (&lat1);
00042 adjust_lat (&lat2);
00043 if (lon1 > lon2)
00044 {
00045 register double temp;
00046 swap(lon1,lon2);
00047 swap(lat1,lat2);
00048 }
00049 if (lon1 == lon2)
00050 {
00051 A = B = 0.0;
00052 return 0;
00053 }
00054 lon1 = Radians(lon1);
00055 lon2 = Radians(lon2);
00056 lat1 = Radians(lat1);
00057 lat2 = Radians(lat2);
00058
00059 sin21 = sin (lon2-lon1);
00060 tan1 = tan (lat1);
00061 tan2 = tan (lat2);
00062
00063 A = ( tan2 * cos(lon1) - tan1 * cos(lon2) ) /sin21;
00064 B = ( tan2 * sin(lon1) - tan1 * sin(lon2) ) /sin21;
00065
00066 return 1;
00067 }
00068
00069
00070
00071 double G_geodesic_lat_from_lon (double lon)
00072 {
00073 adjust_lon (&lon);
00074 lon = Radians(lon);
00075 return Degrees (atan(A * sin(lon) - B * cos(lon)));
00076 }
00077
00078 static int adjust_lon( double *lon)
00079 {
00080 while (*lon > 180.0)
00081 *lon -= 360.0;
00082 while (*lon < -180.0)
00083 *lon += 360.0;
00084
00085 return 0;
00086 }
00087
00088 static int adjust_lat( double *lat)
00089 {
00090 if (*lat > 90.0) *lat = 90.0;
00091 if (*lat < -90.0) *lat = -90.0;
00092
00093 return 0;
00094 }