gfc.tar_CPolarCoordinate.cpp_GFC_earth_地球_地球两点距离

#include <GFC.h> #pragma hdrstop #include <stdio.h> /* ** 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 =