cearth.cpp

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

   coordinate = ( distance_from_center_to_surface_of_the_ellipsoid + polar_coordinate.GetDistanceFromSurfaceInMeters() ) * cosine_of_up_down_radians * CMath::Sine( left_right_radians );
   cartesian_coordinate.SetYCoordinateInMeters( coordinate );

   coordinate = ( distance_from_center_to_surface_of_the_ellipsoid * ( 1.0 - m_EccentricitySquared ) + polar_coordinate.GetDistanceFromSurfaceInMeters() ) * sine_of_up_down_radians;
   cartesian_coordinate.SetZCoordinateInMeters( coordinate );
}

double CEarth::GetDistanceToHorizon( const CEarthCoordinate& point_1 ) const
{
   CPolarCoordinate polar_coordinate;

   Convert( point_1, polar_coordinate );

   return( GetDistanceToHorizon( polar_coordinate ) );
}

double CEarth::GetDistanceToHorizon( const CPolarCoordinate& point_1 ) const
{
   double distance_to_horizon = 0.0;

   // d = ::sqrt( feet ) *  1.144 for nmi
   // optical horizon is 1.317 * sqrt( h );
   // d= ::sqrt( 17 * height_in_meters ); d is in meters

   distance_to_horizon = CMath::SquareRoot( 17.0 * point_1.GetDistanceFromSurfaceInMeters() );

   return( distance_to_horizon );
}

double CEarth::GetEquatorialRadiusInMeters( void ) const
{
   return( m_EquatorialRadiusInMeters );
}

double CEarth::GetPolarRadiusInMeters( void ) const
{
   return( m_PolarRadiusInMeters );
}

double CEarth::GetLineOfSightDistanceFromCourse( const CEarthCoordinate& current_location, const CEarthCoordinate& point_a, const CEarthCoordinate& point_b ) const
{
   // This function tells you how far off course you are from a straight line between
   // point_a and point_b.

   /*
   ** References:
   ** I got the formula from:
   ** Engineering Mathematics Handbook
   ** Jan J. Tuma, Ph.D.
   ** McGraw-Hill Book Company
   ** 1970
   ** Library of Congress Catalog Number 78-101174
   ** page 19, (a) Oblique triangle
   **
   ** Teach Yourself Trigonometry
   ** P. Abbott, B.A.
   ** English Universities Press Ltd.
   ** 102 Newgate Street
   ** London, E.C.I
   ** Originally published 1940
   ** I used the 1964 printing.
   ** Page 22, Figure 12 calls this "the altitude from the vertex A"
   */ 

   double distance_from_current_location_to_point_a = 0.0;
   double distance_from_current_location_to_point_b = 0.0;
   double distance_from_point_a_to_point_b          = 0.0;

   distance_from_current_location_to_point_a = GetLineOfSightDistance( current_location, point_a );
   distance_from_current_location_to_point_b = GetLineOfSightDistance( current_location, point_b );
   distance_from_point_a_to_point_b          = GetLineOfSightDistance( point_a,          point_b );

   double p = 0.0;

   p  = distance_from_current_location_to_point_a;
   p += distance_from_current_location_to_point_b;
   p += distance_from_point_a_to_point_b;
   p /= 2.0;

   double temp_double = 0.0;

   temp_double  = p;
   temp_double *= (double) ( p - distance_from_current_location_to_point_a );
   temp_double *= (double) ( p - distance_from_current_location_to_point_b );
   temp_double *= (double) ( p - distance_from_point_a_to_point_b );

   double area = 0.0;

   area = CMath::SquareRoot( temp_double );

   double distance_from_course = 0.0;

   // The altitude from the vertex A is two times the area of the triangle divided by the baseline

   distance_from_course = ( 2.0 * area ) / distance_from_point_a_to_point_b;

   return( distance_from_course );
}

double CEarth::GetLineOfSightDistance( const CEarthCoordinate& point_1, const CEarthCoordinate& point_2 ) const
{
   // This function implements the Pythagoras method of computing the distance
   // between two points.
   // This is a line-of-sight algorithm. It does not take into acccount the 
   // curvature of the Earth. It is not a distance on the surface algorithm.
   // If you had a laser and connected the two points, this algorithm tells
   // you how long the laser beam is. 

   double distance     = 0.0;
   double x_coordinate = 0.0;
   double y_coordinate = 0.0;
   double z_coordinate = 0.0;

   x_coordinate = point_1.GetXCoordinateInMeters() - point_2.GetXCoordinateInMeters();
   y_coordinate = point_1.GetYCoordinateInMeters() - point_2.GetYCoordinateInMeters();
   z_coordinate = point_1.GetZCoordinateInMeters() - point_2.GetZCoordinateInMeters();

   // Square the coordinates
   x_coordinate *= x_coordinate;
   y_coordinate *= y_coordinate;
   z_coordinate *= z_coordinate;

   distance = CMath::SquareRoot( x_coordinate + y_coordinate + z_coordinate );

   return( distance );
}

double CEarth::GetLineOfSightDistance( const CEarthCoordinate& point_1, const CPolarCoordinate& point_2 ) const
{
   CEarthCoordinate earth_center_earth_fixed_point_2;

   Convert( point_2, earth_center_earth_fixed_point_2 );

   return( GetLineOfSightDistance( point_1, earth_center_earth_fixed_point_2 ) );
}

double CEarth::GetLineOfSightDistance( const CPolarCoordinate& point_1, const CEarthCoordinate& point_2 ) const
{
   CEarthCoordinate earth_center_earth_fixed_point_1;

   Convert( point_1, earth_center_earth_fixed_point_1 );

   return( GetLineOfSightDistance( earth_center_earth_fixed_point_1, point_2 ) );
}

double CEarth::GetLineOfSightDistance( const CPolarCoordinate& point_1, const CPolarCoordinate& point_2 ) const
{
   CEarthCoordinate earth_center_earth_fixed_point_1;
   CEarthCoordinate earth_center_earth_fixed_point_2;

   Convert( point_1, earth_center_earth_fixed_point_1 );
   Convert( point_2, earth_center_earth_fixed_point_2 );

   return( GetLineOfSightDistance( earth_center_earth_fixed_point_1, earth_center_earth_fixed_point_2 ) );
}

double CEarth::GetSurfaceDistance( const CEarthCoordinate& point_1, const CEarthCoordinate& point_2, double * heading_from_point_1_to_point_2_p, double * heading_from_point_2_to_point_1_p ) const
{
   CPolarCoordinate polar_point_1;
   CPolarCoordinate polar_point_2;

   Convert( point_1, polar_point_1 );
   Convert( point_2, polar_point_2 );

   return( GetSurfaceDistance( polar_point_1, polar_point_2, heading_from_point_1_to_point_2_p, heading_from_point_2_to_point_1_p ) );
}

double CEarth::GetSurfaceDistance( const CEarthCoordinate& point_1, const CPolarCoordinate& point_2, double * heading_from_point_1_to_point_2_p, double * heading_from_point_2_to_point_1_p ) const
{
   CPolarCoordinate polar_point_1;

   Convert( point_1, polar_point_1 );

   return( GetSurfaceDistance( polar_point_1, point_2, heading_from_point_1_to_point_2_p, heading_from_point_2_to_point_1_p ) );
}

double CEarth::GetSurfaceDistance( const CPolarCoordinate& point_1, const CEarthCoordinate& point_2, double * heading_from_point_1_to_point_2_p, double * heading_from_point_2_to_point_1_p ) const
{
   CPolarCoordinate polar_point_2;

   Convert( point_2, polar_point_2 );

   return( GetSurfaceDistance( point_1, polar_point_2, heading_from_point_1_to_point_2_p, heading_from_point_2_to_point_1_p ) );
}

double CEarth::GetSurfaceDistance( const CPolarCoordinate& point_1, const CPolarCoordinate& point_2, double * heading_from_point_1_to_point_2_p, double * heading_from_point_2_to_point_1_p ) const
{
   // This is a translation of the Fortran routine INVER1 found in the
   // INVERS3D program at:
   // ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/invers3d.for

   // The ton of variables used...

   double c           = 0.0;
   double c_value_1   = 0.0;
   double c_value_2   = 0.0;
   double c2a         = 0.0;
   double cosine_of_x = 0.0;
   double cy          = 0.0;
   double cz          = 0.0;
   double d           = 0.0;
   double e           = 0.0;
   double r_value     = 0.0;
   double s           = 0.0;
   double s_value_1   = 0.0;
   double sa          = 0.0;
   double sine_of_x   = 0.0;
   double sy          = 0.0;
   double tangent_1   = 0.0;
   double tangent_2   = 0.0;
   double x           = 0.0;
   double y           = 0.0;

   int loop_count = 0;

   double heading_from_point_1_to_point_2 = 0.0;
   double heading_from_point_2_to_point_1 = 0.0;

   // UpDown    == Latitude
   // LeftRight == Longitude

   double point_1_latitude_in_radians  = CMath::ConvertDegreesToRadians( point_1.GetUpDownAngleInDegrees()    );
   double point_1_longitude_in_radians = CMath::ConvertDegreesToRadians( point_1.GetLeftRightAngleInDegrees() );
   double point_2_latitude_in_radians  = CMath::ConvertDegreesToRadians( point_2.GetUpDownAngleInDegrees()    );
   double point_2_longitude_in_radians = CMath::ConvertDegreesToRadians( point_2.GetLeftRightAngleInDegrees() );

   r_value = 1.0 - m_Flattening;
   tangent_1 = ( r_value * CMath::Sine( point_1_latitude_in_radians ) ) / CMath::Cosine( point_1_latitude_in_radians );
   tangent_2 = ( r_value * CMath::Sine( point_2_latitude_in_radians ) ) / CMath::Cosine( point_2_latitude_in_radians );
   c_value_1 = 1.0 / CMath::SquareRoot( ( tangent_1 * tangent_1 ) + 1.0 );
   s_value_1 = c_value_1 * tangent_1;
   c_value_2 = 1.0 / CMath::SquareRoot( ( tangent_2 * tangent_2 ) + 1.0 );
   s = c_value_1 * c_value_2;

   heading_from_point_2_to_point_1 = s * tangent_2; // backward_azimuth
   heading_from_point_1_to_point_2 = heading_from_point_2_to_point_1 * tangent_1;

   x = point_2_longitude_in_radians - point_1_longitude_in_radians;

   bool exit_loop = false;

   while( exit_loop != true )
   {
      sine_of_x   = CMath::Sine( x );
      cosine_of_x = CMath::Cosine( x );
      tangent_1 = c_value_2 * sine_of_x;
      tangent_2 = heading_from_point_2_to_point_1 - ( s_value_1 * c_value_2 * cosine_of_x );
      sy = CMath::SquareRoot( ( tangent_1 * tangent_1 ) + ( tangent_2 * tangent_2 ) );
      cy = ( s * cosine_of_x ) + heading_from_point_1_to_point_2;
      y = CMath::ArcTangentOfYOverX( sy, cy );

      // Thanks to John Werner ([e-mail]) for
      // finding a bug where sy could be zero. Here's his fix:

      if ( ( s * sine_of_x ) == 0.0 && ( sy == 0.0 ) )
      {
         sa = 1.0;
      }
      else
      {
         sa = ( s * sine_of_x ) / sy;
      }

      c2a = ( (-sa) * sa ) + 1.0;
      cz = heading_from_point_1_to_point_2 + heading_from_point_1_to_point_2;

      if ( c2a > 0.0 )
      {
         cz = ( (-cz) / c2a ) + cy;
      }

      e = ( cz * cz * 2.0 ) - 1.0;
      c = ( ( ( ( ( -3.0 * c2a ) + 4.0 ) * m_Flattening ) + 4.0 ) * c2a * m_Flattening ) / 16.0;
      d = x;
      x = ( ( ( ( e * cy * c ) + cz ) * sy * c ) + y ) * sa;
      x = ( ( 1.0 - c ) * x * m_Flattening ) + point_2_longitude_in_radians - point_1_longitude_in_radians;

      if ( CMath::AbsoluteValue( d - x ) > 0.00000000000000000000005 )
      {
         exit_loop = false;
      }
      else
      {
         exit_loop = true;
      }
   }

   heading_from_point_1_to_point_2 = CMath::ArcTangentOfYOverX( tangent_1, tangent_2 );

   double temp_degrees         = 0.0;
   double temp_minutes         = 0.0;
   double temp_seconds         = 0.0;
   double temp_decimal_degrees = 0.0;

   temp_decimal_degrees = CMath::ConvertRadiansToDegrees( heading_from_point_1_to_point_2 );

   if ( temp_decimal_degrees < 0.0 )
   {
      temp_decimal_degrees += 360.0;
   }

   if ( heading_from_point_1_to_point_2_p != NULL )
   {
      // The user passed us a pointer, don't trust it.
      // If you are using Visual C++ on Windows NT, the following
      // try/catch block will ensure you won't blow up when random
      // pointers are passed to you. If you are on a legacy operating
      // system like Unix, you are screwed.

      try
      {
         *heading_from_point_1_to_point_2_p = temp_decimal_degrees;
      }
      catch( ... )
      {
         // Do Nothing
      }
   }

   heading_from_point_2_to_point_1 = CMath::ArcTangentOfYOverX( c_value_1 * sine_of_x, ( (heading_from_point_2_to_point_1 * cosine_of_x ) - ( s_value_1 * c_value_2 ) ) ) + CMath::Pi();

   temp_decimal_degrees = CMath::ConvertRadiansToDegrees( heading_from_point_2_to_point_1 );

   if ( temp_decimal_degrees < 0 )
   {
      temp_decimal_degrees += 360.0;
   }

   if ( heading_from_point_2_to_point_1_p != NULL )
   {
      // The user passed us a pointer, don't trust it.
      // If you are using Visual C++ on Windows NT, the following
      // try/catch block will ensure you won't blow up when random
      // pointers are passed to you. If you are on a legacy operating
      // system like Unix, you are screwed.

      try
      {
         *heading_from_point_2_to_point_1_p = temp_decimal_degrees;
      }
      catch( ... )
      {
         // Do Nothing
      }
   }

   x = CMath::SquareRoot( ( ( ( 1.0 / r_value / r_value ) - 1 ) * c2a ) + 1.0 ) + 1.0;
   x = ( x - 2.0 ) / x;
   c = 1.0 - x;
   c = ( ( ( x * x ) / 4.0 ) + 1.0 ) / c;
   d = ( ( 0.375 * ( x * x ) ) - 1.0 ) * x;

   // 1998-09-01
   // Thanks go to Gerard Murphy ([e-mail]) for finding a typo here.

   x = e * cy;

   s = ( 1.0 - e ) - e;

   double term_1 = 0.0;
   double term_2 = 0.0;
   double term_3 = 0.0;
   double term_4 = 0.0;
   double term_5 = 0.0;

   term_1 = ( sy * sy * 4.0 ) - 3.0;
   term_2 = ( ( s * cz * d ) / 6.0 ) - x;
   term_3 = term_1 * term_2;
   term_4 = ( ( term_3 * d ) / 4.0 ) + cz;
   term_5 = ( term_4 * sy * d ) + y;

   s = term_5 * c * m_EquatorialRadiusInMeters * r_value;

   return( s );
}

void CEarth::m_ComputeEccentricitySquared( void )
{
   if ( m_Flattening == 0.0 )
   {
      m_EccentricitySquared = 0.0;
      return;
   }

   m_EccentricitySquared = ( 2.0 * m_Flattening ) - ( m_Flattening * m_Flattening );
   //printf( "Eccentricity Squared = %.23lf\n", m_EccentricitySquared );
}