[Reprint]: Latitude and WGS84 coordinate conversion

This code is implemented in the WGS84 system geodetic Coordinates (BLH) and space rectangular coordinates (XYZ) of the mutual transformation, conforms to the standard syntax, can be used directly

The following code, the output is:
WGS84:-2175790.73969891 4461032.11207734 3992337.79032463
blh:38.9999999999998 116.000000000000 33.0000069718808

Module Corrtrans!//WGS84 system BLH coordinates and spatial coordinate conversion!//Fortran coder http://fcode.cn!//adapted from Acheng ' s C + + code IMPLI CIT None Integer, parameter:: DP = Selected_real_kind (a) Real (KIND=DP), parameter:: PI = 3.14159265358979323846_d P Real (KIND=DP), parameter, private:: A = 6378137._DP Real (KIND=DP), parameter, private:: F = 1.0_dp/298.2572235 63_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 R EAL (KIND=DP), Intent (out):: B, L, H Real (KIND=DP):: N integer:: I L = Atan (ABS (y/x)) B = ATA N (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=d P), 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 Blh2xyzend Module corrtransprogram 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, HEnd program WWW_FCODE_CN

Transferred from: http://fcode.cn/code_prof-89-1.html

[Reprint]: Latitude and WGS84 coordinate conversion