C DRIVE 01 C ******************************************************************DRIVE 02 C * *DRIVE 03 C * D R I V E R P R O G R A M F O R L S A O M P *DRIVE 04 C * *DRIVE 05 C ******************************************************************DRIVE 06 C * *DRIVE 07 C * 'LSAOMP' LEAST SQUARES ADJUSTMENT ON THE MAPPING PLANE *DRIVE 08 C * *DRIVE 09 C * THE PROGRAM REQUIRES ONE SEQUENTIAL DATA FILE NSEQF TO STORE *DRIVE 10 C * THE COMPACTED DESIGN MATRIX , COLUMN NUMBERS, CONSTANT VECTOR*DRIVE 11 C * AND WEIGHTS, SO THAT THEY CAN BE RECALLED LATER FOR THE *DRIVE 12 C * RESIDUAL COMPUTATION. *DRIVE 13 C * *DRIVE 14 C ******************************************************************DRIVE 15 C DRIVE 16 IMPLICIT REAL*8(A-H,O-Z) DRIVE 17 DIMENSION Y(81),X(81),XVEC(154),WEIGHT(10),CONST(10),STOR(10,20) DIMENSION TA(86,154),BL(154),AL(10),UP(10,10),AA(20,20),NAME(81) DIMENSION KK(4),ELEV( 81),NDIR(10),JCOL(20),ICOL(10),B(154) DRIVE 20 DIMENSION W(20,10),DEG(10),AMIN(10),SEC(10),AZ(10) DRIVE 21 DIMENSION STORAZ(10,4),STORDR(10,4),STORDS(10,4),MAZ(10,4) DRIVE 22 DIMENSION MDIR(10,4),MDIS(10,4),WTAZ(10),WTDR(10),WTDS(10) DRIVE 23 DIMENSION WVAZ(10),WVDR(10),WVDS(10),DIST(10),HTARGT(10),NDIS(10) DRIVE 24 C DRIVE 25 C THE ABOVE DIMENSION STATEMENTS MUST HAVE THE FOLLOWING MINIMUM DRIVE 26 C DIMENSIONS DRIVE 27 C DRIVE 28 C Y(NK1),X(NK1),XVEC(NK2),WEIGHT(NK3),CONST(NK3),STOR(NK3,NK4) DRIVE 29 C TA(NK2,NK5),BL(NK2),AL(NK3),UP(NK3,NK3),AA(NK4,NK4),NAME(NK1) DRIVE 30 C KK(4),ELEV(NK1),NDIR(NK3),JCOL(NK4),ICOL(NK3),B(NK2) DRIVE 31 C W(NK4,NK3),DEG(NK3),AMIN(NK3),SEC(NK3),AZ(NK3) DRIVE 32 C STORAZ(NK3,4),STORDR(NK3,4),STORDS(NK3,4),MAZ(NK3,4) DRIVE 33 C MDIR(NK3,4),MDIS(NK3,4),WTAZ(NK3),WTDR(NK3),WTDS(NK3) DRIVE 34 C WVAZ(NK3),WVDR(NK3),WVDS(NK3),DIST(NK3),HTARGT(NK3),NDIS(NK3) DRIVE 35 C DRIVE 36 C WHERE DRIVE 37 C NK1 = NUMBER OF STATIONS IN THE NETWORK DRIVE 38 C NK2 = NUMBER OF UNKNOWN COORDINATES DRIVE 39 C NK3 = MAX. NO. OF DIRECTIONS OR AZIMUTHS OR DISTANCES MEASURED DRIVE 40 C FROM ANY ONE STATION PLUS ONE DRIVE 41 C NK4 = 2*NK3 DRIVE 42 C NW = BANDWIDTH DRIVE 43 C NB = BORDERWIDTH DRIVE 44 C NK5 = NW + NB DRIVE 45 C DRIVE 46 NK1=81 NK2=154 NK3 = 10 DRIVE 49 NK4 = 20 DRIVE 50 NW=30 NB = 56 NF = 0 C DRIVE 54 C THE ABOVE SEVEN NUMBERS MUST BE COMPUTED FOR EACH PARTICULAR DRIVE 55 C ADJUSTMENT AND ENTERED AS THESE ARE DRIVE 56 C DRIVE 57 C NF = THE NUMBER OF FIXED STATIONS THAT ARE NOT OCCUPIED FOR DRIVE 58 C DIRECTIONS OR AZIMUTHS DRIVE 59 C DRIVE 60 NK5 = NW + NB DRIVE 61 N = NK2 - NB DRIVE 62 NBB = (NB**2 + NB)/2 DRIVE 63 NQQ = NK2 DRIVE 64 C DRIVE 65 C ENTER THE FOLLOWING DRIVE 66 C NSTDIR = NUMBER OF SETS OF DIRECTIONS TO BE ENTERED DRIVE 67 C NSTAZ = NUMBER OF SETS OF AZIMUTHS TO BE ENTERED DRIVE 68 C NSTDI = NUMBER OF SETS OF DISTANCES TO BE ENTERED DRIVE 69 C DRIVE 70 NSTDIR = 36 NSTAZ = 0 NSTDI = 85 NOSETS = NSTDIR + NSTAZ + NSTDI DRIVE 74 C DRIVE 75 C ENTER THE DEGREES OF FREEDOM 'DF' AND C FACTOR 'CFACT' DRIVE 76 C DRIVE 77 NDF = 200 DRIVE 79 CFACT = 2.474D0 DRIVE 80 C DRIVE 81 C ENTER THE FILE NUMBER OF NSEQF DRIVE 82 C DRIVE 83 NSEQF = 25 DRIVE 84 REWIND NSEQF DRIVE 85 C DRIVE 86 C PROGRAM OPTION DRIVE 87 C IF IPREAN = 0 THE SOLUTION IS COMPLETE DRIVE 88 C IF IPREAN = 1 THE X-VECTOR AND RESIDUALS ARE NOT COMPUTED DRIVE 89 C DRIVE 90 IPREAN=1 C DRIVE 92 C CALL LSAOMP DRIVE 93 C DRIVE 94 CALL LSAOMP(Y,X,XVEC,WEIGHT,CONST,STOR,TA,BL,AL,UP,AA,NAME,KK,ELEVDRIVE 95 1,NDIR,JCOL,ICOL,B,W,DEG,AMIN,SEC,AZ,NK1,NK2,NK3,NK4,NK5,NW,NB,NBB,DRIVE 96 2NQQ,NF,NSTDIR,NSTAZ,NSTDI,N,STORAZ,STORDR,STORDS,MAZ,MDIR,MDIS,WTADRIVE 97 3Z,WTDR,WTDS,WVAZ,WVDR,WVDS,DIST,HTARGT,NDIS,NSEQF,CFACT,NDF,NOSETSDRIVE 98 4,IPREAN) DRIVE 99 STOP DRIVE100 END DRIVE101 SUBROUTINE LSAOMP(Y,X,XVEC,WEIGHT,CONST,STOR,TA,BL,AL,UP,AA,NAME,KLSAOMP01 1K,ELEV,NDIR,JCOL,ICOL,B,W,DEG,AMIN,SEC,AZ,NK1,NK2,NK3,NK4,NK5,NW,NLSAOMP02 2B,NBB,NQQ,NF,NSTDIR,NSTAZ,NSTDI,N,STORAZ,STORDR,STORDS,MAZ,MDIR,MDLSAOMP03 3IS,WTAZ,WTDR,WTDS,WVAZ,WVDR,WVDS,DIST,HTARGT,NDIS,NSEQF,CFACT,NDF,LSAOMP04 4NOSETS,IPREAN) LSAOMP05 C LSAOMP06 C ******************************************************************LSAOMP06 C * *LSAOMP07 C * S U B R O U T I N E L S A O M P *LSAOMP08 C * *LSAOMP09 C * THIS SUBROUTINE READS IN THE FIXED AND APPROXIMATE *LSAOMP10 C * COORDINATES OF THE STATIONS IN THE NETWORK, PRINTS OUT *LSAOMP11 C * HEADINGS AND CALLS THE SUBROUTINES FOR THE SOLUTION. *LSAOMP12 C * *LSAOMP13 C ******************************************************************LSAOMP14 C LSAOMP15 IMPLICIT REAL*8(A-H,O-Z) LSAOMP16 DIMENSION Y(NK1),X(NK1),XVEC(NK2),WEIGHT(NK3),CONST(NK3) LSAOMP17 DIMENSION STOR(NK3,NK4),TA(NK5,NK2),BL(NK2),AL(NK3),UP(NK3,NK3) LSAOMP18 DIMENSION AA(NK4,NK4),NAME(NK1),KK(4),ELEV(NK1),NDIR(NK3) LSAOMP19 DIMENSION JCOL(NK4),ICOL(NK3),B(NK2),W(NK4,NK3),DEG(NK3) LSAOMP20 DIMENSION AMIN(NK3),SEC(NK3),AZ(NK3),STORAZ(NK3,4),STORDR(NK3,4) LSAOMP21 DIMENSION STORDS(NK3,4),MAZ(NK3,4),MDIR(NK3,4),MDIS(NK3,4) LSAOMP22 DIMENSION WTAZ(NK3),WTDS(NK3),WVAZ(NK3),WVDR(NK3),WVDS(NK3) LSAOMP23 DIMENSION DIST(NK3),HTARGT(NK3),NDIS(NK3),WTDR(NK3) LSAOMP24 NSPACE = N - NW + 1 LSAOMP25 DX = 0.0D0 LSAOMP26 NMZ = 0 LSAOMP27 NDEV = 0 LSAOMP28 C LSAOMP29 C READ IN COORDINATES OF STATIONS LSAOMP30 C LSAOMP31 PRINT 1 LSAOMP32 1 FORMAT('1',/////,39X,'COORDINATES OF FIXED STATIONS') LSAOMP33 PRINT 2 LSAOMP34 2 FORMAT(' ',///,19X,'STN. NO. SEQ. NO.',13X,'NORTHING',17X,'EASTLSAOMP35 #ING',//) LSAOMP36 DO 3 I = 1,NK1 LSAOMP37 3 READ 4,NAME(I),X(I),Y(I),ELEV(I) LSAOMP38 4 FORMAT(5X,I10,F15.2,F15.2,F15.2) LSAOMP39 MMN = NK1 - NK2/2 LSAOMP40 DO 5 I = 1,MMN LSAOMP41 IF((I/25*25).EQ.I) PRINT 2 LSAOMP42 IF((I/25*25).EQ.I) PRINT 1 LSAOMP4A 5 PRINT 6,NAME(I),I,Y(I),X(I) LSAOMP43 6 FORMAT(20X,I8,9X,I3,6X,F15.2,10X,F15.2,/) LSAOMP44 PRINT 7 LSAOMP45 7 FORMAT('1', /////,30X,'APPROXIMATE COORDINATES OF STATIONS TO BE ALSAOMP46 #DJUSTED') LSAOMP47 PRINT 2 LSAOMP48 MNNN = MMN + 1 LSAOMP49 DO 8 I = MNNN,NK1 LSAOMP50 IF(((I-MNNN+1)/25*25).EQ.(I-MNNN+1)) PRINT 7 LSAOMP5B IF(((I-MNNN+1)/25*25).EQ.(I-MNNN+1)) PRINT 2 LSAOMP51 8 PRINT 6,NAME(I),I,Y(I),X(I) LSAOMP52 PRINT 9 LSAOMP53 9 FORMAT('1',/////,23X,'DIRECTION OBSERVATIONS AND THEIR STANDARD DELSAOMP54 #VIATIONS',///) LSAOMP55 PRINT 10 LSAOMP56 10 FORMAT(18X,'FROM',6X,'TO',7X,'ST. DEV. SECS.',5X,'DIRECTION OBSERVLSAOMP57 1ATION',7X,'MISSCLOSURE VECTOR',//,58X,'DEG',3X,'MIN',4X,'SEC',16X,LSAOMP58 2'SECS',////) LSAOMP59 C LSAOMP60 C CLEAR NORMAL EQUATION MATRIX AND CONSTANT VECTOR LSAOMP61 C LSAOMP62 DO 11 I = 1,NK2 LSAOMP63 DO 11 J = 1,NK5 LSAOMP64 11 TA(J,I) = 0.0D0 LSAOMP65 DO 12 I = 1,NK2 LSAOMP66 12 BL(I) = 0.0D0 LSAOMP67 IAZ = 1 LSAOMP68 INT = 0 LSAOMP69 JNT = 1 LSAOMP70 IF(NSTAZ.EQ.0) GO TO 14 LSAOMP71 DO 13 INST = 1,NSTAZ LSAOMP72 CALL SETUP(Y,X,STOR,WEIGHT,DEG,AMIN,SEC,CONST,NOBS,NODIR,II,ICOL,NLSAOMP73 1DIR,NAME,KK,NF,IJK,IKJ,NK1,NK2,NK3,NK4,INT,JNT,HINST,DIST,HTARGT) LSAOMP74 CALL DIRAZ(Y,X,STOR,WEIGHT,DEG,AMIN,SEC,CONST,NOBS,NODIR,II,ICOL,NLSAOMP75 #DIR,NAME,KK,NF,IJK,IKJ,NK1,NK2,NK3,NK4,INT,AZ,IAZ,STORAZ,MAZ,NSEQFLSAOMP76 @) LSAOMP77 CALL DIRNOR(STOR,WEIGHT,CONST,AL,AA,NODIR,II,UP,W,NK3,NK4,INST,IAZLSAOMP78 #) LSAOMP79 CALL COLUMN(II,ICOL,JCOL,NK3,NK4) LSAOMP80 CALL ADDNOR(TA,AA,II,JCOL,NK2,NK4,NK5,NSPACE,NW) LSAOMP81 CALL ADDVEC(BL,AL,II,JCOL,NK2,NK3,NK4) LSAOMP82 13 CONTINUE LSAOMP83 14 IF(NSTDIR.EQ.0) GO TO 16 LSAOMP84 INT = 0 LSAOMP85 JNT = 1 LSAOMP86 IAZ = 0 LSAOMP87 DO 15 INST = 1,NSTDIR LSAOMP88 CALL SETUP(Y,X,STOR,WEIGHT,DEG,AMIN,SEC,CONST,NOBS,NODIR,II,ICOL,NLSAOMP89 1DIR,NAME,KK,NF,IJK,IKJ,NK1,NK2,NK3,NK4,INT,JNT,HINST,DIST,HTARGT) LSAOMP90 CALL DIRAZ(Y,X,STOR,WEIGHT,DEG,AMIN,SEC,CONST,NOBS,NODIR,II,ICOL,NLSAOMP91 #DIR,NAME,KK,NF,IJK,IKJ,NK1,NK2,NK3,NK4,INT,AZ,IAZ,STORAZ,MAZ,NSEQFLSAOMP92 @) LSAOMP93 CALL DIRNOR(STOR,WEIGHT,CONST,AL,AA,NODIR,II,UP,W,NK3,NK4,INST,IAZLSAOMP94 #) LSAOMP95 CALL COLUMN(II,ICOL,JCOL,NK3,NK4) LSAOMP96 CALL ADDNOR(TA,AA,II,JCOL,NK2,NK4,NK5,NSPACE,NW) LSAOMP97 CALL ADDVEC(BL,AL,II,JCOL,NK2,NK3,NK4) LSAOMP98 15 CONTINUE LSAOMP99 16 IF(NSTDI.EQ.0) GO TO 20 LSAOM100 INT = 1 LSAOM101 JNT = 2 LSAOM102 PRINT 17 LSAOM103 17 FORMAT('1',/////,23X,'DISTANCE OBSERVATIONS AND THEIR STANDARD DEILSAOM104 #ATIONS',///) LSAOM105 PRINT 18 LSAOM106 18 FORMAT(18X,'FROM',6X,'TO',8X,'ST. DEV.',4X,'DISTANCE OBSERVATION'LSAOM107 #,14X,'MISSCLOSURE VECTOR',//,40X,'FEET',15X,'FEET',28X,'FEET',///)LSAOM108 DO 19 INST = 1,NSTDI LSAOM109 CALL SETUP(Y,X,STOR,WEIGHT,DEG,AMIN,SEC,CONST,NOBS,NODIST,II,ICOL,LSAOM110 1NDIS,NAME,KK,NF,IJK,IKJ,NK1,NK2,NK3,NK4,INT,JNT,HINST,DIST,HTARGT)LSAOM111 CALL DISTA(Y,X,STOR,NOBS,NODIST,II,ICOL,NDIS,KK,NF,IJK,IKJ,NK1,NK2LSAOM112 #,NK3,NK4,DIST,CONST,WEIGHT,HINST,HTARGT,ELEV,NSEQF) LSAOM113 CALL DISNOR(STOR,WEIGHT,CONST,AL,AA,NODIST,II,NK3,NK4) LSAOM114 CALL COLUMN(II,ICOL,JCOL,NK3,NK4) LSAOM115 CALL ADDNOR(TA,AA,II,JCOL,NK2,NK4,NK5,NSPACE,NW) LSAOM116 CALL ADDVEC(BL,AL,II,JCOL,NK2,NK3,NK4) LSAOM117 19 CONTINUE LSAOM118 20 CONTINUE LSAOM119 CALL ARROW(TA,B,XVEC,BL,NK2,NK5,NB,NW,N,NBB,NQQ,IPREAN) LSAOM120 IF(IPREAN.EQ.1) GO TO 26 LSAOM121 VARFAC = 0.0D0 LSAOM12* CALL RESIDL(XVEC,STORAZ,STORDR,STORDS,MAZ,MDIR,MDIS,WTAZ,WTDR,WTDSLSAOM122 #,WVAZ,WVDR,WVDS,VARFAC,NK2,NK3,NSEQF,NOSETS,NK1) LSAOM123 DF = DFLOAT(NDF) LSAOM124 VARFAC = VARFAC/DF LSAOM125 PRINT 27,VARFAC,NDF LSAOM126 27 FORMAT('1',//////////,10X,'V A R I A N C E F A C T O R =',F10.LSAOM127 #6,5X,'D E G R E E S O F F R E E D O M =',I6) LSAOM128 L = 0 LSAOM129 I = -1 LSAOM130 J = NK1 - NK2/2 LSAOM131 MI = NK2/2 LSAOM132 PRINT 22 LSAOM133 22 FORMAT('1',/////,16X,'SEQ. NO.',7X,'STATION NO.',6X,'DIFF. IN NORTLSAOM134 #HING',5X,'DIFF. IN EASTING',//,54X,'(DY)',18X,'(DX)',///) LSAOM135 DO 25 K = 1,MI LSAOM136 J = J + 1 LSAOM137 I = I + 2 LSAOM138 L = L + 2 LSAOM139 PRINT 23,J,NAME(J),XVEC(I),XVEC(L) LSAOM140 23 FORMAT(' ',17X,I3,11X,I5,13X,F9.3,10X,F11.3,/) LSAOM141 IF((K/25*25).EQ.K) PRINT 22 LSAOM142 25 CONTINUE LSAOM143 DO 21 I = 1,NK2 LSAOM144 DO 21 J = 1,NK5 LSAOM145 21 TA(J,I) = VARFAC*TA(J,I) LSAOM146 26 CONTINUE LSAOM147 CALL ABSOLU(TA,NK1,NK2,NK5,CFACT,NAME) LSAOM148 CALL RELATI(TA,NK1,NK2,NK5,CFACT,NSPACE,NW,NAME) LSAOM149 RETURN LSAOM150 END LSAOM151 SUBROUTINE READAZ(S,T1,T2,T3,T4,I1,I2,I3,I4,MAZ,WT,WV,E,V,I,NK3) READAZ01 C READAZ02 C ******************************************************************READAZ03 C * *READAZ04 C * S U B R O U T I N E R E A D A Z *READAZ05 C * *READAZ06 C * THIS SUBROUTINE TAKES THE DESIGN MATRIX COEFFICIENTS THAT *READAZ07 C * CAME FROM THE STROAGE DEVICE AND STORES THEM IN A MATRIX *READAZ08 C * SO THAT MATRIX MULTIPLICATION YIELD THE RESIDUALS. *READAZ09 C * *READAZ10 C ******************************************************************READAZ11 C READAZ12 IMPLICIT REAL*8(A-H,O-Z) READAZ13 DIMENSION S(NK3,4),MAZ(NK3,4),WT(NK3),WV(NK3) READAZ14 C READAZ15 S(I,1) = T1 READAZ16 S(I,2) = T2 READAZ17 S(I,3) = T3 READAZ18 S(I,4) = T4 READAZ19 MAZ(I,1) = I1 READAZ20 MAZ(I,2) = I2 READAZ21 MAZ(I,3) = I3 READAZ22 MAZ(I,4) = I4 READAZ23 WV(I) = V READAZ24 WT(I) = E READAZ25 RETURN READAZ26 END READAZ27 SUBROUTINE SETUP(Y,X,STOR,WEIGHT,DEG,AMIN,SEC,CONST,NOBS,NODIR,II,SETUP 01 1ICOL,NDIR,NAME,KK,NF,IJK,IKJ,NK1,NK2,NK3,NK4,INT,JNT,HINST,DIST,HTSETUP 02 2ARGT) SETUP 03 C SETUP 04 C ******************************************************************SETUP 05 C * *SETUP 06 C * S U B R O U T I N E S E T U P *SETUP 07 C * *SETUP 08 C * THIS SUBROUTINE READS IN THE DATA FOR EACH OBSERVING STATION *SETUP 09 C * ONE AT A TIME. IF THE INPUT PARAMETER JNT = 1,THE *SETUP 10 C * OBSERVATIONS ARE EITHER DIRECTIONS OR AZIMUTHS - THESE WILL *SETUP 11 C * BE SORTED IN SUBROUTINE DIRNOR. IF THE INPUT PARAMETER *SETUP 12 C * JNT = 2, THE OBSERVATIONS ARE DISTANCES. *SETUP 13 C * *SETUP 14 C ******************************************************************SETUP 15 C SETUP 16 IMPLICIT REAL*8(A-H,O-Z) SETUP AA DIMENSION Y(NK1),X(NK1),STOR(NK3,NK4),WEIGHT(NK3),CONST(NK3) SETUP 17 DIMENSION NDIR(NK3),KK(4),NAME(NK1),PTOR(4),ICOL(NK3) SETUP 18 DIMENSION DEG(NK3),AMIN(NK3),SEC(NK3) SETUP 19 DIMENSION DIST(NK3),HTARGT(NK3) SETUP 20 C SETUP 21 C READ IN DIRECTION OR DISTANCE OBSERVATIONS ONE STATION AT A TIME SETUP 22 C SETUP 23 IF(JNT.EQ.1) READ 1,NOBS,NODIR,INT SETUP 24 1 FORMAT(3I5) SETUP 25 IF(JNT.EQ.2) READ 2,NOBS,NODIR,HINST SETUP 26 2 FORMAT(2I5,F10.2) SETUP 27 IF(JNT.EQ.2) GO TO 5 SETUP 28 DO 3 I = 1,NODIR SETUP 29 READ 4,NDIR(I),WEIGHT(I),MDEG,MIN,SEC(I) SETUP 30 DEG(I) = DFLOAT(MDEG) SETUP 31 3 AMIN(I) = DFLOAT(MIN) SETUP 32 4 FORMAT(I5,F5.2,12X,I3,3X,I2,3X,F4.1) SETUP 33 GO TO 8 SETUP 34 5 DO 6 I = 1,NODIR SETUP 35 6 READ 7,NDIR(I),WEIGHT(I),DIST(I),HTARGT(I) SETUP 36 7 FORMAT(I5,F5.3,12X,F13.3,F10.2) SETUP 37 GO TO 10 SETUP 38 8 DO 9 I = 1,NODIR SETUP 39 IF(INT.EQ.NDIR(I)) INT = I SETUP 40 9 CONTINUE SETUP 41 10 CONTINUE SETUP 42 C SETUP 43 C SET UP MATRIX TO CONTAIN PARTIAL A MATRIX SETUP 44 C NUMBER OF COLUMNS IN PARTIAL A MATRIX = II SETUP 45 C NUMBER OF ROWS IN PARTIAL A MATRIX = NO OF DIRECTIONS SETUP 46 C IF THE OBSERVATION STATION IS FIXED IT WILL NOT BE PART OF THE SETUP 47 C PARTIAL A MATRIX SETUP 48 C IF ONE OF THE DIRECTION STATIONS IS FIXED IT WILL NOT BE SETUP 49 C PART OF THE PARTIAL A MATRIX EITHER. SETUP 50 C THE OBSERVATION STATION IS REPRESENTED IN THE STOR MATRIX BY SETUP 51 C COLUMNS 2*K+1 AND2*K+2, THE REMAINDER OF THE COLUMNS SETUP 52 C REPRESENT THE DIRECTION STATIONS SETUP 53 C SETUP 54 II = 0 SETUP 55 MMN = 2*NK1 - NK2 SETUP 56 NTEST = NOBS - MMN/2 SETUP 57 IF(NTEST.GT.0) II = 1 SETUP 58 DO 11 I = 1,NODIR SETUP 59 MTEST = NDIR(I) - MMN/2 SETUP 60 IF(MTEST.GT.0) II = II + 1 SETUP 61 11 CONTINUE SETUP 62 II = 2*II SETUP 63 IJK = 0 SETUP 64 K = 0 SETUP 65 IF(NTEST.LE.0) GO TO 13 SETUP 66 DO 12 I = 1,NODIR SETUP 67 MTEST = NDIR(I) - MMN/2 SETUP 68 IF(MTEST.LE.0) GO TO 12 SETUP 69 IF(MTEST.LT.NTEST) K = K + 1 SETUP 70 12 CONTINUE SETUP 71 IJK = 2*K + 1 SETUP 72 IKJ = IJK + 1 SETUP 73 GO TO 14 SETUP 74 13 CONTINUE SETUP 75 14 CONTINUE SETUP 76 C SETUP 77 C CLEAR THE PARTIAL A MATRIX SETUP 78 C SETUP 79 DO 15 I = 1,II SETUP 80 DO 15 J = 1,NODIR SETUP 81 15 STOR(J,I) = 0.0D0 SETUP 82 RETURN SETUP 83 END SETUP 84 SUBROUTINE DIRAZ(Y,X,STOR,WEIGHT,DEG,AMIN,SEC,CONST,NOBS,NODIR,II,DIRAZ 01 %ICOL,NDIR,NAME,KK,NF,IJK,IKJ,NK1,NK2,NK3,NK4,INT,AZ,IAZ,STORAZ,MAZDIRAZ 02 #,NSEQF) DIRAZ02A C DIRAZ 03 C ******************************************************************DIRAZ 04 C * *DIRAZ 05 C * S U B R O U T I N E D I R A Z *DIRAZ 06 C * *DIRAZ 07 C * THIS SUBROUTINE TAKES IN ALL STATION NUMBERS PERTAINING *DIRAZ 08 C * TO ALL THE DIRECTIONS MEASURED AT ONE STATION (OR AZIMUTHS) *DIRAZ 09 C * AND FORMS THE COEFFICIENTS OF THE DESIGN MATRIX. THE *DIRAZ 10 C * COEFFICIENTS ARE COMPACTED INTO A RECTANGULAR MATRIX (STOR), *DIRAZ 11 C * THE NUMBER OF ROWS EQUAL TO THE NUMBER OF DIRECTIONS *DIRAZ 12 C * OBSERVED AT THAT STATION (I), THE NUMBER OF COLUMNS ARE *DIRAZ 13 C * EQUAL TO 2(I+1). A COLUMN INDICATOR IS ALSO GENERATED FROM *DIRAZ 14 C * THE STATION NUMBERS, ONLY THE ODD COLUMN NUMBERS ARE *DIRAZ 15 C * INDICATED :(ICOL) *DIRAZ 16 C * *DIRAZ 17 C ******************************************************************DIRAZ 18 C DIRAZ 19 IMPLICIT REAL*8(A-H,O-Z) DIRAZ 20 DIMENSION AZ(NK3) DIRAZ 21 DIMENSION Y(NK1),X(NK2),STOR(NK3,NK4),WEIGHT(NK3),CONST(NK3) DIRAZ 22 DIMENSION NDIR(NK3),KK(4),NAME(NK1),PTOR(4),ICOL(NK3) DIRAZ 23 DIMENSION DEG(NK3),AMIN(NK3),SEC(NK3),STORAZ(NK3,4),MAZ(NK3,4) DIRAZ 24 MMN = 2*NK1 - NK2 DIRAZ 25 RAD = 206264.806D0 DIRAZ 26 IZ = 0 DIRAZ 27 M1 = NOBS DIRAZ 28 DO 10 J = 1,NODIR DIRAZ 29 M2 = NDIR(J) DIRAZ 30 C DIRAZ 31 C M1= STATION NUMBER OF OBSERVING STATION DIRAZ 32 C M2= STATION NUMBER OF TARGET STATION DIRAZ 33 C DIRAZ 34 C FORM COORDINATE DIFFERENCES DIRAZ 35 C DIRAZ 36 A = Y(M1) DIRAZ 37 B = Y(M2) DIRAZ 38 C = X(M1) DIRAZ 39 D = X(M2) DIRAZ 40 AB = B - A DIRAZ 41 BA = DABS(AB) DIRAZ 42 IF(BA.LT.0.000001) AB = 0.000001D0 DIRAZ 43 CD = D - C DIRAZ 44 C DIRAZ 45 C FORM CONSTANT VECTOR 'L' DIRAZ 46 C DIRAZ 47 CONST(J) = DEG(J) + AMIN(J)/60.0D0 + SEC(J)/3600.0D0 DIRAZ 48 C DIRAZ 49 C FORM APPROXIMATE CONSTANT VECTOR 'LA' DIRAZ 50 C DIRAZ 51 ABCD = DABS(CD/AB) DIRAZ 52 IF(AB.LT.0.0D0) GO TO 2 DIRAZ 53 IF(CD.LT.0.0D0) GO TO 1 DIRAZ 54 AZ(J) = DATAN(ABCD)*57.2957795D0 DIRAZ 55 GO TO 4 DIRAZ 56 1 AZ(J) = 360.0D0 - DATAN(ABCD)*57.2957795D0 DIRAZ 57 GO TO 4 DIRAZ 58 2 IF(CD.LT.0.0D0) GO TO 3 DIRAZ 59 AZ(J) = 180.0D0 - DATAN(ABCD)*57.2957795D0 DIRAZ 60 GO TO 4 DIRAZ 61 3 AZ(J) = 180.0D0 + DATAN(ABCD)*57.2957795D0 DIRAZ 62 4 CONTINUE DIRAZ 63 C DIRMAP01 C DIRMAP02 C ******************************************************************DIRMAP03 C * *DIRMAP04 C * S U B R O U T I N E D I R M A P *DIRMAP05 C * *DIRMAP06 C * THIS SUBROUTINE REDUCES THE MEASURED DIRECTIONS TO THE *DIRMAP07 C * MAPPING PLANE. THE SUBROUTINE HAS TO BE WRITTEN BY THE *DIRMAP08 C * USER FOR HIS PARTICULAR MAP PROJECTION. *DIRMAP09 C * THE INPUT IS AS FOLLOWS *DIRMAP10 C * CONST(J) = DIRECTION NUMBER(J) IN DECIMAL DEGREES *DIRMAP11 C * AZ = THE APPROXIMATE AZIMUTH OF THE DIRECTION J *DIRMAP12 C * Y(M1) = NORTHING OF INSTRUMENT STATION *DIRMAP13 C * Y(M2) = NORTHING OF TARGET STATION FOR DIRECTION (J) *DIRMAP14 C * X(M1) = EASTING OF INSTRUMENT STATION *DIRMAP15 C * X(M1) = EASTING OF INSTRUMENT STATION *DIRMAP16 C * X(M2) = EASTING OF TARGET STATION FOR DIRECTION (J) *DIRMAP17 C * M1 = STATION NUMBER OF INSTRUMENT STATION *DIRMAP18 C * M2 = STATION NUMBER OF TARGET STATION FOR DIRECTION(J) *DIRMAP19 C * *DIRMAP20 C * IF IAZ = 0 (A DIRECTION WAS MEASURED ) *DIRMAP21 C * IF IAZ = 1 (AN AZIMUTH WAS MEASURED) *DIRMAP22 C * *DIRMAP23 C ******************************************************************DIRMAP24 C DIRMAP25 C DIRMAP28 C DIRAZ 64 C FORM COEFFICIENTS OF DESIGN MATRIX, INDICATOR VECTOR AND STOR THEMDIRAZ 65 C DIRAZ 66 ABSQ = AB*AB DIRAZ 67 CDSQ = CD*CD DIRAZ 68 DSQ = ABSQ + CDSQ DIRAZ 69 M = 2*M1 - 1 DIRAZ 70 I = 0 DIRAZ 71 KKK = 0 DIRAZ 72 KKKK = NK2 + M1 - NF DIRAZ 73 IF(M.LE.MMN) GO TO 5 DIRAZ 74 I = I + 1 DIRAZ 75 STOR(J,IJK) = RAD/DSQ*CD DIRAZ 76 PTOR(I) = STOR(J,IJK) DIRAZ 77 KK(I) = M - MMN DIRAZ 78 IF(J.EQ.1) IZ = IZ + 1 DIRAZ 79 IF(J.EQ.1) ICOL(1) = KK(I) DIRAZ 80 M = M + 1 DIRAZ 81 I = I + 1 DIRAZ 82 STOR(J,IKJ) = -RAD/DSQ*AB DIRAZ 83 PTOR(I) = STOR(J,IKJ) DIRAZ 84 KK(I) = M - MMN DIRAZ 85 5 M = M2*2 - 1 DIRAZ 86 IF(M.LE.MMN) GO TO 6 DIRAZ 87 IZ = IZ + 1 DIRAZ 88 I = I + 1 DIRAZ 89 KKL = 2*(J+II/2-NODIR) - 1 DIRAZ 90 IF((IJK.GT.1).AND.(M2.LT.M1)) KKL = KKL - 2 DIRAZ 91 STOR(J,KKL) = -RAD/DSQ*CD DIRAZ 92 PTOR(I) = STOR(J,KKL) DIRAZ 93 KK(I) = M -MMN DIRAZ 94 ICOL(IZ) = KK(I) DIRAZ 95 M = M + 1 DIRAZ 96 I = I + 1 DIRAZ 97 LKL = KKL + 1 DIRAZ 98 STOR(J,LKL) = RAD/DSQ*AB DIRAZ 99 PTOR(I) = STOR(J,LKL) DIRAZ100 KK(I) = M - MMN DIRAZ101 6 CONTINUE DIRAZ102 IF(I.EQ.2) GO TO 7 DIRAZ103 IF(I.EQ.4) GO TO 8 DIRAZ104 STORAZ(J,1) = 0.0D0 DIRAZ105 STORAZ(J,2) = 0.0D0 DIRAZ106 STORAZ(J,3) = 0.0D0 DIRAZ107 STORAZ(J,4) = 0.0D0 DIRAZ108 MAZ(J,1) = 1 DIRAZ109 MAZ(J,2) = 1 DIRAZ110 MAZ(J,3) = 1 DIRAZ111 MAZ(J,4) = 1 DIRAZ112 GO TO 9 DIRAZ113 7 STORAZ(J,1) = 0.0D0 DIRAZ114 STORAZ(J,2) = 0.0D0 DIRAZ115 STORAZ(J,3) = PTOR(1) DIRAZ116 STORAZ(J,4) = PTOR(2) DIRAZ117 MAZ(J,1) = 1 DIRAZ118 MAZ(J,2) = 1 DIRAZ119 MAZ(J,3) = KK(1) DIRAZ120 MAZ(J,4) = KK(2) DIRAZ121 GO TO 9 DIRAZ122 8 STORAZ(J,1) = PTOR(1) DIRAZ123 STORAZ(J,2) = PTOR(2) DIRAZ124 STORAZ(J,3) = PTOR(3) DIRAZ125 STORAZ(J,4) = PTOR(4) DIRAZ126 MAZ(J,1) = KK(1) DIRAZ127 MAZ(J,2) = KK(2) DIRAZ128 MAZ(J,3) = KK(3) DIRAZ129 MAZ(J,4) = KK(4) DIRAZ130 9 CONTINUE DIRAZ131 10 CONTINUE DIRAZ132 C DIRAZ133 C FORM CONSTANT VECTOR 'W' EQUALS APPROXIMATE MINUS OBSERVED DIRAZ134 C DIRAZ135 DO 11 J = 1,NODIR DIRAZ136 BZ = AZ(J) - AZ(INT) DIRAZ137 IF(IAZ.EQ.1) BZ = AZ(J) DIRAZ138 IF(BZ.LT.0.0D0) BZ = BZ + 360.0D0 DIRAZ139 IF((BZ - CONST(J)).GT.359) BZ = BZ - 360.0D0 DIRAZ140 CONST(J) = (BZ - CONST(J))*3600.0D0 DIRAZ141 11 PRINT 12,NOBS,NDIR(J),WEIGHT(J),DEG(J),AMIN(J),SEC(J),CONST(J) DIRAZ142 12 FORMAT(17X,I4,5X,I4,10X,F5.1,10X,F6.0,F6.0,F8.1,10X,F10.1,/) DIRAZ143 DO 13 J = 1,NODIR DIRAZ144 IF(IAZ.EQ.1) NCODE = 50 DIRAZ145 IF(IAZ.EQ.0) NCODE = 30 DIRAZ146 IF(J.EQ.1) WRITE(NSEQF) NODIR,NOBS,NDIR(J),MAZ(J,1),MAZ(J,2),MAZ(JDIRAZ147 1,3),MAZ(J,4),STORAZ(J,1),STORAZ(J,2),STORAZ(J,3),STORAZ(J,4),CONSTDIRAZ148 2(J),WEIGHT(J),NCODE DIRAZ149 IF(J.GT.1) WRITE(NSEQF) NOBS,NDIR(J),MAZ(J,1),MAZ(J,2),MAZ(J,3),MADIRAZ150 1Z(J,4),STORAZ(J,1),STORAZ(J,2),STORAZ(J,3),STORAZ(J,4),CONST(J),WEDIRAZ151 3IGHT(J),NCODE DIRAZ152 13 CONTINUE DIRAZ153 RETURN DIRAZ154 END DIRAZ155 SUBROUTINE DIRNOR(STOR,WEIGHT,CONST,AL,AA,NODIR,II,UP,W,NK3,NK4,INDIRNOR01 1ST,IAZ) DIRNOR02 C DIRNOR03 C ******************************************************************DIRNOR04 C * *DIRNOR05 C * S U B R O U T I N E D I R N O R *DIRNOR06 C * *DIRNOR07 C * THIS SUBROUTINE ACCEPTS THE COMPACTED DESIGN MATRIX FROM *DIRNOR08 C * THE PARTICULAR OBSERVATION STATION RESULTING FROM EITHER *DIRNOR09 C * DIRECTIONS OR AZIMUTHS. IF THE OBSERVATION EQUATIONS ARE *DIRNOR10 C * FROM DIRECTIONS, THE ORIENTATION UNKNOWNS ARE ELIMINATED. *DIRNOR11 C * THE PARTIAL NORMAL EQUATION MATRIX AND CONSTANT VECTOR *DIRNOR12 C * ARE FORMED. *DIRNOR13 C * *DIRNOR14 C ******************************************************************DIRNOR15 C DIRNOR16 IMPLICIT REAL*8(A-H,O-Z) DIRNOR17 DIMENSION STOR(NK3,NK4),WEIGHT(NK3),CONST(NK3),AL(NK3) DIRNOR18 DIMENSION AA(NK4,NK4),UP(NK3,NK3),W(NK4,NK3) DIRNOR19 C DIRNOR20 C WEIGHT EQUALS ONE DIVIDED BY VARIANCE DIRNOR21 C DIRNOR22 DO 1 I = 1,NODIR DIRNOR23 1 WEIGHT(I) = 1.0D0/WEIGHT(I)/WEIGHT(I) DIRNOR24 C DIRNOR25 C MULTIPLY ATUP BY CONSTANT VECTOR 'W' DIRNOR26 C DIRNOR27 SUM = 0.0D0 DIRNOR28 DO 2 I = 1,NODIR DIRNOR29 DO 2 J = 1,NODIR DIRNOR30 2 UP(I,J) = 0.0D0 DIRNOR31 IF(IAZ.EQ.1) GO TO 5 DIRNOR32 DO 3 I = 1,NODIR DIRNOR33 3 SUM = SUM + WEIGHT(I) DIRNOR34 DO 4 I = 1,NODIR DIRNOR35 DO 4 J = I,NODIR DIRNOR36 UP(I,J) = -WEIGHT(I)*WEIGHT(J)/SUM DIRNOR37 4 UP(J,I) = UP(I,J) DIRNOR38 5 DO 6 I = 1,NODIR DIRNOR39 6 UP(I,I) = UP(I,I) + WEIGHT(I) DIRNOR40 C DIRNOR41 C MULTIPLY A TRANSPOSED BY MATRIX UP DIRNOR42 C DIRNOR43 SUM = 0.0D0 DIRNOR44 DO 8 I = 1,II DIRNOR45 DO 8 J = 1,NODIR DIRNOR46 DO 7 K = 1,NODIR DIRNOR47 7 SUM = SUM + STOR(K,I)*UP(K,J) DIRNOR48 W(I,J) = SUM DIRNOR49 8 SUM = 0.0D0 DIRNOR50 C DIRNOR51 C MULTIPLY ATUP BY A MATRIX(PARTIAL) DIRNOR52 C DIRNOR53 SUM = 0.0D0 DIRNOR54 DO 10 I = 1,II DIRNOR55 DO 10 J = I,II DIRNOR56 DO 9 K = 1,NODIR DIRNOR57 9 SUM = SUM + W(I,K)*STOR(K,J) DIRNOR58 AA(I,J) = SUM DIRNOR59 10 SUM = 0.0D0 DIRNOR60 C DIRNOR61 C MULTIPLY ATUP BY CONSTANT VECTOR 'W' DIRNOR62 C DIRNOR63 SUM = 0.0D0 DIRNOR64 DO 12 I = 1,II DIRNOR65 DO 11 K = 1,NODIR DIRNOR66 11 SUM = SUM + W(I,K)*CONST(K) DIRNOR67 AL(I) = SUM DIRNOR68 12 SUM = 0.0D0 DIRNOR69 RETURN DIRNOR70 END DIRNOR71 SUBROUTINE DISTA(Y,X,STOR,NOBS,NODIST,II,ICOL,NDIS,KK,NF,IJK,IKJ,NDISTA 01 %K1,NK2,NK3,NK4,DIST,CONST,WEIGHT,HINST,HTARGT,ELEV,NSEQF) DISTA 02 C DISTA 03 C ******************************************************************DISTA 04 C * *DISTA 05 C * S U B R O U T I N E D I S T A *DISTA 06 C * *DISTA 07 C * THIS SUBROUTINE FORMS THE DESIGN MATRIX ASSOCIATED WITH *DISTA 08 C * THE DISTANCES FROM ONE STATION. IT ALSO FORMS THE CONSTANT *DISTA 09 C * VECTOR AND THE INDICIE VECTOR (WHICH INDICATES WHICH COLUMN *DISTA 10 C * NUMBERS THE COEFFICIENTS OF THE PARTIAL DESIGN MATRIX BELONG *DISTA 11 C * TO. *DISTA 12 C * *DISTA 13 C ******************************************************************DISTA 14 C DISTA 15 IMPLICIT REAL*8(A-H,O-Z) DISTA 16 DIMENSION Y(NK1),X(NK1),STOR(NK3,NK4),NDIS(NK3),KK(4),PTOR(4) DISTA 17 DIMENSION ICOL(NK3),CONST(NK3),HTARGT(NK3),WEIGHT(NK3) DISTA 18 DIMENSION DIST(NK3),ELEV(NK1) DISTA 19 MMN = 2*NK1 - NK2 DISTA 20 IZ = 0 DISTA 21 M1 = NOBS DISTA 22 DO 15 J = 1,NODIST DISTA 23 M2 = NDIS(J) DISTA 24 C DISTA 25 C DISTA 26 C DISMAP01 C DISMAP02 C ******************************************************************DISMAP03 C * *DISMAP04 C * S U B R O U T I N E D I S M A P *DISMAP05 C * *DISMAP06 C * THIS SUBROUTINE ACCEPTS THE HORIZONTAL MEASURED DISTANCE *DISMAP07 C * AND REDUCES IT TO THE MAPPING PLANE. YOU HAVE TO WRITE *DISMAP08 C * YOUR OWN EQUATIONS IN THIS SECTION ASSOCIATED WITH THE *DISMAP09 C * MAP PROJECTION YOU ARE USING. THE EQUATIONS HEREIN ARE *DISMAP10 C * FOR THE NEW BRUNSWICK STEREOGRAPHIC PROJECTION. *DISMAP11 C * *DISMAP12 C ******************************************************************DISMAP13 C DISMAP14 C DISMAP15 C DISMAP16 C DISMAP17 C ELEV(M1) = ELEVATION OF OBSERVING STATION DISMAP18 C ELEV(M2) = ELEVATION OF TARGET STATION DISMAP19 C HINST = HEIGHT OF INSTRUMENT ABOVE OBSERVING STATION DISMAP20 C HTARGT(J) = HEIGHT OF TARGET ABOVE TARGET STATION DISMAP21 C DIST(J) = MEASURED HORIZONTAL DISTANCE (INPUT) DISMAP22 C DIST(J) = CORRECTED HORIZONTAL DISTANCE (OUTPUT) DISMAP23 C DISMAP24 HMEAN = (ELEV(M1) + ELEV(M2) + HINST + HTARGT(J))/2.0D0 DISMAP25 DIST(J) = DIST(J)*(1.0D0 - 4.776D-08*HMEAN) DISMAP26 A = ((X(M1) + X(M2))/2.0D0 - 1.0D06)**2 DISMAP27 B = ((Y(M1) + Y(M2))/2.0D0 - 1.0D06)**2 DISMAP28 C = A + B DISMAP29 DIST(J) = DIST(J)*(1.0D0 + (-8.846D0 + 5.707D0*C*1.0D-11)*1.0D-05)DISMAP30 C DISTA 27 C DISTA 28 C DISTA 29 C FIND COORDINATE DIFFERENCES BETWEEN STATIONS DISTA 30 C DISTA 31 A = Y(M1) DISTA 32 B = Y(M2) DISTA 33 C = X(M1) DISTA 34 D = X(M2) DISTA 35 AB = B - A DISTA 36 CA = DABS(AB) DISTA 37 CD = D - C DISTA 38 DB = DABS(CD) DISTA 39 G = DSQRT(AB*AB + CD*CD) DISTA 40 M = 2*M1 - 1 DISTA 41 I = 0 DISTA 42 C DISTA 43 C COMPUTE DESIGN MATRIX COEFFICIENTS DISTA 44 C DISTA 45 IF(M.LE.MMN) GO TO 5 DISTA 46 IF(CA.LT.0.001D0) GO TO 1 DISTA 47 I = I + 1 DISTA 48 STOR(J,IJK) = -(B-A)/G DISTA 49 PTOR(I) = STOR(J,IJK) DISTA 50 GO TO 2 DISTA 51 1 I = I + 1 DISTA 52 STOR(J,IJK) = 0.0D0 DISTA 53 PTOR(I) = 0.0D0 DISTA 54 2 KK(I) = M - MMN DISTA 55 IF(J.EQ.1) IZ = IZ + 1 DISTA 56 IF(J.EQ.1) ICOL(1) = KK(I) DISTA 57 M = M + 1 DISTA 58 IF(DB.LT.0.001D0) GO TO 3 DISTA 59 I = I + 1 DISTA 60 STOR(J,IKJ) = -(D-C)/G DISTA 61 PTOR(I) = STOR(J,IKJ) DISTA 62 GO TO 4 DISTA 63 3 I = I + 1 DISTA 64 STOR(J,IKJ) = 0.0D0 DISTA 65 PTOR(I) = 0.0D0 DISTA 66 4 KK(I) = M - MMN DISTA 67 5 M = M2*2 - 1 DISTA 68 IF(M.LE.MMN) GO TO 10 DISTA 69 IZ = IZ + 1 DISTA 70 KKL = 2*(J+II/2-NODIST) - 1 DISTA 71 IF((IJK.GT.1).AND.(M2.LT.M1)) KKL = KKL - 2 DISTA 72 IF(CA.LT.0.001D0) GO TO 6 DISTA 73 I = I + 1 DISTA 74 STOR(J,KKL) = (B-A)/G DISTA 75 PTOR(I) = STOR(J,KKL) DISTA 76 GO TO 7 DISTA 77 6 I = I + 1 DISTA 78 STOR(J,KKL) = 0.0D0 DISTA 79 PTOR(I) = 0.0D0 DISTA 80 7 KK(I) = M - MMN DISTA 81 ICOL(IZ) = KK(I) DISTA 82 M = M + 1 DISTA 83 LKL = KKL + 1 DISTA 84 IF(DB.LT.0.001D0) GO TO 8 DISTA 85 I = I + 1 DISTA 86 STOR(J,LKL) = (D-C)/G DISTA 87 PTOR(I) = STOR(J,LKL) DISTA 88 GO TO 9 DISTA 89 8 I = I + 1 DISTA 90 STOR(J,LKL) = 0.0D0 DISTA 91 PTOR(I) = STOR(J,LKL) DISTA 92 9 KK(I) = M - MMN DISTA 93 10 CONTINUE DISTA 94 CONST(J) = G - DIST(J) DISTA 95 PRINT 11,NOBS,NDIS(J),WEIGHT(J),DIST(J),CONST(J) DISTA 96 11 FORMAT(17X,I4,5X,I4,10X,F5.3,10X,F10.3,F31.3,/) DISTA 97 C DISTA 98 C WRITE DESIGN MATRIX COEFFICIENTS AND COLUMN NUMBERS ON A DEVICE DISTA 99 C DISTA100 C DISTA101 IF(I.EQ.4) GO TO 13 DISTA102 IF(I.EQ.2) GO TO 12 DISTA103 T1 = 0.0D0 DISTA104 T2 = 0.0D0 DISTA105 T3 = 0.0D0 DISTA106 T4 = 0.0D0 DISTA107 I1 = 1 DISTA108 I2 = 1 DISTA109 I3 = 1 DISTA110 I4 = 1 DISTA111 GO TO 14 DISTA112 12 T1 = 0.0D0 DISTA113 T2 = 0.0D0 DISTA114 T3 = PTOR(1) DISTA115 T4 = PTOR(2) DISTA116 I1 = 1 DISTA117 I2 = 1 DISTA118 I3 = KK(1) DISTA119 I4 = KK(2) DISTA120 GO TO 14 DISTA121 13 T1 = PTOR(1) DISTA122 T2 = PTOR(2) DISTA123 T3 = PTOR(3) DISTA124 T4 = PTOR(4) DISTA125 I1 = KK(1) DISTA126 I2 = KK(2) DISTA127 I3 = KK(3) DISTA128 I4 = KK(4) DISTA129 14 CONTINUE DISTA130 NCODE = 40 DISTA131 IF(J.EQ.1) WRITE(NSEQF) NODIST,NOBS,NDIS(J),I1,I2,I3,I4,T1,T2,T3,TDISTA133 14,CONST(J),WEIGHT(J),NCODE DISTA133 IF(J.GT.1) WRITE(NSEQF) NOBS,NDIS(J),I1,I2,I3,I4,T1,T2,T3,T4,CONSTDISTA134 1(J),WEIGHT(J),NCODE DISTA135 15 CONTINUE DISTA136 RETURN DISTA137 END DISTA138 SUBROUTINE DISNOR(STOR,WEIGHT,CONST,AL,AA,NODIST,II,NK3,NK4) DISNOR01 C DISNOR02 C ******************************************************************DISNOR03 C * *DISNOR04 C * S U B R O U T I N E D I S N O R *DISNOR05 C * *DISNOR06 C * THIS SUBROUTINE ACCEPTS THE COMPACTED DESIGN MATRIX FROM *DISNOR07 C * THE PARTICULAR OBSERVATION STATION, RESULTING FROM *DISNOR08 C * DISTANCES. THE PARTIAL NORMAL EQUATION MATRIX AND CONSTANT *DISNOR09 C * VECTOR ARE FORMED. *DISNOR10 C * *DISNOR11 C ******************************************************************DISNOR12 C DISNOR13 IMPLICIT REAL*8(A-H,O-Z) DISNOR14 DIMENSION STOR(NK3,NK4),WEIGHT(NK3),CONST(NK3),AL(NK3),AA(NK4,NK4)DISNOR15 C DISNOR16 C WEIGHT EQUALS ONE DIVIDED BY VARIANCE DISNOR17 C DISNOR18 DO 1 I = 1,NODIST DISNOR19 1 WEIGHT(I) = 1.0D0/WEIGHT(I)/WEIGHT(I) DISNOR20 C DISNOR21 C MULTIPLY A TRANSPOSED BY P BY A MATRIX (PARTIAL) DISNOR22 C DISNOR23 SUM = 0.0D0 DISNOR24 DO 3 I = 1,II DISNOR25 DO 3 J = I,II DISNOR26 DO 2 K = 1,NODIST DISNOR27 2 SUM = SUM + STOR(K,I)*WEIGHT(K)*STOR(K,J) DISNOR28 AA(I,J) = SUM DISNOR29 3 SUM = 0.0D0 DISNOR30 C DISNOR31 C MULTIPLY A TRANSPOSED BY P BY CONSTANT VECTOR 'W' DISNOR32 C DISNOR33 SUM = 0.0D0 DISNOR34 DO 5 I = 1,II DISNOR35 DO 4 K = 1,NODIST DISNOR36 4 SUM = SUM + STOR(K,I)*WEIGHT(K)*CONST(K) DISNOR37 AL(I) = SUM DISNOR38 5 SUM = 0.0D0 DISNOR39 RETURN DISNOR40 END DISNOR41 SUBROUTINE COLUMN(II,ICOL,JCOL,NK3,NK4) COLUMN01 C COLUMN02 C ******************************************************************COLUMN03 C * *COLUMN04 C * S U B R O U T I N E C O L U M N *COLUMN05 C * *COLUMN06 C * THIS SUBROUTINE COMPUTES THE COLUMN INDICIES FOR THE *COLUMN07 C * ADDITION OF THE PARTIAL NORMAL EQUATION MATRIX TO THE *COLUMN08 C * NORMAL EQUATION MATRIX. THE COMPUTED COLUMN INDICIES ARE *COLUMN09 C * REPRESENTIVE OF A FULL MATRIX NOT THE COMPACTED MATRIX. *COLUMN10 C * REDUCTION TO THE COMPACTED MATRIX COMES LATER. *COLUMN11 C * *COLUMN12 C ******************************************************************COLUMN13 C COLUMN14 DIMENSION ICOL(NK3),JCOL(NK4) COLUMN15 IT = II/2 - 1 COLUMN16 IF(IT.EQ.0) GO TO 2 COLUMN17 DO 1 I = 1,IT COLUMN18 IF(ICOL(I).LT.ICOL(I+1)) GO TO 1 COLUMN19 IA = ICOL(I) COLUMN20 IB = ICOL(I+1) COLUMN21 ICOL(I) = IB COLUMN22 ICOL(I+1) = IA COLUMN23 1 CONTINUE COLUMN24 2 IT = II/2 COLUMN25 K = 0 COLUMN26 IF(IT.EQ.0) GO TO 4 COLUMN27 DO 3 I = 1,IT COLUMN28 K = K + 1 COLUMN29 JCOL(K) = ICOL(I) COLUMN30 K = K + 1 COLUMN31 3 JCOL(K) = ICOL(I) + 1 COLUMN32 GO TO 5 COLUMN33 4 JCOL(1) = ICOL(1) COLUMN34 JCOL(2) = ICOL(1) + 1 COLUMN35 5 CONTINUE COLUMN36 RETURN COLUMN37 END COLUMN38 SUBROUTINE ADDNOR(TA,AA,II,JCOL,NK2,NK4,NK5,NSPACE,NW) ADDNOR01 C ADDNOR02 C ******************************************************************ADDNOR03 C * *ADDNOR04 C * S U B R O U T I N E A D D N O R *ADDNOR05 C * *ADDNOR06 C * THIS SUBROUTINE ADDS THE PARTIAL NORMAL EQUATION MATRIX *ADDNOR07 C * TO THE NORMAL EQUATION MATRIX. THE NORMAL EQUATION MATRIX *ADDNOR08 C * IS COMPACTED INTO THE ARROWHEAD PATTERN. THIS IS *ADDNOR09 C * ACCOMPLISHED WITH THE USE OF NSPACE, NSPACE EQUALS THE *ADDNOR10 C * NUMBER OF ZEROS BETWEEN THE EDGE OF THE BAND AND THE EDGE *ADDNOR11 C * OF THE BORDER IN THE FIRST ROW OF THE MATRIX. NSPACE IS *ADDNOR12 C * SUBSTRACTED FROM THE COLUMN INDICIE FOR THE FIRST ROW, FOR *ADDNOR13 C * THE SECOND ROW WE MUST ADD ONE, THIRD ROW ADD TWO, ETC. *ADDNOR14 C * *ADDNOR15 C ******************************************************************ADDNOR16 C ADDNOR17 IMPLICIT REAL*8(A-H,O-Z) ADDNOR18 DIMENSION TA(NK5,NK2),AA(NK4,NK4),JCOL(NK4) ADDNOR19 DO 2 I = 1,II ADDNOR20 IJC = JCOL(I) ADDNOR21 DO 2 K = I,II ADDNOR22 KJC = JCOL(K) ADDNOR23 IF(KJC.LE.(NW+IJC-1)) GO TO 1 ADDNOR24 NSKIP = NSPACE - IJC ADDNOR25 IF(NSKIP.LT.0) NSKIP = 0 ADDNOR26 KJC = KJC - NSKIP ADDNOR27 1 IKJC = KJC - IJC + 1 ADDNOR28 2 TA(IKJC,IJC) = TA(IKJC,IJC) + AA(I,K) ADDNOR29 RETURN ADDNOR30 END ADDNOR31 SUBROUTINE ADDVEC(BL,AL,II,JCOL,NK2,NK3,NK4) ADDVEC01 C ADDVEC02 C ******************************************************************ADDVEC03 C * *ADDVEC04 C * S U B R O U T I N E A D D V E C *ADDVEC05 C * *ADDVEC06 C * THIS SUBROUTINE ADDS THE PARTIAL CONSTANT VECTOR ATPW TO *ADDVEC07 C * THE CONSTANT VECTOR ATPW. *ADDVEC08 C * *ADDVEC09 C ******************************************************************ADDVEC10 C ADDVEC11 IMPLICIT REAL*8(A-H,O-Z) ADDVEC12 DIMENSION BL(NK2),AL(NK3),JCOL(NK4) ADDVEC13 DO 1 I = 1,II ADDVEC14 IJC = JCOL(I) ADDVEC15 1 BL(IJC) = BL(IJC) + AL(I) ADDVEC16 RETURN ADDVEC17 END ADDVEC18 SUBROUTINE ARROW(A,B,XVEC,W,NNB,NWB,NB,NW,N,NBB,NQQ,IPREAN) ARROW 01 C ARROW 02 C ******************************************************************ARROW 03 C * *ARROW 04 C * S U B R O U T I N E A R R O W *ARROW 05 C * *ARROW 06 C * THIS SUBROUTINE SOLVES A BANDED AND BORDERED SYMMETRIC *ARROW 07 C * POSITIVE DEFINITE SET OF NORMAL EQUATIONS. THE BANDED AND *ARROW 08 C * BORDERED PORTION ONLY OF THE VARIANCE CO-VARIANCE MATRIX *ARROW 09 C * IS SOLVED FOR,TOGETHER WITH THE X-VECTOR. THE METHOD USED *ARROW 10 C * IS THAT OF CHOLESKI. *ARROW 11 C * *ARROW 12 C ******************************************************************ARROW 13 C ARROW 14 REAL*8 A(NWB,NNB),B(NNB),XVEC(NNB),W(NNB) ARROW 15 C ARROW 16 C THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES ARROW 17 C 1) TRIANG ARROW 18 C 2) FORWAR AND BACK ARROW 19 C 3) ARRINV ARROW 20 C ARROW 21 C FIND THE CHOLESKI SQUARE ROOT ARROW 22 C ARROW 23 CALL TRIANG(A,NNB,NWB,NB,NW,N,NBB) ARROW 24 IF(IPREAN.EQ.1) GO TO 4 ARROR 2A C ARROW 25 C CLEAR VECTOR B ARROW 26 C ARROW 27 DO 2 I = 1,NNB ARROW 28 2 B(I) = 0.0D0 ARROW 29 C ARROW 30 C COMPUTE FORWARD SOLUTION OF X-VECTOR ARROW 31 C ARROW 32 CALL FORWAR(A,W,B,NNB,NWB,NB,NW,N) ARROW 33 C ARROW 34 C CLEAR X-VECTOR ARROW 35 C ARROW 36 DO 3 I = 1,NNB ARROW 37 3 XVEC(I) = 0.0D0 ARROW 38 C ARROW 39 C COMPUTE BACK SOLUTION FOR X-VECTOR ARROW 40 C ARROW 41 CALL BACK(A,B,XVEC,NNB,NWB,NB,NW,N) ARROW 42 C ARROW 43 C ELEMENT (N,N) OF VARIANCE CO-VARIANCE IS ARROW 44 C ARROW 45 4 A(1,NNB) = 1.0D0/A(1,NNB)/A(1,NNB) ARROW 46 MW = NWB - 1 ARROW 47 I = NNB ARROW 48 J1 = 1 ARROW 49 I2 = NWB - 1 ARROW 50 C ARROW 51 C COMPUTE LOWER BLOCK OF VARIANCE CO-VARIANCE MATRIX ARROW 52 C ARROW 53 CALL ARRINV(A,B,NNB,NWB,NB,NW,N,MW,I,J1,I2) ARROW 54 J1 = NWB ARROW 55 MW = NNB - 1 ARROW 56 C ARROW 57 C COMPUTE REMAINING PORTION OF THE VARIANCE CO-VARIANCE PERTAINING ARROW 58 C TO THE BANDED AND BORDERED PORTION OF THE NORMAL EQUATION MATRIX ARROW 59 C ARROW 60 CALL ARRINV(A,B,NNB,NWB,NB,NW,N,MW,I,J1,I2) ARROW 61 RETURN ARROW 62 END ARROW 63 SUBROUTINE TRIANG(A,NNB,NWB,NB,NW,N,NBB) TRIANG01 C TRIANG02 C ******************************************************************TRIANG03 C * *TRIANG04 C * S U B R O U T I N E T R I A N G *TRIANG05 C * *TRIANG06 C * THIS SUBROUTINE TAKES THE BANDED AND BORDERED NORMAL EQUATION*TRIANG07 C * MATRIX AND PERFORMS THE CHOLESKI SQUARE ROOT (THAT IS *TRIANG08 C * DECOMPOSES THE MATRIX INTO TWO TRIANGULAR MATRICES L AND LT; *TRIANG09 C * WHERE LT IS THE TRANSPOSE OF L AND LT*L = N) *TRIANG10 C * *TRIANG11 C ******************************************************************TRIANG12 C TRIANG13 C TRIANG14 C FIND SQRT OF FIRST DIAGONAL COEFFICIENT TRIANG15 C TRIANG16 REAL*8 A(NWB,NNB),SUM,DSQRT TRIANGAA A(1,1) = DSQRT(A(1,1)) TRIANG17 KWNB = NW TRIANG18 KW = NW + NB TRIANG19 II = 0 TRIANG20 C TRIANG21 C DIVIDE REMAINING COEFFICIENTS IN FIRST ROW BY A(1,1) TRIANG22 C TRIANG23 DO 1 J = 2,NWB TRIANG24 1 A(J,1) = A(J,1)/A(1,1) TRIANG25 MW = N - NW + 1 TRIANG26 C TRIANG27 C FIND THE CHOLESKI SQUARE ROOT FOR THE REMAINDER OF THE BAND AND TRIANG28 C BORDER DOWN TO WHERE THE BAND BEGINS TO DECREASE IN SIZE. TRIANG29 C TRIANG30 DO 2 I = 2,MW TRIANG31 I1 = I - 1 TRIANG32 IF(I.GT.NW) I1 = NW - 1 TRIANG33 CALL TRI1(A,NW,I,NB,NWB,NBB,I1,II,KW,NNB) TRIANG34 2 CONTINUE TRIANG35 JW = NW + 1 TRIANG36 L1 = NW - 1 TRIANG37 C TRIANG38 C FIND THE CHOLESKI SQUARE ROOT FOR THE REMAINDER OF THE BAND AND TRIANG39 C BORDER DOWN TO WHERE THE BAND VANISHES. TRIANG40 C TRIANG41 DO 3 II = 1,L1 TRIANG42 JW = JW - 1 TRIANG43 I = II + MW TRIANG44 KW = KW - 1 TRIANG45 CALL TRI1(A,NW,I,NB,NWB,NBB,L1,II,KW,NNB) TRIANG46 3 CONTINUE TRIANG47 MW = N + 1 TRIANG48 LL = 0 TRIANG49 ND = NB + 1 TRIANG50 C TRIANG51 C UPDATE N22 TO EQUAL N22 - L12'*L12 TRIANG52 C TRIANG53 DO 5 I = 1,NB TRIANG54 INW = I + NW TRIANG55 DO 5 J = I,NB TRIANG56 JNW = J + NW TRIANG57 SUM = 0.0D0 TRIANG58 KNW = 0 TRIANG59 DO 4 K = 1,N TRIANG60 IF(K.GT.(N-NW+1)) KNW = KNW + 1 TRIANG61 4 SUM = SUM + A(INW - KNW,K)*A(JNW - KNW,K) TRIANG62 5 A(J-I+1,I+N) = A(J-I+1,I+N) - SUM TRIANG63 C TRIANG64 C FIND CHOLESKI SQUARE ROOT FOR N22 TRIANG65 C TRIANG66 DO 8 I = MW,NNB TRIANG67 ND = ND - 1 TRIANG68 DO 6 J = 1,ND TRIANG69 IF(J.EQ.1) A(J,I) = DSQRT(A(J,I)) TRIANG70 IF(J.NE.1) A(J,I) = A(J,I)/A(1,I) TRIANG71 6 CONTINUE TRIANG72 IF(I.EQ.NNB) GO TO 8 TRIANG73 DO 7 K = 2,ND TRIANG74 DO 7 J = K,ND TRIANG75 7 A(J-K+1,I+K-1) = A(J-K+1,I+K-1) - A(K,I)*A(J,I) TRIANG76 8 CONTINUE TRIANG77 RETURN TRIANG78 END TRIANG79 SUBROUTINE TRI1(A,NW,I,NB,NWB,NBB,I1,II,KW,NNB) TRI1 01 C TRI1 02 C ******************************************************************TRI1 03 C * *TRI1 04 C * S U B R O U T I N E T R I 1 *TRI1 05 C * *TRI1 06 C * THIS SUBROUTINE IS CALLED BY TRIANG AND IS USED TO FIND *TRI1 07 C * THE CHOLESKI SQUARE ROOT. THE SUBROUTINE DOES A DIFFERENT *TRI1 08 C * JOB EACH TIME IT IS CALLED SO CHECK THE COMMENTS IN TRIANG *TRI1 09 C * FOR THE BASIC IDEA BEHIND IT. *TRI1 10 C * *TRI1 11 C ******************************************************************TRI1 12 C TRI1 13 REAL*8 A(NWB,NNB),DSQRT TRI1 14 IF(II.LE.1) J1 = NW - 1 TRI1 15 IF(II.GT.1) J1 = J1 - 1 TRI1 16 DO 1 K = 1,I1 TRI1 17 IK = I - K TRI1 18 K2 = K + 1 TRI1 19 1 A(1,I) = A(1,I) - A(K2,IK)*A(K2,IK) TRI1 20 A(1,I) = DSQRT(A(1,I)) TRI1 21 IF(II.EQ.I1) GO TO 4 TRI1 22 IF(II.EQ.0) A(NW,I) = A(NW,I)/A(1,I) TRI1 23 IF(II.EQ.0) K1 = 1 TRI1 24 IF(II.NE.0) K1 = II TRI1 25 DO 3 J = 2,J1 TRI1 26 JJ = KW - J + 2 - NB TRI1 27 IF(II.EQ.0) JJ = NW - J + 1 TRI1 28 DO 2 K = 1,K1 TRI1 29 IK = I - K TRI1 30 2 A(JJ,I) = A(JJ,I) - A(K+1,IK)*A(JJ+K,IK) TRI1 31 IF(K1.LT.I1) K1 = K1 + 1 TRI1 32 3 A(JJ,I) = A(JJ,I)/A(1,I) TRI1 33 4 CONTINUE TRI1 34 DO 8 J = 1,NB TRI1 35 JJ = J + KW - NB TRI1 36 IF(JJ.EQ.(NW+J)) GO TO 6 TRI1 37 DO 5 K = 1,I1 TRI1 38 IK = I - K TRI1 39 IF(K.LE.II) K3 = K TRI1 40 5 A(JJ,I) = A(JJ,I) - A(JJ+K3,IK)*A(K+1,IK) TRI1 41 GO TO 8 TRI1 42 6 DO 7 K = 1,I1 TRI1 43 IK = I - K TRI1 44 7 A(JJ,I) = A(JJ,I) - A(JJ,IK)*A(K+1,IK) TRI1 45 8 A(JJ,I) = A(JJ,I)/A(1,I) TRI1 46 RETURN TRI1 47 END TRI1 48 SUBROUTINE FORWAR(A,W,B,NNB,NWB,NB,NW,N) FORWAR01 C FORWAR02 C ******************************************************************FORWAR03 C * *FORWAR04 C * S U B R O U T I N E F O R W A R *FORWAR05 C * *FORWAR06 C * THIS SUBROUTINE COMPUTES THE FORWARD SOLUTION OF THE *FORWAR07 C * X-VECTOR FOR THE BANDED AND BORDERED NORMAL EQUATIONS, *FORWAR08 C * METHOD OF CHOLESKI. *FORWAR09 C * *FORWAR10 C ******************************************************************FORWAR11 C FORWAR12 REAL*8 A(NWB,NNB),W(NNB),B(NNB) FORWAR13 B(1) = W(1)/A(1,1) FORWAR14 DO 1 K = 1,NB FORWAR15 IK = K + N FORWAR16 1 B(IK) = B(IK) + B(1)*A(NW+K,1) FORWAR17 KW = NW FORWAR18 MW = N - NW + 1 FORWAR19 DO 4 I = 2,N FORWAR20 IF(I.GT.MW) KW = KW - 1 FORWAR21 I1 = I - 1 FORWAR22 IF(I.GT.NW) I1 = NW - 1 FORWAR23 DO 2 K = 1,I1 FORWAR24 IK = I - K FORWAR25 2 B(I) = B(I) + B(IK)*A(K+1,IK) FORWAR26 B(I) = (W(I) - B(I))/A(1,I) FORWAR27 DO 3 K = 1,NB FORWAR28 IK = K + N FORWAR29 3 B(IK) = B(IK) + B(I)*A(KW+K,I) FORWAR30 4 CONTINUE FORWAR31 NBB = NB FORWAR32 MW = N + 1 FORWAR33 DO 6 I = MW,NNB FORWAR34 NBB =NBB - 1 FORWAR35 B(I) = (W(I) - B(I))/A(1,I) FORWAR36 IF(NBB.EQ.0) GO TO 6 FORWAR37 DO 5 K = 1,NBB FORWAR38 IK = I + K FORWAR39 5 B(IK) = B(IK) + B(I)*A(K+1,I) FORWAR40 6 CONTINUE FORWAR41 RETURN FORWAR42 END FORWAR43 SUBROUTINE BACK(A,B,XVEC,NNB,NWB,NB,NW,N) BACK01 C BACK02 C ****************************************************************** BACK03 C * * BACK04 C * S U B R O U T I N E B A C K * BACK05 C * * BACK06 C * THIS SUBROUTINE COMPUTES THE BACK SOLUTION OF THE X-VECTOR * BACK07 C * FOR THE BANDED AND BORDERED NORMAL EQUATIONS, METHOD OF * BACK08 C * CHOLESKI * BACK09 C * * BACK10 C ****************************************************************** BACK11 C BACK12 REAL*8 A(NWB,NNB),B(NNB),XVEC(NNB) BACK13 XVEC(NNB) = B(NNB)/A(1,NNB) BACK14 DO 3 I = 2,NB BACK15 I1 = I -1 BACK16 I2 = NNB - I + 1 BACK17 DO 2 K = 1,I1 BACK18 2 XVEC(I2) = XVEC(I2) + XVEC(I2+K)*A(K+1,I2) BACK19 3 XVEC(I2) = (B(I2)-XVEC(I2))/A(1,I2) BACK20 DO 5 I = 1,NW BACK21 I2 = N - I + 1 BACK22 I1 = I1 + 1 BACK23 DO 4 K = 1,I1 BACK24 4 XVEC(I2) = XVEC(I2) + XVEC(I2+K)*A(K+1,I2) BACK25 5 XVEC(I2) = (B(I2) - XVEC(I2))/A(1,I2) BACK26 MW = N - NW BACK27 I1 = NW - 1 BACK28 DO 8 I = 1,MW BACK29 I2 = I2 - 1 BACK30 DO 6 K = 1,NB BACK31 6 XVEC(I2) = XVEC(I2) + XVEC(NNB+1-K)*A(NWB+1-K,I2) BACK32 DO 7 K = 1,I1 BACK33 7 XVEC(I2) = XVEC(I2) + XVEC(I2+K)*A(K+1,I2) BACK34 8 XVEC(I2) = (B(I2) - XVEC(I2))/A(1,I2) BACK35 RETURN BACK36 END BACK37 SUBROUTINE ARRINV(A,B,NNB,NWB,NB,NW,N,MW,I,J1,I2) ARRINV01 C ARRINV02 C ******************************************************************ARRINV03 C * *ARRINV04 C * S U B R O U T I N E A R R I N V *ARRINV05 C * *ARRINV06 C * THIS SUBROUTINE TAKES THE CHOLESKI DECOMPOSED MATRIX AND *ARRINV07 C * DOES A BACK SOLUTION FOR THE BANDED AND BORDERED PORTION *ARRINV08 C * OF THE VARIANCE CO-VARIANCE MATRIX. *ARRINV09 C * *ARRINV10 C ******************************************************************ARRINV11 C ARRINV12 REAL*8 A(NWB,NNB),B(NNB),P ARRINV13 DO 13 II = J1,MW ARRINV14 I = I - 1 ARRINV15 IK = 0 ARRINV16 I1 = II ARRINV17 IF(II.GE.NWB) I1 = NW - 1 ARRINV18 DO 5 J = 1,I1 ARRINV19 B(J) = 0.0D0 ARRINV20 DO 1 K = 1,J ARRINV21 1 B(J) = B(J) + A(K+1,I)*A(J-K+1,I+K) ARRINV22 M = I + J ARRINV23 IF(J.EQ.I1) GO TO 3 ARRINV24 L = J + 1 ARRINV25 DO 2 K = L,I1 ARRINV26 KI = K + 1 ARRINV27 2 B(J) = B(J) + A(KI,I)*A(KI-J,M) ARRINV28 3 IF(II.LT.NWB) GO TO 5 ARRINV29 L = I1 + 1 ARRINV30 IF((II-J).LT.I2) IK = IK - 1 ARRINV31 DO 4 K = L,I2 ARRINV32 KI = K + 1 ARRINV33 4 B(J) = B(J) + A(KI,I)*A(KI+IK,M) ARRINV34 5 B(J) = -B(J)/A(1,I) ARRINV35 IF(II.LT.NWB) GO TO 10 ARRINV36 I3 = I1 + 1 ARRINV37 I4 = NWB - 1 ARRINV38 DO 7 J = I3,I4 ARRINV39 IK = 0 ARRINV40 B(J) = 0.0D0 ARRINV41 DO 6 K = 1,I1 ARRINV42 IF((II-K).LT.I2) IK = IK - 1 ARRINV43 6 B(J) = B(J) + A(K+1,I)*A(J+1+IK,I+K) ARRINV44 7 CONTINUE ARRINV45 DO 9 J = I3,I4 ARRINV46 L = 1 ARRINV47 JJ = J - I3 + 1 ARRINV48 DO 8 K = I3,I4 ARRINV49 B(J) = B(J) + A(K+1,I)*A(JJ,N+L) ARRINV50 L = L + 1 ARRINV51 JJ = JJ - 1 ARRINV52 IF(L.GT.(J-I3+1)) JJ = JJ + 2 ARRINV53 IF(L.GT.(J-I3+1)) L = L - 1 ARRINV54 8 CONTINUE ARRINV55 9 B(J) = -B(J)/A(1,I) ARRINV56 10 P = 0.0D0 ARRINV57 K = I2 ARRINV58 IF(II.LT.NWB) K = I1 ARRINV59 DO 11 J = 1,K ARRINV60 11 P = P + B(J)*A(J+1,I) ARRINV61 A(1,I) = (1.0D0/A(1,I)-P)/A(1,I) ARRINV62 DO 12 J = 1,K ARRINV63 12 A(J+1,I) = B(J) ARRINV64 13 CONTINUE ARRINV65 RETURN ARRINV66 END ARRINV67 SUBROUTINE RESIDL(XVEC,STORAZ,STORDR,STORDS,MAZ,MDIR,MDIS,WTAZ,WTDRESID 01 #R,WTDS,WVAZ,WVDR,WVDS,VARFAC,NK2,NK3,NSEQF,NOSETS,NK1) RESID 02 C RESID 03 C ******************************************************************RESID 04 C * *RESID 05 C * S U B R O U T I N E R E S I D L *RESID 06 C * *RESID 07 C * THIS SUBROUTINE READS OFF FILE NSEQF THE COMPACTED DESIGN *RESID 08 C * MATRIX ROW BY ROW , MULTIPLIES IT BY THE X-VECTOR TO COMPUTE *RESID 09 C * THE RESIDUALS. THE SUBROUTINE COMPUTES THE RESIDUALS FOR *RESID 10 C * THE AZIMUTHS, THEN THE DIRECTIONS AND THEN THE DISTANCES. *RESID 11 C * *RESID 12 C ******************************************************************RESID 13 C RESID 14 IMPLICIT REAL*8(A-H,O-Z) RESID 15 DIMENSION XVEC(NK2),STORAZ(NK3,4),STORDR(NK3,4),STORDS(NK3,4) RESID 16 DIMENSION MAZ(NK3,4),MDIR(NK3,4),MDIS(NK3,4),WTAZ(NK3) RESID 17 DIMENSION WTDR(NK3),WTDS(NK3),WVAZ(NK3),WVDR(NK3) RESID 18 DIMENSION WVDS(NK3),V(20),ISTA(20),JSTA(20) RESID 19 PRINT 1 RESID 20 1 FORMAT('1',/////,46X,'R E S I D U A L S') RESID 21 PRINT 2 RESID 22 2 FORMAT(' ',//,30X,'INST',4X,'TARGET',4X,'MISSCLOSURE',7X,'RESIDUALRESID 23 #',8X,'CODE') RESID 24 VARFAC = 0.0D0 RESID 25 REWIND NSEQF RESID 26 IR = 0 RESID 27 DO 20 IJK = 1,NOSETS RESID 28 MEQ = 1 RESID 29 I = 0 RESID 30 J = 0 RESID 31 READ(NSEQF) NEQ,ISTA(1),JSTA(1),I1,I2,I3,I4,T1,T2,T3,T4,WV,EV,NCODRESID 32 1E RESID32A EV = 1.0D0/EV/EV RESID32B GO TO 4 RESID 33 3 MEQ = MEQ + 1 RESID 34 READ(NSEQF) ISTA(MEQ),JSTA(MEQ),I1,I2,I3,I4,T1,T2,T3,T4,WV,EV,NCODRESID 35 1E RESID35A EV = 1.0D0/EV/EV RESID35B 4 IF(NCODE.EQ.30) GO TO 5 RESID 36 IF(NCODE.EQ.40) GO TO 6 RESID 37 I = I + 1 RESID 38 CALL READAZ(STORAZ,T1,T2,T3,T4,I1,I2,I3,I4,MAZ,WTAZ,WVAZ,EV,WV,I,NRESID 39 #K3) RESID 40 GO TO 6 RESID 41 5 J = J + 1 RESID 42 CALL READAZ(STORDR,T1,T2,T3,T4,I1,I2,I3,I4,MDIR,WTDR,WVDR,EV,WV,J,RESID 43 #NK3) RESID 44 6 CONTINUE RESID 45 IF(MEQ.LT.NEQ) GO TO 3 RESID 46 C RESID 47 C FORM RESIDUALS OF AZIMUTHS AND DIRECTIONS RESID 48 C RESID 49 IF(I.EQ.0) GO TO 11 RESID 50 DO 10 IN = 1,I RESID 51 IR = IR + 1 RESID 52 V(IN) = 0.0D0 RESID 53 DO 7 INN = 1,4 RESID 54 7 V(IN) = STORAZ(IN,INN)*XVEC(MAZ(IN,INN)) + V(IN) RESID 55 V(IN) = V(IN) - WVAZ(IN) RESID 56 C RESID 57 C RESID 58 IF((IR/25*25).NE.IR) GO TO 8 RESID 59 PRINT 1 RESID 60 PRINT 2 RESID 61 8 PRINT 9,IN,ISTA(IN),JSTA(IN),WVAZ(IN),V(IN) RESID 62 9 FORMAT(' ',/,I28,I6,I10,2F15.1,' AZIMUTH') RESID 63 10 VARFAC = VARFAC + V(IN)*WTAZ(IN)*V(IN) RESID 64 11 IF(J.EQ.0) GO TO 20 RESID 65 SUM = 0.0D0 RESID 66 DO 12 IN = 1,J RESID 67 12 SUM = SUM + WTDR(IN) RESID 68 SUM = -1.0D0/SUM RESID 69 DO 14 IN = 1,J RESID 70 C RESID 71 V(IN) = 0.0D0 RESID 72 DO 13 INN = 1,4 RESID 73 13 V(IN) = STORDR(IN,INN)*XVEC(MDIR(IN,INN)) + V(IN) RESID 74 14 V(IN) = V(IN) - WVDR(IN) RESID 75 ASUM = 0.0D0 RESID 76 DO 15 IN = 1,J RESID 77 15 ASUM = ASUM + WTDR(IN)*V(IN) RESID 78 DO 16 IN = 1,J RESID 79 16 V(IN) = V(IN) + SUM*ASUM RESID 80 DO 19 IN = 1,J RESID 81 IR = IR + 1 RESID81A C RESID 82 C RESID 83 IF((IR/25*25).NE.IR) GO TO 17 RESID 84 PRINT 1 RESID 85 PRINT 2 RESID 86 17 PRINT 18,IN,ISTA(IN),JSTA(IN),WVDR(IN),V(IN) RESID 87 18 FORMAT(' ',/,I28,I6,I10,2F15.1,' DIRECTION') RESID 88 19 VARFAC = VARFAC + V(IN)*WTDR(IN)*V(IN) RESID 89 20 CONTINUE RESID 90 REWIND NSEQF RESID 91 C RESID 92 DO 28 IJK = 1,NOSETS RESID 93 MEQ = 1 RESID 94 K = 0 RESID 95 READ(NSEQF) NEQ,ISTA(1),JSTA(1),I1,I2,I3,I4,T1,T2,T3,T4,WV,EV,NCODRESID 96 1E RESID96A EV = 1.0D0/EV/EV RESID96B GO TO 22 RESID 97 21 MEQ = MEQ + 1 RESID 98 READ(NSEQF) ISTA(MEQ),JSTA(MEQ),I1,I2,I3,I4,T1,T2,T3,T4,WV,EV,NCODRESID 99 1E RESID99A EV = 1.0D0/EV/EV RESID99B 22 IF(NCODE.EQ.30) GO TO 23 RESID100 IF(NCODE.EQ.50) GO TO 23 RESID101 K = K + 1 RESID102 CALL READAZ(STORDS,T1,T2,T3,T4,I1,I2,I3,I4,MDIS,WTDS,WVDS,EV,WV,K,RESID102 #NK3) RESID103 23 CONTINUE RESID104 IF(MEQ.LT.NEQ) GO TO 21 RESID105 C RESID106 C COMPUTE DISTANCE RESIDUALS RESID107 C RESID108 IF(K.EQ.0) GO TO 28 RESID109 DO 27 IN = 1,K RESID110 IR = IR + 1 RESID111 V(IN) = 0.0D0 RESID112 DO 24 INN = 1,4 RESID113 24 V(IN) = V(IN) + STORDS(IN,INN)*XVEC(MDIS(IN,INN)) RESID114 V(IN) = V(IN) - WVDS(IN) RESID115 C RESID116 C RESID117 IF((IR/25*25).NE.IR) GO TO 25 RESID118 PRINT 1 RESID119 PRINT 2 RESID120 25 PRINT 26,IN,ISTA(IN),JSTA(IN),WVDS(IN),V(IN) RESID121 26 FORMAT(' ',/,I28,I6,I10,2F15.3,' DISTANCE') RESID122 27 VARFAC = VARFAC + V(IN)*WTDS(IN)*V(IN) RESID123 28 CONTINUE RESID124 RETURN RESID125 END RESID126 SUBROUTINE ABSOLU(TA,NK1,NK2,NK5,CFACT,NAME) ABSOLU01 C ABSOLU02 C ******************************************************************ABSOLU03 C * *ABSOLU04 C * S U B R O U T I N E A B S O L U *ABSOLU05 C * *ABSOLU06 C * THIS SUBROUTINE COMPUTES 95 % ABSOLUTE ERROR ELLIPSES FOR *ABSOLU07 C * EACH STATION ADJUSTED. A EQUALS THE SEMI-MAJOR AXIS, B *ABSOLU08 C * EQUALS THE SEMI-MINOR AXIS, BOTH HAVING SAME UNITS AS THE *ABSOLU09 C * INPUT COORDINATES. ANGLE IS THE ORIENTATION FROM NORTH OF *ABSOLU10 C * THE SEMI-MAJOR AXIS,POSITIVE IS CLOCKWISE FROM NORTH, *ABSOLU11 C * UNITS OF MEASURE IS DEGREES. *ABSOLU12 C * *ABSOLU13 C ******************************************************************ABSOLU14 C ABSOLU15 IMPLICIT REAL*8(A-H,O-Z) ABSOLU16 DIMENSION TA(NK5,NK2),NAME(NK1) ABSOLU17 PRINT 1 ABSOLU18 1 FORMAT('1',/////,39X,'ABSOLUTE 95 % ERROR ELLIPSES',//) ABSOLU19 PRINT 2 ABSOLU20 2 FORMAT(15X,'STATION NO.',10X,'SEMI MAJOR AXIS',10X,'SEMI MINOR AXIABSOLU21 #S',15X,'ANGLE',///) ABSOLU22 I = 1 ABSOLU23 K = NK1 - NK2/2 ABSOLU24 3 QYY = TA(1,I) ABSOLU25 QXY = TA(2,I) ABSOLU26 I = I + 1 ABSOLU27 QXX = TA(1,I) ABSOLU28 A = DSQRT((QXX+QYY)/2.0D0 + DSQRT((QYY-QXX)**2/4.0D0 + QXY**2)) ABSOLU29 B = DSQRT((QXX+QYY)/2.0D0 - DSQRT((QYY-QXX)**2/4.0D0 + QXY**2)) ABSOLU30 C = -2.0*QXY ABSOLU31 D = QYY - QXX ABSOLU32 C = -DATAN2(C,D)*57.2957795D0/2.0D0 ABSOLU33 A = A*CFACT ABSOLU34 B = B*CFACT ABSOLU35 K = K + 1 ABSOLU36 PRINT 4,NAME(K),A,B,C ABSOLU37 4 FORMAT(I22,F24.3,F25.3,F25.2,/) ABSOLU38 I = I + 1 ABSOLU39 M = K - NK1 + NK2/2 ABSOLU40 IF((M/26*26).EQ.M) PRINT 1 ABSOLU41 IF((M/26*26).EQ.M) PRINT 2 ABSOLU42 IF(I.LT.NK2) GO TO 3 ABSOLU43 RETURN ABSOLU44 END ABSOLU45 SUBROUTINE RELATI(TA,NK1,NK2,NK5,CFACT,NSPACE,NW,NAME) RELATI01 C RELATI02 C ******************************************************************RELATI03 C * *RELATI04 C * S U B R O U T I N E R E L A T I *RELATI05 C * *RELATI06 C * THIS SUBROUTINE COMPUTES 95% RELATIVE ERROR ELLIPSES *RELATI07 C * BETWEEN THOSE STATIONS REQUESTED. A EQUALS THE SEMI-MAJOR *RELATI08 C * AXIS, B EQUALS THE SEMI-MINOR AXIS, BOTH HAVING SAME UNITS *RELATI09 C * AS THE INPUT COORDINATES. ANGLE IS THE ORIENTATION FROM *RELATI10 C * NORTH OF THE SEMI-MINOR AXIS, POSITIVE IS CLOCKWISE FROM *RELATI11 C * NORTH, UNITS OF MEASURE IS DEGREES *RELATI12 C * *RELATI13 C ******************************************************************RELATI14 C RELATI15 IMPLICIT REAL*8(A-H,O-Z) RELATI16 DIMENSION TA(NK5,NK2),NAME(NK1) RELATI17 PRINT 1 RELATI18 1 FORMAT('1',/////,39X,'RELATIVE 95 % ERROR ELLIPSES',//) RELATI19 PRINT 2 RELATI20 2 FORMAT(10X,'STATION TO STATION',10X,'SEMI MAJOR AXIS',10X,'SEMI MIRELATI21 #NOR AXIS',15X,'ANGLE',//) RELATI22 MMN = 2*NK1 - NK2 RELATI23 C RELATI24 C READ NUMBER OF RELATIVE ERROR ELLIPSES RELATI25 C RELATI26 READ 3,NNNN RELATI27 3 FORMAT(I4) RELATI28 C RELATI29 C READ STATION NUMBERS BETWEEN WHICH THE RELATIVE ERROR ELLIPSE IS RELATI30 C REQUIRED RELATI31 C RELATI32 C THE SMALLEST STATION NUMBER MUST PRECEDE THE LARGER NUMBER RELATI33 C RELATI34 C RELATI35 K = 0 RELATI36 4 READ 5,NSTA1,NSTA2 RELATI37 5 FORMAT(2I5) RELATI38 IF(NSTA1.LT.25) GO TO 4 RELATI39 I = 2*NSTA1 - 1 - MMN RELATI40 A = TA(1,I) RELATI41 B = TA(2,I) RELATI42 I = I + 1 RELATI43 C = TA(1,I) RELATI44 I = I - 1 RELATI45 J = (NSTA2 - NSTA1)*2 + 1 RELATI46 IF(J.GT.NW) J = J - NSPACE + I RELATI47 D = TA(J,I) RELATI48 J = J + 1 RELATI49 E = TA(J,I) RELATI50 I = I + 1 RELATI51 IF(J.LE.NW) J = J - 1 RELATI52 F = TA(J,I) RELATI53 J = J - 1 RELATI54 G = TA(J,I) RELATI55 I = 2*NSTA2 - 1 - MMN RELATI56 H = TA(1,I) RELATI57 U = TA(2,I) RELATI58 I = I + 1 RELATI59 Q = TA(1,I) RELATI60 QDXDX = H - 2.0D0*D + A RELATI61 QDYDY = Q - 2.0D0*F + C RELATI62 QDXDY = U + B - E - G RELATI63 A = DSQRT((QDXDX + QDYDY)/2.0D0 + DSQRT((QDYDY - QDXDX)**2/4.0D0 +RELATI64 # QDXDY**2)) RELATI65 B = DSQRT((QDXDX + QDYDY)/2.0D0 - DSQRT((QDYDY - QDXDX)**2/4.0D0 +RELATI66 # QDXDY**2)) RELATI67 C = -2.0D0*QDXDY RELATI68 D = QDYDY - QDXDX RELATI69 C = -DATAN2(C,D)*57.2957795D0/2.0D0 RELATI70 A = A*CFACT RELATI71 B = B*CFACT RELATI72 PRINT 6,NAME(NSTA1),NAME(NSTA2),A,B,C RELATI73 6 FORMAT(I17,I13,F20.3,F25.3,F25.2,/) RELATI74 K = K + 1 RELATI75 IF((K/25*25).EQ.K) PRINT 1 RELATI76 IF((K/25*25).EQ.K) PRINT 2 RELATI77 IF(K.LT.NNNN) GO TO 4 RELATI78 RETURN RELATI79 END RELATI80