cearth.cpp

计算地球或椭圆体中两点距离的程序

void CEarth::m_ComputeFlattening( void )
{
   if ( m_EquatorialRadiusInMeters == 0.0 || m_PolarRadiusInMeters == 0.0 )
   {
      return;
   }

   m_Flattening = CMath::AbsoluteValue( m_EquatorialRadiusInMeters - m_PolarRadiusInMeters ) / m_EquatorialRadiusInMeters;
   //printf( "Flattening = %.23lf\n", m_Flattening );
}

void CEarth::m_Initialize( void )
{
   m_EllipsoidID              = 0;
   m_PolarRadiusInMeters      = 0.0;
   m_EquatorialRadiusInMeters = 0.0;
   m_Flattening               = 0.0;
   m_EccentricitySquared      = 0.0;
}

void CEarth::SetEllipsoid( int ellipsoid_identifier )
{
   m_EllipsoidID = ellipsoid_identifier;

   switch( ellipsoid_identifier )
   {
      case Perfect_Sphere:

         m_EquatorialRadiusInMeters = 6378137.0;
         m_PolarRadiusInMeters      = 6378137.0;

         break;

      case Airy:

         m_EquatorialRadiusInMeters = 6377563.396;
         m_PolarRadiusInMeters      = 6356256.909237;

         break;

      case Austrailian_National:

         m_EquatorialRadiusInMeters = 6378160.0;
         m_PolarRadiusInMeters      = 6356774.719195;

         break;

      case Bessell_1841:

         m_EquatorialRadiusInMeters = 6377397.155;
         m_PolarRadiusInMeters      = 6356078.962818;

         break;

      case Bessel_1841_Nambia:

         m_EquatorialRadiusInMeters = 6377483.865;
         m_PolarRadiusInMeters      = 6356165.382966;

         break;

      case Clarke_1866:

         m_EquatorialRadiusInMeters = 6378206.4;
         m_PolarRadiusInMeters      = 6356583.799999;

         break;

      case Clarke_1880:

         m_EquatorialRadiusInMeters = 6378249.145;
         m_PolarRadiusInMeters      = 6356514.86955;

         break;

      case Everest:

         m_EquatorialRadiusInMeters = 6377276.345;
         m_PolarRadiusInMeters      = 6356075.41314;

         break;

      case Fischer_1960_Mercury:

         m_EquatorialRadiusInMeters = 6378166.0;
         m_PolarRadiusInMeters      = 6356784.283607;

         break;

      case Fischer_1968:

         m_EquatorialRadiusInMeters = 6378150.0;
         m_PolarRadiusInMeters      = 6356768.337244;

         break;

      case GRS_1967:

         m_EquatorialRadiusInMeters = 6378160.0;
         m_PolarRadiusInMeters      = 6356774.516091;

         break;

      case GRS_1980:

         m_EquatorialRadiusInMeters = 6378137.0;
         m_PolarRadiusInMeters      = 6356752.31414;

         break;

      case Helmert_1906:

         m_EquatorialRadiusInMeters = 6378200.0;
         m_PolarRadiusInMeters      = 6356818.169628;

         break;

      case Hough:

         m_EquatorialRadiusInMeters = 6378270.0;
         m_PolarRadiusInMeters      = 6356794.343434;

         break;

      case International:

         m_EquatorialRadiusInMeters = 6378388.0;
         m_PolarRadiusInMeters      = 6356911.946128;

         break;

      case Krassovsky:

         m_EquatorialRadiusInMeters = 6378245.0;
         m_PolarRadiusInMeters      = 6356863.018773;

         break;

      case Modified_Airy:

         m_EquatorialRadiusInMeters = 6377340.189;
         m_PolarRadiusInMeters      = 6356034.447939;

         break;

      case Modified_Everest:

         m_EquatorialRadiusInMeters = 6377304.063;
         m_PolarRadiusInMeters      = 6356103.038993;

         break;

      case Modified_Fischer_1960:

         m_EquatorialRadiusInMeters = 6378155.0;
         m_PolarRadiusInMeters      = 6356773.320483;

         break;

      case South_American_1969:

         m_EquatorialRadiusInMeters = 6378160.0;
         m_PolarRadiusInMeters      = 6356774.719195;

         break;

      case Topex_Poseidon_Pathfinder_ITRF: // Source is http://neptune.gsfc.nasa.gov/~krachlin/corr/refframe.html

         m_EquatorialRadiusInMeters = 6378136.3;
         m_PolarRadiusInMeters      = 6356751.6005629376;

         break;

      case WGS_60:

         m_EquatorialRadiusInMeters = 6378165.0;
         m_PolarRadiusInMeters      = 6356783.286959;

         break;

      case WGS_66:

         m_EquatorialRadiusInMeters = 6378145.0;
         m_PolarRadiusInMeters      = 6356759.769489;

         break;

      case WGS_72:

         m_EquatorialRadiusInMeters = 6378135.0;
         m_PolarRadiusInMeters      = 6356750.520016;

         break;

      case WGS_84:

         // Computed polar radius from the flattening value specified at
         // http://acro.harvard.edu/SSA/BGA/wg84figs.html
         // because it had the most digits after the decimal point.

         m_EquatorialRadiusInMeters = 6378137.0;
         m_PolarRadiusInMeters      = 6356752.3142451793;

         break;

      case Unknown:
      default:

         m_EllipsoidID = Unknown;

         m_Initialize();
         return;
   }

   m_ComputeFlattening();
   m_ComputeEccentricitySquared();
}

void CEarth::SetEllipsoidByRadii( double equatorial_radius, double polar_radius )
{
   m_EquatorialRadiusInMeters = equatorial_radius;
   m_PolarRadiusInMeters      = polar_radius;
   m_EllipsoidID              = Custom;

   m_ComputeFlattening();
   m_ComputeEccentricitySquared();
}

void CEarth::SetEllipsoidByEquatorialRadiusAndFlattening( double equatorial_radius, double flattening )
{
   m_EquatorialRadiusInMeters = equatorial_radius;
   m_Flattening               = flattening;
   m_EllipsoidID              = Custom;

   // We must compute the polar radius

   double temp_double = m_Flattening * m_EquatorialRadiusInMeters;

   m_PolarRadiusInMeters = m_EquatorialRadiusInMeters - temp_double;

   m_ComputeEccentricitySquared();
}
#if 0
<WFC_DOCUMENTATION>
<HTML>
<HEAD><TITLE>GFC - CEarth</TITLE></HEAD>
<BODY>
<H1>CEarth</H1>
$Revision: 11 $
<HR>
<H2>Description</H2>
This class encapsulates the Earth. It holds the data necessary to perform the
calculations of distance and direction. The Earth is not a perfect sphere.
It is an ellipsoid (flattened at the top and bottom). All angles are expressed
in degrees and all distances are expressed in meters.
<H2>Methods</H2>
<DL COMPACT>
<DT><PRE>void <B>AddLineOfSightDistanceAndDirectionToCoordinate</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_1, double distance, double direction, <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_2, double height_above_surface_of_point_2 = 0.0 )</PRE><DD>
If you were to shine a laser from <CODE>point_1</CODE> pointing towards <CODE>direction</CODE>
for <CODE>distance</CODE> meters, this function will tell you what location you would
be at. It will also let you specify a point <CODE>height_above_surface_of_point_2</CODE>
meters above the surface. This method does not take into account the curvature of
the Earth.
<DT><PRE>void <B>AddSurfaceDistanceAndDirectionToCoordinate</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_1, double distance, double direction, <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_2 )
void <B>AddSurfaceDistanceAndDirectionToCoordinate</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_1, double distance, double direction, <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_2 )
void <B>AddSurfaceDistanceAndDirectionToCoordinate</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_1, double distance, double direction, <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_2 )
void <B>AddSurfaceDistanceAndDirectionToCoordinate</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_1, double distance, double direction, <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_2 )</PRE><DD>
This allows you to add a distance over the surface of the Earth to a location and
get a new location. It answers the question &quot;If I head out in this direction
for that amout of meters, where will I be?&quot;
<DT><PRE>void <B>Convert</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& cartesian_coordinate, <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& polar_coordinate ) const
void <B>Convert</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& polar_coordinate, <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& cartesian_coordinate ) const</PRE><DD>
This method allows you to convert from polar and cartestian coordinates.
<DT><PRE>double <B>GetDistanceToHorizon</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_1 ) const
double <B>GetDistanceToHorizon</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_1 ) const</PRE><DD>
This tells you how far (in meters) from the horizon <CODE>point_1</CODE> is.
<DT><PRE>double <B>GetEquatorialRadiusInMeters</B>( void ) const</PRE><DD>
This tells you what the equatorial radius is in meters for the selected ellipsoid.
<DT><PRE>double <B>GetPolarRadiusInMeters</B>( void ) const</PRE><DD>
This tells you what the polar radius is in meters for the selected ellipsoid.
<DT><PRE>double <B>GetLineOfSightDistanceFromCourse</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& current_location, const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_a, const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_b ) const</PRE><DD>
Draw a line from <CODE>point_a</CODE> to <CODE>point_b</CODE>. This function will tell
you how far <CODE>current_location</CODE> is from that line.
<DT><PRE>double <B>GetLineOfSightDistance</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_1, const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_2 ) const
double <B>GetLineOfSightDistance</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_1, const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_2 ) const
double <B>GetLineOfSightDistance</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_1, const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_2 ) const
double <B>GetLineOfSightDistance</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_1, const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_2 ) const</PRE><DD>
This will tell you how many meters it is between two points. It answers the question, &quot;If I
pointed a laser from <CODE>point_1</CODE> to <CODE>point_2</CODE>, how far would the
laser beam travel?&quot;
<DT><PRE>double <B>GetSurfaceDistance</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_1, const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_2, double * heading_from_point_1_to_point_2 = 0, double * heading_from_point_2_to_point_1 = 0 ) const
double <B>GetSurfaceDistance</B>( const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_1, const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_2, double * heading_from_point_1_to_point_2 = 0, double * heading_from_point_2_to_point_1 = 0 ) const
double <B>GetSurfaceDistance</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_1, const <A HREF="CEarthCoordinate.htm">CEarthCoordinate</A>& point_2, double * heading_from_point_1_to_point_2 = 0, double * heading_from_point_2_to_point_1 = 0 ) const
double <B>GetSurfaceDistance</B>( const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_1, const <A HREF="CPolarCoordinate.htm">CPolarCoordinate</A>& point_2, double * heading_from_point_1_to_point_2 = 0, double * heading_from_point_2_to_point_1 = 0 ) const</PRE><DD>
This will tell you how many meters it is between two points. It answers the question, &quot;If I
were to walk from <CODE>point_1</CODE> to <CODE>point_2</CODE>, how far would I walk?&quot;
<DT><PRE>void <B>SetEllipsoid</B>( int ellipsoid )</PRE><DD>
This allows you to set the ellipsoid used by <B>CEarth</B> in its calculations. The
default is <CODE>WGS84</CODE> which is generally accepted as being the closest
approximation of the Earth's ellipsoid. The <CODE>ellipsoid</CODE> parameter
may be one of the following:
<UL>
<LI>Perfect_Sphere
<LI>Airy,
<LI>Austrailian_National,
<LI>Bessell_1841,
<LI>Bessel_1841_Nambia,
<LI>Clarke_1866,
<LI>Clarke_1880,
<LI>Everest,
<LI>Fischer_1960_Mercury,
<LI>Fischer_1968,
<LI>GRS_1967,
<LI>GRS_1980,
<LI>Helmert_1906,
<LI>Hough,
<LI>International,
<LI>Krassovsky,
<LI>Modified_Airy,
<LI>Modified_Everest,
<LI>Modified_Fischer_1960,
<LI>South_American_1969,
<LI>Topex_Poseidon_Pathfinder_ITRF,
<LI>WGS_60,
<LI>WGS_66,
<LI>WGS_72,
<LI>WGS_84
<LI>NAD_27
<LI>Tokyo
</UL>
<DT><PRE>void <B>SetEllipsoidByRadii</B>( double equatorial_radius, double polar_radius )</PRE><DD>
This let's you use your own (custom) values to describe the ellipsoid of the Earth.
<DT><PRE>void <B>SetEllipsoidByEquatorialRadiusAndFlattening</B>( double equatorial_radius, double flattening )</PRE><DD>
This let's you use your own (custom) values to describe the ellipsoid of the Earth.
</DL>
<H2>Example</H2>
<PRE><CODE>#include &lt;stdio.h&gt;
#include &lt;GFC.h&gt;
#pragma hdrstop

void main( void )
{
   // Let's figure out how far it is from here to there

   CPolarCoordinate here;
   CPolarCoordinate there;

   // Convert from Latitude/Longitude to coordinates our system understands

   // here is 39 degrees 12.152 minutes North Latitude, 76 degrees 46.795 minutes West Longitude
   here.SetUpDownAngleInDegrees(     CMath::ConvertDegreesMinutesSecondsCoordinateToDecimalDegrees(  39.0, 12.152, 0.0 ) );
   here.SetLeftRightAngleInDegrees(  CMath::ConvertDegreesMinutesSecondsCoordinateToDecimalDegrees( -76.0, 46.795, 0.0 ) );

   // there is 12 degrees 8.535 minutes North Latitude, 68 degrees 16.547 West Longitude
   there.SetUpDownAngleInDegrees(    CMath::ConvertDegreesMinutesSecondsCoordinateToDecimalDegrees(  12.0,  8.535, 0.0 ) );
   there.SetLeftRightAngleInDegrees( CMath::ConvertDegreesMinutesSecondsCoordinateToDecimalDegrees( -68.0, 16.547, 0.0 ) );

   <B>CEarth</B> earth; // We are talking about the earth...

   double distance_in_meters         = 0.0;
   double heading_from_here_to_there = 0.0;
   double heading_from_there_to_here = 0.0;

   distance_in_meters = earth.GetSurfaceDistance( here, there, &heading_from_here_to_there, &heading_from_there_to_here );

   printf( &quot;Distance between here and there: %.23lf meters\nHeading from here to there:      %.19lf degrees\nHeading from there to here:      %.19lf degrees\n&quot;,
           distance_in_meters,
           heading_from_here_to_there,
           heading_from_there_to_here );

   double degrees = 0.0;
   double minutes = 0.0;
   double seconds = 0.0;

   CMath::ConvertDecimalDegreesToDegreesMinutesSeconds( heading_from_here_to_there, degrees, minutes, seconds );
   printf( &quot;Heading %lf degrees, %lf minutes, %lf seconds\n&quot;, degrees, minutes, seconds );

   CMath::ConvertDecimalDegreesToDegreesMinutesSeconds( heading_from_there_to_here, degrees, minutes, seconds );
   printf( &quot;Heading %lf degrees, %lf minutes, %lf seconds\n&quot;, degrees, minutes, seconds );
}</CODE></PRE>
<I>Copyright, 1998, Samuel R. Blackburn</I><BR>
$Workfile: CEarth.cpp $<BR>
$Modtime: 9/01/98 9:56p $
</BODY>
</HTML>
</WFC_DOCUMENTATION>
ToolTipFormatLine=CEarth=m_EllipsoidID=<m_EllipsoidID>
#endif