rhumbline.c

Go to the documentation of this file.
00001 /*
00002  * This code is preliminary. I don't know if it is even
00003  * correct.
00004  */
00005 
00006 /*
00007  * From "Map Projections" by Peter Richardus and Ron K. Alder, 1972
00008  * (526.8 R39m in Map & Geography Library)
00009  * page  20,21, formulas 2.21,2.22
00010  *
00011  * Formula is the equation of a rhumbline from (lat1,lon1) to (lat2,lon2)
00012  * Input is lon, output is lat (all in degrees)
00013  *
00014  * Note formula only works if 0 < abs(lon2-lon1) < 180
00015  * If lon1 == lon2 then rhumbline is the merdian lon1 
00016  * (and the formula will fail)
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;   /* a lie */
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 /* only works if lon1 < lon < lon2 */
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 }

Generated on Fri Nov 21 11:02:18 2008 for GRASS by  doxygen 1.5.1