GRASS Programmer's Manual
6.4.1(2011)
|
00001 #include <math.h> 00002 #include <grass/gis.h> 00003 #include "pi.h" 00004 00005 00006 /* 00007 * This code is preliminary. I don't know if it is even 00008 * correct. 00009 */ 00010 00011 /* 00012 * From "Map Projections" by Peter Richardus and Ron K. Alder, 1972 00013 * (526.8 R39m in Map & Geography Library) 00014 * page 19, formula 2.17 00015 * 00016 * Formula is the equation of a geodesic from (lat1,lon1) to (lat2,lon2) 00017 * Input is lon, output is lat (all in degrees) 00018 * 00019 * Note formula only works if 0 < abs(lon2-lon1) < 180 00020 * If lon1 == lon2 then geodesic is the merdian lon1 00021 * (and the formula will fail) 00022 * if lon2-lon1=180 then the geodesic is either meridian lon1 or lon2 00023 */ 00024 00025 /* TODO: 00026 * 00027 * integrate code from raster/r.surf.idw/ll.c 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 /* only works if lon1 < lon < lon2 */ 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 }