C LTC134 C C INSTRUCTIONS FOR THE USE OF THIS PROGRAMME AVAILABLE FROM SURVEYING C DEPARTMENT C C LTC134 GEODETIC TRAVERSE PROGRAM 13400000 C 13400010 IMPLICIT REAL*8 (A-H,O-Z) 13400020 DIMENSION Y(101),X(101),YK(500),XK(500),AZ(100),DIS(100),NUM(101) 13400030 DIMENSION KA(500),KB(500),KC(500),KN(500),NA(101),NB(101),NC(101) 13400031 INTEGER TITLE(10),ER,BLK,ICC(15),FT,SER,ERX,PROC,Z 13400040 DATA ER/' ER'/,BLK/' '/,NN7/' 997'/,FT/'F'/,NN8/' 998'/,M/'M'/13400050 DATA SER/'SERI'/,ERX/'AL X'/,NOT/'NOT '/,PROC/'PROC'/,KN/500*0/ 13400060 COMMON ALAT,ALON,BLAT,BLON,AZAB,AZBA,S,IZ 13400070 DEFINE FILE 4(200,12,U,IR) 13400080 RADS(A,B,C)=.0174532925199D0*A+.290888208666D-3*B+.4848136811D-5*C13400090 CALL REREAD 13470010 C 13470940 C READ LEAD CARDS 13470950 C 13470960 Z=1 13470961 10 DO11 I=1,101 13470962 11 NUM(I)=0 13470963 WRITE(6,12) 13470970 12 FORMAT(1H1) 13470980 390 IER=1 13470990 KER=BLK 13470991 LER=BLK 13470992 READ(5,391,END=1100)ISER 13470993 391 FORMAT(I3) 13470994 IF(ISER.EQ.997)GO TO 392 13470995 LER=NN7 13470996 KER=ER 13470997 IER=2 13470998 392 READ(99,402)ICC 13470999 WRITE(6,405)ICC,KER,LER 13471000 IF(IER.EQ.2)GO TO 390 13471001 400 READ(99,401)ISER,IND1,IND2,NL,NSS,AEANG,AELAT,AELON,MFO,TITLE 13471002 401 FORMAT(I3,2I1,2I3,F3.0,2F3.2,A1,10A4) 13471010 402 FORMAT(15A4) 13471030 405 FORMAT(1H ,17A4) 13471090 406 IND1=IND1+1 13471120 IND3=0 13471121 IND4=0 13471122 IND2=IND2+1 13471130 IF(NL.EQ.0)GO TO 1200 13471131 IF(MFO.EQ.M)GO TO 408 13471140 407 MFO=2 13471150 GO TO 410 13471160 408 MFO=1 13471170 410 L=NL+1 13471180 READ (5,411)ISER,NA(1),NB(1),NC(1),NUM(1),A,B,C,D,E,F,P,Q,R,S,T,U 13471190 411 FORMAT (I3,3A3,I6,2(2F2.0,F6.4,F3.0,F2.0,F6.4)) 13471200 READ(99,402)ICC 13471210 IF(ISER.NE.998)GO TO 426 13471215 Y(1)=RADS(A,B,C) 13471220 X(1)=RADS(D,E,F) 13471225 IF(Y(1).NE.0.D0)GO TO 416 13471230 NSTA=NUM(1) 13471235 IF(NSTA.EQ.0)GO TO 426 13471240 DO412 K=1,500 13471245 IF(NSTA.EQ.KN(K))GO TO 414 13471250 412 CONTINUE 13471255 413 WRITE(6,926)NSTA 13471260 GO TO 426 13471265 414 Y(1)=YK(K) 13471270 X(1)=XK(K) 13471275 IND3=1 13471276 416 BLAT=RADS(P,Q,R) 13471280 BLON=RADS(S,T,U) 13471285 IF(BLAT.NE.0.D0)GO TO 418 13471290 IF(BLON.EQ.0.D0)GO TO 426 13471295 BACKAZ=BLON 13471300 GO TO 430 13471305 418 IF(BLON.NE.0.D0)GO TO 424 13471310 NSTA=(R+.00001D0)*10000.D0 13471315 DO420 K=1,500 13471320 IF(NSTA.EQ.KN(K))GO TO 422 13471325 420 CONTINUE 13471330 GO TO 413 13471335 422 BLAT=YK(K) 13471340 BLON=XK(K) 13471345 424 ALAT=Y(1) 13471350 ALON=X(1) 13471355 IZ=1 13471360 CALL LTCS04 13471365 BACKAZ=AZAB 13471370 GO TO 430 13471375 426 LER=NN8 13471380 KER=ER 13471385 IER=2 13471390 430 WRITE(6,405)ICC,KER,LER 13471392 IF(IER.EQ.2)GO TO 390 13471394 440 IF(IND1.EQ.1.AND.IND2.EQ.1)GO TO 490 13471400 C 13471420 C READ CLOSING STATION DATA 13471430 C 13471440 450 READ (5,411)ISER,NA(L),NB(L),NC(L),NUM(L),A,B,C,D,E,F,P,Q,R,S,T,U 13471450 READ(99,402)ICC 13471455 IF(ISER.NE.999)GO TO 480 13471460 CLOLAT=RADS(A,B,C) 13471465 CLOLON=RADS(D,E,F) 13471470 IND4=2 13471471 IF(CLOLAT.NE.0.D0)GO TO 462 13471475 NSTA=NUM(L) 13471480 IF(NSTA.EQ.0)GO TO 480 13471485 DO455 K=1,500 13471490 IF(NSTA.EQ.KN(K))GO TO 460 13471495 455 CONTINUE 13471500 GO TO 472 13471510 460 CLOLAT=YK(K) 13471515 CLOLON=XK(K) 13471520 IND4=1 13471521 462 IF(IND1.EQ.1)GO TO 485 13471525 BLAT=RADS(P,Q,R) 13471530 BLON=RADS(S,T,U) 13471535 IF(BLAT.NE.0.D0)GO TO 465 13471540 IF(BLON.EQ.0.D0)GO TO 480 13471545 FOREAZ=BLON 13471550 GO TO 485 13471555 465 IF(BLON.NE.0.D0)GO TO 477 13471560 NSTA=(R+.00001D0)*10000.D0 13471565 DO470 K=1,500 13471570 IF(NSTA.EQ.KN(K))GO TO 475 13471575 470 CONTINUE 13471580 472 WRITE(6,926)NSTA 13471585 GO TO 480 13471590 475 BLAT=YK(K) 13471595 BLON=XK(K) 13471600 477 ALAT=CLOLAT 13471605 ALON=CLOLON 13471610 IZ=1 13471615 CALL LTCS04 13471620 FOREAZ=AZAB 13471625 GO TO 485 13471630 480 LER=NNN 13471635 KER=ER 13471640 IER=2 13471645 485 WRITE(6,405)ICC,KER,LER 13471650 IF(IER.EQ.2)GO TO 390 13471655 490 AZBA=BACKAZ 13471660 500 IZ=0 13471670 CUDIS=0. 13471690 C 13471691 C READ AND COMPUTE TRAVERSE CARDS 13471692 C 13471693 DO 560 I=1,NL 13471700 J=I+1 13471710 READ (5,501)ISER,A,B,C,S,MF,NA(J),NB(J),NC(J),NUM(J) 13471720 501 FORMAT (I3,F3.0,F2.0,F3.1,F9.3,A1,3A3,5I6) 13471730 READ(99,402)ICC 13471735 IF(ISER.EQ.I)GO TO 505 13471740 IER=2 13471750 KER=SER 13471760 LER=ERX 13471770 505 WRITE(6,405)ICC,KER,LER 13471780 IF(IER.EQ.2)GO TO 390 13471790 510 AZAB=RADS(A,B,C)+AZBA 13471800 IF(AZAB.GE.6.28318530718D0)AZAB=AZAB-6.28318530718D0 13471820 525 AZ(I)=AZAB 13471840 GO TO (530,535),MFO 13471850 530 IF(MF-FT)545,540,545 13471860 535 IF(MF- M)540,545,540 13471870 540 S=S/3.28083333333 13471880 545 DIS(I)=S 13471890 548 CUDIS=CUDIS+S 13471920 ALAT=Y(I) 13471930 ALON=X(I) 13471940 CALL LTCS04 13471950 550 Y(J)=BLAT 13471960 X(J)=BLON 13471970 560 CONTINUE 13471980 IF(IND1.EQ.1)GO TO 1200 13471990 580 READ (5,501)ISER,A,B,C,S,MF,M1,M2,M3,M4,N1,N2,JF,JT 13472000 READ(99,402)ICC 13472005 IF(ISER.EQ.L)GO TO 585 13472010 KER=SER 13472015 LER=ERX 13472020 IER=2 13472025 585 WRITE(6,405)ICC,KER,LER 13472030 KER=BLK 13472035 LER=BLK 13472040 IF(S.EQ.0.D0.AND.JF.EQ.0)GO TO 1200 13472050 590 WRITE(3,591) 13472060 591 FORMAT(22H CLOSING ANGLE MISSING) 13472065 IND1=1 13472066 IND2=1 13472067 C 13472080 C READ, EDIT, PRINT AND RECORD SIDE OBSERVATIONS 13472090 C 13472100 1200 IF (NSS)1290,1290,1205 13472110 1205 K=NL+IND1 13472120 LL=K+NSS-1 13472130 JER=1 13472140 IR=1 13472150 DO 1280 I=K,LL 13472160 READ(5,501)ISER,D,E,F,S,MF,M1,M2,M3,M4,N1,N2,JF,JT 13472170 READ(99,402)ICC 13472180 IF(ISER.EQ.I)GO TO 1210 13472190 KER=SER 13472200 LER=ERX 13472205 1210 WRITE(6,405)ICC,KER,LER 13472210 KER=BLK 13472215 LER=BLK 13472220 IF(ISER.NE.0)GO TO 1230 13472222 NSS=I-K+1 13472224 GO TO 1290 13472226 1230 D=RADS(D,E,F) 13472230 GO TO (1240,1250),MFO 13472240 1240 IF(MF-FT)1270,1260,1270 13472250 1250 IF(MF- M)1260,1270,1260 13472260 1260 S=S/3.28083333333D0 13472270 1270 WRITE(8'IR)D,S,M1,M2,M3,M4,N1,N2,JF,JT 13472275 1280 CONTINUE 13472280 1290 IF(IER.EQ.2)GO TO 10 13472290 IF(NL.EQ.0)GO TO 890 13472300 600 WRITE (6,601)TITLE 13472310 601 FORMAT(1H1,10A4/1H ) 13472320 LINE=3 13472330 GO TO (710,610),IND1 13472340 C 13472341 C COMPUTE ANGULAR MISCLOSURE 13472342 C 13472343 610 AMIS=RADS(A,B,C) 13472350 AMIS=FOREAZ-AMIS-AZBA 13472360 IF(AMIS+6.28318530718D0)615,615,620 13472370 615 AMIS=AMIS+6.28318530718D0 13472380 620 IF(AMIS+3.14159265359D0)621,625,622 13472390 621 AMIS=AMIS+6.28318530718D0 13472400 GO TO 625 13472410 622 IF(AMIS-3.14159265359D0)625,625,623 13472420 623 AMIS=AMIS-6.28318530718D0 13472430 625 A=AMIS*206264.806247D0 13472440 WRITE (6,630)A 13472450 630 FORMAT (19H ANGULAR MISCLOSURE,F11.2,8H SECONDS) 13472460 LINE=LINE+3 13472470 IF(DABS(A).LE.AEANG)GO TO 650 13472480 WRITE (6,631) 13472490 631 FORMAT (20H ANGLES NOT BALANCED/) 13472500 640 IND2=1 13472510 GO TO 790 13472520 C 13472530 C ADJUST AZIMUTHS AND RECOMPUTE TRAVERSE 13472540 C 13472550 650 WRITE (6,651) 13472560 651 FORMAT (16H ANGLES BALANCED/) 13472570 D=L 13472580 AMIS=AMIS/D 13472590 IZ=0 13472600 DO 700 I=1,NL 13472620 J=I+1 13472630 F=I 13472640 AZAB=AZ(I)+AMIS*F 13472650 IF (AZAB)660,675,665 13472660 660 AZAB=AZAB+6.28318530718D0 13472670 GO TO 675 13472680 665 IF(AZAB-6.28318530718D0)675,670,670 13472690 670 AZAB=AZAB+6.28318530718D0 13472700 675 AZ(I)=AZAB 13472710 685 S=DIS(I) 13472760 ALAT=Y(I) 13472770 ALON=X(I) 13472780 CALL LTCS04 13472790 690 Y(J)=BLAT 13472800 X(J)=BLON 13472810 700 CONTINUE 13472820 710 GO TO (790,720),IND2 13472830 720 CFLAT=CLOLAT-Y(J) 13472840 CFLON=CLOLON-X(J) 13472850 DISMIS=DSQRT((CFLAT*6378200.D0)**2+(CFLON*6392719.D0*DCOS(CLOLAT))13472860 1**2) 13472861 FRAC=CUDIS/DISMIS 13472870 A=CFLAT*206264.806247D0 13472880 B=CFLON*206264.806247D0 13472890 WRITE (6,730)A,B,DISMIS,FRAC 13472900 730 FORMAT (28H MISCLOSURES IN SECONDS LAT,F12.4,6H LONG,F12.4/7H LI13472910 1NEAR,F11.3,13H METRES OR 1/,F9.0/) 13472920 LINE=LINE+3 13472930 IF(DABS(A).GT.AELAT.OR.DABS(B).GT.AELON)IND3=1 13472940 IF(IND3.EQ.1)GO TO 640 13472950 740 CFLAT=CFLAT/CUDIS 13472960 CFLON=(6392719.D0*DCOS(CLOLAT)*CFLON)/CUDIS 13472970 C 13473010 C BALANCE COORDINATES 13473020 C 13473030 760 CORY=0.D0 13473040 CORX=0.D0 13473050 DO 770 I=2,NL 13473060 J=I-1 13473070 CORY=CORY+DIS(J)*CFLAT 13473080 Y(I)=Y(I)+CORY 13473090 CORX=CORX+DIS(J)*CFLON 13473100 770 X(I)=X(I)+CORX/(6392719.D0*DCOS(Y(I))) 13473110 Y(NL+1)=CLOLAT 13473120 X(NL+1)=CLOLON 13473130 C 13473140 C PRINT HEADINGS 13473150 790 IGO=1 13473160 LINE=LINE+3 13473170 GO TO 802 13473180 800 WRITE (6,601)TITLE 13473190 LINE=6 13473200 802 GO TO (805,810,815),IND2 13473210 805 WRITE (6,806) 13473220 806 FORMAT (1H ,5X,7HSTATION,58X,10HUNADJUSTED) 13473230 GO TO 815 13473240 810 WRITE (6,811) 13473250 811 FORMAT (1H ,5X,7HSTATION,59X,8HADJUSTED) 13473260 815 WRITE (6,816) 13473270 816 FORMAT (1H ,2X,4HNAME,5X,6HNUMBER,5X,7HAZIMUTH,5X,8HDISTANCE,5X, 13473280 124HREV AZIMUTH LATITUDE,7X,9HLONGITUDE/) 13473290 GO TO (830,845,952,1005),IGO 13473300 C 13473310 830 IGO=2 13473320 IZ=1 13473330 DO 880 I=1,L 13473350 832 YI=Y(I) 13473360 XI=X(I) 13473370 CALL LTCS03(YI,M1,M2,M3,M4,SM) 13473380 CALL LTCS03(XI,N1,N2,N3,N4,SN) 13473390 835 WRITE(6,836)NA(I),NB(I),NC(I),NUM(I),M1,M2,M3,M4,SM,N1,N2,N3,N4,SN13473400 836 FORMAT(1H ,3A3,I8,41X,2(2X,I3,I2,I1,I2,F6.4)) 13473420 837 IF(NUM(I).EQ.0)GO TO 838 13473421 IF(I.EQ.1.AND.IND3.EQ.1.OR.I.EQ.L.AND.IND4.EQ.1)GO TO 838 13473422 IF(I.NE.L.OR.IND4.NE.2)GO TO 1837 13473423 YI=CLOLAT 13473430 XI=CLOLON 13473431 1837 YK(Z)=YI 13473432 XK(Z)=XI 13473433 KA(Z)=NA(I) 13473434 KB(Z)=NB(I) 13473435 KC(Z)=NC(I) 13473436 KN(Z)=NUM(I) 13473437 Z=Z+1 13473438 838 IF (I-L)840,890,890 13473440 840 LINE=LINE+2 13473450 IF(LINE.GE.60)GO TO 800 13473460 845 J=I+1 13473470 ALAT=Y(I) 13473480 ALON=X(I) 13473490 BLAT=Y(J) 13473500 BLON=X(J) 13473510 CALL LTCS04 13473520 850 CALL LTCS03(AZAB,M1,M2,M3,M4,SM) 13473530 CALL LTCS03(AZBA,N1,N2,N3,N4,SN) 13473540 GO TO (870,860),MFO 13473550 860 S=S*3.28083333333D0 13473560 WRITE(6,861)M1,M2,M3,M4,SM,S,N1,N2,N3,N4,SN 13473570 861 FORMAT(1H ,19X,I3,I2,I1,I2,F4.2,F11.2,5H FT ,I3,I2,I1,I2,F4.2) 13473580 GO TO 872 13473590 870 WRITE(6,871)M1,M2,M3,M4,SM,S,N1,N2,N3,N4,SN 13473610 871 FORMAT(1H ,19X,I3,I2,I1,I2,F4.2,F12.3,4H M ,I3,I2,I1,I2,F4.2) 13473620 872 GO TO (880,880,960),IND2 13473630 880 CONTINUE 13473640 890 IF(NSS.EQ.0)GO TO 10 13473650 C 13473660 C COMPUTE SIDE OBSERVATIONS AND INVERSES 13473670 C 13473680 900 IR=1 13473700 LINE=60 13473710 DO 970 IS=1,NSS 13473730 READ(8'IR)A,SM,NA(101),NB(101),NC(101),NUM(101),NOS,NBS,JF,JT 13473745 IF(JF+JT.NE.0)GO TO 1000 13473746 IGO=3 13473747 C 13473750 C SEARCH FOR OCCUPIED STATION 13473760 905 DO910 K=1,500 13473770 IF(NOS.EQ.KN(K))GO TO 915 13473780 910 CONTINUE 13473790 NBS=NOS 13473800 GO TO 925 13473810 915 NOS=K 13473820 C SEARCH FOR BACKSIGHT STATION 13473830 DO920 K=1,500 13473840 IF(NBS.EQ.KN(K))GO TO 930 13473850 920 CONTINUE 13473860 925 WRITE (6,926)NBS 13473870 926 FORMAT (4H STA,I7,16H NOT COORDINATED) 13473880 GO TO 965 13473890 930 NBS=K 13473900 C 13473910 C COORDINATE SIDE STATION 13473920 C 13473930 ALAT=YK(NOS) 13473940 ALON=XK(NOS) 13473950 BLAT=YK(NBS) 13473960 BLON=XK(NBS) 13473970 IZ=1 13473980 CALL LTCS04 13473990 940 AZAB=AZAB+A 13474010 IF(AZAB.GT.6.28318530718D0)AZAB=AZAB-6.28318530718D0 13474020 945 IZ=0 13474040 S=SM 13474050 CALL LTCS04 13474060 950 IND2=3 13474070 IF(LINE.GT.58)GO TO 800 13474080 952 YI=YK(NOS) 13474090 XI=XK(NOS) 13474100 955 CALL LTCS03(YI,M1,M2,M3,M4,SM) 13474110 CALL LTCS03(XI,N1,N2,N3,N4,SN) 13474120 WRITE(6,836)KA(NOS),KB(NOS),KC(NOS),KN(NOS),M1,M2,M3,M4,SM,N1,N2, 13474130 1N3,N4,SN 13474135 GO TO 850 13474140 960 IF(NUM(101).EQ.0)GO TO 962 13474141 YK(Z)=BLAT 13474142 XK(Z)=BLON 13474143 KA(Z)=NA(101) 13474144 KB(Z)=NB(101) 13474145 KC(Z)=NC(101) 13474146 KN(Z)=NUM(101) 13474147 Z=Z+1 13474148 962 YI=BLAT 13474150 XI=BLON 13474151 CALL LTCS03(YI,M1,M2,M3,M4,SM) 13474152 CALL LTCS03(XI,N1,N2,N3,N4,SN) 13474153 WRITE(6,836)NA(101),NB(101),NC(101),NUM(101),M1,M2,M3,M4,SM,N1,N2,13474160 1N3,N4,SN 13474161 965 WRITE (6,966) 13474170 966 FORMAT (//) 13474180 LINE=LINE+5 13474190 970 CONTINUE 13474200 GO TO 10 13474210 C 13474220 C COMPUTE INVERSE 13474230 C 13474240 1000 IND2=3 13474250 IGO=4 13474251 IF(LINE.GT.58)GO TO 800 13474252 1005 DO1010 K=1,500 13474250 IF(JF.EQ.KN(K))GO TO 1015 13474260 1010 CONTINUE 13474270 NBS=JF 13474280 GO TO 925 13474290 1015 NBS=K 13474300 DO1020 K=1,500 13474310 IF(JT.EQ.KN(K))GO TO 1030 13474320 1020 CONTINUE 13474330 NBS=JT 13474340 GO TO 925 13474350 1030 ALAT=YK(NBS) 13474360 ALON=XK(NBS) 13474370 BLAT=YK(K) 13474380 BLON=XK(K) 13474390 IZ=1 13474400 CALL LTCS04 13474410 CALL LTCS03(AZAB,M1,M2,M3,M4,SM) 13474420 CALL LTCS03(AZBA,N1,N2,N3,N4,SN) 13474430 IF(MFO.EQ.2)S=S*3.28083333333D0 13474440 WRITE(6,1041)KA(NBS),KB(NBS),KC(NBS),JF 13474450 1041 FORMAT(1H ,3A3,I8) 13474460 IF(MFO.EQ.2)GO TO 1050 13474470 WRITE(6,871)M1,M2,M3,M4,SM,S,N1,N2,N3,N4,SN 13474480 GO TO 1060 13474490 1050 WRITE(6,861)M1,M2,M3,M4,SM,S,N1,N2,N3,N4,SN 13474500 1060 WRITE(6,1041)KA(K),KB(K),KC(K),JT 13474510 GO TO 965 13474520 1100 STOP 13474530 END 13474540 C LTCS03 - CONVERT RADIANS TO DEGS, MINS AND SECS FOR OUTPUT 10 SUBROUTINE LTCS03 (RADS,ID,IMT,IM,IST,SECS) 20 IMPLICIT REAL*8 (A-H,O-Z) 30 SECS=RADS*206264.806247D0 + .00005D0 40 ID=SECS/3600.D0 50 SECS=SECS-ID*3600 60 IM=SECS/60.D0 70 SECS=SECS-IM*60 80 IMT=IM/10 90 IM=IM-IMT*10 100 IST=SECS/10.D0 110 SECS=SECS-IST*10-.00005D0 120 IF(SECS.LE.0.D0)SECS=.00001D0 130 RETURN 140 END 150 C LTCS04 GEOGRAPHICAL COORDS FROM AZIMUTH AND DISTANCE AND INVERSE 100 SUBROUTINE LTCS04 110 IMPLICIT REAL*8 (A-H,O-Z) 120 COMMON ALAT,ALON,BLAT,BLON,AZAB,AZBA,S,IZ 130 EQUIVALENCE (SINTH,SINDL,DL,R2,RK),(R3,RK2),(D,F),(D3,P,SINP) 140 EQUIVALENCE (TH2,COSP,SIN1) 150 C 160 SINA=DSIN(ALAT) 170 COSA=DCOS(ALAT) 180 AN=6378206.4D0/DSQRT(1.D0-.676865799781D-2*SINA*SINA) 190 IF(IZ.NE.0)GO TO 200 200 C 210 C DIRECT COMPUTATION 220 C 230 66 COSAB=DCOS(AZAB) 240 R2=-.681478494638D-2*COSA*COSA*COSAB*COSAB 250 R3=.204443548391D-1*(1.D0-R2)*COSA*SINA*COSAB 260 D=S/AN 270 D3=D*D*D 280 THETA=D-R2*(1.D0+R2)*D3/6.D0-R3*(1.D0+R2+R2+R2)*D3*D/24.D0 290 TH2=THETA*THETA 300 F=1.D0-R2*TH2/2.D0-R3*THETA*TH2/6.D0 310 SINTH=DSIN(THETA) 320 SINP=SINA*DCOS(THETA)+COSA*COSAB*SINTH 330 COSP=DSQRT(1.D0-SINP*SINP) 340 SINDL=DSIN(AZAB)*SINTH/COSP 350 DL=DATAN(SINDL/(DSQRT(1.D0-SINDL*SINDL))) 360 BLON=ALON-DL 370 BLAT=DATAN((1.00681478494D0)*(SINP-.676865799781D-2*F*SINA)/COSP) 380 C 400 C INDIRECT COMPUTATION 410 C 420 200 SINB=DSIN(BLAT) 430 COSB=DCOS(BLAT) 440 BN=6378206.4D0/DSQRT(1.D0-.676865799781D-2*SINB*SINB) 450 205 SIN1=SINB 470 SIN2=SINA 480 COS1=COSB 490 COS2=COSA 500 AN1=BN 510 AN2=AN 520 D12=BLON-ALON 530 IX=1 540 GO TO 300 550 210 AZBA=AZ 560 IF(IZ.EQ.0)RETURN 565 215 SIN1=SINA 570 SIN2=SINB 580 COS1=COSA 590 COS2=COSB 600 AN1=AN 610 AN2=BN 620 D12=ALON-BLON 630 IX=2 640 GO TO 300 650 220 AZAB=AZ 660 230 BAN=BN/AN 680 COSAB=DCOS(AZAB)*COS1 690 P=DSQRT((BAN*C-COS1)**2+(BAN*AA)**2+(.993231342002D0* 700 1(BAN*SIN2-SIN1))**2) 710 F=1.D0+.00681478494638D0*COSAB*COSAB 720 RK=P*F 730 RK2=RK*RK 740 F1=.825517107417D-1*SIN1 750 H=.825517107417D-1*COSAB 751 H2=H*H 752 H21=1.D0+H2 753 F=F1*H/H21 754 H=(F1*F1-H2)/H21 755 S=(1.D0+RK2*.0416666666667D0-RK2*RK*F*.125D0)*P*AN1 760 IF(S.LT.100000.D0)RETURN 761 S=S+(.46875D-2+.375D-1*H+.25D0*F*F)*RK2*RK2*P*AN1 762 RETURN 770 C 780 C AZIMUTH FROM COORDINATES 790 C 800 300 AA=COS2*DSIN(D12) 810 C=COS2*DCOS(D12) 820 BB=.993231342002D0*COS1*SIN2-C*SIN1+AN1/AN2*.00676865799781D0* 830 1COS1*SIN1 840 AZ=DATAN(AA/BB) 850 IF(D12)310,360,320 860 310 IF(AZ)330,360,340 870 320 IF(AZ)340,360,350 880 330 AZ=AZ+3.14159265359D0 890 340 AZ=AZ+3.14159265359D0 900 350 GO TO (210,220),IX 910 360 IF(SIN1-SIN2)350,350,340 920 END 930 SUBROUTINE COMBO1 (A,AT,IROWA,ICOLA,B,BT,IROWB,ICOLB,W,P,N,M,XCAP) C ********************************************************************* C * SUBROUTINE 'COMBO1' DOES A LEAST SQUARES (COMBINED * C * CASE) ADJUSTMENT USING ADDITIONAL CONSTRAINTS BETWEEN THE * C * UNKNOWN PARAMETERS TO OBTAIN UNIQUE VALUES FOR THE UNKNOWNS * C * PLUS THEIR VARIANCE-COVARIANCE MATRIX. COMBO1 CALCULATES THE * C * XCAP VECTOR. * C * INPUT: * C * A - IS THE DESIGN MATRIX FOR THE UNKNOWNS * C * AT - IS THE TRANSPOSE OF A * C * IROWA - IS THE ROW DIMENSION OF A * C * ICOLA - IS THE COLUMN DIMENSION OF A * C * B - IS THE DESIGN MATRIX FOR THE OBSERVABLES * C * BT - IS THE TRANSPOSE OF B * C * IROWB - IS THE ROW DIMENSION OF B * C * ICOLB - IS THE COLUMN DIMENSION OF B * C * W - IS THE VECTOR OF MISCLOSURE * C * P - IS THE WEIGHT MATRIX * C * N - IS AT * M INV * A ALL INVERSE * C * M - IS B * P INV * BT ALL INVERSE * C * XCAP - IS A VECTOR WHICH CONTAINS THE APPROXIMATIONS * C * FOR THE UNKNOWNS * C * * C * OUTPUT: * C * XCAP - THE SOLUTION VECTOR * C * N INV - THE INVERSE OF AT * M INV * A ALL INVERSE * C * * C * NOTE: THIS ROUTINE USES SYSTEM SUBROUTINES EXCEPT FOR CHOLD * C * AND MATMUL. * C * * C * BY: R. CALDWELL (DECEMBER, 1976) * C * * C ********************************************************************* C IMPLICIT REAL*8 (A-H,O-Z), INTEGER*4 (I-N ) REAL*8 M,N DIMENSION A(IROWA,ICOLA),AT(ICOLA,IROWA),B(IROWB,ICOLB),BT(ICOLB, @ IROWB),W(IROWA),M(IROWB,IROWB),N(ICOLA,ICOLA), @ XCAP(ICOLA),P(ICOLB,ICOLB) C C MULTIPLYING B * P INV BT C CALL TRNSD (BT,ICOLB,B,IROWB,IROWB,ICOLB) CALL TRNSD (AT,ICOLA,A,IROWA,IROWA,ICOLA) CALL MATMUL (B,IROWB,ICOLB,P,ICOLB,ICOLB,B,IROWB,ICOLB,IROWB, @ ICOLB,ICOLB,2) CALL MATMUL (B,IROWB,ICOLB,BT,ICOLB,IROWB,M,IROWB,IROWB,IROWB, @ ICOLB,IROWB,1) C C INVERTING M C CALL CHOLD (M,IROWB,IROWB,DETA,&100) C C MULTIPLYING AT * M INV * A C CALL MATMUL (AT,ICOLA,IROWA,M,IROWB,IROWB,AT,ICOLA,IROWA,ICOLA, @ IROWA,IROWA,2) CALL MATMUL (AT,ICOLA,IROWA,A,IROWA,ICOLA,N,ICOLA,ICOLA,ICOLA, @ IROWA,ICOLA,1) C C INVERTING N C CALL CHOLD (N,ICOLA,ICOLA,DETA,&100) C C MULTIPLYING AT * M INV * W C CALL MATMUL (AT,ICOLA,IROWA,W,IROWA,1,XCAP,ICOLA,1,ICOLA,IROWA,1, @ 1) C C DETERMINING THE SOLUTION VECTOR C CALL MATMUL (N,ICOLA,ICOLA,XCAP,ICOLA,1,XCAP,ICOLA,1,ICOLA,ICOLA, @ 1,3) C C RETURNING TO THE MAIN PROGRAM / SUBROUTINE C 100 RETURN END