C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**     LFN_CALFUN      **
C*************************
C
C
      SUBROUTINE  LFN_CALFUN (M, N, FVECC, XC)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Calculates function values for refinement routine for LFN_RFN
C
C **** Service Routine Only ****
C
C Author:  Based on code from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER M
      INTEGER N
      DOUBLE PRECISION XC(N)
      DOUBLE PRECISION FVECC(M)
C
C   M        (R)  Number of observations
C   N        (R)  Number of parameters being refined
C   XC(N)    (R) Vector of variables
C   FVECC(M) (W) Vector of function values
C
C====== GLOBALS:
C
      INCLUDE 'refine_limits.finc'
      INCLUDE 'refobs_list.finc'
      INCLUDE 'refine_data.finc'
C
C====== FUNCTION TYPES:
C
      LOGICAL LDM_RTOD_DP
C
C====== LOCALS:
C
      LOGICAL          RECORDED
      INTEGER          I, J 
      INTEGER          II
      INTEGER          ISYST
      INTEGER          DISTOR_TYPE
      INTEGER          IROT, IBEAM
      INTEGER          LFLAG, ND, LERR
      INTEGER          H, K, L
      INTEGER          SI, SJ, SK, ILOOP, JLOOP
      INTEGER          IANG
      INTEGER          NUM_CELL_PARAMS
      INTEGER          NCP(7)
      INTEGER          NSPDX
      INTEGER          NSPDY
      REAL             SPDX_MIN
      REAL             SPDX_MAX
      REAL             SPDY_MIN
      REAL             SPDY_MAX
      REAL             CELL(6)
      REAL             RCELL(6)
      REAL             CELL_PARAMS(6)
      REAL             MM_RAST_X
      REAL             MM_RAST_Y
      REAL             SPINDLE
      DOUBLE PRECISION A(3,3), SPIN(3,3)
      DOUBLE PRECISION X, Y, Z
      DOUBLE PRECISION XF, YF
      DOUBLE PRECISION XM, YM
      DOUBLE PRECISION X1, Y1, Z1
      DOUBLE PRECISION X_CEN, Y_CEN, OMEGA, Y_SCRAST
      DOUBLE PRECISION COS_OMEGA, SIN_OMEGA
      DOUBLE PRECISION RF, K_FACTOR, DEG_TO_RAD
      DOUBLE PRECISION PC(3,3), C(3,3), CX(3,3), CY(3,3), CZ(3,3)
      DOUBLE PRECISION CPX, CPY, CPZ, SPX, SPY, SPZ
      DOUBLE PRECISION RL(3), CC(3), SS(3)
      DOUBLE PRECISION COS_PSI, COS_RHO, A_AXIS
      DOUBLE PRECISION COS_ALPHA, SIN_ALPHA, FUN
      DOUBLE PRECISION RQ, XQ, YQ
      DOUBLE PRECISION CPSI, SPSI
      DOUBLE PRECISION DSTAR2
      DOUBLE PRECISION TORAD, FACT
      DOUBLE PRECISION XFD, YFD, XD, YD, XP, YP, XP2, YP2, D, DD, SV
      CHARACTER*2      VALSTR
      CHARACTER*2      ERRSTR
C
C====== Local double precision copies of the params being refined.
C
      DOUBLE PRECISION     LOC_CELL_PARAMS(6)
      DOUBLE PRECISION     LOC_PHIS(3)
      DOUBLE PRECISION     LOC_CTOF
      DOUBLE PRECISION     LOC_X_C
      DOUBLE PRECISION     LOC_Y_C
      DOUBLE PRECISION     LOC_W_C
      DOUBLE PRECISION     LOC_TWIST
      DOUBLE PRECISION     LOC_TILT
      DOUBLE PRECISION     LOC_BULGE
      DOUBLE PRECISION     LOC_ROFF
      DOUBLE PRECISION     LOC_TOFF
      DOUBLE PRECISION     LOC_Y_SCALE
      DOUBLE PRECISION     LOC_SPDXY
      DOUBLE PRECISION     LOC_SPDX(MAX_SPD)
      DOUBLE PRECISION     LOC_SPDY(MAX_SPD)

C
C====== Data statements
C
      DATA NCP/6,4,3,2,2,2,1/
      DATA TORAD /1.745329D-04/
C
C------------------------------------------------------------------------------
C
      DEG_TO_RAD = 3.1415926536/180.0
      CALL LDM_GET_SYST (VALSTR, ISYST, LFLAG)
      CALL LDM_GET_DSTOR (VALSTR, DISTOR_TYPE, LFLAG)
      CALL LDM_GET_ROT (VALSTR, IROT, LFLAG)
      CALL LDM_GET_BEAM (VALSTR, IBEAM, LFLAG)
      CALL LDM_GET_SPIN (KRPACK, SPINDLE, LFLAG, ND, LERR, ERRSTR)
      CALL LDM_MMRAST (KRPACK, KRPLATE, MM_RAST_X, MM_RAST_Y)
      NUM_CELL_PARAMS = NCP(ISYST)
      IF (DISTOR_TYPE.EQ.3) THEN
         CALL LDM_GET_SPDXN (KRPACK, KRPLATE, NSPDX, LFLAG, 
     +                       LERR, ERRSTR)
         CALL LDM_GET_SXMIN (KRPACK, KRPLATE, SPDX_MIN, LFLAG, ND,
     +                       LERR, ERRSTR)
         CALL LDM_GET_SXMAX (KRPACK, KRPLATE, SPDX_MAX, LFLAG, ND,
     +                       LERR, ERRSTR)
         CALL LDM_GET_SPDYN (KRPACK, KRPLATE, NSPDY, LFLAG, 
     +                       LERR, ERRSTR)
         CALL LDM_GET_SYMIN (KRPACK, KRPLATE, SPDY_MIN, LFLAG, ND,
     +                       LERR, ERRSTR)
         CALL LDM_GET_SYMAX (KRPACK, KRPLATE, SPDY_MAX, LFLAG, ND,
     +                       LERR, ERRSTR)
      END IF
C
C====== Unpack parameters and correct for scaling
C
C
C====== General setting parameters
C
      DO 20 I = 1 , 3
         IF (REF_FLAG(I).GT.0) THEN
            LOC_PHIS(I) = XC(REF_FLAG(I))*SCALING(REF_FLAG(I)) 
     +                    + X_START(I)
         ELSE
            LOC_PHIS(I) = X_START(I)
         END IF
 20   CONTINUE
      IF (REF_FLAG(4).GT.0) THEN
         LOC_CTOF = XC(REF_FLAG(4))*SCALING(REF_FLAG(4)) 
     +              + X_START(4)
      ELSE
         LOC_CTOF = X_START(4)
      END IF
      IF (REF_FLAG(5).GT.0) THEN
         LOC_X_C = XC(REF_FLAG(5))*SCALING(REF_FLAG(5)) 
     +             + X_START(5)
      ELSE
         LOC_X_C = X_START(5)
      END IF
      IF (REF_FLAG(6).GT.0) THEN
         LOC_Y_C = XC(REF_FLAG(6))*SCALING(REF_FLAG(6)) 
     +             + X_START(6)
      ELSE
         LOC_Y_C = X_START(6)
      END IF
      IF (REF_FLAG(7).GT.0) THEN
         LOC_W_C = XC(REF_FLAG(7))*SCALING(REF_FLAG(7)) 
     +             + X_START(7)
      ELSE
         LOC_W_C = X_START(7)
      END IF
C
C====== Cell parameters
C
      DO 30 I = 1, NUM_CELL_PARAMS
         IF (RFCL_FLAG(I).GT.0) THEN
            LOC_CELL_PARAMS(I) = XC(RFCL_FLAG(I))*SCALING(RFCL_FLAG(I)) 
     +                           + XCELL_START(I)
         ELSE
            LOC_CELL_PARAMS(I) = XCELL_START(I)
         END IF
 30   CONTINUE
C
C====== Distortion parameters
C
      IF (DSTOR_FLAG(1).GT.0) THEN
         LOC_Y_SCALE = XC(DSTOR_FLAG(1))*SCALING(DSTOR_FLAG(1)) 
     +                 + XDSTOR_START(1)
      ELSE
         LOC_Y_SCALE = XDSTOR_START(1)
      END IF
      IF (DSTOR_FLAG(2).GT.0) THEN
         LOC_TWIST = XC(DSTOR_FLAG(2))*SCALING(DSTOR_FLAG(2))
     +               + XDSTOR_START(2)
      ELSE
         LOC_TWIST = XDSTOR_START(2)
      END IF
      IF (DSTOR_FLAG(3).GT.0) THEN
         LOC_TILT = XC(DSTOR_FLAG(3))*SCALING(DSTOR_FLAG(3)) 
     +               + XDSTOR_START(3)
      ELSE
         LOC_TILT = XDSTOR_START(3)
      END IF
      IF (DISTOR_TYPE.EQ.1) THEN
         IF (DSTOR_FLAG(4).GT.0) THEN
            LOC_BULGE = XC(DSTOR_FLAG(4))*SCALING(DSTOR_FLAG(4)) 
     +                  + XDSTOR_START(4)
         ELSE
            LOC_BULGE = XDSTOR_START(4)
         END IF
      ELSE IF (DISTOR_TYPE.EQ.2) THEN
         IF (DSTOR_FLAG(4).GT.0) THEN
            LOC_ROFF = XC(DSTOR_FLAG(4))*SCALING(DSTOR_FLAG(4)) 
     +                 + XDSTOR_START(4)
         ELSE
            LOC_ROFF = XDSTOR_START(4)
         END IF
         IF (DSTOR_FLAG(5).GT.0) THEN
            LOC_TOFF = XC(DSTOR_FLAG(5))*SCALING(DSTOR_FLAG(5)) 
     +                 + XDSTOR_START(5)
         ELSE
            LOC_TOFF = XDSTOR_START(5)
         END IF
      ELSE IF (DISTOR_TYPE.EQ.3) THEN
         IF (DSTOR_FLAG(4).GT.0) THEN
            LOC_SPDXY = XC(DSTOR_FLAG(4))*SCALING(DSTOR_FLAG(4)) 
     +                  + XDSTOR_START(4)
         ELSE
            LOC_SPDXY = XDSTOR_START(4)
         END IF
         DO 60 I = 1, NSPDX
            IF (DSTOR_FLAG(I+4).GT.0) THEN
               N = DSTOR_FLAG(I+4)
               LOC_SPDX(I) = XC(N)*SCALING(N) + XDSTOR_START(I+4)
            ELSE
               LOC_SPDX(I) = XDSTOR_START(I+4)
            END IF
60       CONTINUE
         DO 70 I = 1, NSPDY
            IF (DSTOR_FLAG(I+NSPDX+4).GT.0) THEN
               N = DSTOR_FLAG(I+NSPDX+4)
               LOC_SPDY(I) = XC(N)*SCALING(N) 
     +                        + XDSTOR_START(NSPDX+I+4)
            ELSE
               LOC_SPDY(I) = XDSTOR_START(NSPDX+I+4)
            END IF
70       CONTINUE
      END IF
 
C
C====== Now prepare to evaluate residual
C
      X_CEN = LOC_XCENF + LOC_X_C*MM_RAST_X
      Y_CEN = LOC_YCENF + LOC_Y_C*MM_RAST_Y

      OMEGA = LOC_W_C
      COS_OMEGA = COS(OMEGA*DEG_TO_RAD)
      SIN_OMEGA = SIN(OMEGA*DEG_TO_RAD)
 
      IF (REFINE_CELL .OR. FIRST) THEN
C
C====== In-line, R*8 version of MAKE_ORIENT--------------------
C====== i, j, k have special meanings in this piece of code.
C====== i refers to the crystal rotation axis, Z
C
         FIRST = .FALSE.
         I = IABS(IROT)
         SI = 1
         IF (IROT.LT.0) SI = -1
C
C====== j refers to X_axis, the x-ray beam axis, X
C
         J = IABS(IBEAM)
         SJ = 1
         IF (IBEAM.LT.0) SJ = -1
C
C====== k is therefore the other axis, Y
C
         K = 6 - I - J
C
C====== The sign of k relies on i,j,k and X,Y,Z being right-handed
C
         IF (K.EQ.J+1 .OR. K.EQ.J-2) THEN
            SK = SI*SJ
         ELSE
            SK = -SI*SJ
         END IF
C
C====== Zero orientation matrix
C
         DO 90 ILOOP = 1 , 3
            DO 80 JLOOP = 1 , 3
               LOC_ORIENT(ILOOP,JLOOP) = 0.
 80         CONTINUE
 90      CONTINUE
C
C====== Convert cell angles to radians
C
         IANG = 4
C
C====== Check if rhombohedral
C
         IF (ISYST.EQ.6) IANG = 2
C
         DEG_TO_RAD = 3.1415926536/180.0
         LOC_CELL_PARAMS(IANG) = LOC_CELL_PARAMS(IANG)*DEG_TO_RAD
         LOC_CELL_PARAMS(IANG+1) = LOC_CELL_PARAMS(IANG+1)*DEG_TO_RAD
         LOC_CELL_PARAMS(IANG+2) = LOC_CELL_PARAMS(IANG+2)*DEG_TO_RAD
C
C====== Now go and do specific stuff for each crystal system
C
         IF (ISYST.EQ.2) THEN
C
C====== Monoclinic
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(2)
            RL(3) = LOC_CELL_PARAMS(3)
 
            CC(1) = 0.0
            CC(2) = COS(LOC_CELL_PARAMS(4))
            CC(3) = 0.0
 
            SS(1) = 1.0
            SS(2) = SIN(LOC_CELL_PARAMS(4))
            SS(3) = 1.0
 
            LOC_ORIENT(1,J) = SJ*RL(J)*SS(K)
            LOC_ORIENT(1,K) = SJ*RL(K)*CC(I)
            LOC_ORIENT(2,K) = SK*RL(K)*SS(I)*SS(J)
            LOC_ORIENT(3,I) = SI*RL(I)
            LOC_ORIENT(3,J) = SI*RL(J)*CC(K)
            LOC_ORIENT(3,K) = SI*RL(K)*CC(J)
         ELSE IF (ISYST.EQ.3) THEN
C
C====== Orthorhombic
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(2)
            RL(3) = LOC_CELL_PARAMS(3)
 
            LOC_ORIENT(1,J) = SJ*RL(J)
            LOC_ORIENT(2,K) = SK*RL(K)
            LOC_ORIENT(3,I) = SI*RL(I)
         ELSE IF (ISYST.EQ.4) THEN
C
C====== Tetragonal
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(1)
            RL(3) = LOC_CELL_PARAMS(2)
 
            LOC_ORIENT(1,J) = SJ*RL(J)
            LOC_ORIENT(2,K) = SK*RL(K)
            LOC_ORIENT(3,I) = SI*RL(I)
         ELSE IF (ISYST.EQ.5) THEN
C
C====== Hexagonal
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(1)
            RL(3) = LOC_CELL_PARAMS(2)
 
            CC(1) = 0.0
            CC(2) = 0.0
            CC(3) = 0.5
 
            SS(1) = 1.0
            SS(2) = 1.0
            SS(3) = SQRT(3.0)/2.0
 
            LOC_ORIENT(1,J) = SJ*RL(J)*SS(K)
            LOC_ORIENT(1,K) = SJ*RL(K)*CC(I)
            LOC_ORIENT(2,K) = SK*RL(K)*SS(I)*SS(J)
            LOC_ORIENT(3,I) = SI*RL(I)
            LOC_ORIENT(3,J) = SI*RL(J)*CC(K)
            LOC_ORIENT(3,K) = SI*RL(K)*CC(J)
         ELSE IF (ISYST.EQ.6) THEN
C
C====== Rhombohedral
C
            A_AXIS = LOC_CELL_PARAMS(1)
            COS_ALPHA = COS(LOC_CELL_PARAMS(2))
            SIN_ALPHA = SIN(LOC_CELL_PARAMS(2))
            FUN = SQRT(1.-3.*COS_ALPHA**2+2.*COS_ALPHA**3)
 
            LOC_ORIENT(1,J) = SJ*A_AXIS*SIN_ALPHA
            LOC_ORIENT(1,K) = SJ*A_AXIS*COS_ALPHA*(1.-COS_ALPHA)
     +                        /SIN_ALPHA
            LOC_ORIENT(2,K) = SK*A_AXIS*FUN/SIN_ALPHA
            LOC_ORIENT(3,I) = SI*A_AXIS
            LOC_ORIENT(3,J) = SI*A_AXIS*COS_ALPHA
            LOC_ORIENT(3,K) = SI*A_AXIS*COS_ALPHA
         ELSE IF (ISYST.EQ.7) THEN
C
C====== Cubic
C
            A_AXIS = LOC_CELL_PARAMS(1)
 
            LOC_ORIENT(1,J) = SJ*A_AXIS
            LOC_ORIENT(2,K) = SK*A_AXIS
            LOC_ORIENT(3,I) = SI*A_AXIS
         ELSE
C
C====== Triclinic
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(2)
            RL(3) = LOC_CELL_PARAMS(3)
 
            CC(1) = COS(LOC_CELL_PARAMS(4))
            CC(2) = COS(LOC_CELL_PARAMS(5))
            CC(3) = COS(LOC_CELL_PARAMS(6))
 
            SS(1) = SIN(LOC_CELL_PARAMS(4))
            SS(2) = SIN(LOC_CELL_PARAMS(5))
            SS(3) = SIN(LOC_CELL_PARAMS(6))
C
C====== See Buerger's Xray Crystallography for the next formulae, p348
C
            COS_PSI = (CC(I)-CC(K)*CC(J))/SS(K)
            COS_RHO = SQRT(1.0-CC(I)**2-CC(J)**2-CC(K)**2+2.0*CC(I)
     +                *CC(J)*CC(K))/SS(K)
 
            LOC_ORIENT(1,J) = SJ*RL(J)*SS(K)
            LOC_ORIENT(1,K) = SJ*RL(K)*COS_PSI
            LOC_ORIENT(2,K) = SK*RL(K)*COS_RHO
            LOC_ORIENT(3,I) = SI*RL(I)
            LOC_ORIENT(3,J) = SI*RL(J)*CC(K)
            LOC_ORIENT(3,K) = SI*RL(K)*CC(J)
         END IF
C
C====== End of code from MAKE_ORIENT--------------------
C
      END IF
C
C====== In-line version of A_CALC - made DOUBLE PRECISION throughout.---------
C
C
C====== Define phi rotation matrix
C
      CPX = COS(LOC_PHIS(1)*DEG_TO_RAD)
      SPX = SIN(LOC_PHIS(1)*DEG_TO_RAD)
      CPY = COS(LOC_PHIS(2)*DEG_TO_RAD)
      SPY = SIN(LOC_PHIS(2)*DEG_TO_RAD)
      CPZ = COS(LOC_PHIS(3)*DEG_TO_RAD)
      SPZ = SIN(LOC_PHIS(3)*DEG_TO_RAD)
 
      DO 120 I = 1 , 3
         DO 110 J = 1 , 3
            CX(I,J) = 0.
110      CONTINUE
120   CONTINUE
      CX(1,1) = 1.0
      CX(2,2) = CPX
      CX(2,3) = -SPX
      CX(3,2) = SPX
      CX(3,3) = CPX
 
      DO 140 I = 1 , 3
         DO 130 J = 1 , 3
            CY(I,J) = 0.
130      CONTINUE
140   CONTINUE
      CY(1,1) = CPY
      CY(1,3) = SPY
      CY(2,2) = 1.0
      CY(3,1) = -SPY
      CY(3,3) = CPY
 
      DO 160 I = 1 , 3
         DO 150 J = 1 , 3
            CZ(I,J) = 0.
 150     CONTINUE
 160  CONTINUE
      CZ(1,1) = CPZ
      CZ(1,2) = -SPZ
      CZ(2,1) = SPZ
      CZ(2,2) = CPZ
      CZ(3,3) = 1.0
 
      DO 190 I = 1 , 3
         DO 180 J = 1 , 3
            C(I,J) = 0.0
            DO 170 K = 1 , 3
               C(I,J) = C(I,J) + CY(I,K)*CX(K,J)
 170        CONTINUE
 180     CONTINUE
 190  CONTINUE
 
      DO 230 I = 1 , 3
         DO 220 J = 1 , 3
            PC(I,J) = 0.0
            DO 210 K = 1 , 3
               PC(I,J) = PC(I,J) + CZ(I,K)*C(K,J)
 210        CONTINUE
 220     CONTINUE
 230  CONTINUE
C
C====== Apply phi rotations
C
      DO 260 I = 1 , 3
         DO 250 J = 1 , 3
            A(I,J) = 0.0
            DO 240 K = 1 , 3
               A(I,J) = A(I,J) + PC(I,K)*LOC_ORIENT(K,J)
 240        CONTINUE
 250     CONTINUE
 260  CONTINUE
C
C====== End of text included from a_calc --------------------
C
C
C====== Set up spindle rotation matrix
C
      SPIN(1,1) = COS(SPINDLE*DEG_TO_RAD)
      SPIN(1,2) = -SIN(SPINDLE*DEG_TO_RAD)
      SPIN(2,1) = -SPIN(1,2)
      SPIN(2,2) = SPIN(1,1)
 
      DO 500 I = 1 , M
         H = HOBS(I)
         K = KOBS(I)
         L = LOBS(I)
 
         X1 = A(1,1)*H + A(1,2)*K + A(1,3)*L
         Y1 = A(2,1)*H + A(2,2)*K + A(2,3)*L
         Z1 = A(3,1)*H + A(3,2)*K + A(3,3)*L
C
C====== Apply spindle rotation
C
         X = SPIN(1,1)*X1 + SPIN(1,2)*Y1
         Y = SPIN(2,1)*X1 + SPIN(2,2)*Y1
         Z = Z1

         DSTAR2 = X*X + Y*Y + Z*Z
         RECORDED = LDM_RTOD_DP (X, Y, Z, DSTAR2, LOC_CTOF, XF, YF)
         Y_SCRAST = MM_RAST_X*LOC_Y_SCALE
 
         XQ = XF*COS_OMEGA - YF*SIN_OMEGA
         YQ = XF*SIN_OMEGA + YF*COS_OMEGA
         IF (DISTOR_TYPE.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
         END IF
C
C====== Make distortion corrections & add to centre position
C
      FACT = TORAD/LOC_CTOF
      IF (DISTOR_TYPE.EQ.1) THEN
C
C====== K_factor = 1.0 + tilt*xq + twist*yq + bulge*sqrt(xq*xq + yq*yq)
C
         K_FACTOR = 1.0 + FACT*LOC_TILT*XQ + FACT*LOC_TWIST*YQ
     +              + FACT*LOC_BULGE*SQRT(XQ*XQ+YQ*YQ)
 
         XM = X_CEN + K_FACTOR*XQ*MM_RAST_X
         YM = Y_CEN + K_FACTOR*YQ*Y_SCRAST

      ELSE IF (DISTOR_TYPE.EQ.2) THEN
C
C====== K_factor = 1.0 + tilt*xq + twist*yq 
C====== roff/toff correction on x = -toff*sin(psi) + roff*cos(psi)
C======                      on y =  toff*cos(psi) + roff*sin(psi)
C
         K_FACTOR = 1.0 + FACT*LOC_TILT*XQ + FACT*LOC_TWIST*YQ 
 
         XM = X_CEN + K_FACTOR*XQ*MM_RAST_X
     +        + 0.01*MM_RAST_X*(-LOC_TOFF*SPSI+LOC_ROFF*CPSI)
         YM = Y_CEN + K_FACTOR*YQ*Y_SCRAST
     +        + 0.01*Y_SCRAST*(LOC_TOFF*CPSI+LOC_ROFF*SPSI)
C
C====== K_factor = 1.0 + tilt*xq + twist*yq 
C====== Polynomial distortion correction spdx1*xd + spdx2*xd**2 ...
C======                                  spdy1*yd + spdy2*yd**2 ...
C====== (xd, yd, distances from distortion centre)
C
      ELSE IF (DISTOR_TYPE.EQ.3) THEN
         K_FACTOR = 1.0 + FACT*LOC_TILT*XQ + FACT*LOC_TWIST*YQ 
         XFD = K_FACTOR*XQ
         YFD = K_FACTOR*YQ
         XFD = XFD + LOC_SPDXY*YFD
         XD = XFD
         IF (XD.LT.SPDX_MIN) THEN
             XD = SPDX_MIN
             CHEB_OK = .FALSE.
         END IF
         IF (XD.GT.SPDX_MAX) THEN
             XD = SPDX_MAX
             CHEB_OK = .FALSE.
         END IF
         D = 0.0
         DD = 0.0
         XP = (2.0*XD-SPDX_MIN-SPDX_MAX)/(SPDX_MAX-SPDX_MIN)
         XP2 = 2.0*XP
         DO 310 II = NSPDX,2,-1
            SV = D
            D = XP2*D - DD + LOC_SPDX(II)
            DD = SV
310      CONTINUE
         XFD = XFD + XP*D - DD + 0.5*LOC_SPDX(1)
         YD = YFD
         IF (YD.LT.SPDY_MIN) THEN
            YD = SPDY_MIN
            CHEB_OK = .FALSE.
         END IF
         IF (YD.GT.SPDY_MAX) THEN
            YD = SPDY_MAX
            CHEB_OK = .FALSE.
         END IF
         D = 0.0
         DD = 0.0
         YP = (2.0*YD-SPDY_MIN-SPDY_MAX)/(SPDY_MAX-SPDY_MIN)
         YP2 = 2.0*YP
         DO 320 II = NSPDY,2,-1
            SV = D
            D = YP2*D - DD + LOC_SPDY(II) 
            DD = SV
320      CONTINUE
         YFD = YFD + YP*D - DD + 0.5*LOC_SPDY(1)
         XM = X_CEN + XFD*MM_RAST_X
         YM = Y_CEN + YFD*Y_SCRAST
      END IF
C
C====== Calculate new positional error
C
         FVECC(I) = SQRT(((XM_OBS(I)-XM)/MM_RAST_X)**2+((YM_OBS(I)-YM)/
     +              MM_RAST_Y)**2)
 
 500  CONTINUE
 
      RETURN
      END
