本代码实现在WGS84系统的大地坐标(BLH)和空间直角坐标(XYZ)的互相转换,符合标准语法,可直接使用
如下代码,输出为:
WGS84: -2175790.73969891 4461032.11207734 3992337.79032463
BLH: 38.9999999999998 116.000000000000 33.0000069718808
Module CorrTrans !// WGS84 系统BLH坐标与空间直角坐标转换 !// Fortran Coder http://fcode.cn !// Adapted from Acheng\'s C++ code Implicit None Integer , parameter :: DP = Selected_Real_Kind( 12 ) Real(kind=DP) , parameter :: PI = 3.14159265358979323846_DP Real(kind=DP) , parameter , private :: a = 6378137._DP Real(kind=DP) , parameter , private :: f = 1.0_DP / 298.257223563_DP Real(kind=DP) , parameter , private :: t = a * ( 1.0_DP - f ) Real(kind=DP) , parameter , private :: e = sqrt( (a*a - t*t) / (a*a) ) contains Subroutine XYZ2BLH( x , y , z , B , L , H ) Real(Kind=DP) , Intent( IN ) :: x , y , z Real(Kind=DP) , Intent( OUT ) :: B , L , H real(kind=DP) :: N integer :: i L = atan( abs(y/x) ) B = atan( abs(z/sqrt(x*x+y*y) ) ) do i = 1 , 5 N = a / sqrt( 1.0_DP - e*e*sin(B)*sin(B) ) H = z / sin(B) - N * ( 1.0_DP - e*e ) B = atan( z*(N+H) / ( sqrt(x*x+y*y)*(N*(1.0_DP-e*e)+H) ) ) end do if ( x < 0._DP .and. y > 0._DP ) L = PI - L if ( x > 0._DP .and. y < 0._DP ) L = -1.0_DP * L if ( x < 0._DP .and. y < 0._DP ) L = -1.0 * PI + L B = B * 180._DP / PI L = L * 180._DP / PI End Subroutine XYZ2BLH Subroutine BLH2XYZ( B , L , H , x , y , z ) Real(Kind=DP) , Intent( IN ) :: B , L , H Real(Kind=DP) , Intent( OUT ) :: x , y , z real(kind=DP) :: N , Br , Lr Br = B * PI / 180._DP Lr = L * PI / 180._DP N = a / sqrt( 1.0_DP - e*e*sin(Br)*sin(Br) ) x = ( N + H ) * cos(Br)*cos(Lr) y = ( N + H ) * cos(Br)*sin(Lr) z = ( N*(1.0_DP - e*e ) + H ) * sin(Br); End Subroutine BLH2XYZ End Module CorrTrans Program www_fcode_cn Use CorrTrans Implicit None Real(Kind=DP) :: x , y , z Real(Kind=DP) :: B , L , H B = 39._DP ; L = 116._DP ; H = 33._DP call BLH2XYZ( B , L , H , x , y , z ) write(*,*) \'WGS84:\' , x , y , z call XYZ2BLH( x , y , z , B , L , H ) write(*,*) \'BLH:\' , B , L , H End Program www_fcode_cn
转自:http://fcode.cn/code_prof-89-1.html