www.pudn.com > gfc.tar > CEarth.cpp
#include#pragma hdrstop #include /* ** Author: Samuel R. Blackburn ** Internet: sblackbu@erols.com ** ** You can use it any way you like as long as you don't try to sell it. ** ** Any attempt to sell GFC in source code form must have the permission ** of the original author. You can produce commercial executables with ** GFC but you can't sell GFC. ** ** Copyright, 1998, Samuel R. Blackburn ** ** $Workfile: CEarth.cpp $ ** $Revision: 11 $ ** $Modtime: 9/01/98 9:56p $ */ CEarth::CEarth( int ellipsoid_identifier ) { m_Initialize(); SetEllipsoid( ellipsoid_identifier ); } CEarth::~CEarth() { m_Initialize(); } void CEarth::AddLineOfSightDistanceAndDirectionToCoordinate( const CPolarCoordinate& point_1, double distance, double direction, CPolarCoordinate& point_2, double height_above_surface_of_point_2 ) { // The method used here is to convert the straight (line-of-sight) distance to // a surface distance and then find out the position using the surface distance. // This is a translation of the MMDIST routine found in the FORWRD3D program at // ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/forwrd3d.for double c = 0.0; double cosine_of_latitude_squared = 0.0; double cosine_of_direction_squared = 0.0; double cosine_of_point_1_latitude = 0.0; double difference_in_height = 0.0; double direction_in_radians = 0.0; double distance_adjusted_for_differences_in_height = 0.0; double height_above_surface_of_point_1 = 0.0; double n = 0.0; double point_1_latitude_in_radians = 0.0; double polar_eccentricity_squared = 0.0; double polar_flattening = 0.0; double r = 0.0; double sine_of_point_1_latitude = 0.0; double surface_distance = 0.0; double term_1 = 0.0; double term_2 = 0.0; double term_3 = 0.0; double two_r = 0.0; double x = 0.0; // Many thanks to Peter Dana (pdana@mail.utexas.edu) for educating me // on the finer points of Geodesy, one of which was how to compute // "second eccentricity squared" polar_eccentricity_squared = ( ( m_EquatorialRadiusInMeters * m_EquatorialRadiusInMeters ) - ( m_PolarRadiusInMeters * m_PolarRadiusInMeters ) ) / ( m_PolarRadiusInMeters * m_PolarRadiusInMeters ); //printf( "polar_eccentricity_squared is %.23lf\n", polar_eccentricity_squared ); point_1_latitude_in_radians = CMath::ConvertDegreesToRadians( point_1.GetUpDownAngleInDegrees() ); direction_in_radians = CMath::ConvertDegreesToRadians( direction ); cosine_of_point_1_latitude = CMath::Cosine( point_1_latitude_in_radians ); cosine_of_latitude_squared = cosine_of_point_1_latitude * cosine_of_point_1_latitude; cosine_of_direction_squared = CMath::Cosine( direction_in_radians ) * CMath::Cosine( direction_in_radians ); c = ( m_EquatorialRadiusInMeters * m_EquatorialRadiusInMeters ) / m_PolarRadiusInMeters; n = c / CMath::SquareRoot( 1.0 + ( polar_eccentricity_squared * cosine_of_latitude_squared ) ); r = n / ( 1.0 + ( polar_eccentricity_squared * cosine_of_latitude_squared * cosine_of_direction_squared ) ); height_above_surface_of_point_1 = point_1.GetDistanceFromSurfaceInMeters(); difference_in_height = height_above_surface_of_point_2 - height_above_surface_of_point_1; term_1 = ( distance * distance ) - ( difference_in_height * difference_in_height ); term_2 = 1.0 + ( height_above_surface_of_point_1 / r ); term_3 = 1.0 + ( height_above_surface_of_point_2 / r ); distance_adjusted_for_differences_in_height = CMath::SquareRoot( term_1 / ( term_2 * term_3 ) ); // printf( "distance_adjusted_for_differences_in_height is %.11lf\n", distance_adjusted_for_differences_in_height ); two_r = 2.0 * r; surface_distance = two_r * CMath::ArcSine( distance_adjusted_for_differences_in_height / two_r ); // printf( "surface_distance is %.11lf\n", surface_distance ); AddSurfaceDistanceAndDirectionToCoordinate( point_1, surface_distance, direction, point_2 ); } void CEarth::AddSurfaceDistanceAndDirectionToCoordinate( const CEarthCoordinate& point_1, double distance, double direction, CPolarCoordinate& point_2 ) { CPolarCoordinate polar_point_1; Convert( point_1, polar_point_1 ); AddSurfaceDistanceAndDirectionToCoordinate( polar_point_1, distance, direction, point_2 ); } void CEarth::AddSurfaceDistanceAndDirectionToCoordinate( const CEarthCoordinate& point_1, double distance, double direction, CEarthCoordinate& point_2 ) { CPolarCoordinate polar_point_1; CPolarCoordinate polar_point_2; Convert( point_1, polar_point_1 ); AddSurfaceDistanceAndDirectionToCoordinate( polar_point_1, distance, direction, polar_point_2 ); Convert( polar_point_2, point_2 ); } void CEarth::AddSurfaceDistanceAndDirectionToCoordinate( const CPolarCoordinate& point_1, double distance, double direction, CEarthCoordinate& point_2 ) { CPolarCoordinate polar_coordinate; AddSurfaceDistanceAndDirectionToCoordinate( point_1, distance, direction, polar_coordinate ); Convert( polar_coordinate, point_2 ); } void CEarth::AddSurfaceDistanceAndDirectionToCoordinate( const CPolarCoordinate& point_1, double distance, double direction, CPolarCoordinate& point_2 ) { // This is a translation of the Fortran routine DIRCT1 found in the // FORWRD3D program at: // ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/forwrd3d.for double c = 0.0; double c2a = 0.0; double cosine_of_direction = 0.0; double cosine_of_y = 0.0; double cu = 0.0; double cz = 0.0; double d = 0.0; double e = 0.0; double direction_in_radians = 0.0; double eps = 0.0; double heading_from_point_2_to_point_1_in_radians = 0.0; double point_1_latitude_in_radians = 0.0; double point_1_longitude_in_radians = 0.0; double point_2_latitude_in_radians = 0.0; double point_2_longitude_in_radians = 0.0; double r = 0.0; double sa = 0.0; double sine_of_direction = 0.0; double sine_of_y = 0.0; double su = 0.0; double tangent_u = 0.0; double term_1 = 0.0; double term_2 = 0.0; double term_3 = 0.0; double x = 0.0; double y = 0.0; direction_in_radians = CMath::ConvertDegreesToRadians( direction ); eps = 0.000000000000005; r = 1.0 - m_Flattening; point_1_latitude_in_radians = CMath::ConvertDegreesToRadians( point_1.GetUpDownAngleInDegrees() ); point_1_longitude_in_radians = CMath::ConvertDegreesToRadians( point_1.GetLeftRightAngleInDegrees() ); tangent_u = ( r * CMath::Sine( point_1_latitude_in_radians ) ) / CMath::Cosine( point_1_latitude_in_radians ); sine_of_direction = CMath::Sine( direction_in_radians ); cosine_of_direction = CMath::Cosine( direction_in_radians ); heading_from_point_2_to_point_1_in_radians = 0.0; if ( cosine_of_direction != 0.0 ) { heading_from_point_2_to_point_1_in_radians = CMath::ArcTangentOfYOverX( tangent_u, cosine_of_direction ) * 2.0; } cu = 1.0 / CMath::SquareRoot( ( tangent_u * tangent_u ) + 1.0 ); su = tangent_u * cu; sa = cu * sine_of_direction; c2a = ( (-sa) * sa ) + 1.0; x = CMath::SquareRoot( ( ( ( 1.0 / r / r ) - 1.0 ) * 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; tangent_u = distance / r / m_EquatorialRadiusInMeters / c; y = tangent_u; bool exit_loop = false; while( exit_loop != true ) { sine_of_y = CMath::Sine( y ); cosine_of_y = CMath::Cosine( y ); cz = CMath::Cosine( heading_from_point_2_to_point_1_in_radians + y ); e = ( cz * cz * 2.0 ) - 1.0; c = y; x = e * cosine_of_y; y = ( e + e ) - 1.0; term_1 = ( sine_of_y * sine_of_y * 4.0 ) - 3.0; term_2 = ( ( term_1 * y * cz * d ) / 6.0 ) + x; term_3 = ( ( term_2 * d ) / 4.0 ) - cz; y = ( term_3 * sine_of_y * d ) + tangent_u; if ( CMath::AbsoluteValue( y - c ) > eps ) { exit_loop = false; } else { exit_loop = true; } } heading_from_point_2_to_point_1_in_radians = ( cu * cosine_of_y * cosine_of_direction ) - ( su * sine_of_y ); c = r * CMath::SquareRoot( ( sa * sa ) + ( heading_from_point_2_to_point_1_in_radians * heading_from_point_2_to_point_1_in_radians ) ); d = ( su * cosine_of_y ) + ( cu * sine_of_y * cosine_of_direction ); point_2_latitude_in_radians = CMath::ArcTangentOfYOverX( d, c ); c = ( cu * cosine_of_y ) - ( su * sine_of_y * cosine_of_direction ); x = CMath::ArcTangentOfYOverX( sine_of_y * sine_of_direction, c ); c = ( ( ( ( ( -3.0 * c2a ) + 4.0 ) * m_Flattening ) + 4.0 ) * c2a * m_Flattening ) / 16.0; d = ( ( ( ( e * cosine_of_y * c ) + cz ) * sine_of_y * c ) + y ) * sa; point_2_longitude_in_radians = ( point_1_longitude_in_radians + x ) - ( ( 1.0 - c ) * d * m_Flattening ); heading_from_point_2_to_point_1_in_radians = CMath::ArcTangentOfYOverX( sa, heading_from_point_2_to_point_1_in_radians ) + CMath::Pi(); point_2.SetUpDownAngleInDegrees( CMath::ConvertRadiansToDegrees( point_2_latitude_in_radians ) ); point_2.SetLeftRightAngleInDegrees( CMath::ConvertRadiansToDegrees( point_2_longitude_in_radians ) ); } void CEarth::Convert( const CEarthCoordinate& cartesian_coordinate, CPolarCoordinate& polar_coordinate ) const { // convert from cartesian to polar double equatorial_radius_times_eccentricity_squared = 0.0; equatorial_radius_times_eccentricity_squared = m_EquatorialRadiusInMeters * m_EccentricitySquared; double p = 0.0; p = CMath::SquareRoot( ( cartesian_coordinate.GetXCoordinateInMeters() * cartesian_coordinate.GetXCoordinateInMeters() ) + ( cartesian_coordinate.GetYCoordinateInMeters() * cartesian_coordinate.GetYCoordinateInMeters() ) ); double temp_latitude = 0.0; double z_coordinate = cartesian_coordinate.GetZCoordinateInMeters(); // for convienance double one_minus_eccentricity_squared = 1.0 - m_EccentricitySquared; temp_latitude = z_coordinate / p / one_minus_eccentricity_squared; double old_value = 0.0; double temp_value = 0.0; double part_a = 0.0; double part_b = 0.0; double part_c = 0.0; unsigned long loop_index = 0; unsigned long maximum_number_of_tries = 1024; bool convergence_was_acheived = false; while( convergence_was_acheived != true && loop_index < maximum_number_of_tries ) { old_value = temp_latitude; part_a = one_minus_eccentricity_squared * temp_latitude * temp_latitude; part_b = equatorial_radius_times_eccentricity_squared / CMath::SquareRoot( 1.0 + part_a ); part_c = p - part_b; temp_latitude = z_coordinate / part_c; loop_index++; if ( CMath::AbsoluteValue( temp_latitude - old_value ) > 0.000000000000000000001 ) { // Oh well, try again... } else { // YIPEE!! We've reached convergence! convergence_was_acheived = true; } } if ( convergence_was_acheived == true ) { double latitude_angle_in_radians = 0.0; // Save the UpDown angle in degrees latitude_angle_in_radians = CMath::ArcTangent( temp_latitude ); polar_coordinate.SetUpDownAngleInDegrees( CMath::ConvertRadiansToDegrees( latitude_angle_in_radians ) ); // Latitude double sine_of_latitude_in_radians = 0.0; double cosine_of_latitude_in_radians = 0.0; sine_of_latitude_in_radians = CMath::Sine( latitude_angle_in_radians ); cosine_of_latitude_in_radians = CMath::Cosine( latitude_angle_in_radians ); double longitude_in_radians = 0.0; longitude_in_radians = CMath::ArcTangentOfYOverX( cartesian_coordinate.GetYCoordinateInMeters(), cartesian_coordinate.GetXCoordinateInMeters() ); polar_coordinate.SetLeftRightAngleInDegrees( CMath::ConvertRadiansToDegrees( longitude_in_radians ) ); // Longitude double w = 0.0; w = CMath::SquareRoot( 1.0 - ( m_EccentricitySquared * sine_of_latitude_in_radians * sine_of_latitude_in_radians ) ); double distance_from_center_to_surface_of_the_ellipsoid = 0.0; distance_from_center_to_surface_of_the_ellipsoid = m_EquatorialRadiusInMeters / w; double distance_from_surface = 0.0; if ( CMath::AbsoluteValue( latitude_angle_in_radians ) < 0.7854 ) { //printf( "fabs( %14.10lf ) < 0.7854\n", latitude_angle_in_radians ); distance_from_surface = ( p / cosine_of_latitude_in_radians ) - distance_from_center_to_surface_of_the_ellipsoid; } else { //printf( "fabs( %14.10lf ) >= 0.7854\n", latitude_angle_in_radians ); distance_from_surface = ( z_coordinate / sine_of_latitude_in_radians ) - distance_from_center_to_surface_of_the_ellipsoid + ( m_EccentricitySquared * distance_from_center_to_surface_of_the_ellipsoid ); } // printf( "CEarth::Convert() : First method produced a length of %14.10lf\n", distance_from_surface ); polar_coordinate.SetDistanceFromSurfaceInMeters( distance_from_surface ); } else { // Oh well, we gave it a shot.. polar_coordinate.Set( 0.0, 0.0, 0.0 ); } } void CEarth::Convert( const CPolarCoordinate& polar_coordinate, CEarthCoordinate& cartesian_coordinate ) const { // convert from polar to cartesian double up_down_radians = 0.0; // latitude double left_right_radians = 0.0; // longitude angle up_down_radians = CMath::ConvertDegreesToRadians( polar_coordinate.GetUpDownAngleInDegrees() ); left_right_radians = CMath::ConvertDegreesToRadians( polar_coordinate.GetLeftRightAngleInDegrees() ); double sine_of_up_down_radians = 0.0; double cosine_of_left_right_radians = 0.0; // cosine_of_longitude double cosine_of_up_down_radians = 0.0; // cosine_of_latitude sine_of_up_down_radians = CMath::Sine( up_down_radians ); cosine_of_left_right_radians = CMath::Cosine( left_right_radians ); cosine_of_up_down_radians = CMath::Cosine( up_down_radians ); // Now we need to calculate the distance from the center of the ellipsoid to the surface of the ellipsoid double w = 0.0; w = CMath::SquareRoot( 1.0 - ( m_EccentricitySquared * sine_of_up_down_radians * sine_of_up_down_radians ) ); double distance_from_center_to_surface_of_the_ellipsoid = 0.0; distance_from_center_to_surface_of_the_ellipsoid = m_EquatorialRadiusInMeters / w; // printf( "en = %.25lf\n", distance_from_center_to_surface_of_the_ellipsoid ); double coordinate = 0.0; coordinate = ( distance_from_center_to_surface_of_the_ellipsoid + polar_coordinate.GetDistanceFromSurfaceInMeters() ) * cosine_of_up_down_radians * cosine_of_left_right_radians; cartesian_coordinate.SetXCoordinateInMeters( coordinate ); 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 (werner@tij.wb.xerox.com) 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 (bjmillar@dera.gov.uk) 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 ); } 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 ToolTipFormatLine=CEarth=m_EllipsoidID= GFC - CEarth CEarth
$Revision: 11 $
Description
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.Methods
- If you were to shine a laser from
point_1pointing towardsdirectionfordistancemeters, this function will tell you what location you would be at. It will also let you specify a pointheight_above_surface_of_point_2meters above the surface. This method does not take into account the curvature of the Earth.- 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 "If I head out in this direction for that amout of meters, where will I be?"
- This method allows you to convert from polar and cartestian coordinates.
- This tells you how far (in meters) from the horizon
point_1is.- This tells you what the equatorial radius is in meters for the selected ellipsoid.
- This tells you what the polar radius is in meters for the selected ellipsoid.
- Draw a line from
point_atopoint_b. This function will tell you how farcurrent_locationis from that line.- This will tell you how many meters it is between two points. It answers the question, "If I pointed a laser from
point_1topoint_2, how far would the laser beam travel?"- This will tell you how many meters it is between two points. It answers the question, "If I were to walk from
point_1topoint_2, how far would I walk?"- This allows you to set the ellipsoid used by CEarth in its calculations. The default is
WGS84which is generally accepted as being the closest approximation of the Earth's ellipsoid. Theellipsoidparameter may be one of the following:
- Perfect_Sphere
- Airy,
- Austrailian_National,
- Bessell_1841,
- Bessel_1841_Nambia,
- Clarke_1866,
- Clarke_1880,
- Everest,
- Fischer_1960_Mercury,
- Fischer_1968,
- GRS_1967,
- GRS_1980,
- Helmert_1906,
- Hough,
- International,
- Krassovsky,
- Modified_Airy,
- Modified_Everest,
- Modified_Fischer_1960,
- South_American_1969,
- Topex_Poseidon_Pathfinder_ITRF,
- WGS_60,
- WGS_66,
- WGS_72,
- WGS_84
- NAD_27
- Tokyo
- This let's you use your own (custom) values to describe the ellipsoid of the Earth.
- This let's you use your own (custom) values to describe the ellipsoid of the Earth.
Example
Copyright, 1998, Samuel R. Blackburn#include <stdio.h> #include <GFC.h> #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 ) ); CEarth 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( "Distance between here and there: %.23lf meters\nHeading from here to there: %.19lf degrees\nHeading from there to here: %.19lf degrees\n", 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( "Heading %lf degrees, %lf minutes, %lf seconds\n", degrees, minutes, seconds ); CMath::ConvertDecimalDegreesToDegreesMinutesSeconds( heading_from_there_to_here, degrees, minutes, seconds ); printf( "Heading %lf degrees, %lf minutes, %lf seconds\n", degrees, minutes, seconds ); }
$Workfile: CEarth.cpp $
$Modtime: 9/01/98 9:56p $#endif