C 00000010 C ADERIN2 00000020 C 00000030 C COMBINATION OF DOPPLER AND TERRESTRIAL NETWORKS FOR PARAMETERISING C SYSTEMATIC ERRORS IN A TERRESTRIAL GEODETIC NETWORK USING THE C METHOD OF LEAST SQUARES APPROXIMATION. C 00000070 C DATA INPUT 00000080 C 00000090 C 1. JOB HEADING (20A4) 00000100 C 00000110 C 2. PROGRAMME VARIABLES 00000120 C MDI (I5) =0 ALL NETWORK DATA FROM CARDS 00000130 C =1 SYSTEM 1 DATA FROM DISC(10) 00000140 C =2 SYSTEM 2 DATA FROM DISC(12) 00000150 C =3 ALL NETWORK DATA FROM DISC (10&12) 00000160 C 00000170 C 00000200 C =1 SATELLITE NET HAS DIAGONAL VAR-COVAR MATRIX 00000190 C ISIG (I5) =0 SATELLITE NETWORK HAS FULL VAR-COVAR MATRIX 00000180 C N1 (I5) INNER ZONE COMMON POINTS (EXCL. INITIAL POINT WHEN INPS=0 00000210 C 00000220 C N2 (I5) OUTER ZONE COMMON POINTS 00000230 C 00000240 C NUS (I5) =6 SOLVE FOR 3 TRANS,3 ROTATIONS 00000250 C =3 SOLVE FOR 3 TRANS 00000270 C 00000280 C C NUN (I5) =4 SOLVE FOR ALL CO-EFFICIENTS C INPS (I5) =0 INIT. PT. FXD.,NO SAT COORDS 00000350 C =1 INIT. PT. FXD.,SAT. COORDS. 00000360 C =2 INIT. PT. WTD.,NO SAT COORDS 00000370 C =3 INIT. PT. WTD.,SAT. COORDS 00000380 C ** NOTE : IF INPS=1 OR 3 , TERR. INIT. PT. IS INCLUDED IN SYSTEMS 00000390 C 1 AND 2 DATA SETS. 00000400 C IPHAS (I5) =1 ONLY ONE PHASE SPECIFIED 00000410 C =2 TWO PHASES REQUIRED 00000420 IMPLICIT REAL*8 (A-H,O-Z) 00000840 C APVF (F10.3) APRIORI VARIANCE FACTOR (DEFAULT IS 1.00) 00000430 C RA,RB (2F15.5) REFERENCE ELLIPSOID PARAMETERS ( A&B OR A& 1/F) 00000440 C 00000450 C 3. INITIAL POINT DATA (TERRESTRIAL) 00000460 C STA. NO. (I9) COLS 2-10) 00000470 C LAT. (2I4,F7.3) 00000480 C LONG. (2I4,F7.3) (+VE E.) 00000490 C ORTHO. HT. (F10.3) 00000500 C GEOID HT. (F10.3) 00000510 C STA. NAME (A8) 00000520 C 00000530 C 4. SYSTEM 1 (SATELLITE) COORDINATE DATA - ENTER ALL COMMON COORDS 00000540 C (1) HEADING (20A4) 00000550 C (2) STA.NO.,CARTESIAN COORDS,STA. NAME (1X,I9,3F15.3,2X,A8) 00000560 C IF MDI= 1 OR 3, READ FROM DISC DATA SET (10) 00000570 C 00000580 C 5. SYSTEM 1 VAR-COVAR MATRIX 00000590 C IF ISIG=0 ALL ELEMENTS (8F10.6 PER CARD) 00000600 C =1 DIAGONAL ELLEMENTS (8F10.6 PER CARD) 00000610 C IF MDI=1 OR 3 READ FROM DISC DATA SET (10) 00000620 C 00000630 REAL*8 NORMM,NORMN C 6. SYSTEM 2 (TERRESTRIAL) COORDINATE DATA 00000640 C SAME FORMAT AS PER SYSTEM COORDINATE DATA, OR DISC DATA SET (12) 00000650 C 00000660 C 7. SYSTEM 2 VAR-COVAR MATRIX 00000670 C ALL ELLEMENTS, 8F10.6 PER CARD, OR DISC DATA SET (12) 00000680 C 00000690 C **** NOTE : 1. IF INPS EQ. 1 OR 3, THE TERR. INIT. PT. DATA MUST BE 00000700 C INCLUDED IN THE SYSTEM 2 DATA INPUT(COORDS AND VAR- 00000710 C COVAR MATRIX). IF INPS EQ. 0 OR 2, TERR. INIT. PT. 00000720 C NOT TO BE INCLUDED IN SYSTEM 2 DATA INPUT. 00000730 C 2. INITIAL POINT COORDINATES MUST BE INPUT FIRST IN BOTH 00000740 C SYSTEMS WHEN THEY ARE REQUIRED. 00000750 C 3. ALWAYS ENTER FULL TERR. VAR-COVAR MATRIX(INCLUDING 00000760 C INITIAL POINT VALUES). 00000770 C 00000780 C 00000790 C DATA OUTPUT : DATA TO BE USED IN A SECOND PHASE IS PLACED ON 00000800 C 1. DATA SET (11) JOB DESC.,JOB PARAMETERS, INIT.PT.DATA, SYSTEM 00000810 C 1 AND 2 COORDINATE DATA. 00000820 C 2. DATA SET (14) X VECTOR, QX, QL, AND COV(X,L) 00000830 DIMENSION DJOB(20),DJOB1(20),IST1(21),X1(21),Y1(21),Z1(21), 00000860 1 STA1(21),VC1(63,63),DJOB2(20),X2(21),Y2(21),Z2(21),STA2(21), 00000870 2 VC2(63,63),A(12,6 ),B(12,24),W(12),VC(126,126),X(6),V(24), 3 SIGL(126,126) ,IST2(21),XKI(21),YKI(21),ZKI(21),SIGX(6,6) DIMENSION PNX(21),PNY(21),PNZ(21),FX(21),FY(21),FZ(21) DIMENSION B2(51,108),P2(108,108),ABP(51,108),ABT(51,51),PW(51) DIMENSION A2(51,24),VZ(63),PX(51),PPHI(21),PILAM(21),PH(21) DIMENSION SIGG(51), G(51,24),SIGJ(63,63),GST(51) DIMENSION SUMC(24),C(24),STDC(24),SIGMAC(24,24),SC2(24) DIMENSION GSA(51,24),SGPB(51,51),SGPT(51,51),SGPF(51),GSAT(51,51) DIMENSION STDP(24),SX(6), D(24),ALPHA(24,24),SL(24),VCK(24,24) MAX1=63 00000900 MAX2=MAX1*2 00000910 RHO=206264.8062471D0 00000920 C READ JOB HEADING, PROGRAMME VARIABLES, AND INITIAL POINT DATA 00000930 READ(5,10) DJOB 00000940 10 FORMAT(20A4) 00000950 READ(5,20) MDI,ISIG,N1,N2,NUS,NUN,INPS,IPHAS,APVF,RA,RB 00000960 20 FORMAT(8I5,F10.3,2F15.5) 00000970 IF(APVF.EQ.0) APVF=1.D0 00000980 READ(5,30) ISTK,IDK,IMK,SK,IIDK,IIMK,SSK,HK,GNK,STAK 00000990 30 FORMAT(I10,2(2I4,F7.3),2F10.3,A8) 00001000 C PRINT PAGE 1 00001010 WRITE(6,40) DJOB,N1,N2,APVF 00001020 40 FORMAT(' ',/////,6X,'COMBINATION OF 2 GEODETIC NETWORKS USING THE 00001030 1MODIFIED KRAKIWSKY-THOMSON MODEL',///,6X,20A4,////,6X,'NO. COMMON 00001040 2INNER ZONE POINTS : ',I5,//,6X,'NO. COMMON OUTER ZONE POINTS : ',I00001050 35,//,6X,'APRIORI VARIANCE FACTOR : ',F10.3,//) 00001060 WRITE(6,50) 00001070 50 FORMAT(' ',5X,'NO. SYSTEM UNKNOWNS : ') 00001080 IF(NUS.EQ.6) WRITE(6,70) NUS 00001090 IF(NUS.EQ.3) WRITE(6,90) NUS 00001110 70 FORMAT('+',30X,I5,' (3TRANSLATIONS,3ROTATIONS)',//) 00001120 90 FORMAT('+',30X,I5,' (3TRANSLATIONS)',//) 00001140 WRITE(6,110) 00001150 110 FORMAT(' ',5X,'NO. NETWORK UNKNOWNS : ') 00001160 IF(NUN.EQ.4) WRITE(6,133) 00001210 133 FORMAT('+',33X,' 24 CO-EFFICIENTS',//) C COMPUTE CARTESIAN COORDINATES OF THE TERR. INIT. PT. 00001270 IF(RB.LT.5.D6) ESQ=2.D0*(1.D0/RB)-(1.D0/RB)**2 00001280 IF(RB.GT.6.D6) ESQ=(RA**2-RB**2)/RA**2 00001290 PHI=(IABS(IDK)*3600.D0+IABS(IMK)*60.D0+DABS(SK))/RHO IF(IDK.LT.0.OR.IMK.LT.0.OR.SK.LT.0.E0) PHI=-PHI 00001310 DLAM=(IABS(IIDK)*3600.D0+IABS(IIMK)*60.D0+DABS(SSK))/RHO 0000132 IF(IIDK.LT.0.OR.IIMK.LT.0.OR.SSK.LT.0.E0) DLAM=-DLAM 00001330 RN=RA/DSQRT(1.D0-ESQ*(DSIN(PHI)**2)) 00001340 XK=(RN+HK+GNK)*DCOS(PHI)*DCOS(DLAM) 00001350 YK=(RN+HK+GNK)*DCOS(PHI)*DSIN(DLAM) 00001360 ZK=(RN*(1.D0-ESQ)+HK+GNK)*DSIN(PHI) 00001370 C PRINT INITIAL POINT DATA 00001380 WRITE(6,140) STAK,ISTK,IDK,IMK,SK,IIDK,IIMK,SSK,HK,GNK,XK,YK,ZK 00001390 140 FORMAT(' ',////,6X,'TERRESTRIAL INITIAL POINT',5X,A8,3X,I9,//,6X,'00001400 1LAT',2I4,F7.3,4X,'LONG(+VE E)',2I4,F7.3,4X,'ORTHO. HT.(M)',F10.3,400001410 2X,'GEOID HT.(M)',F10.3,///,6X,'CARTESIAN COORDINATES',//,6X,'X(M)'00001420 3,F15.3,4X,'Y(M)',F15.3,4X,'Z(M)',F15.3,///) 00001430 IF(INPS.EQ.0) WRITE(6,150) 00001440 IF(INPS.EQ.1) WRITE(6,160) 00001450 IF(INPS.EQ.2) WRITE(6,170) 00001460 150 FORMAT(' ',5X,'*** TERRESTRIAL INITIAL POINT IS FIXED, NO SATELLIT00001480 1E COORDINATES AVAILABLE ***') 00001490 IF(INPS.EQ.3) WRITE(6,180) 00001470 160 FORMAT(' ',5X,'*** TERRESTRIAL INITIAL POINT IS FIXED, SATELLITE C00001500 1OORDINATES ARE USED ***') 00001510 170 FORMAT(' ',5X,'*** TERRESTRIAL INITIAL POINT IS REDIFINED, NO SATE00001520 1LLITE COORDINATES AVAILABLE ***') 00001530 180 FORMAT(' ',5X,'*** TERRESTRIAL INITIAL POINT IS REDEFINED, SATELLI00001540 1TE COORDINATES ARE USED ***') 00001550 C COMPUTE VALUES OF PARAMETERS FOR THIS JOB 00001560 NCPT=N1+N2 00001570 NT=NCPT*3 00001580 NT6=NT*2 00001590 NVC=NT*2 00001600 IF(INPS.EQ.0.OR.INPS.EQ.2) NVC=NVC+3 00001610 NIP=N1 00001620 NIP3=NIP*3 00001630 NIP6=NIP*6 00001640 IF(INPS.EQ.0.OR.INPS.EQ.2) NIP6=NIP6+3 00001650 NTEMP=NIP+1 00001660 NVC2=NT 00001670 IF(INPS.EQ.0.OR.INPS.EQ.2) NVC2=NVC2+3 00001680 NU=30 NU1=6 NU2=24 NT1=N1*3 NT2=NT1*2 NT3=N2*3 NT4=NT3*2 NT5=NT4+NU1 NU3=NU1+1 C READ SYSTEM 1 DATA 00001700 DO 190 I=1,NT 00001710 DO 190 J=1,NT 00001720 190 VC1(I,J)=0.D0 00001730 IF(MDI.EQ.1.OR.MDI.EQ.3) GO TO 230 00001740 READ(5,10) DJOB1 00001750 DO 200 I=1,NCPT 00001760 200 READ(5,210) IST1(I),X1(I),Y1(I),Z1(I),STA1(I) 00001770 210 FORMAT(I10,3F15.3,2X,A8) 00001780 C READ SYSTEM 1 VAR-COVAR MATRIX 00001790 IF(ISIG.EQ.1) GO TO 225 00001800 READ(5,220) ((VC1(I,J),J=1,NT),I=1,NT) 00001810 220 FORMAT(8F10.6) 00001820 GO TO 250 00001830 225 READ(5,220) (VC1(I,I),I=1,NT) 00001840 GO TO 250 00001850 230 READ(10) DJOB1 00001860 DO 240 I=1,NCPT 00001870 240 READ(10) IST1(I),X1(I),Y1(I),Z1(I),STA1(I) 00001880 READ(10) ((VC1(I,J),J=1,NT),I=1,NT) 00001890 250 CONTINUE 00001900 C READ SYSTEM 2 DATA 00001910 IF(MDI.EQ.2.OR.MDI.EQ.3) GO TO 270 00001920 READ(5,10) DJOB2 00001930 DO 260 I=1,NCPT 00001940 260 READ(5,210) IST2(I),X2(I),Y2(I),Z2(I),STA2(I) 00001950 READ(5,220) ((VC2(I,J),J=1,NVC2),I=1,NVC2) 00001960 GO TO 290 00001970 270 READ(12) DJOB2 00001980 DO 280 I=1,NCPT 00001990 280 READ(12) IST2(I),X2(I),Y2(I),Z2(I),STA2(I) 00002000 READ(12) ((VC2(I,J),J=1,NVC2),I=1,NVC2) 00002010 290 CONTINUE 00002020 C PRINT SYSTEM 1 DATA 00002030 WRITE(6,300) DJOB,DJOB1 00002040 300 FORMAT('1',////,6X,20A4,//,6X,20A4,//,6X,'INNER ZONE',///) 00002050 WRITE(6,310) 00002060 310 FORMAT(' ',3X,'STA. NAME',6X,'STA. NO.',9X,'X(M)',7X,'STD.DEV.',8X00002070 1,'Y(M)',7X,'STD.DEV.',8X,'Z(M)',7X,'STD.DEV.',//) 00002080 DO 320 I=1,NIP 00002090 S1=DSQRT(VC1(3*I-2,3*I-2)) 00002100 S2=DSQRT(VC1(3*I-1,3*I-1)) 00002110 S3=DSQRT(VC1(3*I,3*I)) 00002120 320 WRITE(6,330) STA1(I),IST1(I),X1(I),S1,Y1(I),S2,Z1(I),S3 00002130 330 FORMAT(' ',4X,A8,3X,I10,1X,3(F17.3,F10.3),/) 00002140 IF(IPHAS.EQ.1) GO TO 355 00002150 WRITE(6,340) DJOB,DJOB1 00002160 340 FORMAT('1',////,6X,20A4,//,6X,20A4,//,6X,'OUTER ZONE',///) 00002170 WRITE(6,310) 00002180 DO 350 I=NTEMP,NCPT 00002190 S1=DSQRT(VC1(3*I-2,3*I-2)) 00002200 S2=DSQRT(VC1(3*I-1,3*I-1)) 00002210 S3=DSQRT(VC1(3*I,3*I)) 00002220 350 WRITE(6,330) STA1(I),IST1(I),X1(I),S1,Y1(I),S2,Z1(I),S3 00002230 C PRINT SYSTEM 2 DATA 00002240 355 WRITE(6,300) DJOB,DJOB2 00002250 WRITE(6,310) 00002260 DO 370 I=1,NIP 00002270 IF(INPS.EQ.0.OR.INPS.EQ.2) GO TO 360 00002280 S1=DSQRT(VC2(3*I-2,3*I-2)) 00002290 SI=DSQRT(VC2(3*I-1,3*I-1)) 00002300 S3=DSQRT(VC2(3*I,3*I)) 00002310 GO TO 370 00002320 360 S1=DSQRT(VC2(3*I+1,3*I+1)) 00002330 S2=DSQRT(VC2(3*I+2,3*I+2)) 00002340 S3=DSQRT(VC2(3*I+3,3*I+3)) 00002350 370 WRITE(6,330) STA2(I),IST2(I),X2(I),S1,Y2(I),S2,Z2(I),S3 00002360 WRITE(6,340) DJOB,DJOB2 00002370 WRITE(6,310) 00002380 DO 380 I=NTEMP,NCPT 00002390 IF(INPS.EQ.0.OR.INPS.EQ.2) GO TO 375 00002400 S1=DSQRT(VC2(3*I-2,3*I-2)) 00002410 SI=DSQRT(VC2(3*I-1,3*I-1)) 00002420 S3=DSQRT(VC2(3*I,3*I)) 00002430 GO TO 380 00002440 375 S1=DSQRT(VC2(3*I+1,3*I+1)) 00002450 S2=DSQRT(VC2(3*I+2,3*I+2)) 00002460 S3=DSQRT(VC2(3*I+3,3*I+3)) 00002470 380 WRITE(6,330) STA2(I),IST2(I),X2(I),S1,Y2(I),S2,Z2(I),S3 00002480 C DESIGN MATRIX A 00002490 DO 390 I=1,NT1 DO 390 J=1,NU1 390 A(I,J)=0.D0 00002520 C FORM A1 MATRIX DO 400 I=1,N1 A(3*I-2,1)=1.D0 00002630 A(3*I-1,2)=1.D0 00002640 A(3*I,3)=1.D0 00002650 A(3*I-2,5)=-Z2(I)*1.D-6 00002660 A(3*I-2,6)=Y2(I)*1.D-6 00002670 A(3*I-1,4)=Z2(I)*1.D-6 00002680 A(3*I-1,6)=-X2(I)*1.D-6 00002690 A(3*I,4)=-Y2(I)*1.D-6 00002700 400 A(3*I,5)=X2(I)*1.D-6 00002710 C FORM B MATRIX 00002860 DO 410 I=1,NT1 DO 410 J=1,NT2 410 B(I,J)=0.D0 00002890 NSAT=3 00002900 IF(INPS.EQ.1.OR.INPS.EQ.3) NSAT=NSAT-3 00002910 C INSERT B11 ELEMENTS 00002920 DO 430 I=1,NIP 00002930 B(3*I-2,6*I-5)=1.D0 B(3*I-1,6*I-4)=1.D0 B(3*I,6*I-3)=1.D0 B(3*I-2,6*I-2)=-1.D0 B(3*I-1,6*I-1)=-1.D0 B(3*I,6*I)=-1.D0 430 CONTINUE 00003090 C APPLY SCALE CORRECTION OF MINUS 1:1,000,000 TO ALL DOPPLER COORDINATES DO 11 I=1,NCPT X1(I)=X1(I)*.999999D0 Y1(I)=Y1(I)*.999999D0 11 Z1(I)=Z1(I)*.999999D0 C FORM MISCLOSURE VECTOR 00003180 DO 440 I=1,N1 W(3*I-2)=X2(I)-X1(I) 00003200 W(3*I-1)=Y2(I)-Y1(I) 00003210 440 W(3*I)=Z2(I)-Z1(I) 00003220 C FORM VAR-COVAR MATRIX VC 00003230 DO 450 I=1,NVC 00003240 DO 450 J=1,NVC 00003250 450 VC(I,J)=0.D0 00003260 NSKIP=0 00003270 IF(INPS.EQ.0.OR.INPS.EQ.2) NSKIP=NSKIP+3 00003280 C ENETER VC2 INTO VC 00003290 KK=1 00003300 II=1 00003310 K=1 00003320 L=1 00003330 DO 510 I=1,NVC2 00003340 DO 510 J=1,NVC2 00003350 JJ=J 00003360 IF(I.GT.3.AND.J.GT.3) GO TO 460 00003370 VC(K,L)=VC2(I,J) 00003380 GO TO 470 00003390 460 VC(K,L)=VC2(KK,II)+VC2(I,J)-VC2(KK,J)-VC2(I,II) 00003400 470 CONTINUE 00003410 IF(JJ.EQ.NVC2) GO TO 490 00003420 L=L+1 00003430 II=II+1 00003440 IF(II.GT.3) GO TO 480 00003450 GO TO 510 00003460 480 L=L+3-NSKIP 00003470 IF(NSKIP.EQ.3.AND.J.GE.6) L=L+3 00003480 II=1 00003490 GO TO 510 00003500 490 K=K+1 00003510 KK=KK+1 00003520 L=1 00003530 II=1 00003540 IF(KK.GT.3) GO TO 500 00003550 GO TO 510 00003560 500 K=K+3-NSKIP 00003570 IF(NSKIP.EQ.3.AND.I.GE.6) K=K+3 00003580 KK=1 00003590 510 CONTINUE 00003600 C ENTER VC1 IN VC 00003610 KK=1 00003620 II=1 00003630 K=4+NSKIP 00003640 L=4+NSKIP 00003650 DO 550 I=1,NT 00003660 DO 550 J=1,NT 00003670 JJ=J 00003680 VC(K,L)=VC1(I,J) 00003690 IF(JJ.EQ.NT) GO TO 530 00003700 L=L+1 00003710 II=II+1 00003720 IF(II.GT.3) GO TO 520 00003730 GO TO 550 00003740 520 L=L+3 00003750 II=1 00003760 GO TO 550 00003770 530 K=K+1 00003780 KK=KK+1 00003790 L=4+NSKIP 00003800 II=1 00003810 IF(KK.GT.3) GO TO 540 00003820 GO TO 550 00003830 540 K=K+3 00003840 KK=1 00003850 550 CONTINUE 00003860 C FORM VCK MATRIX DO 551 I=1,NT2 DO 551 J=1,NT2 551 VCK(I,J)=VC(I,J) C LEAST SQUARES ESTIMATION SOLUTION 00003870 WRITE(6,560) 00003880 560 FORMAT('1',////,6X,'SUMMARY ESTIMATION PROCEDURE',///) 00003890 CALL CELS1(A,B,W,VCK,APVF,NT1,NU1,NT2,NT1,NT2,NUS,NUN,X,V,SIGX, 1SIGL,SX,SL,&9999) C PRINT RESULTS 00003920 X(4)=X(4)*RHO*1.D-6 00003930 X(5)=X(5)*RHO*1.D-6 00003940 X(6)=X(6)*RHO*1.D-6 00003950 SX(4)=SX(4)*RHO*1.D-6 SX(5)=SX(5)*RHO*1.D-6 SX(6)=SX(6)*RHO*1.D-6 WRITE(6,570) DJOB 570 FORMAT('-',6X,'RESULTS OF THE GEODETIC NETWORK COMBINATION USING 1LEAST SQUARES ESTIMATION AND LEAST SQUARES ',//,32X,'APPROXIMATION 2PROCEDURES',//,6X,20A4,//) WRITE(6,590) X(1), SX(1),X(2), SX(2),X(3), SX(3) 590 FORMAT(' ',5X,'TRANSLATION COMPONENTS (METRES)',///,6X,'DX : ',F8.00004090 13,3X,'STD. DEV. : ',F8.3,//,6X,'DY : ',F8.3,3X,'STD. DEV. : ',F8.300004100 2,//,6X,'DZ : ',F8.3,3X,'STD. DEV. : ',F8.3,////) 00004110 IF(NUS.GT.3) WRITE(6,600) X(4), SX(4),X(5), SX(5),X(6), SX(6) 00004120 600 FORMAT(' ',5X,'ROTATION ELEMENTS(ARC-SECS)',//,6X,'X-AXIS : ',F8.400004130 1,3X,'STD. DEV. : ',F8.3,//,6X,'Y-AXIS : ',F8.4,3X,'STD. DEV. : ',F00004140 28.3,//,6X,'Z-AXIS : ',F8.4,3X,'STD. DEV. : ',F8.3,////) 00004150 C SOLUTION OF F2 USING LEAST SQUARES APPROXIMATION C FORM A2 ELEMENTS DO 441 I=1,NT3 DO 441 J=1,NU2 441 A2(I,J)=0.D0 DO 402 I=NTEMP,NCPT 00002730 XKI(I)=X2(I)-XK YKI(I)=Y2(I)-YK ZKI(I)=Z2(I)-ZK C SCALE DIFFERENCES IN COORDINATES (TERRESTRIAL ONLY) BY 1:10008000. XKI(I)=XKI(I)*1.D-6 YKI(I)=YKI(I)*1.D-6 ZKI(I)=ZKI(I)*1.D-6 C DEFINE BASE FUNCTIONS FI1=1.0D0 FI2=XKI(I) FI3=YKI(I) FI4=ZKI(I) FI5=XKI(I)*YKI(I) FI6=XKI(I)*ZKI(I) FI7=YKI(I)*ZKI(I) FI8=XKI(I)*YKI(I)*ZKI(I) J=I-N1 A2(3*J-2,1)=FI1 A2(3*J-2,2)=FI2 A2(3*J-2,3)=FI3 A2(3*J-2,4)=FI4 A2(3*J-2,5)=FI5 A2(3*J-2,6)=FI6 A2(3*J-2,7)=FI7 A2(3*J-2,8)=FI8 A2(3*J-1,9)=FI1 A2(3*J-1,10)=FI2 A2(3*J-1,11)=FI3 A2(3*J-1,12)=FI4 A2(3*J-1,13)=FI5 A2(3*J-1,14)=FI6 A2(3*J-1,15)=FI7 A2(3*J-1,16)=FI8 A2(3*J,17)=FI1 A2(3*J,18)=FI2 A2(3*J,19)=FI3 A2(3*J,20)=FI4 A2(3*J,21)=FI5 A2(3*J,22)=FI6 A2(3*J,23)=FI7 A2(3*J,24)=FI8 402 CONTINUE C C FORM B2 MATRIX AND PW VECTOR C XO=X(1) YO=X(2) ZO=X(3) EX=X(4)/RHO EY=X(5)/RHO EZ=X(6)/RHO KO=0 800 CONTINUE DO 431 I=1,NT3 DO 431 J=1,NT5 431 B2(I,J)=0.D0 DO 432 I=1,N2 B2(3*I-2,1)=1.D0 B2(3*I-2,5)=-Z2(I+N1)*1.D-6 B2(3*I-2,6)=Y2(I+N1)*1.D-6 B2(3*I-1,2)=1.D0 B2(3*I-1,4)=Z2(I+N1)*1.D-6 B2(3*I-1,6)=-X2(I+N1)*1.D-6 B2(3*I,3)=1.D0 B2(3*I,4)=-Y2(I+N1)*1.D-6 B2(3*I,5)=X2(I+N1)*1.D-6 B2(3*I-2,6*I+1)=1.D0 B2(3*I-2,6*I+4)=-1.D0 B2(3*I-1,6*I+2)=1.D0 B2(3*I-1,6*I+5)=-1.D0 B2(3*I,6*I+3)=1.D0 B2(3*I,6*I+6)=-1.D0 PW(3*I-2)=X1(I+N1)-X2(I+N1)-EZ*Y2(I+N1)+EY*Z2(I+N1)-XO PW(3*I-1)=Y1(I+N1)-Y2(I+N1)+EZ*X2(I+N1)-EX*Z2(I+N1)-YO PW(3*I)=Z1(I+N1)-Z2(I+N1)-EY*X2(I+N1)+EX*Y2(I+N1)-ZO 432 CONTINUE C FORM P2 MATRIX DO 442 I=1,NT5 DO 442 J=1,NT5 442 P2(I,J)=0.D0 DO 443 I=1,NU1 DO 443 J=1,NU1 443 P2(I,J)=SIGX(I,J) DO 4422 I=NU3,NT5 DO 4422 J=NU3,NT5 4422 P2(I,J)=VC(I+18,J+18) C COMPUTE 'M' MATRIX DO 444 I=1,NT3 DO 444 J=1,NT5 444 ABP(I,J)=0.D0 DO 445 I=1,NT3 DO 445 J=1,NT5 DO 445 K=1,NT5 445 ABP(I,J)=ABP(I,J)+B2(I,K)*P2(K,J) DO 446 I=1,NT3 DO 446 J=1,NT3 446 ABT(I,J)=0.D0 DO 447 I=1,NT3 DO 447 J=1,NT3 DO 447 K=1,NT5 447 ABT(I,J)=ABT(I,J)+ABP(I,K)*B2(J,K) CALL CHOLD(ABT,NT3,NT3,DET,&448) GO TO 449 448 WRITE(6,451) 451 FORMAT(' ',//,10X,'DETERMINANT OF M IS LESS THAN 1.D-10',//) 449 CONTINUE DO 452 I=1,NT3 452 PX(I)=ABT(I,I) INDEX0 = 0 INDEX1 = 1 INDEX2 = 2 INDEX3 = 3 CALL ORTHO(NT3,NU2,1.D0,A2,NU2,SIGMAF,VFAC,NPC,INDEX2,VZ,SIGMAC, 1PW,PX,D,ALPHA,C,SUMC,SC2,STDP,6) C CHECK WEATHER ALL FOURIER COEFFICIENTS ARE DISCARDED IF(NPC.GT.0) GO TO 408 WRITE(6,409) 409 FORMAT(10X,'** WARNING ** ALL FOURIER COEFFICIENTS DISCARDED. TRY 1AGAIN WITH LOWER LEVEL OF PROBABILITY') GO TO 9999 408 CONTINUE C* PRINT RESULTS WRITE(6,418) NPC 418 FORMAT(10X,'NUMBER OF COEF. OF THE ORIGINAL POLYNOMIAL AFTER THE 1TEST IS PERFORMED : ',I4,//) WRITE(6,421) 421 FORMAT(10X,'PARAMETERS OF ORTHORGONALIZED POLYNOMIAL',///,16X, 1'FOURIER COEFFICIENTS',8X,'STANDARD DEVIATION',//) WRITE(6,411) (I,C(I),SUMC(I),I=1,NU2) 411 FORMAT(18X,I2,'. ',F11.6,6X,F11.6) WRITE(6,412) SIGMAF 412 FORMAT(//,10X,'THE FOURIER POLYNOMIAL A POSTEORI VARIANCE FACTOR : 1',F9.5,//) WRITE(6,415) VFAC 415 FORMAT(//,10X,'THE ORIGINAL POLYNOMIAL A POSTEORI VARIANCE FACTOR 1:',F9.5,//) WRITE(6,413) 413 FORMAT(///,40X,'COEFFICIENTS OF ORIGINAL POLYNOMIALS',//,13X, 1 'COEFF.NO.',4X,'CX',8X,'STD.DEV.',9X,'CY',8X,'STD.DEV.',9X, 2 'CZ',8X,'STD.DEV.') DO 417 I=1,NU2 417 STDC(I)=DSQRT(SIGMAC(I,I)) WRITE(6,414) (I,D(I),STDC(I),D(I+8),STDC(I+8),D(I+16),STDC(I+16), 1I=1,8) 414 FORMAT('0',15X,I2,'.',F12.6,F13.6,F14.6,F13.6,F14.6,F13.6) C CX1=D(1) CX2=D(2) CX3=D(3) CX4=D(4) CX5=D(5) CX6=D(6) CX7=D(7) CX8=D(8) CY1=D(9) CY2=D(10) CY3=D(11) CY4=D(12) CY5=D(13) CY6=D(14) CY7=D(15) CY8=D(16) CZ1=D(17) CZ2=D(18) CZ3=D(19) CZ4=D(20) CZ5=D(21) CZ6=D(22) CZ7=D(23) CZ8=D(24) C WRITE(6,1001) 1001 FORMAT('-',10X,'FUNCTIONAL VALUES OF POLYNOMIALS USING THE TERREST 1RIAL COORDINATES',//) 3114 FORMAT('0',8X,A7,3X,I10,3X,F12.4,3X,F10.4,3X,F12.4,3X,F10.4,3X 1,F10.4,F11.4) DO 1006 I=NTEMP,NCPT FI1=1.0D0 FI2=XKI(I) FI3=YKI(I) FI4=ZKI(I) FI5=XKI(I)*YKI(I) FI6=XKI(I)*ZKI(I) FI7=YKI(I)*ZKI(I) FI8=XKI(I)*YKI(I)*ZKI(I) PNX(I)=FI1*CX1+FI2*CX2+FI3*CX3+FI4*CX4+FI5*CX5+FI6*CX6+FI7*CX7 1 +FI8*CX8 PNY(I)=FI1*CY1+FI2*CY2+FI3*CY3+FI4*CY4+FI5*CY5+FI6*CY6+FI7*CY7 1 +FI8*CY8 PNZ(I)=FI1*CZ1+FI2*CZ2+FI3*CZ3+FI4*CZ4+FI5*CZ5+FI6*CZ6+FI7*CZ7 1 +FI8*CZ8 1006 CONTINUE C C FORM COVARIANCE MATRICES WRITE(6,3333) 3333 FORMAT(//,7X,'STA.NAME',9X,'STA.NO',8X,'X(M)',6X,'STD.DEV.',6X, 1 'Y(M)',6X,'STD.DEV.',6X,'Z(M)',6X,'STD.DEV.',//) DO 3101 I=1,NT3 DO 3101 J=1,NU2 3101 G(I,J)=0.0D0 DO 3103 I=1,NT3 DO 3103 J=1,NU2 G(I,J)=A2(I,J) 3103 CONTINUE DO 4111 I=1,NT DO 4111 J=1,NT 4111 SIGJ(I,J)=0.D0 AA=6378206.4D0 BB=6356583.8D0 DO 3001 I=1,NCPT CALL PLHXYZ(X2(I),Y2(I),Z2(I),RP,RL,HK) CALL NORMAL(AA,BB,RP,XM,XN) H1=XM+HK H2=(XN+HK)*DCOS(RP) SIGJ(3*I-2,3*I-2)=-DSIN(RP)*DCOS(RL)/H1 SIGJ(3*I-2,3*I-1)=-DSIN(RP)*DSIN(RL)/H1 SIGJ(3*I-2,3*I)=DCOS(RP)/H1 SIGJ(3*I-1,3*I-2)=-DSIN(RL)/H2 SIGJ(3*I-1,3*I-1)=DCOS(RL)/H2 SIGJ(3*I,3*I-2)=DCOS(RP)*DCOS(RL) SIGJ(3*I,3*I-1)=DCOS(RP)*DSIN(RL) SIGJ(3*I,3*I)=DSIN(RP) 3001 CONTINUE DO 4013 I=NTEMP,NCPT CALL PLHXYZ(X2(I),Y2(I),Z2(I),RP,RL,HK) CALL NORMAL(AA,BB,RP,XM,XN) H1=XM+HK H2=(XN+HK)*DCOS(RP) PPHI(I)=(RHO/H1)*(-DSIN(RP)*DCOS(RL)*PNX(I)-DSIN(RP)*DSIN(RL)*PNY( 1I)+DCOS(RP)*PNZ(I)) PILAM(I)=(RHO/H2)*(-DSIN(RL)*PNX(I)+DCOS(RL)*PNY(I)) PH(I)=DCOS(RP)*DCOS(RL)*PNX(I)+DCOS(RP)*DSIN(RL)*PNY(I)+DSIN(RP)* 1 PNZ(I) 4013 CONTINUE MN=N1*3 NN=N2*3 DO 4115 I=1,NT3 DO 4115 J=1,NU2 4115 GSA(I,J)=0.D0 DO 4112 I=1,NT3 DO 4112 J=1,NT3 4112 GSAT(I,J)=0.D0 DO 4005 I=1,NT3 DO 4005 J=1,NT3 4005 SGPB(I,J)=0.D0 DO 4006 I=1,NT3 DO 4006 J=1,NT3 4006 SGPT(I,J)=0.D0 DO 4113 I=1,NT3 DO 4113 J=1,NU2 DO 4113 K=1,NU2 4113 GSA(I,J)=GSA(I,J)+G(I,K)*SIGMAC(K,J) DO 4114 I=1,NT3 DO 4114 J=1,NT3 DO 4114 K=1,NU2 4114 GSAT(I,J)=GSAT(I,J)+GSA(I,K)*G(J,K) DO 3111 I=1,NT3 SIGG(I)=DSQRT(GSAT(I,I)) 3111 CONTINUE DO 3115 I=NTEMP,NCPT J=I-N1 3115 WRITE(6,3334) STA2(I),IST2(I),PNX(I),SIGG(3*J-2),PNY(I), 1SIGG(3*J-1),PNZ(I),SIGG(3*J) 3334 FORMAT('0',8X,A7,3X,I10,3X,6F12.4) DO 4007 I=1,NT3 DO 4007 J=1,NT3 DO 4007 K=1,NT3 4007 SGPB(I,J)=SGPB(I,J)+SIGJ(I+MN,K+MN)*GSAT(K,J) DO 4008 I=1,NT3 DO 4008 J=1,NT3 DO 4008 K=1,NT3 4008 SGPT(I,J)=SGPT(I,J)+SGPB(I,K)*SIGJ(J+MN,K+MN) DO 4009 I=1,N2 DO 4009 J=1,2 4009 SGPF(3*I-J)=DSQRT(SGPT(3*I-J,3*I-J))*RHO DO 4010 I=1,N2 4010 SGPF(3*I)=DSQRT(SGPT(3*I,3*I)) WRITE(6,4011) WRITE(6,3112) 4011 FORMAT('-',10X,'FUNCTIONAL VALUES OF POLYNOMIALS IN CURVILINEAR CO 1-ORDINATES',/) 3112 FORMAT('-',7X,'STA. NAME',6X,'STA. NO.',6X,'LAT(SEC)',5X,'STD.DEV. 1',6X,'LONG(SEC)',4X,'STD.DEV.',6X,'HT.(M)',4X,'STD.DEV.') DO 4012 I=NTEMP,NCPT WRITE(6,3114) STA2(I),IST2(I),PPHI(I),SGPF(3*(I-N1)-2),PILAM(I), 1SGPF(3*(I-N1)-1),PH(I),SGPF(3*(I-N1)) 4012 CONTINUE C C COMPUTE DISTORTION AT SELECTED POINTS CALL DSTORT(X2(1),Y2(1),Z2(1),NU2,D,SIGMAC,STDC,KO) C 9999 STOP 00005290 END 00005300 SUBROUTINE CELS1(A,B,W,P,APVF,NP,NU,NPP,MAX1,MAX2,NUS,NUN,X,V, 00005310 1 SIGX,SIGL, SX,SL,*) C* CELS PERFORMS A COMBINED LEAST SQUARES ADJUSTMENT , GIVEN THE DESIGN00005330 C* MATRICIES A AND B, THE MISCLOSURE VECTOR W, AND THE INVERSE OF THE *00005340 C* WEIGHT MATRIX OF OBSERVATIONS P. ALLOWANCE IS MADE FOR WEIGHTED *00005350 C* PARAMETERS( SEE NOTES ON THE PARAMETER 'IP') *00005360 IMPLICIT REAL*8 (A-H,O-Z) 00005370 REAL*8 A(MAX1,NU),B(MAX1,NPP),W(MAX1),P(MAX2,NPP),BP(12,24), 1 M(12,12),N( 6 ,6),ATM(6,12),AMW(6),X(6),AX(12),MK(12), 2 PBT(24,12),V(24),SIGX(6,6),SX(6),AQ(12,6),AQM(12,12), 3 QK(12,12),PBQ(24,12),SIGL(24,24),SL(24) RHO=206264.8062471D0 00005420 IPX=0 00005430 C* COMPUTE M AND INVERT 00005440 DO 1 I=1,NP 00005450 DO 1 J=1,NPP 00005460 1 BP(I,J)=0.D0 00005470 DO 2 I=1,NP 00005480 DO 2 J=1,NPP 00005490 DO 2 K=1,NPP 00005500 2 BP(I,J)=BP(I,J)+B(I,K)*P(K,J) 00005510 DO 3 I=1,NP 00005520 DO 3 J=1,NP 00005530 3 M(I,J)=0.D0 00005540 DO 4 I=1,NP 00005550 DO 4 J=1,NP 00005560 DO 4 K=1,NPP 00005570 4 M(I,J)=M(I,J)+BP(I,K)*B(J,K) 00005580 CALL CHOLD (M,MAX1,NP,DET,&5) 00005590 GO TO 20 00005600 5 WRITE(6,10) 00005610 10 FORMAT(' ',//,10X,'DETERMINANT OF M IS LESS THAN 1.D-10') 00005620 RETURN 1 00005630 20 WRITE(6,30) DET 00005640 30 FORMAT(' ',//,5X,'DET OF M IS : ',D16.8,//) 00005650 C* SOLVE FOR THE SOLUTION VECTOR X 00005660 DO 40 I=1,NU 00005670 DO 40 J=1,NP 00005680 40 ATM(I,J)=0.D0 00005690 DO 50 I=1,NU 00005700 DO 50 J=1,NP 00005710 DO 50 K=1,NP 00005720 50 ATM(I,J)=ATM(I,J)+A(K,I)*M(K,J) 00005730 DO 60 I=1,NU 00005740 DO 60 J=1,NU 00005750 60 N(I,J)=0.D0 00005760 DO 70 I=1,NU 00005770 DO 70 J=1,NU 00005780 DO 70 K=1,NP 00005790 70 N(I,J)=N(I,J)+ATM(I,K)*A(K,J) 00005800 CALL CHOLD(N,NU,NU,DET,&90) GO TO 100 00006100 90 WRITE(6,95) 00006110 95 FORMAT(' ',//,10X,'DETERMINANT OF N IS LESS THAN 1.D-10') 00006120 RETURN 1 00006130 100 WRITE(6,110) DET 00006140 110 FORMAT(' ',//,5X,'DET OF N IS : ',D16.8,//) 00006150 DO 120 I=1,NU 00006160 120 AMW(I)=0.D0 00006170 DO 130 I=1,NU 00006180 DO 130 J=1,NP 00006190 130 AMW(I)=AMW(I)+ATM(I,J)*W(J) 00006200 DO 140 I=1,NU 00006210 140 X(I)=0.D0 00006220 DO 150 I=1,NU 00006230 DO 150 J=1,NU 00006240 150 X(I)=X(I)+N(I,J)*AMW(J) 00006250 DO 160 I=1,NU 00006260 160 X(I)=-1.D0*X(I) 00006270 C* COMPUTE CORRELATE MATRIX K 00006280 DO 170 I=1,NP 00006290 170 AX(I)=0.D0 00006300 DO 180 I=1,NP 00006310 DO 180 J=1,NU 00006320 180 AX(I)=AX(I)+A(I,J)*X(J) 00006330 DO 190 I=1,NP 00006340 190 AX(I)=AX(I)+W(I) 00006350 DO 200 I=1,NP 00006360 200 MK(I)=0.D0 00006370 DO 210 I=1,NP 00006380 DO 210 J=1,NP 00006390 210 MK(I)=MK(I)+((-1.D0)*M(I,J)*AX(J)) 00006400 C* COMPUTE RESIDUAL VECTOR V 00006410 DO 220 I=1,NPP 00006420 DO 220 J=1,NP 00006430 220 PBT(I,J)=0.D0 00006440 DO 230 I=1,NPP 00006450 DO 230 J=1,NP 00006460 DO 230 K=1,NPP 00006470 230 PBT(I,J)=PBT(I,J)+P(I,K)*B(J,K) 00006480 DO 240 I=1,NPP 00006490 240 V(I)=0.D0 00006500 DO 250 I=1,NPP 00006510 DO 250 J=1,NP 00006520 250 V(I)=V(I)+PBT(I,J)*MK(J) 00006530 C* COMPUTE VTPV,DF,AND EVF 00006540 VPV=0.D0 00006550 DO 260 I=1,NP 00006560 260 VPV=VPV+((-1.D0)*MK(I)*W(I)) 00006570 IDF=NP-NU EVF=VPV/IDF 00006590 WRITE(6,270) IDF,VPV,EVF 00006600 270 FORMAT(' ',///,5X,'DEGREES OF FREEDOM :',I10,//,10X,'VPV :',F15.5,00006610 1//,5X,'ESTIMATED VARIANCE FACTOR :',F10.6,//) 00006620 DO 280 I=1,NU 00006630 DO 280 J=1,NU 00006640 280 SIGX(I,J)=APVF*N(I,J) 00006650 DO 290 I=1,NU 00006660 290 SX(I)=DSQRT(SIGX(I,I)) 00006670 C* COMPUTE SIGL 00006680 DO 300 I=1,NP 00006690 DO 300 J=1,NU 00006700 300 AQ(I,J)=0.D0 00006710 DO 310 I=1,NP 00006720 DO 310 J=1,NU 00006730 DO 310 K=1,NU 00006740 310 AQ(I,J)=AQ(I,J)+A(I,K)*N(K,J) 00006750 DO 320 I=1,NP 00006760 DO 320 J=1,NP 00006770 320 AQM(I,J)=0.D0 00006780 DO 330 I=1,NP 00006790 DO 330 J=1,NP 00006800 DO 330 K=1,NU 00006810 330 AQM(I,J)=AQM(I,J)+AQ(I,K)*ATM(K,J) 00006820 DO 340 I=1,NP 00006830 DO 340 J=1,NP 00006840 340 QK(I,J)=0.D0 00006850 DO 350 I=1,NP 00006860 DO 350 J=1,NP 00006870 DO 350 K=1,NP 00006880 350 QK(I,J)=QK(I,J)+M(I,K)*AQM(K,J) 00006890 DO 360 I=1,NP 00006900 DO 360 J=1,NP 00006910 360 QK(I,J)=M(I,J)-QK(I,J) 00006920 DO 370 I=1,NPP 00006930 DO 370 J=1,NP 00006940 370 PBQ(I,J)=0.D0 00006950 DO 380 I=1,NPP 00006960 DO 380 J=1,NP 00006970 DO 380 K=1,NP 00006980 380 PBQ(I,J)=PBQ(I,J)+PBT(I,K)*QK(K,J) 00006990 DO 390 I=1,NPP 00007000 DO 390 J=1,NPP 00007010 390 SIGL(I,J)=0.D0 00007020 DO 400 I=1,NPP 00007030 DO 400 J=1,NPP 00007040 DO 400 K=1,NP 00007050 400 SIGL(I,J)=SIGL(I,J)+PBQ(I,K)*BP(K,J) 00007060 DO 405 I=1,NPP 00007070 DO 405 J=1,NPP 00007080 405 SIGL(I,J)=P(I,J)-SIGL(I,J) 00007090 DO 410 I=1,NPP 00007100 DO 410 J=1,NPP 00007110 410 SIGL(I,J)=APVF*SIGL(I,J) 00007120 DO 420 I=1,NPP 00007130 420 SL(I)=DSQRT(SIGL(I,I)) 00007140 C* PRINT X AND V VECTORS 00007150 WRITE(6,430) 00007160 430 FORMAT('1',//,20X,'SOLUTION VECTOR',///) 00007170 WRITE(6,440) (X(I),I=1,NU) 00007180 440 FORMAT(20X,D19.12) 00007190 WRITE(6,450) 00007200 450 FORMAT('1',//,20X,'RESIDUAL VECTOR',///) 00007210 WRITE(6,460) (V(I),I=1,NPP) 00007220 460 FORMAT(20X,F15.5) 00007230 RETURN 00007240 END 00007250 SUBROUTINE CHOLD(A,IRDA,NA,DETA,*) 00007260 C 00007270 C MATRIX INVERSION USING CHOLESKI DECOMPOSITION 00007280 C 00007290 C INPUT ARGUMENTS 00007300 C A = ARRAY CONTAINING POSITIVE DEFINITE SYMMETRIC INPUT MATRIX 00007310 C IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX 00007320 C NA = SIZE OF INPUT MATRIX 00007330 C OUTPUT ARGUMENTS 00007340 C DETA = DETERMINANT OF INPUT MATRIX 00007360 C A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED) 00007350 C * = ERROR RETURN (TAKEN IF NA .LT. 1 OR IF DETA .LT. SING) 00007370 C 00007380 DOUBLE PRECISION A,DETA,SUM,SQRT,DSQRT,ABS,DABS,SING 00007390 DIMENSION A(IRDA,NA) 00007400 SQRT(SUM) = DSQRT(SUM) 00007410 ABS(DETA) = DABS(DETA) 00007420 DATA SING/1D-50/ 00007430 C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 00007440 IF(NA .LT. 1) GO TO 18 00007450 DETA = A(1,1) 00007460 A(1,1) = SQRT(A(1,1)) 00007470 IF(NA .EQ. 1) GO TO 6 00007480 DO 1 I = 2,NA 00007490 1 A(I,1) = A(I,1) / A(1,1) 00007500 DO 5 J = 2,NA 00007510 SUM = 0. 00007520 J1 = J - 1 00007530 DO 2 K = 1,J1 00007540 2 SUM = SUM + A(J,K) ** 2 00007550 A(J,J) = SQRT(A(J,J) - SUM) 00007570 DETA = DETA * (A(J,J) - SUM) 00007560 IF(J .EQ. NA) GO TO 5 00007580 J2 = J + 1 00007590 DO 4 I = J2,NA 00007600 SUM = 0. 00007610 DO 3 K = 1,J1 00007620 3 SUM = SUM + A(I,K) * A(J,K) 00007630 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 00007640 5 CONTINUE 00007650 6 IF(ABS(DETA) .LT. SING) GO TO 16 00007660 C INVERSION OF LOWER TRIANGULAR MATRIX 00007670 DO 7 I = 1,NA 00007680 A(I,I) = 1. / A(I,I) 00007690 7 CONTINUE IF(NA .EQ. 1) GO TO 10 00007700 N1 = NA - 1 00007710 DO 9 J = 1,N1 00007720 J2 = J + 1 00007730 DO 9 I = J2,NA 00007740 SUM = 0. 00007750 I1 = I - 1 00007760 DO 8 K = J,I1 00007770 8 SUM = SUM + A(I,K) * A(K,J) 00007780 9 A(I,J) = - A(I,I) * SUM 00007790 C CONSTRUCTION OF INVERSE OF INPUT MATRIX 00007800 10 DO 15 J = 1,NA 00007810 IF(J .EQ. 1) GO TO 12 00007820 J1 = J - 1 00007830 DO 11 I = 1,J1 00007840 11 A(I,J) = A(J,I) 00007850 12 DO 14 I = J,NA 00007860 SUM = 0. 00007870 DO 13 K = I,NA 00007880 13 SUM = SUM + A(K,I) * A(K,J) 00007890 14 A(I,J) = SUM 00007900 15 CONTINUE 00007910 RETURN 00007920 16 WRITE(6,17) DETA 00007930 17 FORMAT(10X, 'SINGULAR MATRIX IN CHOLD. DET =',E20.5) 00007940 RETURN 1 00007950 18 WRITE(6,19) 00007960 19 FORMAT(10X,'MATRIX OF DIMENSION ZERO IN CHOLD') 00007970 RETURN 1 00007980 END 00007990 SUBROUTINE RADARC (A,I,J,S) 00008000 C * SUBROUTINE 'RADARC' CONVERTS DOUBLE PRECISION RADIANS TO SINGLE*00008010 C * PRECISION DEGREES,MINUTES AND SECONDS *00008020 C ******************************************************************00008030 REAL*8 A,SEC,AD,AJ,RHO/206264.8062470964D0/ 00008040 IF(A)2,1,1 00008050 1 SIGN=1. 00008060 GO TO 3 00008070 2 SIGN=-1. 00008080 A=-A 00008090 3 SEC=A*RHO+0.0005 00008100 I=SEC/3600. 00008110 AD=I 00008120 J=SEC/60.-AD*60. 00008130 AJ=J 00008140 S=SEC-AD*3600.-AJ*60.-0.0005 00008150 IF(I)5,4,5 00008160 4 IF(J)6,7,6 00008170 6 J=J*SIGN 00008180 GO TO 10 00008190 7 S=S*SIGN 00008200 GO TO 10 00008210 5 I=I*SIGN 00008220 10 CONTINUE 00008230 RETURN 00008240 END 00008250 SUBROUTINE PLHXYZ (X,Y,Z,PHI,DAMBDA,H) 00008260 C PLHXYZ COMPUTES GEODETIC COORDINATES FROM CARTESIAN COORDINATES 00008270 C (CLARKE 1866 ELLIPSOID ONLY) 00008280 IMPLICIT REAL*8 (A-H,O-Z) 00008290 A=6378206.4D0 00008300 B=6356583.8D0 00008310 ESQ=(A**2-B**2)/A**2 00008320 P=DSQRT(X**2+Y**2) 00008330 H=0.D0 00008340 PHI1=DATAN((Z/P)*(1.D0/(1.D0-ESQ))) 00008350 1 RN=A/DSQRT(1.D0-ESQ*(DSIN(PHI1)**2)) 00008360 H=(P/DCOS(PHI1))-RN 00008370 PHI=DATAN((Z/P)*(1.D0/(1.D0-ESQ*(RN/(RN+H))))) 00008380 IF((PHI-PHI1).LE.1.D-12) GO TO 2 00008390 PHI1=PHI 00008400 GO TO 1 00008410 2 H=(P/DCOS(PHI))-RN 00008420 DAMBDA=DATAN2(Y,X) 00008430 RETURN 00008440 END 00008450 SUBROUTINE NORMAL(A,B,PHI,XM,XN) C*********************************************************************** C* * C* THIS SUBROUTINE COMPUTES THE RADII OF CURVATURE ALONG THE MERIDIAN * C* (XM) AND ALONG THE PRIME VERTICAL (XN). A & B ARE THE SEMI-MAJOR * C* AND SEMI-MINOR AXES RESPECTIVELY. PHI IS THE GEODETIC LATITUDE. * C* SUPPLY VALUES FOR A , B & PHI IN THE MAIN PROGRAM * C* * C*********************************************************************** IMPLICIT REAL*8 (A-H,O-Z) E2=1.0D0-B**2/A**2 W=DSQRT(1.0D0-E2*DSIN(PHI)**2) XM=A*(1.0D0-E2)/W**3 XN=A/W RETURN END SUBROUTINE ORTHO(N,M,SIGMA,PHI,MRD,SIGMAF,VFC,NPC,INDEX,V,SUMD,F,W00007370 &,D,ALPHA, C,SUMC,SC2,STDP,IW) 00007380 C THIS SUBROUTINE ORTHOGONALIZES THE MATRIX PHI USING THE GRAM-SCHMIDT 00007390 C METHOD, COMPUTES THE FOURIER COEFFICIENTS OF THE ORTHOGONALIZED MATRIX00007400 C DERIVES THE COEFFICIENTS OF PHI,COMPUTES THE VARIANCES OF THE FOURIER 00007410 C COEFFICIENTS AND THE VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 00007420 C INPUTS : 00007430 C 1. PHI(OPTIONAL - COULD BE FUNCTION SUBPROGRAM INSTEAD) - AN N BY M 00007440 C CONTAINING THE BASE FUNCTIONS EVALUATED FOR EACH OBSERVATION 00007450 C 2. N - THE NUMBER OF OBSERVATIONS 00007460 C 3. M - THE NUMBER OF BASE FUNCTIONS (EQUAL OR GREATER THAN 2) 00007470 C 4. W - A VECTOR OF LENGTH N CONTAINING THE COMPUTED WEIGHT FUNCTIONS00007480 C 5. F - FUNCTIONAL VALUES 00007490 C 6. SIGMA - THE A PRIORI VARIANCE FACTOR 00007500 C 7. NRD - THE MAXIMUM ROW DIMENSION OF PHI 00007510 C 8. MRD - THE MAXIMUM COLUMN DIMENSION OF PHI 00007520 C 9. INDEX - PERMITS OPTIONAL TEST FOR STATISTICAL SIGNIFICANCE 00007530 C OF FOURIER COEFFICIENTS.... 00007540 C IF 0,STATISTICAL TEST FOR FOURIER COEFFICIENTS ABANDONED 00007550 C IF 1,TESTS AGAINST ONE TIME ITS STANDARD DEVIATION 00007560 C IF 2,TESTS AGAINST TWICE ITS ST.DEVIATION 00007570 C IF 3,TESTS AGAINST THREE TIMES ITS ST.DEVIATION 00007580 IMPLICIT REAL*8(A-H,O-Z) 00007750 C 10. IW - WRITE CODE OF THE COMPUTER 00007590 C OUTPUTS : 00007600 C 1. ALPHA - AN MRD BY M MATRIX CONTAINING THE ALPHA'S USED IN COMPUTI00007610 C THE ORTHOGONALIZED MATRIX AND IN COMPUTING THE COEFFICIENTS OF PH00007620 C 2. C - THE M FOURIER COEFFICIENTS OF THE ORTHOGONALIZED MATRIX 00007630 C 3. D - THE M COEFFICIENTS OF THE INPUT MATRIX PHI 00007640 C 4. SUMC - THE VARIANCES OF THE FOURIER COEFFICIENTS 00007650 C 5. SUMD - THE VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 00007660 C 6. SC2 - THE SQUARES OF THE NORMS OF THE ORTHOGONALIZED MATRIX 00007670 C 7. SIGMAF - THE FOURIER POLYNOMIAL A POSTEORI VARIANCE FACTOR 00007680 C 8. V - THE N RESIDUALS 00007690 C 9. VFC - THE ORIGINAL POLYNOMIAL A POSTEORI VARIANCE FACTOR 00007700 C 10. NPC - NUMBER OF THE COEFFICIENTS OF THE ORIGINAL POLYNOMIAL 00007710 C AFTER THE STATISTICAL TEST IS PERFORMED 00007720 C 11. STDP - VECTOR AGAINST WHICH THE ABSOLUTE VALUES OF FOURIER 00007730 C COEFFICIENTS ARE TESTED 00007740 DIMENSION ALPHA(MRD,M),W(N),F(N),C(M),D(M) 00007770 DIMENSION SUMD(MRD,M),SUMC(M),SC2(M),V(N) 00007780 DIMENSION STDP(M),PHI(N,M) C TEST FOR NEGATIVE DEGREES OF FREEDOM 00007800 IF (N.LT.M) GO TO 100 00007810 K=1 00007820 ALPHA(M,M)=1.D0 00007830 C DETERMINE THE ALPHA'S FOR COMPUTATION OF ORTHOGONALIZED MATRIX 00007840 10 DO 3 J=K,M 00007850 IF(J.NE.K) GO TO 6 00007860 ALPHA(K,K)=1.D0 00007870 GO TO 3 00007880 6 SC1=0.D0 00007890 SC2(K)=0.D0 00007900 SC3=0.D0 00007910 DO 2 I=1,N 00007920 P=PHI(I,K) 00007940 IF(K.EQ.1) GO TO 4 00007950 K1=K-1 00007960 DO 5 J1=1,K1 00007970 5 P=P+ALPHA(J1,K)*PHI(I,J1) 00007980 4 SC1=SC1+W(I)*PHI(I,J)*P 00007990 SC3=SC3+F(I)*W(I)*P 00008000 2 SC2(K)=SC2(K)+W(I)*P**2 00008010 ALPHA(J,K)=-SC1/SC2(K) 00008020 ALPHA(K,J)=ALPHA(J,K) 00008030 3 CONTINUE 00008040 C DETERMINE THE FOURIER COEFFICIENTS FOR THE ORTHOGONALIZED MATRIX 00008050 C(K)=SC3/SC2(K) 00008060 K=K+1 00008070 IF(M.EQ.2) GO TO 34 00008080 IF(K.LT.3) GO TO 10 00008090 C DETERMINE THE ALPHA'S USED IN COMPUTING THE COEFFICIENTS OF PHI 00008100 JK=K-1 00008110 9 JL=K 00008120 JK=JK-1 00008130 JJ=K-JK-1 00008140 DO 8 LM=1,JJ 00008150 JL=JL-1 00008160 8 ALPHA(JK,K)=ALPHA(JK,K)+ALPHA(JK,JL)*ALPHA(K,JL) 00008170 IF(JK.NE.1) GO TO 9 00008180 IF(K.LT.M) GO TO 10 00008190 C DETERMINE THE LAST FOURIER COEFFICIENT 00008200 34 SC2(K)=0.D0 00008210 SC3=0.D0 00008220 DO 7 I=1,N 00008230 P=PHI(I,K) 00008250 K1=K-1 00008260 DO 1 J=1,K1 00008270 1 P=P+ALPHA(J,K)*PHI(I,J) 00008280 SC2(K)=SC2(K)+W(I)*P**2 00008290 7 SC3=SC3+F(I)*W(I)*P 00008300 C(K)=SC3/SC2(K) 00008310 C DETERMINE THE COEFFICIENTS OF PHI 00008320 IDEKT=1 00008330 ICOUNT=0 00008340 1000 CONTINUE 00008350 DO 13 I=1,M 00008360 D(I)=C(I) 00008370 IF(I.EQ.M) GO TO 13 00008380 II=I+1 00008390 DO 14 J=II,M 00008400 14 D(I)=D(I)+ALPHA(I,J)*C(J) 00008410 13 CONTINUE 00008420 C COMPUTE THE VARIANCE OF THE FOURIER COEFFICIENTS AND THE VARIANCE-COVA00008430 C MATRIX OF THE COEFFICIENTS 00008440 DO 15 I=1,M 00008450 DO 15 J=1,M 00008460 15 SUMD(I,J)=0.D0 00008470 SC4=0.D0 00008480 DO 22 I=1,N 00008490 PN=0.D0 00008510 DO 21 J=1,M 00008520 21 PN=PN+D(J)*PHI(I,J) 00008530 V(I)=F(I)-PN 00008540 V2=V(I)**2 00008550 22 SC4=SC4+V2*W(I) 00008560 SIGMAF=SC4/(N-M+ICOUNT)*SIGMA 00008570 VFC=SIGMAF 00008580 IF(IDEKT.EQ.2) VFC=SC4/(N-NPC)*SIGMA 00008590 IF(INDEX.EQ.0) NPC=M 00008600 DO 28 I=1,M 00008610 SUMC(I)=SIGMAF/SC2(I) 00008620 IF(IDEKT.EQ.1) GO TO 28 00008630 IF(C(I).EQ.0D0) SUMC(I)=0D0 00008640 28 CONTINUE 00008650 DO 23 I=1,M 00008660 DO 23 J=1,I 00008670 DO 23 K=J,I 00008680 23 SUMD(J,K)=SUMD(J,K)+ALPHA(J,I)*ALPHA(K,I)*SUMC(I) 00008690 DO 24 I=1,M 00008700 IT=I+1 00008710 IF(IT.GT.M) GO TO 30 00008720 DO 24 J=IT,M 00008730 24 SUMD(J,I)=SUMD(I,J) 00008740 30 CONTINUE 00008760 C OPTIONAL CHECK FOR STATISTICALLY SIGNIFICANT FOURIER COEFFICIENTS 00008750 IF(INDEX.EQ.0) GO TO 40 00008770 IF(IDEKT.EQ.2) GO TO 40 00008780 PINDEX=DFLOAT(INDEX) 00008790 DO 31 I=1,M 00008800 STDP(I)=PINDEX*DSQRT(SUMC(I)) 00008810 IF(DABS(C(I)).LT.STDP(I)) GO TO 32 00008820 GO TO 31 00008830 32 C(I)=0D0 00008840 ICOUNT=ICOUNT+1 00008850 SUMC(I)=0D0 00008860 31 CONTINUE 00008870 NPC=0 00008880 DO 33 I=1,M 00008890 IF(C(I).NE.0D0) NPC=I 00008900 33 CONTINUE 00008910 IDEKT=2 00008920 GO TO 1000 00008930 40 RETURN 00008940 100 WRITE(IW,102) 00008950 102 FORMAT('0','*ERROR* NEGATIVE DEGREES OF FREEDOM') 00008960 RETURN 00008970 END 00008980 SUBROUTINE DSTORT(X1,Y1,Z1,NU2,C,SIGMAC,STDC,KO) C THIS SUBROUTINE COMPUTES THE DISTORTION IN X Y Z COORDINATES OF POINTS C USING THE COEFFICIENTS DERIVED FROM LEAST SQUARES APPROXIMATION. C THE BASE FUNCTIONS ARE DERIVED FROM CARTESIAN COORDINATES IN MEGA METERS C C CARD 1 : ASSUMED HEIGHT OF STATION (F10.2) C CONTROL CARD: C CARD 2 COL.1-3 ICONT(I4)=1: DISTORTION DUE TO ALL TRANSFORM- C ATION PARAMETERS FROM GEODETIC C TO AVERAGE TERRESTRIAL C =2: DISTORTION DUE TO SYSTEMATIC C ERRORS IN GEODETIC NETWORK ALONE C COL.5-8 NUMB(I4) =NUMBER OF STATIONS (DATA POINTS) C COORDINATES C CARD 3-47: C COL.1-4 STA(A4)=STATION NAME C COL.5-9 PHI(F5.0)=LATITUDE(DEGREES) C COL.10-14 DLAM(F5.0)=LONGITUDE (DEGREES) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION STA(44),PHI(44),DLAM(44),CX(8),CY(8),CZ(8),X(44), 1 Y(44),Z(44),PNX(44),PNY(44),PNZ(44),DPHI(44),DILAM(44),DH(44) DIMENSION C(24),SIGMAC(24,24),STDC(24),G(132,24),SIGJ(132,132) DIMENSION GSA(132,24),GSAT(132,132),SGPB(132,132),SGPT(132,132) DIMENSION SIGG(132),SGPF(132),RP(44),XM(44),XN(44),SDPHI(44), 1SDLAM(44),R(44), DJ(132,132),DJPT(132,132),DJT(132,132),DJTF(132) DIMENSION RL(44) IF(KO.GE.1) GO TO 31 C READ IN HEIGHT READ(5,1) HI 1 FORMAT(F10.2) C READ IN THE CONTROL NUMBERS READ(5,10) ICONT,NUMB 10 FORMAT(2I4) C READ STATION NUMBERS AND THEIR GEODETIC COORDINATES IN DEGREES DO 20 I=1,NUMB 20 READ(5,30) STA(I),PHI(I),DLAM(I) 30 FORMAT(A4,2F5.0) 31 CONTINUE NT3=NUMB*3 DO 40 I=1,8 CX(I)=C(I) CY(I)=C(I+8) CZ(I)=C(I+16) 40 CONTINUE IF(ICONT.EQ.2) GO TO 57 WRITE(6,51) 51 FORMAT(///,10X,'**DISTORTION DUE TO ALL TRANSFORMATION PARAMETERS 1FROM GEODETIC TO AVERAGE TERRESTRIAL SYSTEM**',//) GO TO 59 57 WRITE(6,58) 58 FORMAT(///,20X,'**DISTORTION DUE TO SYSTEMATIC ERRORS IN TERRESTRI 1AL NETWORK ONLY**',///) 59 CONTINUE AA=6378206.4D0 BB=6356583.8D0 E2=1.D0-BB**2/AA**2 PI=DARCOS(-1.D0) RHO=(180.D0/PI)*3600.D0 C COMPUTE CARTESIAN COORDINATES OF STATIONS WRITE(6,52) 52 FORMAT(//,6X,'STA.NAME',6X,'X(M)',9X,'Y(M)',9X,'Z(M)',//) DO 60 I=1,NUMB RP(I)=PHI(I)*3600.D0/RHO RL(I)=DLAM(I)*3600.D0/RHO CALL NORMAL(AA,BB,RP(I),XM(I),XN(I)) X(I)=(XN(I)+HI)*DCOS(RP(I))*DCOS(RL(I)) Y(I)=(XN(I)+HI)*DCOS(RP(I))*DSIN(RL(I)) Z(I)=(XN(I)*(1.D0-E2)+HI)*DSIN(RP(I)) 60 WRITE(6,70) STA(I),X(I),Y(I),Z(I) 70 FORMAT('0',8X,A5,3F15.3) WRITE(6,413) 413 FORMAT(///,40X,'COEFFICIENTS OF ORIGINAL POLYNOMIALS',//,13X, 1 'COEFF.NO.',4X,'CX',8X,'STD.DEV.',9X,'CY',8X,'STD.DEV.',9X, 2 'CZ',8X,'STD.DEV.') WRITE(6,414) (I,C(I),STDC(I),C(I+8),STDC(I+8),C(I+16),STDC(I+16), 1I=1,8) 414 FORMAT('0',15X,I2,'.',F12.6,F13.6,F14.6,F13.6,F14.6,F13.6) WRITE(6,72) 72 FORMAT(//,6X,'STA.NAME',6X,'LAT(DEG)',4X,'LONG(DEG)',4X,'HT(M)',/) DO 64 I=1,NUMB 64 WRITE(6,61) STA(I),PHI(I),DLAM(I),HI 61 FORMAT('0',8X,A5,5F12.3) DO 3101 I=1,NT3 DO 3101 J=1,NU2 3101 G(I,J)=0.0D0 DO 4111 I=1,NT3 DO 4111 J=1,NT3 4111 SIGJ(I,J)=0.D0 C COMPUTE POLYNOMIAL FIT EFFECT DUE TO ALL TRANSFORMATION PARAMETERS IF(ICONT.EQ.2) GO TO 100 DO 80 I=1,NUMB X(I)=X(I)*1.D-6 Y(I)=Y(I)*1.D-6 Z(I)=Z(I)*1.D-6 FI1=1.D0 FI2=X(I) FI3=Y(I) FI4=Z(I) FI5=X(I)*Y(I) FI6=X(I)*Z(I) FI7=Y(I)*Z(I) FI8=X(I)*Y(I)*Z(I) PNX(I)=FI1*CX(1)+FI2*CX(2)+FI3*CX(3)+FI4*CX(4)+FI5*CX(5)+FI6*CX(6) 1 +FI7*CX(7)+FI8*CX(8) PNY(I)=FI1*CY(1)+FI2*CY(2)+FI3*CY(3)+FI4*CY(4)+FI5*CY(5)+FI6*CY(6) 1 +FI7*CY(7)+FI8*CY(8) PNZ(I)=FI1*CZ(1)+FI2*CZ(2)+FI3*CZ(3)+FI4*CZ(4)+FI5*CZ(5)+FI6*CZ(6) 1 +FI7*CZ(7)+FI8*CZ(8) G(3*I-2,1)=FI1 G(3*I-2,2)=FI2 G(3*I-2,3)=FI3 G(3*I-2,4)=FI4 G(3*I-2,5)=FI5 G(3*I-2,6)=FI6 G(3*I-2,7)=FI7 G(3*I-2,8)=FI8 G(3*I-1,9)=FI1 G(3*I-1,10)=FI2 G(3*I-1,11)=FI3 G(3*I-1,12)=FI4 G(3*I-1,13)=FI5 G(3*I-1,14)=FI6 G(3*I-1,15)=FI7 G(3*I-1,16)=FI8 G(3*I,17)=FI1 G(3*I,18)=FI2 G(3*I,19)=FI3 G(3*I,20)=FI4 G(3*I,21)=FI5 G(3*I,22)=FI6 G(3*I,23)=FI7 G(3*I,24)=FI8 X(I)=X(I)*1.D6 Y(I)=Y(I)*1.D6 Z(I)=Z(I)*1.D6 80 CONTINUE GO TO 120 100 CONTINUE C COMPUTE EFFECT DUE TO SYSTEMATIC ERRORS IN NETWORK ALONE DO 110 I=1,NUMB XKI=(X(I)-X1)*1.D-6 YKI=(Y(I)-Y1)*1.D-6 ZKI=(Z(I)-Z1)*1.D-6 FI1=1.D0 FI2=XKI FI3=YKI FI4=ZKI FI5=XKI*YKI FI6=XKI*ZKI FI7=YKI*ZKI FI8=XKI*YKI*ZKI PNX(I)=FI1*CX(1)+FI2*CX(2)+FI3*CX(3)+FI4*CX(4)+FI5*CX(5)+FI6*CX(6) 1 +FI7*CX(7)+FI8*CX(8) PNY(I)=FI1*CY(1)+FI2*CY(2)+FI3*CY(3)+FI4*CY(4)+FI5*CY(5)+FI6*CY(6) 1 +FI7*CY(7)+FI8*CY(8) PNZ(I)=FI1*CZ(1)+FI2*CZ(2)+FI3*CZ(3)+FI4*CZ(4)+FI5*CZ(5)+FI6*CZ(6) G(3*I-2,1)=FI1 G(3*I-2,2)=FI2 G(3*I-2,3)=FI3 G(3*I-2,4)=FI4 G(3*I-2,5)=FI5 G(3*I-2,6)=FI6 G(3*I-2,7)=FI7 G(3*I-2,8)=FI8 G(3*I-1,9)=FI1 G(3*I-1,10)=FI2 G(3*I-1,11)=FI3 G(3*I-1,12)=FI4 G(3*I-1,13)=FI5 G(3*I-1,14)=FI6 G(3*I-1,15)=FI7 G(3*I-1,16)=FI8 G(3*I,17)=FI1 G(3*I,18)=FI2 G(3*I,19)=FI3 G(3*I,20)=FI4 G(3*I,21)=FI5 G(3*I,22)=FI6 G(3*I,23)=FI7 G(3*I,24)=FI8 110 CONTINUE 120 CONTINUE DO130 I=1,NUMB H1=XM(I)+HI H2=(XN(I)+HI)*DCOS(RP(I)) DPHI(I)=(RHO/H1)*(-DSIN(RP(I))*DCOS(RL(I))*PNX(I)-DSIN(RP(I)) 1*DSIN(RL(I))*PNY(I)+DCOS(RP(I))*PNZ(I)) DILAM(I)=(RHO/H2)*(-DSIN(RL(I))*PNX(I)+DCOS(RL(I))*PNY(I)) DH(I)=DCOS(RP(I))*DCOS(RL(I))*PNX(I)+DCOS(RP(I))*DSIN(RL(I)) 1*PNY(I)+DSIN(RP(I))*PNZ(I) SIGJ(3*I-2,3*I-2)=-DSIN(RP(I))*DCOS(RL(I))/H1 SIGJ(3*I-2,3*I-1)=-DSIN(RP(I))*DSIN(RL(I))/H1 SIGJ(3*I-2,3*I)=DCOS(RP(I))/H1 SIGJ(3*I-1,3*I-2)=-DSIN(RL(I))/H2 SIGJ(3*I-1,3*I-1)=DCOS(RL(I))/H2 SIGJ(3*I,3*I-2)=DCOS(RP(I))*DCOS(RL(I)) SIGJ(3*I,3*I-1)=DCOS(RP(I))*DSIN(RL(I)) SIGJ(3*I,3*I)=DSIN(RP(I)) 130 CONTINUE C FORM VARIANCE COVARIANCE MATRICES DO 4115 I=1,NT3 DO 4115 J=1,NU2 4115 GSA(I,J)=0.D0 DO 4112 I=1,NT3 DO 4112 J=1,NT3 4112 GSAT(I,J)=0.D0 DO 4005 I=1,NT3 DO 4005 J=1,NT3 4005 SGPB(I,J)=0.D0 DO 4006 I=1,NT3 DO 4006 J=1,NT3 4006 SGPT(I,J)=0.D0 DO 4113 I=1,NT3 DO 4113 J=1,NU2 DO 4113 K=1,NU2 4113 GSA(I,J)=GSA(I,J)+G(I,K)*SIGMAC(K,J) DO 4114 I=1,NT3 DO 4114 J=1,NT3 DO 4114 K=1,NU2 4114 GSAT(I,J)=GSAT(I,J)+GSA(I,K)*G(J,K) DO 3111 I=1,NT3 SIGG(I)=DSQRT(GSAT(I,I)) 3111 CONTINUE DO 4007 I=1,NT3 DO 4007 J=1,NT3 DO 4007 K=1,NT3 4007 SGPB(I,J)=SGPB(I,J)+SIGJ(I,K)*GSAT(K,J) DO 4008 I=1,NT3 DO 4008 J=1,NT3 DO 4008 K=1,NT3 4008 SGPT(I,J)=SGPT(I,J)+SGPB(I,K)*SIGJ(J,K) DO 4009 I=1,NUMB DO 4009 J=1,2 4009 SGPF(3*I-J)=DSQRT(SGPT(3*I-J,3*I-J))*RHO DO 4010 I=1,NUMB 4010 SGPF(3*I)=DSQRT(SGPT(3*I,3*I)) WRITE(6,4011) WRITE(6,3112) 4011 FORMAT('-',10X,'FUNCTIONAL VALUES OF POLYNOMIALS IN CURVILINEAR CO 1-ORDINATES',/) 3112 FORMAT(///,7X,'STA.NAME',6X,'LAT(SEC)',6X,'STD.DEV.',6X,'LONG(SEC) 1',5X,'STA.DEV.',6X,'HT(M)',4X,'STD.DEV.',//) DO 4012 I=1,NUMB WRITE(6,3114) STA(I),DPHI(I),SGPF(3*I-2),DILAM(I),SGPF(3*I-1), 1DH(I),SGPF(3*I) 3114 FORMAT('0',8X,A5,3X,F12.4,3X,F10.4,3X,F12.4,3X,F10.4,3X,F10.4, 1F11.4) 4012 CONTINUE C APPLY LAMES COEFFICIETS TO DETERMINE THE DIFFERENTIAL CHANGES IN C DISTANCE DS ALONG THE CURVILINEARBCOODINATES LINES DO 200 I=1,NUMB R(I)=DSQRT(XM(I)*XN(I)) SDPHI(I)=R(I)/RHO*DPHI(I) SDLAM(I)=R(I)/RHO*DCOS(RP(I))*DILAM(I) 200 CONTINUE C CALCULATE STANDARD DEVIATION OF THE TRANSFORMED DISTANCES DO 210 I=1,NT3 DO 210 J=1,NT3 210 DJ(I,J)=0.D0 DO 220 I=1,NT3 DO 220 J=1,NT3 220 DJPT(I,J)=0.D0 DO 230 I=1,NT3 DO 230 J=1,NT3 230 DJT(I,J)=0.D0 DO 240 I=1,NUMB DJ(3*I-2,3*I-2)=R(I) DJ(3*I-1,3*I-1)=R(I)*DCOS(RP(I)) 240 DJ(3*I,3*I)=1.D0 DO 250 I=1,NT3 DO 250 J=1,NT3 DO 250 K=1,NT3 250 DJPT(I,J)=DJPT(I,J)+DJ(I,K)*SGPT(K,J) DO 260 I=1,NT3 DO 260 J=1,NT3 DO 260 K=1,NT3 260 DJT(I,J)=DJT(I,J)+DJPT(I,K)*DJ(J,K) DO 270 I=1,NT3 270 DJTF(I)=DSQRT(DJT(I,I)) C PRINT DISTORTIONS IN METERS WRITE(6,280) 280 FORMAT(//,10X,'DISTORTION EXPRESSED IN DISTANCE ( METERS )') WRITE(6,290) 290 FORMAT(///,7X,'STA.NAME',7X,'LAT(M)',6X,'STD.DEV.',7X,'LONG(M)', 16X,'STD.DEV.',6X,'HT.(M)',5X,'STD.DEV.',/) WRITE(6,3114) (STA(I),SDPHI(I),DJTF(3*I-2),SDLAM(I),DJTF(3*I-1), 1 DH(I),DJTF(3*I),I=1,NUMB) RETURN END