00001 /* TODO: 00002 00003 Suggestion: all "lon"s in the file "radii.c" should read as "lat" 00004 00005 Comments: 00006 on page http://www.mentorsoftwareinc.com/cc/gistips/TIPS0899.HTM 00007 down where it says "Meridional Radius of Curvature" is the exact formula 00008 out of "radii.c". 00009 Quote: "essentially, the radius of curvature, at a specific latitude ...". 00010 00011 See also http://williams.best.vwh.net/ellipsoid/node1.html which has a nice 00012 picture showning the parametric latitude and phi, the geodetic latitude. 00013 On the next page, 00014 http://williams.best.vwh.net/ellipsoid/node2.html, in equation 3, the 00015 Meridional Radius of Curvature shows up. 00016 00017 So, it looks like you are calculating the Meridional Radius of Curvature 00018 as a function of GEODETIC LATITUDE. 00019 */ 00020 00021 #include <math.h> 00022 #include <grass/gis.h> 00023 #include "pi.h" 00024 00025 00026 /**************************************************************** 00027 Various formulas for the ellipsoid. 00028 Reference: Map Projections by Peter Richardus and Ron K. Alder 00029 University of Illinois Library Call Number: 526.8 R39m 00030 Parameters are: 00031 lon = longitude of the meridian 00032 a = ellipsoid semi-major axis 00033 e2 = ellipsoid eccentricity squared 00034 00035 00036 meridional radius of curvature (p. 16) 00037 2 00038 a ( 1 - e ) 00039 M = ------------------ 00040 2 2 3/2 00041 (1 - e sin lon) 00042 00043 transverse radius of curvature (p. 16) 00044 00045 a 00046 N = ------------------ 00047 2 2 1/2 00048 (1 - e sin lon) 00049 00050 radius of the tangent sphere onto which angles are mapped 00051 conformally (p. 24) 00052 00053 R = sqrt ( N * M ) 00054 00055 ***************************************************************************/ 00056 00073 double G_meridional_radius_of_curvature(double lon, double a, double e2) 00074 { 00075 double x; 00076 double s; 00077 00078 s = sin(Radians(lon)); 00079 x = 1 - e2 * s * s; 00080 00081 return a * (1 - e2) / (x * sqrt(x)); 00082 } 00083 00084 00085 00102 double G_transverse_radius_of_curvature(double lon, double a, double e2) 00103 { 00104 double x; 00105 double s; 00106 00107 s = sin(Radians(lon)); 00108 x = 1 - e2 * s * s; 00109 00110 return a / sqrt(x); 00111 } 00112 00113 00130 double G_radius_of_conformal_tangent_sphere(double lon, double a, double e2) 00131 { 00132 double x; 00133 double s; 00134 00135 s = sin(Radians(lon)); 00136 x = 1 - e2 * s * s; 00137 00138 return a * sqrt(1 - e2) / x; 00139 }