IMPLICIT REAL*8(A-H,O-Z) REAL*8 KO,K,L1,L2,MC,MPAZ,MPAZAP,MPDIST,LS ,L3,MC2,LS2 REAL*8 VCINIT(4,4)/16*0.D0/ DIMENSION DM(2,2),CM(2,2), CM2(2,2), VCDIR(4,4),VCINV(2,2), + DM2(2,2) 30 FORMAT(2F11.2) 31 FORMAT(2I3,F7.3) 32 FORMAT(F14.4) 34 FORMAT(I16,2(2I4,F9.4)) 36 FORMAT(F7.5) 37 FORMAT(I4,I3,F8.4,F10.4,F5.2,F7.4) 41 FORMAT('0',20X,'OBSERVED GEODETIC DATA') 42 FORMAT('0',10X,'AZIMUTH FROM 1 TO 2',5X,'DISTANCE') 43 FORMAT('0',14X,I3,I3,F6.2,8X,F8.3) 44 FORMAT(10X,'SEMI-MAJOR AXIS =',F9.1,5X,'SEMI-MINOR AXIS =',F9.1,// +) 45 FORMAT('1',37X,'TRANSVERSE MERCATOR PROJECTION'////) 46 FORMAT(10X,'LONGITUDE OF CENTRAL MERIDIAN',2I5,F8.3,///) 47 FORMAT('0',10X,'STATION',5X,'LATITUDE',8X,'LONGITUDE',11X,'X',13X, 1'Y') 48 FORMAT(10X,'SCALE FACTOR=',F6.4,5X,'FALSE EASTING =',F9.1,//) 49 FORMAT('0',11X,I3,5X,I3,I3,F8.4,3X,I3,I3,F8.4,3X,F11.3,3X,F11.3) 50 FORMAT(10X,'COVARIANCE MATRIX OF THE TRANSVERSE MERCATOR COORDINAT 1ES IN METRES SQUARED OF POINT 1') 51 FORMAT('1',10X,'COVARIANCE INFORMATION - DIRECT CASE') 52 FORMAT('1',10X,'COVARIANCE INFORMATION - INVERSE CASE') 53 FORMAT('0',10X,'COVARIANCE MATRIX OF POINTS 1 AND 2 IN METRES SQUA 1RED') 54 FORMAT(10X,'COVARIANCE MARRIX OF POINT 1, THE AZIMUTH FROM 1 TO 2, IN SECO + IN SECONDS SQUARED AND THE DISTANCE, IN METRES SQUARED',/) 55 FORMAT('0',10X,'AZIMUTH FROM 1 TO 2',5X,'AZIMUTH FROM 2 TO 1',5X,' 1DISTANCE',//) 56 FORMAT(13X,I4,I3,F7.3,10X,I4,I3,F7.3,7X,F9.3,//) 57 FORMAT(10X,'COVARIANCE MATRIX OF PHI 2 AND LAMDA 2, IN SECONDS SQU +ARED',//) 58 FORMAT('0',20X,'COMPUTED DATA') 59 FORMAT('0',10X,'T-T CORRECTION',3X,'MERIDIAN CONVERGENCE',3X,'LINE 1 SCALE') 60 FORMAT('0',10X,I3,I3,F8.4,6X,I3,I3,F8.4,8X,F10.8) 61 FORMAT('1',20X,'DIRECT PROBLEM') 62 FORMAT('1',20X,'INVERSE PROBLEM') 63 FORMAT('1') PI=3.141592653589793D0 ICODE=1 C C READ IN DATA - ELLIPSOIDAL PARAMETERS, MAPPING PLANE PARAMETERS, POSITION C DATA AND VARIANCE COVARIANCE INFORMATION OF POINT 1. C READ 30,A,B READ 36,SF READ 32,XO READ 31,JLONGD,JLONGM,OLONGS READ 34,NAME,IPHID,IPHIM,PHIS,LONGD,LONGM,ALONGS READ *,DM(1,1),DM(1,2),DM(2,2) PRINT 45 PRINT 44,A,B PRINT 48,SF,XO PRINT 46,JLONGD,JLONGM,OLONGS C C CONVERT ANGLES TO RADIANS C CALL DMSRAD(JLONGD,JLONGM,OLONGS,CMRAD) CALL DMSRAD(LONGD,LONGM,ALONGS,L1) CALL DMSRAD(IPHID,IPHIM,PHIS,P1) DL1=L1-CMRAD C C TRANSFORM THE GEOGRAPHIC COORDINATES OF POINT 1 TO THE MAPPING PLANE C COORDINATES. C CALL TMPLXY(P1,DL1,A,B,SF,XO,CMRAD,X1,Y1) PRINT 61 PRINT 41 C C READ IN THE OBSERVATIONAL DATA AND THEIR STANDARD DEVIATIONS. C READ 37,IGD12,IGM12,GS12,GDIST,SDAZ,SDD PRINT 47 PRINT 49,NAME,IPHID,IPHIM,PHIS,LONGD,LONGM,ALONGS,X1,Y1 PRINT 42 PRINT 43,IGD12,IGM12,GS12,GDIST C C CONVERT THE OBSERVED AZIMUTH TO RADIANS. C CALL DMSRAD(IGD12,IGM12,GS12,GAZ) C C COMPUTE THE MERIDIAN CONVERGENCE ,T-T CORRECTION AND LINE SCALE FACTOR. C CALL TMSFMC(P1,DL1,SF,A,B,SF1,MC) MPAZAP=GAZ-MC X2A=X1+DSIN(MPAZAP)*GDIST Y2A=Y1+DCOS(MPAZAP)*GDIST CALL TTLS(A,B,P1,GDIST,X1,Y1,X2A,Y2A,XO,SF,TT,LS) CALL RADMS(MC,IMCD,IMCM,SMC) CALL RADMS(TT,ITTD,ITTM,STT) C C SOLVE FOR X2 AND Y2. C CALL MPDIR(X1,Y1,GDIST,GAZ,MC,LS,TT,MPDIST,MPAZ,X2,Y2) PRINT 58 CALL TMXYPL(X2,Y2,A,B,SF,XO,CMRAD,P2,L2) DL2=L2-CMRAD CALL RADMS(P2,IPD2,IPM2,PS2) CALL RADMS(L2,LD2,LM2,SL2) PRINT 59 PRINT 60,ITTD,ITTM,STT,IMCD,IMCM,SMC,LS N2=2 PRINT 47 PRINT 49,N2,IPD2,IPM2,PS2,LD2,LM2,SL2,X2,Y2 C C FORM VARIANCE COVARIANCE MATRIX OF POINT 1 AND THE OBSEVED DATA AND C CONVERT TO THE PROPER UNITS. C PRINT 51 DM(2,1)=DM(1,2) DO 14 L=1,2 DO 14 J=1,2 DP2=DM(L,J)*PI**2/180.D0**2 DM(L,J)=DP2/1.296D7 14 CONTINUE CALL ERRTM(A,B,P1,DL1,SF,ICODE,CM,DM) PRINT 50 CALL MOUTD(CM,2,2,2) DO 10 I=1,2 DO 10 J=1,2 VCINIT(I,J)=CM(I,J) 10 CONTINUE VCINIT(3,3)=(2.96D-2)**2 VCINIT(4,4)=((5.D0/36.D2)*PI/18.D1)**2 C C COMPUTE THE VARIANCE COVARIANCE MATRIX OF POINTS 1 AND 2. C CALL ERRMPD(MPAZ,MPDIST,VCINIT,VCDIR) C C PRINT OUT THE COMPUTED DATA. C PRINT 53 CALL MOUTD(VCDIR,4,4,4) DO 11 I=1,2 DO 11 J=1,2 CM2(I,J)=VCDIR(I+2,J+2) 11 CONTINUE ICODE=-1 CALL ERRTM(A,B,P2,DL2,SF,ICODE,CM2,DM2) DO 12 I=1,2 DO 12 J=1,2 DM2(I,J)=DM2(I,J)*180.D0**2/PI**2*36.D2**2 12 CONTINUE PRINT 57 CALL MOUTD(DM2,2,2,2) C C ***** INVERSE CASE ***** C C C SOLVE FOR THE MAPPING PLANE AZIMUTH AND DISTANCE. C PRINT 62 PRINT 41 PRINT 47 PRINT 49,NAME,IPHID,IPHIM,PHIS,LONGD,LONGM,ALONGS,X1,Y1 PRINT 49,N2,IPD2,IPM2,PS2,LD2,LM2,SL2,X2,Y2 CALL TMSFMC(P2,DL2,SF,A,B,SF2,MC2) CALL TTLS(A,B,P2,GDIST,X2,Y2,X1,Y1,XO,SF,TT2,LS2) CALL MPINV(X1,Y1,X2,Y2,TT,MC,LS,TT2,MC2 ,GDIST,GAZ12,GAZ21) CALL RADMS(GAZ12,IGD,IGM,SG) CALL RADMS(GAZ21,IG2D,IG2M,SG2) C C PRINT OUT THE COMPUTED INFORMATION. C PRINT 58 PRINT 55 PRINT 56,IGD,IGM,SG,IG2D,IG2M,SG2,GDIST C C COMPUTE THE VARIANCE COVARIANCE MATRIX OF THE COMPUTED MAPPING PLANE AZIMUTH C AND DISTANCE. C PRINT 52 CALL ERRMPI(MPDIST,X1,Y1,X2,Y2,VCDIR,VCINV) C C CONVERT TO THE PROPER UNITS. C VCINV(1,2)=VCINV(1,2)*180.D0/PI*3600.D0 VCINV(2,1)=VCINV(1,2) VCINV(1,1)=VCINV(1,1)*180.D0**2/PI**2*3600.D0**2 PRINT 54 CALL MOUTD(VCINV,2,2,2) PRINT 63 STOP END SUBROUTINE MPDIR(X1,Y1,GDIST,GAZ,MC,LSK,TT,MPDIST,MPAZ,X2,Y2) C C THIS SUBROUTINE WILL SOLVE THE DIRECT PROBLEM OF GEODETIC POSITIONING C ON THE MAPPING PLANE. C C INPUT: C C X1,Y1 - X,Y COORDINATES OF THE INITIAL POINT. C GDIST - GEODESIC DISTANCE. C GEODETIC AZIMUTH (IN RADIANS). C MC - MERIDIAN CONVERENCE (IN RADIANS). C LSK - LINE SCALE FACTOR. C TT - T-T CORRECTION (IN RADIANS). C C OUTPUT: C C MPDIST - DISTANCE FROM POINT 1 TO 2 ON THE MAPPING PLANE. C MPAZ - AZIMUTH ON THE MAPPING PLANE (IN RADIANS). C X2 Y2 - X,Y,COORDINATES OF THE OBSERVED POINT. C C WRITTEN BY G.BOWIE, DEC. 1977. C IMPLICIT REAL*8 (A-Z) MPAZ=GAZ-MC-TT MPDIST=GDIST*LSK X2=X1+DSIN(MPAZ)*MPDIST Y2=Y1+DCOS(MPAZ)*MPDIST RETURN END SUBROUTINE MPINV(X1,Y1,X2,Y2,TT1,MC1,LSK,TT2,MC2,GDIST,GAZ12,GAZ21 1) C C THIS SUBROUTINE WILL SOLVE THE INVERSE PROBLEM OF GEODETIC POSITIONING C ON THE MAPPING PLANE. C C C INPUT: C C X1,Y1 - X,Y COORDINATES OF THE INITIAL POINT. C X2 Y2 - X,Y,COORDINATES OF POINT 2. C TT1 - T-TCORRECTION AT POINT 1.(IN RADIANS) C TT2 - T-TCORRECTION AT POINT 2.(IN RADIANS) C MC1 - MERIDIAN CONVERGENCE AT POINT 1.(IN RADIANS) C MC2 - MERIDIAN CONVERGENCE AT POINT 2.(IN RADIANS) C LSK - LINE SCALE FACTOR. C C OUTPUT: C C GDIST - GEODESIC DISTANCE. C GEODETIC AZIMUTH (IN RADIANS). C C WRITTEN BY G.BOWIE, DEC. 1977. C IMPLICIT REAL*8 (A-Z) MPDIST=DSQRT((X2-X1)**2+(Y2-Y1)**2) MPAZ=DATAN2((X2-X1),(Y2-Y1)) MPAZ2=DATAN2((X1-X2),(Y1-Y2)) GDIST=MPDIST/LSK GAZ12=MPAZ+TT1+MC1 GAZ21=MPAZ2+TT2+MC2 RETURN END SUBROUTINE ERRMPD(MPAZ,MPDIST,VCINIT,VCDIR) C C THIS SUBROUTINE WILL COMPUTE THE VARIANCE COVARIANCE MATRIX OF THE C X,Y COORDINATES OF POINTS 1 AND 2 FROM THE DIRECT SOLUTION OF GEODETIC C POSITIONING ON A MAPPING PLANE. C C INPUT: C C MPAZ - AZIMUTH ON THE MAPPING PLANE (IN RADIANS). C MPDIST - DISTANCE FROM POINT 1 TO 2 ON THE MAPPING PLANE. C VCINIT - VARIANCE COVARIANCE MATRIX OF POINT 1 AND THE C OBSERVATIONS. C C OUTPUT: C C VCDIR - VARIANCE COVARIANCE MATRIX OF THE X,Y COORDINATES OF C POINTS 1 AND 2. C C WRITTEN BY G.BOWIE, DEC. 1977. C IMPLICIT REAL*8 (A-Z) REAL*8 B(4,4)/16*0.D0/,BVC(4,4),BT(4,4),VCDIR(4,4),VCINIT(4,4) B(1,1)=1.D0 B(2,2)=1.D0 B(3,1)=1.D0 B(3,3)=DSIN(MPAZ) B(3,4)=DCOS(MPAZ)*MPDIST B(4,2)=1.D0 B(4,3)=DCOS(MPAZ) B(4,4)=-DSIN(MPAZ)*MPDIST CALL TRNSD(BT,4,B,4,4,4) CALL MMULD(BVC,4,B,4,VCINIT,4,4,4,4) CALL MMULD(VCDIR,4,BVC,4,BT,4,4,4,4) PRINT 1 1 FORMAT(' DESIGN MATRIX-DIRECT CASE') CALL MOUTD(B,4,4,4) PRINT 2 2 FORMAT(' VC MATRIX-INITIAL') CALL MOUTD(VCINIT,4,4,4) PRINT 3 3 FORMAT(' VC MATRIX-DIRECT CASE') CALL MOUTD(VCDIR,4,4,4) RETURN END SUBROUTINE ERRMPI(MPDIST,XI,YI,XJ,YJ,VCDIR,VCINV) C C THIS SUBROUTINE WILL COMPUTE THE VARIANCE COVARIANCE MATRIX OF THE C MAPPING PLANE AZIMUTH AND DISTANCE SOLVED FOR IN THE INVERSE GEODETIC C POSITIONING PROBLEM ON A MAPPING PLANE. C C INPUT: C C MPDIST - DISTANCE ON THE MAPPING PLANE. C XI,YI,XJ,YJ - X,Y COORDINATES OF POINTS 1 AND 2. C VCDIR - VARIANCE COVARIANCE MATRIX OF THE X,Y COORDINATES OF C POINTS 1 AND 2. C C OUTPUT: C C VCINV - VARIANCE COVARIANCE MATRIX OF MAPPING PLANE AZIMUTH AND C DISTANCE. C C WRITTEN BY G.BOWIE, DEC. 1977. C IMPLICIT REAL*8 (A-Z) REAL*8 B(2,4),BVC(2,4),BT(4,2),VCDIR(4,4),VCINV(2,2) B(1,1)=(YI-YJ)/MPDIST**2 B(1,2)=(XJ-XI)/MPDIST**2 B(1,3)=(YJ-YI)/MPDIST**2 B(1,4)=(XI-XJ)/MPDIST**2 B(2,1)=(XI-XJ)/MPDIST B(2,2)=(YI-YJ)/MPDIST B(2,3)=(XJ-XI)/MPDIST B(2,4)=(YJ-YI)/MPDIST CALL TRNSD(BT,4,B,2,2,4) CALL MMULD(BVC,2,B,2,VCDIR,4,2,4,4) CALL MMULD(VCINV,2,BVC,2,BT,4,2,4,2) PRINT 1 1 FORMAT(' DESIGN MATRIX-INVERSE CASE') CALL MOUTD(B,2,2,4) PRINT 2 2 FORMAT(' VC MATRIX-INVERSE CASE') CALL MOUTD(VCINV,2,2,2) RETURN END SUBROUTINE RADMS(RAD,IDEG,IMIN,SEC) C C SUBROUTINE TO CONVERT FROM RADIANS TO DEGREES C C ***INPUT*** C C RAD----ANGLE (RADIANS) C C ***OUTPUT*** C C IDEG----DEGREES(INTEGER) C IMIN----MINUTES (INTEGER) C SEC----SECONDS (REAL*8) C IMPLICIT REAL*8(A-H,O-Z) PI=DARCOS(-1.0D0) A=RAD*180.D0/PI IDEG=A B=A-IDEG A=B*60.0D0 IMIN=A B=A-IMIN SEC=B*60.0D0 RETURN END SUBROUTINE DMSRAD(IDEG,IMIN,SEC,RAD) C C SUBROUTINE TO CONVERT FROM DEGREES TO RADIANS C C ***INPUT*** C C IDEG----DEGREES(INTEGER) C IMIN----MINUTES (INTEGER) C SEC----SECONDS (REAL*8) C C ***OUTPUT*** C C RAD----ANGLE (RADIANS) IMPLICIT REAL*8(A-H,O-Z) PI=DARCOS(-1.0D0) SIGN=1.0D0 IF(IDEG.LT.0)SIGN=-1.0D0 RAD=(IABS(IDEG)+IMIN/DFLOAT(60)+SEC/3600.D0)*PI/180.0D0 RAD=RAD*SIGN RETURN END SUBROUTINE TMPLXY(PHI,DLAM,A,B,SF,XO,CMRAD,X,Y) C C C C SUBROUTINE TMPLXY COMPUTES THE X,Y COORDINATES FOR THE C TRANSVERSE MERCATOR PROJECTION GIVEN THE GEOGRAPHIC COORDINATES C -LATITUDE AND LONGITUDE. THE EQUATIONS USED TO COMPUTE X AND Y C ARE FROM THOMAS (1952). SUBROUTINE MERARC IS USED TO COMPUTE C THE MERIDIAN ARC LENGTH. C C C INPUT: C C PHI -LATITUDE IN RADIANS C C DLAM-LONGITUDE OF POINT MINUS LONGITUDE OF CENTRAL C MERIDIAN (IN RADIANS) FOR LONGITUDE POSITIVE EAST. C C A -SEMI-MAJOR AXES OF THE REFERENCE ELLIPSOID. C C B -SEMI-MINOR AXES OF THE REFERENCE ELLIPSOID. C C XO - FALSE EASTING OF THE CENTRAL MERIDIAN. C C SF - SCALE OF THE CENTRAL MERIDIAN. C C CMRAD - THE CENTRAL MERIDIAN,IN RADIANS. C C OUTPUT: C C X -EASTING COORDINATE OF THE TRANSVERSE MERCATOR C PROJECTION. C C Y -NORTHING COORDINATE OF THE TRANSVERSE MERCATOR C PROJECTION. C C WRITTEN BY R.R.STEEVES. C MAY,1977 C C IMPLICIT REAL*8(A-H,O-Z) SP=DSIN(PHI) CP=DCOS(PHI) T=DTAN(PHI) E=DSQRT((A**2-B**2)/A**2) ETA=DSQRT((A**2-B**2)/B**2*CP**2) CALL MERARC(PHI,A,B,SPHI) DN=A/DSQRT(1.0D0-E**2*SP**2) X=DN*(DLAM*CP+DLAM**3*CP**3/6.0D0*(1.0D0-T**2+ETA**2)+DLAM**5*CP 1 **5/120.0D0*(5.0D0-18.0D0*T**2+T**4+14.0D0*ETA**2-58.0D0*T**2* 2 ETA**2+13.D0*ETA**4+4.D0*ETA**6-64.D0*ETA**4*T**2-24.D0*ETA**6 3 *T**2) +DLAM**7/5040.0D0*CP**7*(61.0D0-479.0D0*T**2+179.0D0* 4 T**4-T**6)) Y=SPHI+DN* (DLAM**2/2.0D0*SP*CP+DLAM**4/24.0D0*SP*CP**3*(5.0D0- 1 T**2+9.0D0*ETA**2+4.0D0*ETA**4)+DLAM**6/720.0D0*SP*CP**5* 2 (61.0D0-58.0D0*T**2+T**4+270.0D0*ETA**2-330.0D0*T**2*ETA**2 3 +445.D0*ETA**4+324.D0*ETA**6-680.D0*ETA**4*T**2+88.D0*ETA**8 4 -600.D0*ETA**6*T**2-192.D0*ETA**8*T**2)+DLAM**8/40320.D0*SP*CP* 5 *7*(1385.D0-3111.D0*T**2+543.D0*T**4-T**6)) X=SF*X+XO Y=SF*Y RETURN END SUBROUTINE TMXYPL(X,Y,A,B,SF,XO,CMRAD,PHI,OLAM) C C C SUBROUTINE TMXYPL COMPUTES THE GEOGRAPHIC COORDINATES- LATITUDE C AND LONGITUDE - GIVEN THE X,Y COORDINATES OF THE TRANSVERSE C MERCATOR PROJECTION. THE EQUATIONS USED TO COMPUTE THE LONGITUDE C AND LATITUDE ARE FROM THOMAS (1952). SUBROUTINE FPLAT IS USED C TO COMPUTE THE FOOT-POINT LATITUDE. C C C INPUT; C C X -EASTING COORDINATE OF THE TRANSVERSE MERCATOR C PROJECTION. C C Y -NORTHING COORDINATE OF THE TRANSVERSE MERCATOR C PROJECTION. C C A -SEMI-MAJOR AXES OF THE REFERENCE ELLIPSOID. C C B -SEMI-MINOR AXES OF THE REFERENCE ELLIPSOID. C C SF - SCALE OF THE CENTRAL MERIDIAN. C C XO - FALSE EASTING OF THE CENTRAL MERIDIAN. C C CMRAD - THE CENTRAL MERIDIAN,IN RADIANS. C C OUTPUT: C C PHI -LATITUDE OF THE POINT IN RADIANS C C OLAM-LONGITUDE OF THE POINT IN RADIANS C C WRITTEN BY R.R.STEEVES. C MAY,1977 C C IMPLICIT REAL*8(A-H,O-Z) X=(X-XO)/SF Y=Y/SF E=DSQRT((A**2-B**2)/A**2) CALL FPLAT(A,B,Y,PHI1) T=DTAN(PHI1) SP=DSIN(PHI1) CP=DCOS(PHI1) ETA=DSQRT((A**2-B**2)/B**2*CP**2) DN=A/DSQRT(1.0D0-E**2*SP**2) DM=A*(1.0D0-E**2)/DSQRT((1.0D0-E**2*SP**2)**3) PHI=PHI1-T*X**2/2.0D0/DM/DN+T*X**4/24.0D0/DM/DN**3*(5.0D0+3.0D0* 1 T**2+ETA**2-4.0D0*ETA**4-9.0D0*ETA**2*T**2)-T*X**6/720.0D0/DM/ 2 DN**5*(61.0D0+90.0D0*T**2+46.0D0*ETA**2+45.0D0*T**4-252.0D0*T** 3 2*ETA**2-3.0D0*ETA**4+100.0D0*ETA**6-66.0D0*T**2*ETA**4-90.0D0 4 *T**4*ETA**2+88.0D0*ETA**8+225.0D0*T**4*ETA**4+84.0D0*T**2* 5 ETA**6-192.0D0*T**2*ETA**8) PHI=PHI+T*X**8/40320.0D0/DM/DN**7*(1385.0D0+3633.0D0*T**2+4095.0D0 1 *T**4+1575.0D0*T**6) DLAM=(X/DN-(X/DN)**3/6.0D0*(1.0D0+2.0D0*T**2+ETA**2)+(X/DN)**5/ 1 120.0D0*(5.0D0+6.0D0*ETA**2+28.0D0*T**2-3.0D0*ETA**4+8.0D0*T**2 2 *ETA**2+24.0D0*T**4-4.0D0*ETA**6+4.0D0*T**2*ETA**4+24.0D0*T**2* 3 ETA**6)-(X/DN)**7/5040.0D0*(61.0D0+662.0D0*T**2+1320.0D0*T**4+ 4 720.0D0*T**6))/CP OLAM=CMRAD+DLAM X=X*SF+XO Y=Y*SF RETURN END SUBROUTINE TTLS(A,B,PHI1,DIST,X1,Y1,X2,Y2,XO,SF,TT,LS) C C THIS SUBROUTINE WILL COMPUTE THE T-T AND LINE SCALE CORRECTIONS IN C A TRANSVERSE MERCATOR PROJECTION. C C INPUT: C A,B - SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE ELLIPSOID. C PHI1 - GEODETIC LATITUDE OF POINT 1. C DIST - GEODESIC DISTANCE FROM POINT 1 TO 2. C X1,Y1,X2,Y2 - X,Y,COORDINATES OF POINTS 1 AND 2. C X1,Y1,X2,Y2 - X,Y COORDINATES OF POINTS 1 AND 2. C SF - SCALE FACTOR AT THE CENTRAL MERIDIAN. C C OUTPUT: C C TT - T-T CORRECTION FROM POINT 1 TO 2. C LS - LINE SCALE FOR THE DISTANCE FROM 1 TO 2. C C WRITTEN BY G.BOWIE, NOV. 1977. C IMPLICIT REAL*8(A-Z) PI=3.141592653589793D0 X1=X1-XO X2=X2-XO CP=DCOS(PHI1) SP=DSIN(PHI1) E=DSQRT((A**2-B**2)/A**2) DN=A/DSQRT(1.0D0-E**2*SP**2) DM=A*(1.0D0-E**2)/DSQRT((1.0D0-E**2*SP**2)**3) N1=DSQRT(DM*DN) DEN=6.D0*(N1**2) TT=(Y2-Y1)*(X2+2*X1)/DEN*(1.0D0-(2.0D0*X1+X2)**2/(3.5D0*DEN)) P=X1**2+X1*X2+X2**2 LS=1.0D0+P/DEN*(1.0D0+P/(6.0D0*DEN)) LS=LS*SF X1=X1+XO X2=X2+XO RETURN END SUBROUTINE TMSFMC(PHI,DLAM,SFO,A,B,SF,C) C C C THIS ROUTINE COMPUTES THE POINT SCALE FACTOR AND MERIDIAN C CONVERGENCE (FOR A POINT DEFINED BY PHI,DLAM) FOR A TRANSVERSE C MERCATOR PROJECTION DEFINED BY THE SCALE FACTOR SFO AT THE C CENTRAL MERIDIAN. C C INPUT: C C PHI - ELLIPSOIDAL LATITUDE OF THE POINT, IN RADIANS. C DLAM - ELLIPSOIDAL LONITUDE OF THE POINT MINUS THE C ELLIPSOIDAL LONGITUDE OF THE CENTRAL MERIDIAN OF C PROJECTION, (LONGITUDE POSITIVE EAST), IN RADIANS. C SFO - SCALE AT THE CENTRAL MERIDIAN. C A,B - SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE C ELLIPSOID RESPECTIVELY, IN METRES. C C OUTPUT: C C SF - POINT SCALE AT THE POINT. C C - MERIDIAN CONVERGENCE AT THE POINT, IN RADIANS. C C IMPLICIT REAL*8(A-Z) CP=DCOS(PHI) T=DTAN(PHI) ETA=DSQRT((A*A-B*B)/(B*B)*CP**2) C=DLAM*DSIN(PHI)*(1.D0+DLAM**2*CP**2/3.D0*(1.D0+3.D0*ETA**2+2.D0* 1 ETA**4)+DLAM**4*CP**4/15.D0*(2.D0-T**2)) SF=1.D0+DLAM**2*CP**2/2.D0*(1.D0+ETA**2)+DLAM**4*CP**4/24.D0*(5.D0 1 -4.D0*T**2) RETURN END SUBROUTINE MERARC(PHI,A,B,S) C C C THIS ROUTINE COMPUTES THE MERIDIAN ARC LENGTH FROM THE EQUATOR C TO LATITUDE PHI ON AN ELLIPSOID OF REVOLUTION DEFINED BY ITS C SEMI-MAJOR AXIS A AND ITS SEMI-MINOR AXIS B. THE COMPUTED ARC C LENGTH IS ACCURATE TO APPROXIMATELY 10 MICROMETRES OVER THE C ENTIRE RANGE: EQUATOR TO POLE. C C INPUT: PHI-ELLIPSOIDAL LATITUDE IN RADIANS.(MAY BE POSITIVE C (NORTH) OR NEGATIVE(SOUTH OF EQUATOR)) C C A,B-SEMI-MAJOR AND SEMI-MINOR AXES OF THE ELLIPSOID. C (THE COMPUTED ARC LENGTH WILL BE OUTPUT IN THE SAME C UNITS AS A AND B). C C OUTPUT: S -MERIDIAN ARC LENGTH C C C WRITTEN BY R.R.STEEVES AND G.BOWIE C MAY,1977 C IMPLICIT REAL*8(A-Z) E2=(A*A-B*B)/(A*A) E4=E2*E2 E6=E4*E2 E8=E6*E2 A0=1.0D0-E2/4.D0-3.D0*E4/64.D0-5.D0*E6/256.D0-175.D0*E8/16384.D0 A2=3.D0/8.D0*(E2+E4/4.D0+15.D0*E6/128.D0-455.D0*E8/4096.D0) A4=15.D0/256.D0*(E4+3.D0*E6/4.D0-77.D0*E8/128.D0) A6=35.D0/3072.D0*(E6-41.D0*E8/32.D0) A8=-315.D0*E8/131072.D0 S=A*(A0*PHI-A2*DSIN(2.D0*PHI)+A4*DSIN(4.D0*PHI)-A6*DSIN(6.D0*PHI) 1 +A8*DSIN(8.D0*PHI)) RETURN END SUBROUTINE FPLAT(A,B,Y,PHI1) C C C THIS ROUTINE COMPUTES THE FOOT-POINT LATITUDE REQUIRED IN C TRANSFORMING TRANSVERSE MERCATOR PLANE COORDINATES X,Y TO C ELLIPSOIDAL COORDINATES. C C C INPUT: C C A - SEMI-MAJOR AXES OF THE REFERENCE ELLIPSOID. C C B - SEMI-MINOR AXES OF THE REFERENCE ELLIPSOID. C C Y - NORTHING OF THE POINT FOR WHICH THE FOOT-POINT C LATITUDE IS TO BE COMPUTED. C C OUTPUT: C C PHI1 - FOOT-POINT LATITUDE IN RADIANS. C C WRITTEN BY R.R.STEEVES C JUNE,1977 C C IMPLICIT REAL*8(A-Z) F(PHI)=A*(A0*PHI-A2*DSIN(2.D0*PHI)+A4*DSIN(4.D0*PHI)-A6*DSIN(6.D0* 1 PHI)+A8*DSIN(8.D0*PHI))-Y FP(PHI)=A*(A0-2.D0*A2*DCOS(2.D0*PHI)+4.D0*A4*DCOS(4.D0*PHI)-6.D0* 1 A6*DCOS(6.D0*PHI)+8.D0*A8*DCOS(8.D0*PHI)) E2=(A*A-B*B)/(A*A) E4=E2*E2 E6=E4*E2 E8=E6*E2 A0=1.D0-E2/4.D0-3.D0*E4/64.D0-5.D0*E6/256.D0-175.D0*E8/16384.D0 A2=3.D0/8.D0*(E2+E4/4.D0+15.D0*E6/128.D0-455.D0*E8/4096.D0) A4=15.D0/256.D0*(E4+3.D0*E6/4.D0-77.D0*E8/128.D0) A6=35.D0/3072.D0*(E6-41.D0*E8/32.D0) A8=-315.D0*E8/131072.D0 PHI1=Y/A 1 DPHI=F(PHI1)/FP(PHI1) PHI1=PHI1-DPHI IF(DABS(DPHI).LT.1.D-11)GOTO 2 GO TO 1 2 CONTINUE RETURN END SUBROUTINE ERRTM(A,B,PHI,DLAM,KNOT,ICODE,CM,DM) C C C THIS SUBROUTINE COMPUTES THE COVARIANCE MATRIX CM OF THE C TRANSVERSE MERCATOR COORDINATES X,Y GIVEN THE COVARIANCE MATRIX C DM OF THE ELLIPSOIDAL COORDINATES PHI,LAMDA. (ICODE=1) C IT WILL ALSO COMPUTE THE COVARIANCE MATRIX DM OF PHI,LAMDA C GIVEN THE COVARIANCE MATRIX CM OF X,Y. (IF ICODE=-1) C C C INPUT: C C PHI - ELLIPSOIDAL LATITUDE. (IN RADIANS) C C DLAM - LONGITUDE OF THE POINT MINUS THE LONGITUDE OF THE C CENTRAL MERIDIAN (IN RADIANS) FOR LONGITUDE C POSITIVE EAST. C C C KNOT - SCALE FACTOR AT THE CENTRAL MERIDIAN. C C ICODE - 1 PHI,LAMDA TO X,Y C C -1 X,Y TO PHI,LAMDA. C C CM - COVARIANCE MATRIX OF THE TRANSVERSE MERCATOR C COORDINATES . (IN METRES SQUARED) C C OR C C DM - COVARIANCE MATRIX OF THE ELLIPSOIDAL C COORDINATES. (IN RADIANS SQUARED) C C OUTPUT: C C CM - COVARIANCE MATRIX OF THE TRANSVERSE MERCATOR C COORDINATES . (IN METRES SQUARED) C C OR C C DM - COVARIANCE MATRIX OF THE ELLIPSOIDAL C COORDINATES. (IN RADIANS SQUARED) C C C WRITTEN BY R.R.STEEVES C JUNE,1977 C IMPLICIT REAL*8(A-Z) INTEGER ICODE DIMENSION CM(2,2),DM(2,2) CP=DCOS(PHI) SP=DSIN(PHI) C C COMPUTE THE JACOBIAN MATRIX C J11=21670.D0*DLAM*CP-639.D4*DLAM*SP J12=639.D4*CP J21=637.D4 J22=639.D4*DLAM*SP*CP C C IF ICODE=1 CM IS COMPUTED; IF ICODE=-1 DM IS COMPUTED C IF(ICODE.LT.0)GO TO 1 C C COMPUTE THE COVARIANCE MATRIX OF X,Y C CM(1,1)=J11**2*DM(1,1)+2.D0*J11*J12*DM(1,2)+J12**2*DM(2,2) CM(1,2)=J11*J21*DM(1,1)+J12*J21*DM(1,2)+J11*J22*DM(1,2)+J12*J22*DM 1 (2,2) CM(2,1)=CM(1,2) CM(2,2)=J21**2*DM(1,1)+2.D0*J21*J22*DM(1,2)+J22**2*DM(2,2) GO TO 2 1 CONTINUE C C INVERT THE JACOBIAN MATRIX C J1=J11 DET=J11*J22-J21*J12 J11=J22/DET J22=J1/DET J12=-J12/DET J21=-J21/DET C C COMPUTE THE COVARIANCE MATRIX OF PHI,LAMDA C DM(1,1)=J11**2*CM(1,1)+2.D0*J11*J12*CM(1,2)+J12**2*CM(2,2) DM(1,2)=J11*J21*CM(1,1)+J12*J21*CM(1,2)+J11*J22*CM(1,2)+J12*J22*CM 1 (2,2) DM(2,1)=DM(1,2) DM(2,2)=J21**2*CM(1,1)+2.D0*J21*J22*CM(1,2)+J22**2*CM(2,2) 2 RETURN END