
CD-Section: Make Distortion Corrections to Coordinates
C A routine is avaiable to make the appropriate distortion corrections
C to an ideal predicted spot position. The routine takes account of the
C current distortion type and the various distortion parameters.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**********************
C**     LDM_CORR     **
C**********************
C
CD-Routine: Make Distortion Corrections to Ideal Coordinates - LDM_CORR
C This routine will make the appropriate distortion corrections
C to an ideal predicted spot position. The routine takes account of the
C current distortion type and the various distortion parameters. cos(w_c)
C and sin(w_c) are passed as parameters to the routine so that, in
C suitable circumstances, they need not be re-calculated each time that
C LDM_CORR is called.
CD-end
CD-Fortran:
      SUBROUTINE LDM_CORR (IPACK, IPLATE, COS_WC, SIN_WC, CF, XF, YF, 
     +                     XFD, YFD)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Apply distortion corrections to xf, yf spot coordinates
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER IPACK
      INTEGER IPLATE
      REAL COS_WC
      REAL SIN_WC
      REAL CF
      REAL XF
      REAL YF
      REAL XFD
      REAL YFD

CD-Parameters:
C IPACK     i (R)   Pack number 
C IPLATE    i (R)   Plate number
C COS_WC    r (R)   cos(w_c) for the pack and plate
C SIN_WC    r (R)   sin(w_c) for the pack and plate
C CF        r (R)   crystal-film distance for the pack and plate
C XF        r (R)   Predicted reflection 'xf' coordinate in mm with 
C                   respect to the predicted pattern centre.
C YF        r (R)   Predicted reflection 'yf' coordinate in mm with 
C                   respect to the predicted pattern centre.
C XFD       r (W)   Distortion corrected 'xfd' parameter in mm with 
C                   respect to the predicted centre. The coordinate is 
C                   now along the direction of the actual image 'xd' axis
C YFD       r (W)   Distortion corrected 'yfd' parameter in mm with 
C                   respect to the predicted centre. The coordinate is 
C                   now along the direction of the actual image 'yd' axis
CD-end
C
C
C====== GLOBALS:
C
       INCLUDE 'ldm_data.finc'
C
C====== LOCALS:
C
      INTEGER I
      INTEGER NSPDX
      INTEGER NSPDY
      INTEGER IERR
      REAL SPD_XY
      REAL SPDX_MIN
      REAL SPDX_MAX
      REAL SPDY_MIN
      REAL SPDY_MAX
      REAL SXVAL(LDM_MAXSPD)
      REAL SYVAL(LDM_MAXSPD)
      REAL K_FACTOR
      REAL RQ
      REAL XQ
      REAL YQ
      REAL CPSI
      REAL SPSI
      REAL TORAD
      REAL FACT
      REAL XD
      REAL XP
      REAL XP2
      REAL YD
      REAL YP
      REAL YP2
      REAL D
      REAL DD
      REAL SV
      CHARACTER*2 ERRSTR

C     INTEGER I         !Loop variable
C     INTEGER NSPDX     !Spatial distortion x polynomial order
C     INTEGER NSPDY     !Spatial distortion y polynomial order
C     INTEGER IERR      !Error flag
C     REAL SPD_XY       !Spatial distortion xy cross term
C     REAL SPDX_MIN     !Spatial distortion Chebyshev x lower bound
C     REAL SPDX_MAX     !Spatial distortion Chebyshev x upper bound
C     REAL SPDY_MIN     !Spatial distortion Chebyshev y lower bound
C     REAL SPDY_MAX     !Spatial distortion Chebyshev y upper bound
C     REAL SXVAL(LDM_MAXSPD)      !Spatial distortion x polynomial coeffts
C     REAL SYVAL(LDM_MAXSPD)      !Spatial distortion y polynomial coeffts
C     REAL K_FACTOR     !Distortion correction scale factor
C     REAL RQ           !Spot distance from centre of normal flat film (mm)
C     REAL XQ           !X coordinate of spot corrected for omega (mm)
C     REAL YQ           !Y coordinate of spot corrected for omega (mm)
C     REAL CPSI         !cos(psi) for roff/toff correction
C     REAL SPSI         !sin(psi) for roff/toff correction
C     REAL TORAD        !Conversion factor from 1/100deg to radians
C     REAL FACT         !Scaling factor TORAD/CF
C     REAL XD           !x-offset from type 3 distortion correction centre
C     REAL XP           !Normalised coordinate
C     REAL XP2          !2*XP
C     REAL YD           !y-offset from type 3 distortion correction centre
C     REAL YP           !Normalised coordinate
C     REAL YP2          !2*YP
C     REAL D            !Temp for chebyshev evaluation
C     REAL DD           !Temp for chebyshev evaluation
C     REAL SV           !Temp for chebyshev evaluation
C     CHARACTER*2 ERRSTR !Error string (dummy)
C
      DATA TORAD/1.745329E-04/
C
C-------------------------------------------------------------------------------
C
      FACT = TORAD/CF
      XFD = XF
      YFD = YF
      IF (IPACK.LT.1.OR.IPACK.GT.LDM_MAXPK) RETURN
      IF (IPLATE.LT.1.OR.IPLATE.GT.LDM_MAXPL) RETURN
C
C====== Twist, tilt, bulge distortion
C
      XQ = XF*COS_WC - YF*SIN_WC
      YQ = XF*SIN_WC + YF*COS_WC
      IF (DSTOR.EQ.1) THEN
         K_FACTOR = 1.0 + FACT*(TILT(IPACK,IPLATE)*XQ
     +              + TWIST(IPACK,IPLATE)*YQ
     +              + BULGE(IPACK,IPLATE)*SQRT(XQ*XQ+YQ*YQ))
         XFD = K_FACTOR*XQ
         YFD = K_FACTOR*YQ
C
C====== Twist, tilt, roff, toff (Note: assumes predicted centre very close
C       to scanner centre in roff/toff correction
C
      ELSE IF (DSTOR.EQ.2) THEN
         RQ = SQRT(XF*XF+YF*YF)
         IF (RQ.LT.0.0001) THEN
            CPSI = 1.0
            SPSI = 0.0
         ELSE
            CPSI = XQ/RQ
            SPSI = YQ/RQ
         END IF
         K_FACTOR = 1.0 + FACT*(TILT(IPACK,IPLATE)*XQ
     +              + TWIST(IPACK,IPLATE)*YQ)
         XFD = K_FACTOR*XQ
     +        + 0.01*(-TOFF(IPACK,IPLATE)*SPSI+ROFF(IPACK,IPLATE)*CPSI)
         YFD = K_FACTOR*YQ
     +        +0.01*(TOFF(IPACK,IPLATE)*CPSI+ROFF(IPACK,IPLATE)*SPSI)
      ELSE IF (DSTOR.EQ.3) THEN
         CALL LDM_SPDXY (IPACK, IPLATE, SPD_XY,
     +                   NSPDX, SPDX_MIN, SPDX_MAX, SXVAL,
     +                   NSPDY, SPDY_MIN, SPDY_MAX, SYVAL,
     +                   LDM_MAXSPD, IERR, ERRSTR)
         K_FACTOR = 1.0 + FACT*(TILT(IPACK,IPLATE)*XQ
     +              + TWIST(IPACK,IPLATE)*YQ)
         XFD = K_FACTOR*XQ
         YFD = K_FACTOR*YQ
         XFD = XFD + SPD_XY*YFD
         XD = XFD
         IF (XD.LT.SPDX_MIN) XD = SPDX_MIN
         IF (XD.GT.SPDX_MAX) XD = SPDX_MAX
         D = 0.0
         DD = 0.0
         XP = (2.0*XD-SPDX_MIN-SPDX_MAX)/(SPDX_MAX-SPDX_MIN)
         XP2 = 2.0*XP
         DO 100 I=NSPDX,2,-1
            SV = D
            D = XP2*D - DD + SXVAL(I)
            DD = SV
100      CONTINUE
         XFD = XFD + XP*D - DD + 0.5*SXVAL(1)
         YD = YFD
         IF (YD.LT.SPDY_MIN) YD = SPDY_MIN
         IF (YD.GT.SPDY_MAX) YD = SPDY_MAX
         D = 0.0
         DD = 0.0
         YP = (2.0*YD-SPDY_MIN-SPDY_MAX)/(SPDY_MAX-SPDY_MIN)
         YP2 = 2.0*YP
         DO 110 I = NSPDY,2,-1
            SV = D
            D = YP2*D - DD + SYVAL(I) 
            DD = SV
110      CONTINUE
         YFD = YFD + YP*D - DD + 0.5*SYVAL(1)
      END IF
      RETURN
      END
