• Main Page
  • Related Pages
  • Data Structures
  • Files
  • File List
  • Globals

geodesic.c

Go to the documentation of this file.
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 }

Generated on Thu Dec 9 2010 20:46:05 for GRASS Programmer's Manual by  doxygen 1.7.2