home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / MISCELLA.PAK / GEODESY.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  5.4 KB  |  181 lines
  1.  
  2. (* :Title: Geodesy *)
  3.  
  4. (* :Author: John M. Novak *)
  5.  
  6. (* :Summary:
  7.     This package contains functions useful for or derived from geodesy, 
  8.     the science of measuring and mapping the surface of the Earth.
  9.     For instance, the package includes functions for finding
  10.     the distance between two points on the surface of the planet,
  11.     using different models.
  12. *)
  13.  
  14. (* :Context: Miscellaneous`Geodesy` *)
  15.  
  16. (* :Package Version: 1.0 *)
  17.  
  18. (* :History:
  19.     V.1.0, April 1991, by John M. Novak
  20. *)
  21.  
  22. (* :Keywords:
  23.     geodesy, geography, distance, Earth
  24. *)
  25.  
  26. (* :Sources:
  27.     Griffin, Frank, An Introduction to Spherical Trigonometry,
  28.         (Houghton Mifflin Co., 1943)
  29.     Pearson, Frederick II, Map Projections: Theory and
  30.         Applications, (CRC Press, 1990)
  31. *)
  32.  
  33. (* :Mathematica Version: 2.0 *)
  34.  
  35. BeginPackage["Miscellaneous`Geodesy`"]
  36.  
  37. SphericalDistance::usage =
  38.     "SphericalDistance[pt1,pt2] gives the distance between two
  39.     points on the earth using the spherical model of the planet.
  40.     The points are expressed as {lat,long}, where lat, long can
  41.     be in degrees, or {d,m} or {d,m,s} form.";
  42.  
  43. SpheroidalDistance::usage =
  44.     "SpheroidalDistance[pt1,pt2] gives the distance between two
  45.     points on the Earth in km, using the spheroidal model of
  46.     the planet. Coordinates are expressed as in SphericalDistance.";
  47.  
  48. ToAuthalicRadius::usage =
  49.     "ToAuthalicRadius[a,e] gives the radius of the authalic
  50.     sphere based on the spheroid with semimajor axis a and
  51.     eccentricity e."
  52.  
  53. GeodeticToAuthalic::usage =
  54.     "GeodeticToAuthalic[{lat,long},e] returns the latitude and
  55.     longitude of a point on an authalic sphere based on a
  56.     spheroid of eccentricity e, where lat and long are the
  57.     geodetic coordinates of that point on the spheroid."
  58.  
  59. SemimajorAxis::usage =
  60.     "SemimajorAxis is an option for SpheroidalDistance, that
  61.     sets the semimajor axis for the spheroid defining the
  62.     planet, in km.  Default is for the Earth, from the WGS-84
  63.     standard.";
  64.  
  65. Eccentricity::usage =
  66.     "Eccentricity is an option for SpheroidalDistance, which
  67.     sets the eccentricity of the ellipsoid defining the planet.  
  68.     The default is for the Earth, from the WGS-84 standard.";
  69.  
  70. Radius::usage =
  71.     "Radius is an option for SphericalDistance, which sets the
  72.     radius of the sphere defining the planet, in km.  The default
  73.     is for the Earth, for the authalic sphere based on the
  74.     WGS-84 standard spheroid.";
  75.  
  76. ToDegrees::usage =
  77.     "ToDegrees[{d,m}] or ToDegrees[{d,m,s}] converts degree-minute-
  78.     second forms of coordinates to degrees.  The coordinates are
  79.     adjusted to stay within [-180,180] degrees.  Note that the
  80.     sign of d is enforced on m and s; so {-34,3,2} means 34 deg,
  81.     3 min, and 2 sec W (or S), as does {-34,-3,-2}."
  82.  
  83. ToDMS::usage =
  84.     "ToDMS[deg] takes a coordinate in degrees and converts it
  85.     to a degree-minute-second format, to the nearest second.
  86.     The coordinates are within [-180,180] degrees."
  87.  
  88. Begin["`Private`"]
  89.  
  90. (* Ground distance between points on the surface of the planet.  *)
  91.  
  92. Options[SpheroidalDistance] =
  93.     {SemimajorAxis->6378137/1000,
  94.     Eccentricity->.081819}
  95.  
  96. SpheroidalDistance[{p1_,l1_},{p2_,l2_},opts___] :=
  97.     Module[{a,e,phi = ToDegrees[p1] Degree,
  98.             dlam = deltalambda[ToDegrees[l1] Degree,
  99.                 ToDegrees[l2] Degree]},
  100.         {a,e} = {SemimajorAxis,Eccentricity}/.{opts}/.
  101.                 Options[SpheroidalDistance];
  102.         a dlam Cos[phi]/Sqrt[1 - e^2 Sin[phi]^2]
  103.     ]/; p1 == p2
  104.  
  105. deltalambda[l1_,l2_] :=
  106.     Module[{diff},
  107.         If[((diff = Abs[l1 - l2]) > Pi)//N,
  108.             2 Pi - diff,
  109.             diff]
  110.     ]
  111.  
  112. SpheroidalDistance[{p1_,l1_},{p2_,l2_},opts___] :=
  113.     Module[{a,e,phi,phi1,phi2,lam1,lam2},
  114.         {a,e} = {SemimajorAxis,Eccentricity}/.{opts}/.
  115.                 Options[SpheroidalDistance];
  116.         {phi1,phi2,lam1,lam2} = Map[ToDegrees,{p1,p2,l1,l2}] Degree;
  117.         Abs[a NIntegrate[Evaluate[Sqrt[
  118.             (1-e^2)^2/(1-e^2 Sin[phi]^2)^3 +
  119.             (lam2-lam1)^2 (1-e^2 Sin[phi]^2)/
  120.             (((1-e^2)(Tan[phi2] - Tan[phi1]) -
  121.                 e^2(phi2 - phi1))^2 Cos[phi]^2)]],
  122.             Evaluate[{phi,phi1,phi2}]]]
  123.     ]
  124.  
  125. Options[SphericalDistance] =
  126.     {Radius->6371007/1000}
  127.  
  128. SphericalDistance[{p1_,l1_},{p2_,l2_},opts___] :=
  129.     Module[{r,phi1,phi2,lam1,lam2},
  130.         {r} = {Radius}/.{opts}/.Options[SphericalDistance];
  131.         {phi1,phi2,lam1,lam2} = Map[ToDegrees,{p1,p2,l1,l2}] Degree;
  132.         r ArcCos[Sin[phi1] Sin[phi2] +
  133.             Cos[phi1] Cos[phi2] Cos[Abs[lam1 - lam2]]]
  134.     ]
  135.  
  136. (* Conversions between models. *)
  137.  
  138. ToAuthalicRadius[a_?NumberQ,0] := a
  139.  
  140. ToAuthalicRadius[a_?NumberQ,1] := a Sqrt[1/2]
  141.  
  142. ToAuthalicRadius[a_?NumberQ,e_?NumberQ] :=
  143.     a Sqrt[(2 e + (1 - e^2) (Log[1 + e] - Log[1 - e]))/(4e)]
  144.  
  145. GeodeticToAuthalic[{lt_,ln_},0] := {ToDegrees[lt], ToDegrees[ln]}
  146.  
  147. GeodeticToAuthalic[{lt_,ln_},1] := {0,ToDegrees[ln]}
  148.  
  149. GeodeticToAuthalic[{lt_,ln_},e_?NumberQ] :=
  150.     Module[{lat = ToDegrees[lt] Degree,long = ToDegrees[ln] Degree//N,
  151.             arg},
  152.         arg = ((1 - e^2) (2 e Sin[lat] -
  153.                 (1 - e^2 Sin[lat]^2) (Log[1 - e Sin[lat]] -
  154.                 Log[1 + e Sin[lat]])))/
  155.             (Sin[lat] (1 - e^2 Sin[lat]^2) (2 e -
  156.                 (1 - e^2) (Log[1 - e] - Log[1 + e])));
  157.         {Re[ArcSin[Sin[lat] arg]],long}/Degree
  158.     ]
  159.  
  160. (* Auxilliary Functions *)
  161.  
  162. ToDegrees[deg_?NumberQ] :=
  163.     If[Abs[deg] > 180,
  164.         deg - 360 Ceiling[Quotient[deg,180]/2],
  165.         deg]
  166.  
  167. ToDegrees[{d_,m_:0,s_:0}] :=
  168.     ToDegrees[d + (Sign[d] Abs[m])/60 + (Sign[d] Abs[s])/3600]
  169.  
  170. ToDMS[deg_?NumberQ] :=
  171.     Module[{tmp,d,m,s},
  172.         tmp = ToDegrees[deg]; (* Make sure deg is in valid range *)
  173.         d = Floor[Abs[tmp]];
  174.         m = Floor[(Abs[tmp] - d) 60];
  175.         s = Round[(Abs[tmp] - d - m/60) 3600];
  176.         Sign[tmp]{d,m,s}]
  177.  
  178. End[]
  179.  
  180. EndPackage[]
  181.