SUBROUTINE AFTRAN(COVECT, XARG, YARG, NEWX, NEWY)
C**********************************************************************
C* *
C* THIS ROUTINE COMPUTES TRANSFORMED COORDINATES ACCORDING TO THE *
C* AFFINE MODEL. THE VECTOR COVECT CONTAINS THE AFFINE COEFFIC- *
C* IENTS, X AND Y ARE COORDINATES IN THE OLD SYSTEM. *
C* *
C* WRITTEN BY CLIFF PETERSOHN, OCT., 1981 *
C* *
C**********************************************************************
C
DOUBLE PRECISION COVECT(6), X, Y, NEWX, NEWY
REAL XARG,YARG
C
C
X=DBLE(XARG)
Y=DBLE(YARG)
NEWX = COVECT(1) * X + COVECT(2) * Y + COVECT(3)
NEWY = COVECT(4) * X + COVECT(5) * Y + COVECT(6)
RETURN
END
SUBROUTINE AFFINE(DEVICE, IDSIZE, NPASS, PLANE, TABLET, DTABLT,
@ W, V, UNITS)
C**********************************************************************
C* *
C* THIS ROUTINE PERFORMS A LEAST-SQUARES AFFINE TRANSFORMATION *
C* OF COORDINATES FROM THE ARRAY (TABLET) TO THE ARRAY (PLANE) *
C* INPUTS ARE : *
C* - IDSIZE: THE DECLARED DIMENSION OF THE ARRAYS *
C* (PLANE) AND (TABLET) *
C* *
C* - NPASS : THE NUMBER OF COORDINATES IN THE *
C* ARRAYS (PLANE) AND (TABLET) *
C* *
C* - PLANE, *
C* TABLET: THE (X, Y) COORDINATES IN BOTH SYSTEMS *
C* *
C* DTABLT: A SERVICE ARRAY OF DIMENSION (IDSIZE). *
C* OUTPUTS ARE : *
C* - W : THE COEFFICIENT VECTOR *
C* *
C* - V : RESIDUAL VECTOR *
C* *
C* *
C* WRITTEN BY CLIFF PETERSOHN, OCT., 1981 *
C* *
C**********************************************************************
C
DOUBLE PRECISION PLANE(IDSIZE), DTABLT(IDSIZE), V(IDSIZE),
@ EN(21), W(6), X, Y, MK, ML, DETEN, DPASS,
@ DSQRT, DBLE
C
REAL TABLET(IDSIZE)
INTEGER XIND, YIND, LIMIT, KOUNT
INTEGER UNITS, DEVICE, LUNIT
C
C
C INITIALIZE SOME STUFF
C
LUNIT = 6
C IF(DEVICE.EQ.'AED512')LUNIT=9
DO 1 XIND = 1, 21
EN(XIND) = 0.0D0
1 CONTINUE
C
DO 2 XIND = 1, 6
W(XIND) = 0.0D0
2 CONTINUE
C
DO 3 XIND = 1, NPASS
DTABLT(XIND) = DBLE(TABLET(XIND))
3 CONTINUE
C
MK = 0.0D0
ML = 0.0D0
DPASS = DBLE(FLOAT(NPASS / 2))
LIMIT = NPASS - 1
WRITE(LUNIT, 999)
WRITE(LUNIT, 1000) DPASS
WRITE(LUNIT, 1001)
IF (UNITS .EQ. 70) WRITE(LUNIT, 1002)
IF (UNITS .EQ. 77) WRITE(LUNIT, 1003)
C
C FORM NORMAL EQUATIONS
C
DO 4 XIND = 1, LIMIT, 2
YIND = XIND + 1
EN(1) = EN(1) + DTABLT(XIND)**2
EN(2) = EN(2) + DTABLT(XIND) * DTABLT(YIND)
EN(3) = EN(3) + DTABLT(YIND)**2
EN(4) = EN(4) + DTABLT(XIND)
EN(5) = EN(5) + DTABLT(YIND)
EN(6) = EN(6) + 1.0D0
C
C ALSO FORM CONSTANT VECTOR IN THIS LOOP
C
W(1) = W(1) + DTABLT(XIND) * PLANE(XIND)
W(2) = W(2) + DTABLT(YIND) * PLANE(XIND)
W(3) = W(3) + PLANE(XIND)
W(4) = W(4) + DTABLT(XIND) * PLANE(YIND)
W(5) = W(5) + DTABLT(YIND) * PLANE(YIND)
W(6) = W(6) + PLANE(YIND)
4 CONTINUE
C
EN(10) = EN(1)
EN(14) = EN(2)
EN(15) = EN(3)
EN(19) = EN(4)
EN(20) = EN(5)
EN(21) = EN(6)
C
C SOLVE THE SYSTEM
C
CALL CHLSKY(EN, 21, W, DETEN, 6, 1)
C
C FORM RESIDUALS
C
KOUNT = 0
DO 5 XIND = 1, LIMIT, 2
YIND = XIND + 1
KOUNT = KOUNT + 1
XTAB= TABLET(XIND)
YTAB= TABLET(YIND)
CALL AFTRAN(W, XTAB, YTAB, X, Y)
V(XIND) = PLANE(XIND) - X
V(YIND) = PLANE(YIND) - Y
DTABLT (XIND) = X
DTABLT (YIND) = Y
MK = MK + V(XIND)**2 + V(YIND)**2
WRITE(LUNIT, 1004) KOUNT, Y, X, V(YIND), V(XIND)
5 CONTINUE
C
C COMPUTE MEAN STANDARD COORDINATE AND POSITIONAL POINT DEVIATION
C
MK = DSQRT(MK / (2.0D0 * DPASS - 6.0D0))
ML = MK * 1.4142135D0
WRITE(LUNIT, 1005) MK, ML
WRITE(LUNIT, 1006) (W(XIND), XIND = 1, 6)
C
999 FORMAT('-', 23X, 'AFFINE TRANSFORMATION')
1000 FORMAT(' ', 10X, 'ADJUSTMENT RESULTS', 10X, F4.0,
@ ' CONTROL POINTS USED')
1001 FORMAT('-', 3X, 'POINT', 5X, 'NORTHING', 10X, 'EASTING',
@ 12X, 'V(N)', 10X, 'V(E)' )
1002 FORMAT(' ', 13X, '(FEET)', 12X, '(FEET)', 9X, '(FEET)',
@ 8X, '(FEET)' )
1003 FORMAT(' ', 13X, '(METERS)', 10X, '(METERS)', 9X, '(METERS)',
@ 6X, '(METERS)' )
1004 FORMAT(' ', 4X, I2, 4X, F12.3, 4X, F12.3, 2F16.5)
1005 FORMAT('-', 2X, 'MEAN COORDINATE STANDARD DEVIATION', F10.2,
@ /, 3X, 'MEAN POSITIONAL STANDARD DEVIATION', F10.2)
1006 FORMAT('-', 'COEFFICIENTS ARE :', 3(1X, F15.5)/19X,3(1X, F15.5))
C
RETURN
END
SUBROUTINE CHLSKY(A, IRDA, B, DETA, N, ICODE)
C**********************************************************************
C* *
C* THIS ROUTINE SOLVES A SYSTEM OF NORMAL EQUATIONS GIVEN *
C* *
C* -A : THE LOWER TRIANGLE OF THE COEFFICIENT MATRIX *
C* STORED IN A VECTOR SUCH THAT A(I, J) OF THE *
C* COEFFICIENT MATRIX RESIDES IN A((I**2 - I)/2 *
C* + J). *
C* *
C* -B : THE CONSTANT VECTOR. *
C* *
C* MATRIX A IS ASSUMED POSITIVE DEFINITE SYMMETRIC, WHICH THIS *
C* ROUTINE TAKES FULL ADVANTAGE OF IN SOLVING THE SYSTEM BY THE *
C* METHOD OF CHOLESKY. OTHER INPUT QUANTITIES ARE: *
C* *
C* -IRDA : THE DECLARED SIZE OF VECTOR A *
C* *
C* -N : THE ACTUAL MATRIX DIMENSION OF A *
C* *
C* -ICODE : AN INTEGER FLAG OF THE FOLLOWING VALUES: *
C* *
C* 1 : IF ICODE = 1, THE SYSTEM IS SOLVED AND *
C* THE SOLUTION VECTOR IS RETURNED IN B . A *
C* CONTAINS THE LOWER TRIANGULAR DECOMPO- *
C* SITION OF A . *
C* *
C* 0 : IF ICODE = 0, B IS IGNORED AND A IS *
C* RETURNED AS THE INVERSE OF A . *
C* *
C* -1 : IF ICODE = -1, B IS IGNORED AND A IS *
C* THE LOWER TRIANGULAR DECOMPOSITION *
C* OF A . A IS RETURNED AS THE INVERSE OF *
C* ITSELF. *
C* *
C* WRITTEN BY J. A. MILLER, 1977 *
C* *
C* MODIFIED BY C. R. PETERSOHN, OCT.1981*
C* (MATRICIES REMOVED) *
C* LAST UPDATE: OCTOBER 14, 1981 *
C* *
C**********************************************************************
C
DOUBLE PRECISION A (IRDA), B (N), DETA, SUM1, SUM2,
@ DSQRT, DLOG10
C
INTEGER I, J, K, IB, KB, KC, KS, JS, NA
C
IDXX(I, J) = (I**2 - I) / 2 + J
C
C SEE IF ALL WE DO IS INVERT.
C
IF (ICODE. LT. 0) GO TO 99
C
C OTHERWISE FORM THE LOWER TRIANGULAR DECOMPOSITION OF A.
C
A(1) = DSQRT(A(1))
DETA = DLOG10(A(1))
DO 1 J = 2, N
A(IDXX(J, 1)) = A(IDXX(J, 1)) / A(1)
1 CONTINUE
DO 2 I = 2, N
SUM1 = 0.0D0
NK = I - 1
DO 3 K = 1, NK
SUM1 = SUM1 + A(IDXX(I, K))**2
3 CONTINUE
A(IDXX(I, I)) = DSQRT(A(IDXX(I, I)) - SUM1)
DETA = DETA + DLOG10(A(IDXX(I, I)))
IF (I. EQ. N) GO TO 2
JS = I + 1
DO 4 J = JS, N
SUM2 = 0.0D0
DO 5 K = 1, NK
SUM2 = SUM2 + A(IDXX(I, K)) * A(IDXX(J, K))
5 CONTINUE
A(IDXX(J, I)) = (A(IDXX(J, I)) - SUM2) /
@ A(IDXX(I, I))
4 CONTINUE
2 CONTINUE
C
C ANTILOG BASE TEN IS THE DETERMINANT
C
DETA = 10.0D0**DETA
C
C INVERT THE MATRIX IF NO SOLUTION IS DESIRED.
C
IF(ICODE .EQ. 0) GO TO 99
C
C OTHERWISE, PERFORM A FORWARD SUBSTITUTION.
C
NA = N - 1
DO 6 I = 1, NA
B(I) = B(I) / A(IDXX(I, I))
DO 7 J = 1, I
IB = I + 1
B(IB) = B(IB) - A(IDXX(IB, J)) * B(J)
7 CONTINUE
6 CONTINUE
B(N) = B(N) / A(IDXX(N, N))
C
C BACK SUBSTITUTION
C
DO 8 I = 1, NA
KB = N + 1 - I
B(KB) = B(KB) / A(IDXX(KB, KB))
DO 9 J = 1, I
KC = N + 1 - J
IB = KB - 1
B(IB) = B(IB) - A(IDXX(KC, IB)) * B(KC)
9 CONTINUE
8 CONTINUE
B(1) = B(1) / A(1)
C
C SYSTEM SOLVED. QUIT WITH RESULTS IN A AND B.
C
RETURN
C
C THE FOLLOWING CODE INVERTS A DECOMPOSED LOWER TRIANGULAR MATRIX
C
99 CONTINUE
DO 10 I = 1, N
A(IDXX(I, I)) = 1.0D0 / A(IDXX(I, I))
10 CONTINUE
NA = N - 1
DO 11 J = 1, NA
IB = J + 1
DO 12 I = IB, N
SUM1 = 0.0D0
KS = I - 1
DO 13 K = J, KS
SUM1 = SUM1 + A(IDXX(I, K)) *
@ A(IDXX(K,J))
13 CONTINUE
A(IDXX(I, J)) = -SUM1 * A(IDXX(I, I))
12 CONTINUE
11 CONTINUE
C
C CONSTRUCT THE INVERSE OF A
C
DO 14 I = 1, N
DO 15 J = I, N
SUM1 = 0.0D0
DO 16 K = J, N
SUM1 = SUM1 + A(IDXX(K, I)) / A(IDXX(K, J))
16 CONTINUE
A(IDXX(I, J)) = SUM1
15 CONTINUE
14 CONTINUE
RETURN
END
SUBROUTINE PRTRAN(COVECT, IX, IY, NEWX, NEWY)
C***********************************************************************
C* *
C* THIS ROUTINE COMPUTES TRANSFORMED COORDINATES OF POINTS ACCORDING *
C* TO THE PROJECTIVE MODEL, GIVEN: *
C* *
C* - COVECT : A VECTOR OF PROJECTIVE COEFFICIENTS. *
C* - X, Y : A COORDINATE PAIR IN THE OLD SYSTEM. *
C* *
C* OUTPUT IS THE (NEWX, NEWY) COORDINATE PAIR IN THE TRANSFORMED *
C* SYSTEM (IN DOUBLE PRECISION) *
C* *
C* -WRITTEN BY CLIFF PETERSOHN, OCT. 1981. *
C* *
C***********************************************************************
C
DOUBLE PRECISION COVECT(8), X, Y, NEWX, NEWY, DENOM
C
REAL IX, IY
C
C
X = DBLE(IX)
Y = DBLE(IY)
DENOM = COVECT(7) * X + COVECT(8) * Y + 1.0D0
NEWX = (COVECT(1) * X + COVECT(2) * Y + COVECT(3)) / DENOM
NEWY = (COVECT(4) * X + COVECT(5) * Y + COVECT(6)) / DENOM
NEWX = NEWX * 1.0D3
NEWY = NEWY * 1.0D3
RETURN
END
SUBROUTINE PROJEC(DEVICE, IDSIZE, NPASS, PLANE, SPLANE, TABLET,
@ DTABLT, EPS, DELTA, V, UNITS)
C**********************************************************************
C* *
C* THIS ROUTINE PERFORMS A LEAST-SQUARES PROJECTIVE TRANSFORMATION *
C* OF COORDINATES FROM THE ARRAY (TABLET) TO THE ARRAY (PLANE) *
C* INPUTS ARE : *
C* - IDSIZE: THE DECLARED DIMENSION OF THE ARRAYS *
C* (PLANE) AND (TABLET) *
C* *
C* - NPASS : THE NUMBER OF COORDINATES IN THE *
C* ARRAYS (PLANE) AND (TABLET) *
C* *
C* - PLANE, *
C* TABLET: THE (X, Y) COORDINATES IN BOTH SYSTEMS *
C* *
C* - EPS : THE CONVERGENCE CRITERION. *
C* *
C* OUTPUTS ARE : *
C* - DELTA : THE COEFFICIENT VECTOR *
C* *
C* - V : RESIDUAL VECTOR *
C* *
C* *
C* WRITTEN BY CLIFF PETERSOHN, OCT., 1981 *
C* *
C**********************************************************************
C
C
DOUBLE PRECISION PLANE(IDSIZE), SPLANE(IDSIZE), DTABLT(IDSIZE),
@ V(IDSIZE),
@ EN(36), W(8), DELTA(8), X, Y, MK, ML, XY, X2,
@ Y2, XNUMI, YNUMI, XYNUM2, DENOMI, DENOM2,
@ DENOM3, DENOM4, DPASS, DETEN, EPS,
@ DABS, DBLE, DSQRT, PDENOM, PXNUM, PYNUM
C
REAL TABLET(IDSIZE)
INTEGER XIND, YIND, LIMIT, KOUNT, IKOUNT
C
INTEGER UNITS, DEVICE, LUNIT
C
C PDENOM(C1, C2, X, Y) = C1 * X + C2 * Y + 1.0D0
C PXNUM (A1, A2, A3, X, Y) = A1 * X + A2 * Y + A3
C PYNUM (B1, B2, B3, X, Y) = B1 * X + B2 * Y + B3
C
C INITIALIZE SOME STUFF
C
LUNIT=6
C IF(DEVICE.EQ.'AED512')LUNIT=9
DPASS = DBLE(FLOAT(NPASS / 2))
LIMIT = NPASS - 1
IKOUNT = 0
WRITE(LUNIT, 999)
WRITE(LUNIT, 1000) DPASS
WRITE(LUNIT, 1001)
IF(UNITS .EQ. 70) WRITE(LUNIT, 1002)
IF(UNITS .EQ. 77) WRITE(LUNIT, 1003)
MK = 0.0D0
ML = 0.0D0
C
C CONVERT INTEGER TABLE COORDINATES TO DOUBLE PRECISION
C SCALE GROUND COORDINATES FOR NUMERICAL STABILITY.
C
DO 1 KOUNT = 1, NPASS
DTABLT(KOUNT) = DBLE(TABLET(KOUNT))
SPLANE(KOUNT) = PLANE(KOUNT) * 1.0D-3
1 CONTINUE
C
C APPROXIMATE VALUES FOR UNKNOWNS
C
DO 2 KOUNT = 1, 8
DELTA(KOUNT) = 0.0D0
2 CONTINUE
DELTA(3) = SPLANE(1)
DELTA(6) = SPLANE(2)
DELTA(7) = -1.0D-6
DELTA(8) = -DELTA(7)
C
C BEGIN ITERATIVE BODY. ZERO NORMAL EQUATION AND MISCLOSURE VECTORS.
C
99 IKOUNT = IKOUNT + 1
DO 3 KOUNT = 1, 36
EN(KOUNT) = 0.0D0
3 CONTINUE
DO 4 KOUNT = 1, 8
W(KOUNT) = 0.0D0
4 CONTINUE
C
C FORM NORMAL EQUATIONS
C
DO 5 XIND = 1, LIMIT, 2
YIND = XIND + 1
C
C CARRY SOME ITERATION CONSTANTS
C
C XNUMI = PXNUM(DELTA(1), DELTA(2), DELTA(3), DTABLT(XIND),
C @ DTABLT(YIND))
C YNUMI = PYNUM(DELTA(4), DELTA(5), DELTA(6), DTABLT(XIND),
C @ DTABLT(YIND))
C DENOMI = PDENOM(DELTA(7), DELTA(8), DTABLT(XIND),
C @ DTABLT(YIND))
XNUMI = DELTA(1) * DTABLT(XIND) + DELTA(2)*DTABLT(YIND) +DELTA(3)
YNUMI = DELTA(4) * DTABLT(XIND) + DELTA(5)*DTABLT(YIND) +DELTA(6)
DENOMI= DELTA(7) * DTABLT(XIND) + DELTA(8)*DTABLT(YIND) +1.0D0
XYNUM2 = XNUMI**2 + YNUMI**2
DENOM2 = DENOMI**2
DENOM3 = DENOM2 * DENOMI
DENOM4 = DENOM2**2
XY = DTABLT(XIND) * DTABLT(YIND)
X2 = DTABLT(XIND)**2
Y2 = DTABLT(YIND)**2
C
C
EN(1) = EN(1) + X2 / DENOM2
EN(2) = EN(2) + XY / DENOM2
EN(3) = EN(3) + Y2 / DENOM2
EN(4) = EN(4) + DTABLT(XIND) / DENOM2
EN(5) = EN(5) + DTABLT(YIND) / DENOM2
EN(6) = EN(6) + 1.0D0 / DENOM2
EN(22) = EN(22) - X2 * XNUMI / DENOM3
EN(23) = EN(23) - XY * XNUMI / DENOM3
EN(24) = EN(24) - DTABLT(XIND) * XNUMI / DENOM3
EN(25) = EN(25) - X2 * YNUMI / DENOM3
EN(26) = EN(26) - XY * YNUMI / DENOM3
EN(27) = EN(27) - DTABLT(XIND) * YNUMI / DENOM3
EN(28) = EN(28) + X2 * XYNUM2 / DENOM4
EN(30) = EN(30) - Y2 * XNUMI /DENOM3
EN(31) = EN(31) - DTABLT(YIND) * XNUMI / DENOM3
EN(33) = EN(33) - Y2 * YNUMI / DENOM3
EN(34) = EN(34) - DTABLT(YIND) * YNUMI / DENOM3
EN(35) = EN(35) + XY * XYNUM2 / DENOM4
EN(36) = EN(36) + Y2 * XYNUM2 / DENOM4
C
C ALSO FORM MISCLOSURE VECTOR IN THIS LOOP. X2 AND Y2 ARE MISCLOSURES
C
X2 = SPLANE(XIND) - XNUMI / DENOMI
Y2 = SPLANE(YIND) - YNUMI / DENOMI
W(1) = W(1) + DTABLT(XIND) * X2 / DENOMI
W(2) = W(2) + DTABLT(YIND) * X2 / DENOMI
W(3) = W(3) + X2 / DENOMI
W(4) = W(4) + DTABLT(XIND) * Y2 / DENOMI
W(5) = W(5) + DTABLT(YIND) * Y2 / DENOMI
W(6) = W(6) + Y2 / DENOMI
W(7) = W(7) + ((-DTABLT(XIND) * X2 * XNUMI) -
@ (DTABLT(XIND) * Y2 * YNUMI)) /
@ DENOM2
W(8) = W(8) + ((-DTABLT(YIND) * X2 * XNUMI) -
@ (DTABLT(YIND) * Y2 * YNUMI)) /
@ DENOM2
5 CONTINUE
EN(10) = EN(1)
EN(14) = EN(2)
EN(15) = EN(3)
EN(19) = EN(4)
EN(20) = EN(5)
EN(21) = EN(6)
EN(29) = EN(23)
EN(32) = EN(26)
C
C SOLVE THE SYSTEM
C
CALL CHLSKY(EN, 36, W, DETEN, 8, 1)
C
C UPDATE SOLUTION
C
DO 6 KOUNT = 1, 8
DELTA(KOUNT) = DELTA(KOUNT) + W(KOUNT)
6 CONTINUE
C
C SEE IF WE CAN QUIT THIS FOOLISHNESS
C
DO 7 KOUNT = 1, 8
IF(DABS(W(KOUNT)) .GT. EPS) GO TO 99
7 CONTINUE
C
C FINISHED. ALMOST. COMPUTE RESIDUALS, PRINT SOME STUFF AND QUIT.
C
KOUNT = 0
DO 8 XIND = 1, LIMIT, 2
YIND = XIND + 1
KOUNT = KOUNT + 1
CALL PRTRAN(DELTA, TABLET(XIND), TABLET(YIND), X, Y)
V(XIND) = PLANE(XIND) - X
V(YIND) = PLANE(YIND) - Y
MK = MK + V(XIND)**2 + V(YIND)**2
WRITE(LUNIT, 1004) KOUNT, Y, X, V(YIND), V(XIND)
8 CONTINUE
MK = DSQRT(MK / (2 * DPASS - 8.0D0))
ML = 1.4142135D0 * MK
WRITE(LUNIT, 1005) MK, ML
WRITE(LUNIT, 1006) IKOUNT, EPS
WRITE(LUNIT, 1007) (DELTA(KOUNT), KOUNT = 1, 8)
C
999 FORMAT('-', 23X, 'PROJECTIVE TRANSFORMATION')
1000 FORMAT(' ', 10X, 'ADJUSTMENT RESULTS', 10X, F4.0,
@ ' CONTROL POINTS USED')
1001 FORMAT('-', 3X, 'POINT', 5X, 'NORTHING', 10X, 'EASTING',
@ 12X, 'V(N)', 10X, 'V(E)' )
1002 FORMAT(' ', 13X, '(FEET)', 12X, '(FEET)', 9X, '(FEET)',
@ 8X, '(FEET)' )
1003 FORMAT(' ', 13X, '(METERS)', 10X, '(METERS)', 9X, '(METERS)',
@ 6X, '(METERS)' )
1004 FORMAT(' ', 4X, I2, 4X, F12.3, 6X, F12.3, 2(4X, F10.5))
1005 FORMAT('-', 2X, 'MEAN COORDINATE STANDARD DEVIATION', 2X, F8.4,
@ /, 3X, 'MEAN POSITIONAL STANDARD DEVIATION', 2X, F8.4)
1006 FORMAT('-', 2X, I1, ' ITERATIONS REQUIRED TO CONVERGE TO', E15.5)
1007 FORMAT('-', 'COEFFICIENTS ARE :', 4(1X, F10.5)/19X,4(1X, F10.5))
C
RETURN
END
SUBROUTINE SPREAD(IDSIZE, NPASS, PASS, V, X, Y,
@ XHAT, YHAT, SPRSW , WRESID)
C**********************************************************************
C* *
C* THIS ROUTINE ACCEPTS AN (X, Y) COORDINATE PAIR AND SOME *
C* INFORMATION ABOUT THE CONTROL POINTS SURROUNDING THEM. *
C* THIS INFORMATION INCLUDES: *
C* *
C* -PASS : A VECTOR OF REAL*8 COORDINATE *
C* PAIRS OF CONTROL SURROUNDING *
C* POINTS (X, Y) *
C* *
C* -V : A VECTOR OF RESIDUALS CORRES- *
C* PONDING TO THE VECTOR PASS *
C* AND RESULTING FROM A TRANSFOR- *
C* MATION OF THE (X, Y) SYSTEM *
C* TO THE PASS SYSTEM *
C* *
C* -X, Y : A COORDINATE PAIR TRANSFORMED *
C* INTO THE PASS SYSTEM *
C* - WRESID : control Xformed, w. residue *
C* THE ROUTINE DISTRIBUTES, OR SPREADS A FRACTIONAL AMOUNT OF *
C* EACH RESIDUAL INVERSELY PROPORTIONAL TO THE GIVEN POINT'S *
C* DISTANCE FROM THE PASS POINT. *
C* *
C* OUTPUT IS THE SPREADED COORDINATE (IN SINGLE PRECISION) *
C* *
C* WRITTEN BY CLIFF PETERSOHN, OCT., 1981 *
C* LAST UPDATE: GREG MILLS. DATE: 28 FEB 82. *
C* changed 1 October 1985 to treat points near control differently *
C* by Nick Chrisman under direction of AP Vonderohe *
C**********************************************************************
C
DOUBLE PRECISION PASS(IDSIZE), V(IDSIZE), WRESID(IDSIZE), SIGWT,
@ WEIGHT, VXSUM, VYSUM, DSQRT , X, Y
C
REAL XHAT, YHAT
C
INTEGER NPTS, LIMIT, XIND, YIND , NPASS
C
LOGICAL SPRSW
C
C INITIALIZE SOME STUFF.
C
NPTS = NPASS / 2
LIMIT = NPASS - 1
SIGWT = 0.0D0
VXSUM = 0.0D0
VYSUM = 0.0D0
C
IF(.NOT.SPRSW)GO TO 100
C
C FOR EACH PASS POINT, COMPUTE ITS DISTANCE FROM (X, Y) AND
C WEIGHT THE CONTRIBUTION OF ITS RESIDUAL INVERSELY PROPORTIONAL.
C
DO 1 XIND = 1, LIMIT, 2
YIND = XIND + 1
IF ( DSQRT ((X- WRESID(XIND)) **2 +
X (Y- WRESID(YIND)) **2 ) .LT. 1D0 ) GO TO 400
WEIGHT = (DSQRT((X - PASS(XIND)) ** 2 +
@ (Y - PASS(YIND)) ** 2 ))
C - THE OLD BUG - IF(WEIGHT. LT. 1.0D0) WEIGHT = 1.0D0
IF (WEIGHT .LT. 1D-9) WEIGHT = 1D-9
WEIGHT = 1.0D0 / WEIGHT
SIGWT = SIGWT + WEIGHT
VXSUM = VXSUM + V(XIND) * WEIGHT
VYSUM = VYSUM + V(YIND) * WEIGHT
1 CONTINUE
C
VXSUM = VXSUM / SIGWT
VYSUM = VYSUM / SIGWT
100 CONTINUE
XHAT = SNGL(X + VXSUM)
YHAT = SNGL(Y + VYSUM)
C
RETURN
C CASE OF NEARLY EXACT MATCH TO A CONTROL POINT- GIVE TRUE LOC.
400 XHAT= SNGL (PASS(XIND) )
YHAT= SNGL (PASS(YIND))
RETURN
END