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 <grass/gis.h> 00022 #include "pi.h" 00023 00024 extern double sin(), sqrt(); 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 00074 G_meridional_radius_of_curvature (double lon, double a, double e2) 00075 00076 { 00077 double x; 00078 double s; 00079 00080 s = sin (Radians(lon)); 00081 x = 1 - e2 * s*s; 00082 00083 return a * (1 - e2) / (x * sqrt(x)); 00084 } 00085 00086 00087 00104 double 00105 G_transverse_radius_of_curvature (double lon, double a, double e2) 00106 00107 { 00108 double x; 00109 double s; 00110 00111 s = sin (Radians(lon)); 00112 x = 1 - e2 * s*s; 00113 00114 return a / sqrt(x); 00115 } 00116 00117 00134 double 00135 G_radius_of_conformal_tangent_sphere (double lon, double a, double e2) 00136 00137 { 00138 double x; 00139 double s; 00140 00141 s = sin (Radians(lon)); 00142 x = 1 - e2 * s*s; 00143 00144 return a * sqrt (1 - e2) / x; 00145 }