C ********* C *XYZ/PLH* C ********* C CONVERT XYZ TO AND FROM PHI,LAMBDA,H C INPUT DECK CARD 1 CONTAINS A,B,X0,Y0,Z0 FOR ELLIPSOID AND DATUM (5F10.0) C B CAN BE EITHER SEMI-MINOR AXIS OR RECIPROCAL OF FLATTENING CARD 2 CONTAINS SWITCH .LT. 0. FOR XYZ TO PLH (I2) C .GT. 0. FOR PLH TO XYZ C .EQ. FOR END OF DATA C TWO CARDS PER DATA POINT C FIRST CARD IS 80 CHARACTER DESCRIPTION FOR LABELLING (80A1) C SECOND CARD CONTAINS INPUT VALES C FOR XYZ - 3 VALUES IN METRES (3F15.0) C FOR PLH - 8 INPUT VALUES (8F10.0) C LAT IN DEG,MIN,SEC (POSITIVE NORTH) C LONG IN DEG,MIN,SEC (POSITIVE WEST) C HEIGHT = SUM OF ORTH HT(IN FEET) AND GEOID UNDUL(IN M) C IF ALL INPUT VALUES ARE ZERO PROGRAM CYCLES BACK TO CARD 2 C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) 0390 COMMON A,B,X0,Y0,Z0,PI,ICR,IPR PI=3.141592653589793D0 0420 ICR = 5 IPR = 6 C READ AND PRINT ELLIPSOID PARAMETERS AND DATUM TRANSLATION COMPONENTS C IF B .LT. 6000000 ASSUME IT IS 1./F AND COMPUTE B FROM IT READ *,A,B,X0,Y0,Z0 C READ(ICR,1001) A,B,X0,Y0,Z0 WRITE(IPR,1002) A,B,X0,Y0,Z0 IF(B.LT.6D6) B=A*(1D0-1D0/B) 0540 C CYCLE THROUGH CONVERSION ROUTINES AS DIRECTED BY SWITCH 10 CONTINUE READ *,ISW C READ(ICR,1003) ISW IF(ISW .EQ. 0) GO TO 20 IF(ISW .LT. 0) CALL PLH IF(ISW .GT. 0) CALL XYZ GO TO 10 20 STOP C FORMAT STATEMENTS 1001 FORMAT(5F10.0) 1002 FORMAT(1H1,15(/),50X,9(1H*),/,50X,9H*XYZ/PLH*,/,50X,9(1H*),///, * 20X,38HELLIPSOID PARAMETERS USED FOR THIS RUN//25X, * 15HSEMI-MAJOR AXIS,28X,F20.3/25X,28HSEMI-MINOR AXIS OR RECIPROCA * ,15HL OF FLATTENING,F20.3//25X,28HDATUM TRANSLATION COMPONENTS, * 5H - X0,10X,F20.3,/,56X,2HY0,10X,F20.3,/,56X,2HZ0,10X,F20.3) 1003 FORMAT(I2) END 0670 C SUBROUTINE XYZ 0680 C CONVERSION OF LAT/LONG/HT TO XYZ C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) 0740 COMMON A,B,X0,Y0,Z0,PI,ICR,IPR DIMENSION DESC(80) C READ STATION DESCRIPTION CARD AND INPUT VALUES WRITE(IPR,1001) 10 READ(ICR,1002) DESC WRITE(IPR,1003) DESC READ *,DLAT,XMLAT,SLAT,DLONG,XMLONG,SLONG,ORTH,UNDUL C READ(ICR,1004) DLAT,XMLAT,SLAT,DLONG,XMLONG,SLONG,ORTH,UNDUL IF(DLAT .EQ. 0.) GO TO 20 C COMPUTE LAT IN RADIANS AND PRINT LAT IN DEG/MIN/SEC AND RAD RLAT=(DLAT+XMLAT/60D0+SLAT/3600D0)*PI/180D0 0950 WRITE(IPR,1005) DLAT,XMLAT,SLAT,RLAT C CHANGE SIGN OF INPUT LONGITUDE (POS WEST TO USER, POS EAST TO PROGRAM) C COMPUTE LONG IN RADIANS AND PRINT IN DEG/MIN/SEC AND RAD DLONG=-DLONG 1010 XMLONG=-XMLONG 1020 SLONG=-SLONG 1030 RLONG=(DLONG+XMLONG/60D0+SLONG/3600D0)*PI/180D0 1070 WRITE(IPR,1006) DLONG,XMLONG,SLONG,RLONG C COMPUTE AND PRINT ELLIPSOID HEIGHT IN METRES XHT=ORTH*0.3048 + UNDUL 1120 WRITE(IPR,1007) ORTH,UNDUL,XHT C CONVERT FROM PLH TO XYZ AND PRINT XYZ CALL GEOREC(A,B,X0,Y0,Z0,RLAT,RLONG,XHT,X,Y,Z) 1170 WRITE(IPR,1008) X,Y,Z C RECONVERT BACK TO PLH AND PRINT CALL RECGEO(A,B,X0,Y0,Z0,X,Y,Z,RLAT,RLONG,XHT) 1220 DLAT=DEG(RLAT,XMLAT,SLAT) 1230 DLONG = DEG(RLONG,XMLONG,SLONG) 1240 WRITE(IPR,1009) DLAT,XMLAT,SLAT WRITE(IPR,1010) DLONG,XMLONG,SLONG WRITE(IPR,1011) XHT C GET ANOTHER DATA POINT GO TO 10 C END OF DATA 20 RETURN C FORMAT STATEMENTS 1001 FORMAT(1H1/35X,38(1H*)/35X,31H* CONVERSION OF LAT/LONG/HT TO , * 7HX/Y/Z */35X,38(1H*)//) 1002 FORMAT(80A1) 1003 FORMAT(///10X,80A1) 1004 FORMAT(8F10.0) 1005 FORMAT(6H LAT ,6X,2F7.1,F10.4,F15.11) 1006 FORMAT(6H LONG,6X,2F7.1,F10.4,F15.11) 1007 FORMAT(6H HT ,6X,F8.2,7H FEET +,F7.2,9H METRES =,F7.2,7H METRES) 1008 FORMAT(52X,3HXYZ,6X,3F14.4) 1009 FORMAT(12H CHECK LAT ,2F7.1,F10.4) 1010 FORMAT(12H CHECK LONG,2F7.1,F10.4) 1011 FORMAT(12H CHECK HT ,F8.2) END 1500 C SUBROUTINE PLH 1510 C CONVERSION OF XYZ TO LAT/LONG/HT C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) 1570 COMMON A,B,X0,Y0,Z0,PI,ICR,IPR DIMENSION DESC(80) C READ AND PRINT STATION DESCRIPTION AND INPUT VALUES WRITE(IPR,1001) 10 READ(ICR,1002) DESC WRITE(IPR,1003) DESC READ *,X,Y,Z C READ(ICR,1004) X,Y,Z IF(X .EQ. 0.) GO TO 20 WRITE(IPR,1005) X,Y,Z C CONVERT FROM XYZ TO PLH AND PRINT PLH CALL RECGEO(A,B,X0,Y0,Z0,X,Y,Z,RLAT,RLONG,XHT) 1780 DLAT=DEG(RLAT,XMLAT,SLAT) 1790 DLONG=DEG(RLONG,XMLONG,SLONG) 1800 WRITE(IPR,1006) DLAT,XMLAT,SLAT,RLAT WRITE(IPR,1007) DLONG,XMLONG,SLONG,RLONG WRITE(IPR,1008) XHT C RECONVERT BACK TO XYZ AND PRINT CALL GEOREC(A,B,X0,Y0,Z0,RLAT,RLONG,XHT,X,Y,Z) 1910 WRITE(IPR,1009) X,Y,Z C GET ANOTHER DATA POINT GO TO 10 C END OF DATA 20 RETURN C FORMAT STATEMENTS 1001 FORMAT(1H1/35X,38(1H*)/35X,36H* CONVERSION OF X/Y/Z TO LAT/LONG/HT * ,2H */35X,38(1H*)//) 1002 FORMAT(80A1) 1003 FORMAT(///10X,80A1) 1004 FORMAT(3F15.0) 1005 FORMAT(52X,3HXYZ,6X,3F14.4) 1006 FORMAT(6H LAT ,6X,2F7.1,F10.4,F15.11) 1007 FORMAT(6H LONG,6X,2F7.1,F10.4,F15.11) 1008 FORMAT(6H HT ,F8.2) 1009 FORMAT(52X,9HCHECK XYZ,3F14.4) END 2130 C FUNCTION DEG(RAD,XMIN,SEC) 2140 C FUNCTION TO CONVERT FROM RADIANS TO DEGREES/MINUTES/SECONDS C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) 2200 PI=3.141592653589793D0 2210 DEG=RAD*180./PI 2220 IDEG=DEG 2230 MIN=(DEG-IDEG)*60. 2240 XMIN = MIN 2245 SEC=((DEG-IDEG)*60.-MIN)*60. 2250 DEG = IDEG 2255 RETURN 2260 END 2270 C SUBROUTINE GEOREC(A,B,X0,Y0,Z0,PHI,XLAMB,HT,X,Y,Z) C GEODETIC TO CARTESIAN COORDINATE CONVERSION C INPUT GEODETIC COORDINATES PHI,XLAMB,HT C REFERRED TO ELLIPSOID WITH SEMI AXES A,B AND ORIGIN OFFSETS X0,Y0,Z0 C OUTPUT CARTESIAN GEOCENTRIC COORDINATES X,Y,Z C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) 2400 SQRT(PI) = DSQRT(PI) SIN(PI) = DSIN(PI) COS(PI) = DCOS(PI) BOA=(B/A)*(B/A) 2410 CP = COS(PHI) SP = SIN(PHI) CL = COS(XLAMB) SL = SIN(XLAMB) XN = A / SQRT(CP**2 + BOA * SP**2) X = X0 + (XN + HT) * CP * CL Y = Y0 + (XN + HT) * CP * SL Z = Z0 + (XN*BOA + HT) * SP RETURN 2460 END 2470 C SUBROUTINE RECGEO(A,B,X0,Y0,Z0,X,Y,Z,PHI,XLAMB,HT) C CARTESIAN TO GEODETIC COORDINATE CONVERSION C INPUT CARTESIAN GEOCENTRIC COORDINATES X,Y,Z C SEMI-AXES A,B AND ORIGIN OFFSETS X0,Y0,Z0 OF ELLIPSOID TO C WHICH GEODETIC COORDINATES TO BE REFERRED C OUTPUT GEODETIC COORDINATES PHI/XLAMB/HT C ITERATIVE COMPUTATION CONTINUED UNTIL CONVERGENCE TO 1E-10 RADIANS IN C PHI AND A*1E-10 METRES IN HT, OR FOR MAXIMUM OF 20 ITERATIONS C**MACHINE DEPENDENT STATEMENT IMPLICIT REAL*8(A-H,O-Z) 2640 ABS(PI) = DABS(PI) SQRT(PI) = DSQRT(PI) SIN(PI) = DSIN(PI) COS(PI) = DCOS(PI) ATAN(PI) = DATAN(PI) ATAN2(X,Y) = DATAN2(X,Y) C INITIALIZE EPS2 = 1E-10 EPS1 = A * EPS2 ITER = 20 BOA=(B/A)*(B/A) 2680 E2=1.-BOA 2690 C TRANSLATE CARTESIAN ORIGIN TO GEODETIC ORIGIN X1 = X - X0 Y1 = Y - Y0 Z1 = Z - Z0 C COMPUTE LONGITUDE NON-ITERATIVELY P = SQRT(X1**2 + Y1**2) XLAMB = ATAN2(Y1,X1) C INITIAL APPROXIMATIONS FOR ITERATION XN=A 2840 HT = SQRT(X1**2 + Y1**2 + Z1**2) - A PHI = ATAN((Z1/P)/(1.-E2*XN/(XN+HT))) C ITERATION LOOP DO 10 I = 1,ITER HT1=HT 2920 PHI1=PHI 2930 XN = A / SQRT(COS(PHI)**2 + BOA * SIN(PHI)**2) HT = P / COS(PHI) - XN PHI = ATAN((Z1/P)/(1.-E2*XN/(XN+HT))) IF(ABS(HT1-HT) .LT. EPS1 .AND. ABS(PHI1-PHI) .LT. EPS2) GO TO 20 10 CONTINUE C TOO MANY ITERATIONS RETURN C CONVERGENCE ACHIEVED 20 RETURN END 3100