C Program for quadralateral adjustments
C Smith & Varnes, U.S. Geological Survey
C Modified by K.M. Cruikshank, Portland State University
	
     
     
      IMPLICIT REAL*8 (A-H, O-Z, L)
      CHARACTER*1 BEL
      CHARACTER*13 FILEIN, FILEOUT
      CHARACTER*8 LABL(14)
      REAL*8 N(5,5), NINV(5,5)
      DIMENSION O(14), B(14,5), BT(5,14), C(14), F(14), T(5), D(5),
     1          WORK(100), A(8), DIFF(14), DSEC(8), W(14,14),
     2          BTW(5,14)
      BEL = CHAR(7)
      PI = 3.141592653589793D00
      RAD = 180./PI
      TOL = 0.0001
      ISTOP = 0
      MAXIT = 10
      DATA LABL/'ALPHA(1)','ALPHA(2)','ALPHA(3)','ALPHA(4)','ALPHA(5)',
     1          'ALPHA(6)','ALPHA(7)','ALPHA(8)','LENG 1-2','LENG 1-3',
     2          'LENG 1-4','LENG 2-3','LENG 3-4','LENG 2-4'/
      TYPE 498, BEL
  498 FORMAT (' ENTER NAME OF INPUT FILE', A1)
      ACCEPT 497, FILEIN
  497 FORMAT (A)
      TYPE 499, BEL
  499 FORMAT (' ENTER NAME OF OUTPUT FILE', A1)
      ACCEPT 497, FILEOUT
      TYPE 496, BEL
  496 FORMAT (' ENTER WEIGHT FACTOR', A1)
      ACCEPT *, WFACT
      OPEN (5, FILE=FILEIN, STATUS='OLD', FORM='FORMATTED')
      OPEN (6, FILE=FILEOUT, STATUS='NEW', FORM='FORMATTED')
C *** READ OBSERVED ANGLES AND LENGTHS
      READ (5, *) (O(I), I=1,14)
C *** PRINT INITIAL DATA
      WRITE (6, 500) (LABL(I), I=1,8)
      WRITE (6, 501) (O(I), I=1,8)
      WRITE (6, 502) (LABL(I), I=9,14), (O(I), I=9,14)

C *** STORE OBSERVED ANGLES IN DECIMAL DEGREES
      DO 50 I=1,8
		A(I) = O(I)
   50 CONTINUE

C *** CONVERT OBSERVED ANGLES TO RADIANS
      DO 100 I=1,8
		O(I) = O(I)/RAD
  100 CONTINUE

C *** INITIALIZE [B] MATRIX
      DO 150 I=1,14
		DO 150 J=1,5
		B(I,J) = 0.0
  150 CONTINUE

C *** COMPUTE APPROXIMATE COORDINATES FOR FIRST ITERATION
      X2 = O(9) * DCOS(O(1) + O(2))
      Y2 = O(9) * DSIN(O(1) + O(2))
      X3 = O(10) * DCOS(O(1))
      Y3 = O(10) * DSIN(O(1))
      X4 = O(11)

C *** COMPUTE WEIGHT MATRIX
      DO 170 I=1,14
		DO 170 J=1,14
  			W(I,J) = 0.
  170	continue
	DO 180 I=1,8
  		W(I,I) = WFACT
  180	continue
	DO 190 I=9,14
  		W(I,I) = 1.
  190	continue 

C --------------- Start of iteration loop -------
  200 ITER = ITER + 1

C *** COMPUTE ESTIMATED LENGTHS FROM COORDINATES
      L122 = X2*X2 + Y2*Y2
      L132 = X3*X3 + Y3*Y3
      L142 = X4*X4
      L232 = (X3 - X2)**2 + (Y3 -Y2)**2
      L342 = (X4 - X3)**2 + Y3*Y3
      L242 = (X4 -X2)**2 + Y2*Y2
      L12 = DSQRT(L122)
      L13 = DSQRT(L132)
      L14 = X4
      L23 = DSQRT(L232)
      L34 = DSQRT(L342)
      L24 = DSQRT(L242)

C *** COMPUTE [B] MATRIX
      B(1,3) = Y3/L132
      B(1,4) = -X3/L132
      B(2,1) = Y2/L122
      B(2,2) = -X2/L122
      B(2,3) = -B(1,3)
      B(2,4) = -B(1,4)
      B(3,1) = -Y2/L122 + Y2/L242
      B(3,2) = X2/L122 + (X4-X2)/L242
      B(3,5) = -Y2/L242
      B(4,1) = -Y2/L242 + (Y2-Y3)/L232
      B(4,2) = (X2-X4)/L242 + (X3-X2)/L232
      B(4,3) = (Y3-Y2)/L232
      B(4,4) = (X2-X3)/L232
      B(4,5) = -B(3,5)
      B(5,1) = B(4,3)
      B(5,2) = B(4,4)
      B(5,3) = (Y2-Y3)/L232 + Y3/L132
      B(5,4) = (X3-X2)/L232 -X3/L132
      B(6,3) = Y3/L342 - Y3/L132
      B(6,4) = X3/L132 + (X4-X3)/L342
      B(6,5) = -Y3/L342
      B(7,1) = -B(3,5)
      B(7,2) = (X4-X2)/L242
      B(7,3) = B(6,5)
      B(7,4) = (X3-X4)/L342
      B(7,5) = Y3/L342 - Y2/L242
      B(8,1) = B(3,5)
      B(8,2) = (X2-X4)/L242
      B(8,5) = -B(3,5)
      B(9,1) = -X2/L12
      B(9,2) = -Y2/L12
      B(10,3) = -X3/L13
      B(10,4) = -Y3/L13
      B(11,5) = -1.0
      B(12,1) = (X3-X2)/L23
      B(12,2) = (Y3-Y2)/L23
      B(12,3) = -B(12,1)
      B(12,4) = -B(12,2)
      B(13,3) = (X4-X3)/L34
      B(13,4) = -Y3/L34
      B(13,5) = -B(13,3)
      B(14,1) = (X4-X2)/L24
      B(14,2) = -Y2/L24
      B(14,5) = -B(14,1)

C *** COMPUTE MATRIX [BT], THE TRANSPOSE OF [B]
      DO 250 I=1,5
		DO 250 J=1,14
			BT(I,J) = B(J,I)
  250 CONTINUE

C *** COMPUTE THE [C] MATRIX, THE ANGLES AND LENGTHS FROM COORDINATES
      C(1) = DACOS((L132+L142-L342)/(2.*L13*L14))
      C(2) = DACOS((L122+L132-L232)/(2.*L12*L13))
      C(3) = DACOS((L122+L242-L142)/(2.*L12*L24))
      C(4) = DACOS((L232+L242-L342)/(2.*L23*L24))
      C(5) = DACOS((L132+L232-L122)/(2.*L13*L23))
      C(6) = DACOS((L132+L342-L142)/(2.*L13*L34))
      C(7) = DACOS((L242+L342-L232)/(2.*L24*L34))
      C(8) = DACOS((L142+L242-L122)/(2.*L14*L24))
      C(9) = L12
      C(10) = L13
      C(11) = L14
      C(12) = L23
      C(13) = L34
      C(14) = L24

C *** COMPUTE THE [F] MATRIX = [C] - [O]
      DO 300 I=1,14
		F(I) = C(I) - O(I)
  300 CONTINUE

C *** COMPUTE [N] AND [T] MATRICES USING IMSL SUBROUTINE VMULFF
      CALL VMULFF (BT, W, 5, 14, 14, 5, 14, BTW, 5, IER0)
      CALL VMULFF (BTW, B, 5, 14, 5, 5, 14, N, 5, IER1)
      CALL VMULFF (BTW, F, 5, 14, 1, 5, 14, T, 5, IER2)

C *** COMPUTE INVERSE OF [N] MATRIX, USING IMSL SUBROUTINE LINV2F
      CALL LINV1F (N, 5, 5, NINV, 0, WORK, IER3)

C *** COMPUTE DELTA MATRIX [D] USING VMULFF
      CALL VMULFF(NINV, T, 5, 5, 1, 5, 5, D, 5, IER4)

C *** TEST IF CORRECTIONS ARE LESS THAN SPECIFIED TOLERANCE
      IF ((DMAX1(DABS(D(1)),DABS(D(2)),DABS(D(3)),DABS(D(4)),
     1           DABS(D(5)))).LT.TOL) ISTOP = 1

C *** COMPUTE NEW COORDINATES
      X2A = X2 + D(1)
      Y2A = Y2 + D(2)
      X3A = X3 + D(3)
      Y3A = Y3 + D(4)
      X4A = X4 + D(5)

C *** WRITE CALCULATIONS FOR EACH ITERATION
      WRITE (6, 503) ITER
      WRITE (6, 504)
      WRITE (6, 505) X2, D(1), X2A
      WRITE (6, 506) Y2, D(2), Y2A
      WRITE (6, 507) X3, D(3), X3A
      WRITE (6, 508) Y3, D(4), Y3A
      WRITE (6, 509) X4, D(5), X4A
      WRITE (6, 514) IER0, IER1, IER2, IER3, IER4

C *** Test to see if we need to go for another iteration
      IF (ISTOP.EQ.1) GO TO 350
      IF (ITER.EQ.MAXIT) GO TO 350
      X2 = X2A
      Y2 = Y2A
      X3 = X3A
      Y3 = Y3A
      X4 = X4A
      GO TO 200
C -------------- End of iteration loop ------ 

C *** IF FINAL ITERATION, PRINT ANGLES, LENGTHS, AND DIFFERENCES
C *** CONVERT ADJUSTED ANGLES TO DECIMAL DEGREES
  350 DO 360 I=1,8
		C(I) = C(I)*RAD
  360 CONTINUE
C *** COMPUTE DIFFERENCES
      DO 370 I=1,8
		DIFF(I) = C(I) - A(I)
		DSEC (I) = DIFF(I)*3600.
  370 CONTINUE
      DO 380 I=9,14
		DIFF(I) = C(I) - O(I)
  380 CONTINUE
C *** COMPUTE SUMS OF SQUARES OF ANGLES AND LENGTHS
      SSQA = 0.
      SSQL = 0.
      DO 390 I=1,8
		SSQA = SSQA + DSEC(I)*DSEC(I)
  390 continue
        DO 400 I=9,14
		SSQL = SSQL + DIFF(I)*DIFF(I)*(1.D+06)
  400 continue
C *** PRINT FINAL VALUES
      WRITE (6, 510)
      WRITE (6,511) (LABL(I), C(I), A(I), DIFF(I), DSEC(I), I=1,8)
      WRITE (6, 512) 
      WRITE (6, 513) (LABL(I), C(I), O(I), DIFF(I), I=9,14)
      WRITE (6, 515) WFACT, SSQA, SSQL
      CLOSE (5)
      CLOSE (6)
      STOP
  500 FORMAT('1ADJUSTMENT OF QUADRILATERAL USING TRIANGULATION AND TRILA
     1TERATION',/, ' (METHOD OF INDIRECT OBSERVATIONS)', //,
     2' INPUT DATA (OBSERVED VALUES):',//, 8(2X, A8))
  501 FORMAT (8F10.5)
  502 FORMAT (/,6(2X, A8),/,6F10.4)
  503 FORMAT ('0ITERATION NUMBER ', I2)
  504 FORMAT (6X, 'APPROXIMATE', 15X, 'ADJUSTED',/,
     1        6X, 'COORDINATES  CORRECTION  COORDINATES')
  505 FORMAT ('   X2', F12.4, F12.5, F13.4)
  506 FORMAT ('   Y2', F12.4, F12.5, F13.4)
  507 FORMAT ('   X3', F12.4, F12.5, F13.4)
  508 FORMAT ('   Y3', F12.4, F12.5, F13.4)
  509 FORMAT ('   X4', F12.4, F12.5, F13.4)
  510 FORMAT (/, 6X, 'ADJUSTED LENGTHS AND ANGLES', /,
     1        14X, 'ADJUSTED   OBSERVED   DIFF.      "')
  511 FORMAT (3X, A8, 2F11.5, F9.5, F6.1)
  512 FORMAT (' ')
  513 FORMAT (3X, A8, 2F11.4, F9.4)
  514 FORMAT ('0ERROR CODES:', 5I5)
  515 FORMAT ('0WEIGHT FACTOR = ', 1PE10.2, /,
     1        ' SSQ(ANGLES) = ', 0PF8.2, ' SEC**2', /,
     2        ' SSQ(LENGTHS) = ', F7.2, '  MM**2')
      END



