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