C C C THIS PROGRAM WILL SOLVE THE DIRECT AND INVERSE PROBLEMS OF GEODETIC C POSITION COMPUTATIONS IN THE 3-DIMENSIONAL OR AVERAGE TERRESTRIAL C COORDINATE SYSTEM. ALSO IT WILL COMPUTE THE RESULTING VARIANCE COVARIANCE C MATRIX ASSOCIATED WITH EACH CASE. C C NOTE: ALL ANGULAR DATA ENTERED INTO THE MAIN PROGRAM IS IN DEGREES, C MINUTES AND SECONDS. C ALL ANGULAR DATA ENTERED INTO THE SUBROUTINES IS IN RADIANS. C ALL ANGULAR DATA FROM THE SUBROUTINES IS IN RADIANS. C ALL ANGULAR DATA IN THE OUTPUT IS IN DEGREES, MINUTES AND SECONDS C C G.BOWIE, FEB. 1978. C IMPLICIT REAL*8(A-H,O-Z) REAL*8 L1,L,CM(3,3),DM(3,3),VCI(3,3),VC2(6,6),VCP2(3,3),JD(6,6), 1 JI(3,6),DM2(3,3) REAL*8 VCD(6,6)/36*0.D0/ 21 FORMAT(15X,'PHI',15X,'LAMDA',15X,'HEIGHT',//) 22 FORMAT(10X,2I3, F9.5,3X,2I3,F9.5,3X,F14.4) 23 FORMAT(10X,'DELTA AZIMUTH',5X,I3,I3,F6.2/////) 24 FORMAT('1', 'VARIANCE COVARIANCE MATRIX OF POINT 1 IN TERMS OF 1PHI,LAMDA AND HEIGHT. UNITS ARE SEC**2,SEC**2 AND M**2') 25 FORMAT('0', 'VARIANCE COVARIANCE MATRIX OF POINT 2 IN TERMS OF 1PHI,LAMDA AND HEIGHT. UNITS ARE SEC**2, SEC**2 AND M**2') 30 FORMAT('1',30X,'DIRECT PROBLEM OF GEODETIC POSITIONING'/31X,'_____ 1_ _______ __ ________ ___________') 31 FORMAT('0',30X,'OBSERVATIONAL DATA'////10X, 'GEODETIC LATITUDE', 1I6,I3,F5.1,10X, 'GEODETIC LONGITUDE',I6,I3,F5.1// 210X,'ELLIPSOIDAL HEIGHT',3X,F6.2// 310X,'COMPONENTS OF DEFLECTION OF THE VERTICAL XI',F6.1,5X,'ETA 4',F6.1//10X,'SPATIAL DISTANCE',F12.3,5X,'AZIMUTH',I5,I3,F5.1,5X, 5'ZENITH DISTANCE',I5,I3,F5.1/////) 32 FORMAT('0',30X,'COMPUTED INFORMATION'////10X,'GEODETIC COORDINATES 1 OF THE OBSERVED POINT'//21X,'X',16X,'Y',16X,'Z'//10X,3F17.4/////) 33 FORMAT('1',30X,'INVERSE PROBLEM OF GEODETIC POSITIONING'/31X,'____ 1___ _______ __ ________ ___________') 34 FORMAT('0',30X,'OBSERVATIONAL DATA'////10X, 'GEODETIC LATITUDE', 1I6,I3,F5.1,10X, 'GEODETIC LONGITUDE',I6,I3,F5.1// 210X,'ELLIPSOIDAL HEIGHT',3X,F6.2// 310X,'COMPONENTS OF DEFLECTION OF THE VERTICAL XI',F6.1,5X,'ETA 4',F6.1//10X,'GEODETIC COORDINATES OF POINT 1',28X,'GEODETIC COORDI 4NATES OF POINT 2' 5//16X,'X',14X,'Y',14X,'Z',19X,'X',14X,'Y',14X,'Z'//6X,3F15.3,5X, 6 3F15.3/////) 35 FORMAT('0',30X,'COMPUTED INFORMATION',////10X,'SPATIAL DISTANCE',F 112.3, 3X// 2 10X,'AZIMUTH FROM 1 TO 2',I5,I3,F6.2,3X,'AZIMUTH FROM 2 TO 1' 2,I5,I3,F6.2,//10X,'ZENITH DISTANCE FROM 1 TO 2',I5,I3,F6.2,3X,'ZEN 3ITH DISTANCE FROM 2 TO 1',I5,I3,F6.2) 36 FORMAT('0',10X,'DESIGN MATRIX - DIRECT CASE.') 37 FORMAT('0',10X,'DESIGN MATRIX - INVERSE CASE.') 39 FORMAT('1',20X,'VARIANCE COVARIANCE INFORMATION - DIRECT CASE',//) 40 FORMAT('1',20X,'VARIANCE COVARIANCE INFORMATION - INVERSE CASE',/) 41 FORMAT('0',10X,'VARIANCE COVARIANCE MATRIX OF POINT 1 IN METRES**2 1') 42 FORMAT('0',10X,'VARIANCE COVARIANCE MATRIX OF DIRECT PROBLEM - X, 1Y, Z POINT 1'/58X,'DISTANCE, AZIMUTH AND ZENITH ANGLE (GIVEN)') 43 FORMAT('0',10X,'VARIANCE COVARIANCE MATRIX OF POINT 1 AND POINT 2 1IN METRES**2') 44 FORMAT('0',10X,'VARIANCE COVARIANCE MATRIX OF DISTANCE, AZIMUTH AN 1D ZENITH ANGLE IN METRES**2 AND SEC**2') 45 FORMAT('1') PI=3.141592653589793D0 C C READ IN DATA - ELLIPSOIDAL PARAMETERS, MAP PROJETION PARAMETERS, C POSITION COORDINATES AND OBSERVATIONAL DATA. C C C IF GEODETIC X,Y,Z ARE REQUIRED INPUT: XO=0.0, YO=0.0, ZO=0.0 C C C IF AVERAGE TERRESTRIAL X,Y,Z, ARE REQUIRED INPUT TRANSLATION COMPONENTS. C (XO,YO,ZO) C READ(5,*)A,B READ(5,*)XO,YO,ZO READ(5,*)IPD,IPM,PS,LD,LM,SL,H1 READ(5,*)XXII,ETA READ(5,*)R,IAD,IAM,ASEC,IVD,IVM,VSEC C C READ IN VARIANCE COVARIANCE INFORMATION OF POINT I AND STANDARD C DEVIATIONS OF THE OBSERVATIONS. C READ(5,*)DM(1,1),DM(1,2),DM(1,3),DM(2,2),DM(2,3),DM(3,3) READ(5,*)SD,SA,SZ C C CONVERT INFORMATION TO PROPER UNITS (RADIANS FOR THE OBSERVATIONS AND C RADIANS SQUARED OR RADIAN-METRE FOR THE COVARIANCE INFORMATION). C DM(2,1)=DM(1,2) DM(3,1)=DM(1,3) DM(3,2)=DM(2,3) DO 11 I=1,2 DM(I,3)=(DM(I,3)*PI/180.D0)/36.D2 DM(3,I)=(DM(3,I)*PI/180.D0)/36.D2 DO 11 J=1,2 DM(I,J)=(DM(I,J)/(3600.D0**2))*(PI**2)/(180.D0**2) 11 CONTINUE CALL DMSRAD(IPD,IPM,PS,P1) CALL DMSRAD(LD,LM,SL,L1) CALL DMSRAD(IAD,IAM,ASEC,AZ) CALL DMSRAD(IVD,IVM,VSEC,VERT) ZE=PI/2.D0-VERT C C PRINT INITIAL DATA. C PRINT 30 PRINT 31,IPD,IPM,PS,LD,LM,SL,H1,XXII,ETA,R,IAD,IAM,ASEC,IVD,IVM,VS 1EC XXII=(XXII/36.D2)*PI/18.D1 ETA=(ETA/36.D2)*PI/18.D1 C C COMPUTE X,Y,Z FOR POINT J. C CALL GPLAGE(P1,L1,H1,XO,YO,ZO,A,B,ZE,AZ,R,XXII,ETA,XJ,YJ,ZJ,XI,YI, 1 ZI,DAZ) CALL RADMS(DAZ,IDGD,IDGM,DGS) PRINT 23,IDGD,IDGM,DGS C C FORM VARIANCE COVARIANCE MATRIX FOR THE DIRECT CASE. C P=P1 L=L1 H=H1 I=1 CALL ERR3D(P,L,H,A,B,I,CM,DM) DO 10 I=1,3 DO 10 J=1,3 VCD(I,J)=CM(I,J) 10 CONTINUE VCD(4,4)=SD**2 VCD(5,5)=(SA/3600.D0*PI/180.D0)**2 VCD(6,6)=(SZ/3600.D0*PI/180.D0)**2 C C COMPUTE VARIANCE COVARIANCE MATRIX OF POINT J. C CALL ERRD3D(XI,YI,ZI,XJ,YJ,ZJ,P,L,R,DAZ,XXII,ETA,AZ,ZE,VCD,VC2,JD) C C PRINT COMPUTED INFORMATION X,Y,Z POINT J AND VARIANCE COVARIANCE C INFORMATION. C PRINT 32,XJ,YJ,ZJ CALL XYZPLH(XJ,YJ,ZJ,XO,YO,ZO,A,B,PHI,RLAM,H2) CALL RADMS (PHI,IPID,IPIM,SECP) CALL RADMS (RLAM,IRD,IRM,SECR) PRINT 21 PRINT 22,IPID,IPIM,SECP,IRD,IRM,SECR,H2 PRINT 39 PRINT 41 CALL MOUTD(CM,3,3,3) PRINT 36 CALL MOUTD(JD,6,6,6) PRINT 42 VCD(5,5)=VCD(5,5)*18.D1**2/PI**2*36.D2**2 VCD(6,6)=VCD(6,6)*18.D1**2/PI**2*36.D2**2 CALL MOUTD(VCD,6,6,6) PRINT 43 CALL MOUTD(VC2,6,6,6) C C ***** INVERSE CASE ***** C C GIVEN X,Y,Z POINTS I,J COMPUTE DISTANCE AZIMUTH AND ZENITH (OR VERTICAL) C ANGLE. C CALL GPGELA(XI,YI,ZI,XJ,YJ,ZJ,P,L,DAZ,XXII,ETA,DIJ,AIJ,AJI,ZIJ,ZJI 1,DXLG,DYLG,DZLG,DXLA,DYLA,DZLA) C C COMPUTE VARIANCE COVARIANCE MATRIX OF DISTANCE,AZIMUTH AND ZENITH ANGLE. C CALL ERRI3D(XI,YI,ZI,XJ,YJ,ZJ,VC2,DXLG,DYLG,DZLG,DXLA,DYLA,DZLA, 1 DIJ,AIJ,ZIJ,DAZ,XXII,ETA,P,L,VCI,JI) C C CONVERT TO THE PROPER UNITS C CALL RADMS(AIJ,IAD,IAM,ASEC) CALL RADMS(AJI,JAD,JAM,AJSEC) CALL RADMS(ZIJ,IZD,IZM,ZSEC) CALL RADMS(ZJI,JZD,JZM,ZJSEC) XXII=XXII*3600.D0*180.D0/PI ETA = ETA*3600.D0*180.D0/PI C C PRINT OUT THE COMPUTED INFORMATION AND VARIANCE COVARIANCE INFORMATION C PRINT 33 PRINT 34,IPD,IPM,PS,LD,LM,SL,H1,XXII,ETA,XI,YI,ZI,XJ,YJ,ZJ PRINT 35,DIJ,IAD,IAM,ASEC,JAD,JAM,AJSEC,IZD,IZM,ZSEC,JZD,JZM,ZJSEC PRINT 40 PRINT 37 CALL MOUTD(JI,3,3,6) PRINT 43 CALL MOUTD(VC2,6,6,6) PRINT 44 DO 14 I=2,3 VCI(1,I)=VCI(1,I)*18.D1/PI*36.D2 VCI(I,1)=VCI(I,1)*18.D1/PI*36.D2 DO 14 J=2,3 VCI(I,J)=VCI(I,J)*18.D1**2/PI**2*36.D2**2 14 CONTINUE CALL MOUTD(VCI,3,3,3) DO 12 I=1,3 DO 12 J=1,3 VCP2(I,J)=VC2(3+I,3+J) 12 CONTINUE I=-1 CALL ERR3D(PHI,RLAM,H2,A,B,I,VCP2,DM2) DO 13 I=1,2 DM(I,3)= DM(I,3)*18.D1/PI*36.D2 DM2(I,3)=DM2(I,3)*18.D1/PI*36.D2 DM(3,I)= DM(3,I)*18.D1/PI*36.D2 DM2(3,I)=DM2(3,I)*18.D1/PI*36.D2 DO 13 J=1,2 DM(I,J)= DM(I,J)*18.D1**2/PI**2*36.D2**2 DM2(I,J)=DM2(I,J)*18.D1**2/PI**2*36.D2**2 13 CONTINUE PRINT 24 CALL MOUTD(DM,3,3,3) PRINT 25 CALL MOUTD(DM2,3,3,3) PRINT 45 STOP 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 GPLAGE(P,L,H,XO,YO,ZO,A,B,ZE,AZ,R,XI,ETA,X,Y,Z,XP,YP,ZP 1,DA) C C C THIS SUBROUTINE SOLVES THE DIRECT PROBLEM IN 3 DIMENSIONAL GEODETIC C POSITIONING.GIVEN THE GEODETIC LATITUDE AND LONGITUDE, HEIGHT AND C THE COMPONENTS OF THE DEFLECTION OF THE VERTICALFOR POINT 1 AS WELL C AS THE OBSERVED DATA IT WILL COMPUTE X,Y,Z OF POINT 2. C C INPUT: C P,L,H - ELLIPSOIDAL COORDINATES OF POINT 1. C XO,YO,ZO - TRANSLATION COMPONENTS OF THE ORIGIN OF THE C GEODETIC SYSTEM. C A,B - SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE C ELLIPSOID. C ZE,AZ,R - OBSERVED DATA - ZE - ZENITH DISTANCE C AZ - AZIMUTH C R DISTANCE C XI,ETA - DEFLECTION OF THE VERTICAL COMPONENTS. C OUTPUT: C P,L,H - GEOGRAPHICAL COORDINATES OF POINT 1. C XP,YP,ZP - CARTESIAN COORDINATES OF POINT 1. C X,Y,Z - CARTESIAN COORDINATES OF POINT 2. C DA - ASTRONOMIC AZIMUTH MINUS GEODETIC AZIMUTH. C C WRITTEN BY G.BOWIE,DEC.,1977. C IMPLICIT REAL*8(A-Z) PI=3.141592653589793D0 CALL PLHXYZ(P,L,H,XO,YO,ZO,A,B,XP,YP,ZP) CP=DCOS(P) SP=DSIN(P) SL=DSIN(L) CL=DCOS(L) SZ=DSIN(ZE) CZ=DCOS(ZE) SA=DSIN(AZ) CA=DCOS(AZ) DA=ETA*DTAN(P)-(XI*SA-ETA*CA)*DCOTAN(ZE) X=XP-R*(SP*CL*(SZ*CA+DA*SZ*SA+XI*CZ)+SL*(-DA*SZ*CA+SZ*SA+ETA*CZ) 1 -CP*CL*(-XI*SZ*CA-ETA*SZ*SA+CZ)) Y=YP-R*(SP*SL*(SZ*CA+DA*SZ*SA+XI*CZ)-CL*(-DA*SZ*CA+SZ*SA+ETA*CZ) 1 -CP*SL*(-XI*SZ*CA-ETA*SZ*SA+CZ)) Z=ZP+R*(CP*(SZ*CA+DA*SZ*SA+XI*CZ)+SP*(-XI*SZ*CA-ETA*SZ*SA+CZ)) RETURN END SUBROUTINE XYZPLH(X,Y,Z,XO,YO,ZO,A,B,PHI,RLAM,H) C C THIS ROUTINE COMPUTES THE ELLIPSOIDAL COORDINATES PHI,RLAM,H C GIVEN THE CARTESION COORDINATES X,Y,Z. C C INPUT: C X,Y,Z-CARTESION COORDINATES OF THE POINT IN METRES. C XO,YO,ZO-TRANSLATION COMPONENTS FROM THE ORIGIN OF THE C CARTESIAN COORDINATE SYTEM (X,Y,Z)TO THE CENTER C OF THE REFERENCE ELLIPSOID. (IN METRES.) C A,B-SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE C ELLIPSOID IN METRES. C C OUTPUT: C PHI-ELLIPSOIDAL LATITUDE IN RADIANS. C RLAM-ELLIPSOIDAL LONGITUDE IN RADIANS. C (POSITIVE EAST OF GREENWICH) C H-ELLIPSOIDAL HEIGHT IN METRES. C C WRITTEN BY R.R.STEEVES C JUNE,1977 C C IMPLICIT REAL*8(A-Z) E2=(A*A-B*B)/(A*A) XP=X-XO YP=Y-YO ZP=Z-ZO S=DSQRT(XP**2+YP**2) RLAM=DATAN2(YP,XP) ZPS=ZP/S H=DSQRT(XP**2+YP**2+ZP**2)-A PHI=DATAN(ZPS/(1.D0-E2*A/(A+H))) 1 N=A/DSQRT(1.D0-E2*DSIN(PHI)**2) HP=H PHIP=PHI H=S/ DCOS(PHI)-N PHI=DATAN(ZPS/(1.D0-E2*N/(N+H))) IF(DABS(PHIP-PHI).GT.1.D-11.OR.DABS(HP-H).GT.1.D-5)GO TO 1 RETURN END SUBROUTINE PLHXYZ(PHI,RLAM,H,XO,YO,ZO,A,B,X,Y,Z) C C THIS ROUTINE COMPUTES THE CARTESIAN COORDINATES X,Y,Z GIVEN THE C ELLIPSOIDAL COORDINATES PHI,RLAM,H. C C INPUT: C PHI-ELLIPSOIDAL LATITUDE IN RADIANS. C RLAM-ELLIPSOIDAL LONGITUDE IN RADIANS. C (POSITIVE EAST OF GREENWICH) C H-ELLIPSOIDAL HEIGHT IN METRES. C XO,YO,ZO-TRANSLATION COMPONENTS FROM THE ORIGIN OF THE C CARTESIAN COORDINATE SYTEM (X,Y,Z)TO THE CENTER C OF THE REFERENCE ELLIPSOID. (IN METRES.) C A,B-SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE C ELLIPSOID IN METRES. C C OUTPUT: C X,Y,Z-CARTESIAN COORDINATES OF THE POINT IN METRES. C C WRITTEN BY R.R.STEEVES C JUNE,1977 C C IMPLICIT REAL*8(A-Z) E2=(A*A-B*B)/(A*A) SP=DSIN(PHI) CP=DCOS(PHI) N=A/DSQRT(1.D0-E2*SP**2) X=XO+(N+H)*CP*DCOS(RLAM) Y=YO+(N+H)*CP*DSIN(RLAM) Z=ZO+(N*(1.D0-E2)+H)*SP RETURN END SUBROUTINE GPGELA(XI,YI,ZI,XJ,YJ,ZJ,PHII,LAMI,DAZ,XXII,ETA,DIJ,AIJ 1,AJI,ZIJ,ZJI,DXLG,DYLG,DZLG,DXLA,DYLA,DZLA) C C THIS SUBROUTINE WILL COMPUTE THE INVERSE POSITIONING PROBLEM IN C 3-DIMENSIONS. GIVEN THE X,Y,Z COORDINATES OF 2 POINTS, IT WILL C COMPUTE THE SPATIAL DISTANCE AND THE LOCAL ASTRONOMIC AZIMUTH AND C ZENITH ANGLE BETWEEN THE TWO POINTS. C C INPUT: C XI,YI,ZI,XJ,YJ,ZJ - CARTESIAN COORDINATES OF POINTS I AND J. C PHII,LAMI - GEODETIC LATITUDE AND LONGITUDE. C DAZ - ASTRO-AZIMUTH MINUS GEODETIC AZIMUTH. C XXII,ETA - DEFLECTION OF THE VERTICAL COMPONENTS OF POINT I. C OUTPUT: C DIJ - SPATAL DISTANCE. C AIJ,AJI - LOCAL ASTRONOMIC AZIMUTHS FROM I TO J AND FROM J TO I. C ZIJ,ZJI - ZENITH ANGLE FROM I TO J AND FROM J TO I. C DXLA,DYLA,DZLA - X,Y,Z LOCAL ASTRONOMIC POINT J. C DXLG,DYLG,DZLG - X,Y,Z LOCAL GEODETIC POINT J. C C WRITTEN BY G.BOWIE,NOV.1977. C IMPLICIT REAL*8(A-Z) PI=DARCOS(-1.D0) SP=DSIN(PHII) SL=DSIN(LAMI) CP=DCOS(PHII) CL=DCOS(LAMI) DXLG=-(XJ-XI)*SP*CL-(YJ-YI)*SP*SL+(ZJ-ZI)*CP DYLG=-(XJ-XI)*SL+(YJ-YI)*CL DZLG=(XJ-XI)*CP*CL+(YJ-YI)*CP*SL+(ZJ-ZI)*SP DXLA=DXLG-DAZ*DYLG-XXII*DZLG DYLA=DAZ*DXLG+DYLG-ETA*DZLG DZLA=XXII*DXLG+ETA*DYLG+DZLG DIJ=DSQRT((XJ-XI)**2+(YJ-YI)**2+(ZJ-ZI)**2) AIJ=DATAN2(DYLA,DXLA) AJI=AIJ+PI ZIJ=DARCOS(DZLA/DIJ) ZJI=PI-ZIJ DD=DIJ AA=AIJ ZZ=ZIJ RETURN END SUBROUTINE ERRI3D(XI,YI,ZI,XJ,YJ,ZJ,VC1,DXLG,DYLG,DZLG,DXLA,DYLA, 1 DZLA,DIJ,AIJ,ZIJ,DAZ,XXII,ETA,P,L,VC2,J) C C THIS SUBROUTINE WILL COMPUTE THE VARIANCE-COVARIANCE MATRIX OF THE C SPATIAL DISTANCE AND ASTRONOMIC AZIMUTH AND ZENITH ANGLE AS COMPUTED C IN THE INVERSE POSITIONING PROBLEM. C C INPUT: C XI,YI,ZI,XJ,YI,ZJ - CARTESIAN COORDINATES OF POINTS I AND J. C VC1 - VARIANCE COVARIANCE MATRIX OF CARTESIAN COORDINATES C POINTS I AND J. C DIJ,AIJ,ZIJ - COMPUTED INFORMATION - DIJ - SPATIAL DISTANCE. C - AIJ - ASTRONOMIC AZIMUTH. C - ZIJ - ZENITH ANGLE. C DAZ - ASTRONOMIC AZIMUTH MINUS THE GEODETIC AZIMUTH. C XXII,ETA - DEFLECTION OF THE VERTICAL COMPONENTS AT I. C P,L - GEODETIC LATITUDE AND LONGITUDE OF POINT I. C OUTPUT: C VC2 - VARIANCE COVARIANCE MATRIX OF THE SPATIAL DISTANCE, C ASTRONOMIC AZIMUTH AND ZENITH ANGLE. C C WRITTEN BY G.BOWIE,DEC.1977. C IMPLICIT REAL*8(A-Z) DIMENSION J(3,6),JT(6,3),VC1(6,6),VC2(3,3),JVC(3,6) SP=DSIN(P) CP=DCOS(P) SL=DSIN(L) CL=DCOS(L) COTANA=DCOTAN(AIJ) J(1,1)=-(XJ-XI)/DIJ J(1,2)=-(YJ-YI)/DIJ J(1,3)=-(ZJ-ZI)/DIJ J(1,4)=-J(1,1) J(1,5)=-J(1,2) J(1,6)=-J(1,3) J(2,1)=DXLA*((DAZ*SP*CL+SL+ETA*CP*CL)- DTAN(AIJ)*(SP*CL-DAZ*SL+XXI 1I*CP*CL))/(DXLA**2+DYLA**2) J(2,2)=DXLA*((DAZ*SP*SL-CL+ETA*CP*SL)-DTAN(AIJ)*(SP*SL+DAZ*CL+XXII 1 *CP*SL))/(DXLA**2+DYLA**2) J(2,3)=DXLA*((-DAZ*CP+ETA*SP)-DTAN(AIJ)*(-CP+XXII*SP))/(DXLA**2+DY 1LA**2) J(2,4)=-J(2,1) J(2,5)=-J(2,2) J(2,6)=-J(2,3) J(3,1)=-(XXII*SP*CL+ETA*SL-CP*CL+DCOS(ZIJ)*(XJ-XI)/DIJ)/DSQRT(DXLA 1 **2+DYLA**2) J(3,2)=-(XXII*SP*SL-ETA*CL-CP*SL+DCOS(ZIJ)*(YJ-YI)/DIJ)/DSQRT(DXLA 1 **2+DYLA**2) J(3,3)=-(-XXII*CP-SP+DCOS(ZIJ)*(ZJ-ZI)/DIJ)/DSQRT(DXLA**2+DYLA**2) J(3,4)=-J(3,1) J(3,5)=-J(3,2) J(3,6)=-J(3,3) CALL TRNSD(JT,6,J,3,3,6) CALL MMULD(JVC,3,J,3,VC1,6,3,6,6) CALL MMULD(VC2,3,JVC,3,JT,6,3,6,3) RETURN END SUBROUTINE ERRD3D(X,Y,Z,X1,Y1,Z1,P,L,R,DA,XI,ETA,AZ,ZE,VCDIR,VC,J) C C C THIS ROUTINE COMPUTES THE VARIANCE-COVARIANCE MATRIX INVOLVED IN THE C DIRECT PROBLEM OF GEODETIC POSITIONING IN A 3-DIMENSIONAL SYSTEM. C C INPUT: C C X,Y,Z - THE COORDINATES OF THE OBSERVATION STATION(IN METRES) C X1,Y1,Z1 - THE COORDINATES OF THE OBSERVED STATION.(IN METRES) C P,L - PHI,LAMDA OF THE OBSERVATION STATION IN THE GEODETIC C SYSTEM.(IN RADIANS) C R - THE OBSERVED SPATIAL DISTANCE.(IN METRES) C AZ - THE OBSERVED AZIMUTH IN THE LOCAL ASTONOMIC SYSTEM. C (IN RADIANS) C ZE - THE OBSERVED ZENITH ANGLE IN THE LOCAL ASTONOMIC SYSTEM. C (IN RADIANS) C XI,ETA - THE DEFLECTION OF THE VERTICAL COMPONENTS AT THE C OBSERVATION STATION.(IN RADIANS) C VCDIR - THE VARIANCE-COVARIANCE MATRIX OF THE COORDINATES OF THE C OBSERVATION STATION AND OF THE ACTUAL OBSERVATIONS. C C OUTPUT: C C VC - THE VARIANCE-COVARIANCE MATRIX OF THE COORDINATES OF THE C OBSERVATION STATION AND OF THE OBSERVED STATION. C C WRITTEN BY G.BOWIE,DEC.1977. C IMPLICIT REAL*8(A-Z) INTEGER I,K DIMENSION VC(6,6),VCDIR(6,6),JT(6,6),JVC(6,6),J(6,6) DO 1 I=1,6 DO 1 K=1,6 J(I,K)=0.D0 1 CONTINUE PI=3.141592653589793D0 CP=DCOS(P) CL=DCOS(L) CA=DCOS(AZ) CZ=DCOS(ZE) SP=DSIN(P) SL=DSIN(L) SA=DSIN(AZ) SZ=DSIN(ZE) T1=-SZ*SA+DA*SZ*CA T2=DA*SZ*SA+SZ*CA T3=XI*SZ*SA-ETA*SZ*CA T4=CZ*CA+DA*CZ*SA-XI*SZ T5=-DA*CZ*CA+CZ*SA-ETA*SZ T6=-XI*CZ*CA-ETA*CZ*SA-SZ J(1,1)=1.D0 J(2,2)=1.D0 J(3,3)=1.D0 J(4,1)=1.D0 J(5,2)=1.D0 J(6,3)=1.D0 J(4,4)=(X1-X)/R J(4,5)=-R*(SP*CL*T1+SL*T2-CP*CL*T3) J(4,6)=-R*(SP*CL*T4+SL*T5-CP*CL*T6) J(5,4)=(Y1-Y)/R J(5,5)=-R*(SP*SL*T1-CL*T2-CP*SL*T3) J(5,6)=-R*(SP*SL*T4-CL*T5-CP*SL*T6) J(6,4)=(Z1-Z)/R J(6,5)=R*(CP*T1+SP*T3) J(6,6)=R*(CP*T4+SP*T6) CALL TRNSD(JT,6,J,6,6,6) CALL MMULD(JVC,6,J,6,VCDIR,6,6,6,6) CALL MMULD(VC,6,JVC,6,JT,6,6,6,6) RETURN END SUBROUTINE ERR3D(PHI,RLAM,H,A,B,ICODE,CM,DM) C THIS ROUTINE COMPUTES THE COVARIANCE MATRIX CM (3X3) OF THE C CARTESIAN COORDINATES X,Y,Z GIVEN THE COVARIANCE MATRIX DM OF THE C ELLIPSOIDAL COORDINATES PHI,RLAM,H (OF THE SAME POINT) WHEN ICODE C IS 1. C IF ICODE IS -1 THE ROUTINE COMPUTES DM FROM THE GIVEN CM. C INPUT: C PHI-ELLIPSOIDAL LATITUDE OF THE POINT (IN RADIANS) C RLAM-ELLIPSOIDAL LONGITUDE OF THE POINT IN RADIANS C (POSITIVE EAST OF GREENWICH). C H-ELLIPSOIDAL HEIGHT OF THE POINT IN METRES. C A,B-SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE C ELLIPSOID IN METRES. C ICODE- 1 IF CM IS DESIRED FROM THE GIVEN DM C -1 IF DM IS DESIRED FROM THE GIVEN CM C C OUTPUT: C CM-3X3 COVARIANCE MATRIX OF THE CARTESIAN COORDINATES IN C METRES SQUARED. C DM-3X3 COVARIANCE MATRIX OF THE ELLIPSOIDAL COORDINATES; C IN RADIANS SQUARED FOR PHI,RLAM AND METRES SQARED FOR C H;IN RADIAN-METRES FOR THE COVARIANCES BETWEEN C PHI,RLAM AND H. C C SPACE MUST BE PROVIDED FOR CM AND DM IN THE CALLING ROUTINE. C C R.R.STEEVES, JUNE,1977. C IMPLICIT REAL*8(A-Z) INTEGER ICODE DIMENSION CM(3,3),DM(3,3) E2=(A*A-B*B)/(A*A) SP=DSIN(PHI) CP=DCOS(PHI) DEN=DSQRT(1.D0-E2*SP**2) SL=DSIN(RLAM) CL=DCOS(RLAM) N=A/DEN M=A*(1.D0-E2)/DEN**3 J11=-(M+H)*SP*CL J12=-(N+H)*CP*SL J13=CP*CL J21=-(M+H)*SP*SL J22=(N+H)*CP*CL J23=CP*SL J31=(M+H)*CP J32=0.D0 J33=SP IF(ICODE.LT.0)GO TO 1 CM(1,1)=J11**2*DM(1,1)+2.D0*J11*J12*DM(1,2)+2.D0*J11*J13*DM(1,3)+2 1 .D0*J12*J13*DM(2,3)+J13**2*DM(3,3)+J12**2*DM(2,2) CM(1,2)=J11*J21*DM(1,1)+(J11*J22+J12*J21)*DM(1,2)+(J11*J23+J13*J21 1 )*DM(1,3)+J12*J22*DM(2,2)+(J12*J23+J13*J22)*DM(2,3)+J13*J2 2 3*DM(3,3) CM(2,1)=CM(1,2) CM(1,3)=J11*J31*DM(1,1)+(J11*J32+J12*J31)*DM(1,2)+(J11*J33+J13*J31 1 )*DM(1,3)+J12*J32*DM(2,2)+(J12*J33+J13*J32)*DM(2,3)+J13* 2 J33*DM(3,3) CM(3,1)=CM(1,3) CM(2,2)=J21**2*DM(1,1)+2.D0*J21*J22*DM(1,2)+2.D0*J21*J23*DM(1,3)+ 1 J22**2*DM(2,2)+2.D0*J22*J23*DM(2,3)+J23**2*DM(3,3) CM(2,3)=J21*J31*DM(1,1)+(J21*J32+J22*J31)*DM(1,2)+(J21*J33+J23*J31 1 )*DM(1,3)+J22*J32*DM(2,2)+(J22*J33+J23*J32)*DM(2,3)+J23* 2 J33*DM(3,3) CM(3,2)=CM(2,3) CM(3,3)=J31**2*DM(1,1)+2.D0*J31*J32*DM(1,2)+2.D0*J31*J33*DM(1,3)+ 1 J32**2*DM(2,2)+2.D0*J32*J33*DM(2,3)+J33**2*DM(3,3) GO TO 2 1 CONTINUE DET=J11*J22*J33+J12*J23*J31+J13*J21*J32-J13*J22*J31-J12*J21*J33-J1 1 1*J23*J32 J1=J11 J2=J12 J3=J13 J4=J21 J5=J22 J6=J23 J7=J31 J11=(J22*J33-J23*J32)/DET J21=-(J21*J33-J31*J23)/DET J31=(J4*J32-J22*J31)/DET J12=-(J2*J33-J3*J32)/DET J22=(J1*J33-J3*J7)/DET J32=-(J1*J32-J2*J7)/DET J13=(J2*J6-J3*J5)/DET J23=-(J1*J6-J3*J4)/DET J33=(J1*J5-J2*J4)/DET DM(1,1)=J11**2*CM(1,1)+2.D0*J11*J12*CM(1,2)+2.D0*J11*J13*CM(1,3)+2 1 .D0*J12*J13*CM(2,3)+J13**2*CM(3,3)+J12**2*CM(2,2) DM(1,2)=J11*J21*CM(1,1)+(J11*J22+J12*J21)*CM(1,2)+(J11*J23+J13*J21 1 )*CM(1,3)+J12*J22*CM(2,2)+(J12*J23+J13*J22)*CM(2,3)+J13*J2 2 3*CM(3,3) DM(2,1)=DM(1,2) DM(1,3)=J11*J31*CM(1,1)+(J11*J32+J12*J31)*CM(1,2)+(J11*J33+J13* 1 J31)*CM(1,3)+J12*J32*CM(2,2)+(J12*J33+J13*J32)*CM(2,3)+J13 1 *J33*CM(3,3) DM(3,1)=DM(1,3) DM(2,2)=J21**2*CM(1,1)+2.D0*J21*J22*CM(1,2)+2.D0*J21*J23*CM(1,3)+ 1 J22**2*CM(2,2)+2.D0*J22*J23*CM(2,3)+J23**2*CM(3,3) DM(2,3)=J21*J31*CM(1,1)+(J21*J32+J22*J31)*CM(1,2)+(J21*J33+J23*J31 1 )*CM(1,3)+J22*J32*CM(2,2)+(J22*J33+J23*J32)*CM(2,3)+J23* 2 J33*CM(3,3) DM(3,2)=DM(2,3) DM(3,3)=J31**2*CM(1,1)+2.D0*J31*J32*CM(1,2)+2.D0*J31*J33*CM(1,3)+ 1 J32**2*CM(2,2)+2.D0*J32*J33*CM(2,3)+J33**2*CM(3,3) 2 RETURN END