geodesic.c

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

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