SUBROUTINE T_RLP ! File [MAD.SRC]T_RLP.FOR c ================ c cdec$ Ident 'V01A' c------------------------------------------------------------------ c Updates: c V01A 7-May-1996 DM. Put error output to IOLUN, use IMPLICIT NONE and c get the code indented so that it is readable! c------------------------------------------------------------------ c Routines to deal with the reciprocical lattice PB c------------------------------------------------------------------ c Entry points in this file: c c SETRLP : CALCULATION OF S AND INVS , ORIENTATION MATRIX C RL2SPV : TRANSFO FROM RECIP LAT TO SCAT PLANE C SP2RLV : TRANSFO FROM SCAT PLANE TO RECIP LAT C INVS : INVERT MATRIX S, GENERATED BY SETRLP. C ERRESO : DEAL ITH ERROR MESSAGES FOR ALL MODULES C C SUBROUTINE SETRLP(SAM,IER) C SUBROUTINE RL2SPV(QHKL,QT,QM,QS,IER) C SUBROUTINE SP2RLV(QHKL,QT,QM,QS,IER) C SUBROUTINE INVS(S,SINV,IER) C SUBROUTINE ERRESO(MODULE,IER) c------------------------------------------------------------------ implicit none end c SUBROUTINE SETRLP (SAM, IER) ! Part of [MAD.SRC]T_RLP.FOR c ============================ c C SETRLP: Computation of matrix S which transforms (QH,QK,QL) to C vector (Q1,Q2) in scattering plane (defined by vectors A1,A2) C and SINV matrix for the inverse transformation C C INPUT SAM SAMPLE CHARACTERISTICS C SAM(1)=AS LATTICE PARAMETERS C SAM(2)=BS ------------------ C SAM(3)=CS ------------------ C SAM(4)=AA LATTICE ANGLES C SAM(5)=BB -------------- C SAM(6)=CC -------------- C SAM(7)=AX VECTOR A IN SCATTERING PLANE C SAM(8)=AY ---------------------------- C SAM(9)=AZ ---------------------------- C SAM(10)=BX VECTOR B IN SCATTERING PLANE C SAM(11)=BY ---------------------------- C SAM(12)=BZ ---------------------------- C OUTPUT IER ERROR RETURN TO BE TREATED BY ERRESO C IER=1 ERROR ON LATTICE PARAMETERS C IER=2 ERROR ON LATTICE ANGLES C IER=3 ERROR ON VECTORS A1, A2 c------------------------------------------------------------------ IMPLICIT NONE c DOUBLE PRECISION PI PARAMETER (PI=3.14159265358979323846264338327950D0) DOUBLE PRECISION RD PARAMETER (RD=57.29577951308232087679815481410517D0) DOUBLE PRECISION EPS PARAMETER (EPS=1.D-1) DOUBLE PRECISION EPS6 PARAMETER (EPS6=1.D-6) DOUBLE PRECISION EPS8 PARAMETER (EPS8=1.D-8) c------------------------------------------------------------------ C Define the dummy arguments real*4 sam(64) integer*4 ier c------------------------------------------------------------------ C DO NOT EXPORT THE FOLLOWING COMMON ! C IT IS JUST FOR PERMANENT STORAGE USE C DOUBLE PRECISION S, SINV integer*4 iok COMMON /OSOLEM/ S(4,4), SINV(4,4), iok C----------------------------------------------------------------------- double precision a(3), aspv(3,2), rlb(3,2) double precision alfa(3), sina(3), sinb(3), cosa(3), cosb(3) double precision b(3), c(3) double precision vv(3,3), bb(3,3) C double precision zp, cc integer*4 imod integer*4 id, ie integer*4 jd, je, jf integer*4 kg integer*4 lh, lf integer*4 md, me integer*4 ne C----------------------------------------------------------------------- C SOME TESTS AND INIT OF CALCUALTION C ier = 0 IMOD = 1 ZP = 2.D0*PI IOK = 0 DO ID = 1,3 A(ID) = SAM(ID) ALFA(ID) = SAM(ID+3) ASPV(ID,1) = SAM(6+ID) ASPV(ID,2) = SAM(9+ID) ENDDO C DO ID = 1,3 IER = 1 IF(ABS(A(ID)).LE.EPS8) GOTO 999 IER = 0 ENDDO DO ID = 1,3 A(ID) = A(ID)/ZP ALFA(ID) = ALFA(ID)/RD COSA(ID) = COS(ALFA(ID)) SINA(ID) = SIN(ALFA(ID)) ENDDO CC = COSA(1)*COSA(1)+COSA(2)*COSA(2)+COSA(3)*COSA(3) CC = (1.D0+2.D0*COSA(1)*COSA(2)*COSA(3)-CC) IER = 2 IF(CC.LE.EPS) GOTO 999 IER = 0 CC = SQRT(CC) JE = 2 KG = 3 DO ID = 1,3 B(ID) = SINA(ID)/(A(ID)*CC) COSB(ID) = (COSA(JE)*COSA(KG)-COSA(ID))/(SINA(JE)*SINA(KG)) SINB(ID) = SQRT(1.D0-COSB(ID)*COSB(ID)) RLB(ID,2) = ABS(ATAN(SINB(ID)/COSB(ID)))*RD JE = KG KG = ID ENDDO BB(1,1) = B(1) BB(2,1) = 0.D0 BB(3,1) = 0.D0 BB(1,2) = B(2)*COSB(3) BB(2,2) = B(2)*SINB(3) BB(3,2) = 0.D0 BB(1,3) = B(3)*COSB(2) BB(2,3) = -B(3)*SINB(2)*COSA(1) BB(3,3) = 1.D0/A(3) C DO ID = 1,3 RLB(ID,1) = 0.D0 DO JE = 1,3 RLB(ID,1) = RLB(ID,1)+BB(JE,ID)**2 ENDDO IER = 1 IF (ABS(RLB(ID,1)).LE.EPS8) GOTO 999 IER = 0 RLB(ID,1) = SQRT(RLB(ID,1)) ENDDO C----------------------------------------------------------------------- C GENERATION OF S ORIENTATION MATRIX REC. LATTICE TO SCATTERING PLANE C DO KG = 1,2 DO IE = 1,3 VV(KG,IE) = 0.D0 DO JF = 1,3 VV(KG,IE) = VV(KG,IE)+BB(IE,JF)*ASPV(JF,KG) ENDDO ENDDO ENDDO DO MD = 3,2,-1 DO NE = 1,3 ID = MOD(MD,3)+1 JE = MOD(MD+1,3)+1 KG = MOD(NE,3)+1 LH = MOD(NE+1,3)+1 VV(MD,NE) = VV(ID,KG)*VV(JE,LH)-VV(ID,LH)*VV(JE,KG) ENDDO ENDDO C DO ID = 1,3 C(ID) = 0.D0 DO JE = 1,3 C(ID) = C(ID)+VV(ID,JE)**2 ENDDO IER = 3 IF (ABS(C(ID)).LE.EPS6) GOTO 999 IER = 0 C(ID) = SQRT(C(ID)) ENDDO C DO ID = 1,3 DO JE = 1,3 VV(JE,ID) = VV(JE,ID)/C(JE) ENDDO ENDDO DO KG = 1,3 DO ME = 1,3 S(KG,ME) = 0.D0 DO LF = 1,3 S(KG,ME) = S(KG,ME)+VV(KG,LF)*BB(LF,ME) ENDDO ENDDO ENDDO S(4,4) = 1.D0 DO JD = 1,3 S(4,JD) = 0.D0 S(JD,4) = 0.D0 ENDDO C----------------------------------------------------------------------- C INVERT TRANSFORMATION MATRIX S AND PU RESULT IN SINV C IER = 3 CALL INVS(S,SINV,IER) IER = 0 IF (IER.NE.0)GOTO 999 IOK = 123 C--------------------------------------------------------------------------- C SORTIE C 999 CONTINUE IF (IER.NE.0) CALL ERRESO(IMOD,IER) RETURN END c SUBROUTINE RL2SPV (QHKL, QT, QM, QS, IER) ! Part of [MAD.SRC]T_RLP.FOR c ========================================= c c------------------------------------------------------------------ C INPUT QHKL QHKL -> QT C A Q-VECTOR TO BE TRANSFORM FROM RECIP LATTICE TO SCATTERING PLANE C CHECK THAT Q-VECTOR IS IN THE PLANE C C INPUT Q-VECTOR QHKL(3) Q-VECTOR IN RECIPROCICAL LATTICVE C C OUTPUT Q-VECTOR QT(3) Q-VECTOR IN SCATTERING PLANE C OUTPUT QM AND QS QMODULUS AND ITS SQUARE ( TO BE VERIFIED ) C OUTPUT ERROR IER C IER=1 MATRIX S NOT OK C IER=2 Q NOT IN SCATTERING PLANE C IER=3 Q MODULUS TOO SMALL c------------------------------------------------------------------ IMPLICIT NONE c DOUBLE PRECISION DEPS4 PARAMETER (DEPS4=1.D-4) DOUBLE PRECISION DEPS8 PARAMETER (DEPS8=1.D-8) c------------------------------------------------------------------ c Define the dummy arguments DOUBLE PRECISION QHKL(3) DOUBLE PRECISION QT(3) DOUBLE PRECISION QM DOUBLE PRECISION QS integer*4 IER c------------------------------------------------------------------ C DO NOT EXPORT THE FOLLWING COOMON ! C IT IS JUST FOR PERMANENT STORAGE USE C DOUBLE PRECISION S, SINV integer*4 iok COMMON /OSOLEM/ S(4,4), SINV(4,4), iok C--------------------------------------------------------------------------- integer*4 id, je C--------------------------------------------------------------------------- C INIT AND TEST IF TRANSFO MATRICES ARE OK IER = 1 IF (IOK.NE.123) GOTO 999 IER = 0 C----------------------------------------------------------------------- C DO ID = 1,3 QT(ID) = 0.D0 DO JE = 1,3 QT(ID) = QT(ID) + QHKL(JE)*S(ID,JE) ENDDO ENDDO IER = 2 IF(ABS(QT(3)).GT.DEPS4) GOTO 999 IER = 0 QS = 0.D0 DO ID = 1,3 QS = QS+QT(ID)**2 ENDDO IF(QS.LT.DEPS8) THEN IER = 3 ELSE QM = SQRT(QS) ENDIF C--------------------------------------------------------------------------- C 999 CONTINUE RETURN END c SUBROUTINE SP2RLV (QHKL, QT, QM, QS, IER) ! Part of [MAD.SRC]T_RLP.FOR c ========================================= c c------------------------------------------------------------------ C INPUT QT QHKL <- QT C A Q-VECTOR TO BE TRANSFORM FROM SCATTERING PLANE TO RECIP LATTICE C CHECK THAT Q, D & G VECTORS ARE IN THE SCATTERING PLANE C C INPUT Q-VECTOR QT(3) Q-VECTOR IN SCATTERING PLANE C C OUTPUT Q-VECTOR QHKL(3) Q-VECTOR IN RECIPROCICAL LATTICVE C OUTPUT QM AND QS QMODULUS AND ITS SQUARE ( TO BE VERIFIED ) C OUTPUT ERROR IER C IER=1 MATRIX S NOT OK C IER=2 Q NOT IN SCATTERING PLANE C IER=3 Q MODULUS TOO SMALL c------------------------------------------------------------------ IMPLICIT NONE C DOUBLE PRECISION EPS4 PARAMETER (EPS4=1.D-4) DOUBLE PRECISION EPS8 PARAMETER (EPS8=1.D-8) c------------------------------------------------------------------ c Define the dummy arguments DOUBLE PRECISION QHKL(3) DOUBLE PRECISION QT(3) DOUBLE PRECISION QM DOUBLE PRECISION QS integer*4 IER c------------------------------------------------------------------ C DO NOT EXPORT THE FOLLWING COOMON ! C IT IS JUST FOR PERMANENT STORAGE USE C DOUBLE PRECISION S, SINV integer*4 iok COMMON /OSOLEM/ S(4,4), SINV(4,4), iok C--------------------------------------------------------------------------- integer*4 id, je C--------------------------------------------------------------------------- C INIT AND TEST IF TRANSFO MATRICES ARE OK IER = 1 IF (IOK.NE.123) GOTO 999 IER = 2 IF(ABS(QT(3)).GT.EPS4) GOTO 999 IER = 0 C----------------------------------------------------------------------- QS = 0.D0 DO ID = 1,3 QS = QS+QT(ID)**2 ENDDO IF(QS.LT.EPS8) THEN IER = 3 ELSE QM = SQRT(QS) ENDIF C----------------------------------------------------------------------- C DO ID = 1,3 QHKL(ID) = 0.D0 DO JE = 1,3 QHKL(ID) = QHKL(ID)+SINV(ID,JE)*QT(JE) ENDDO ENDDO C--------------------------------------------------------------------------- C 999 CONTINUE RETURN END c SUBROUTINE INVS (S, SINV, IER) ! Part of [MAD.SRC]T_RLP.FOR c ============================== c c------------------------------------------------------------------ C ROUTINE TO INVERT MATRIX S, GENERATED BY SETRLP, WHICH TRANSFORMS C (QH,QK,QL) TO (Q1,Q2) IN THE SCATTERING PLANE C INPUT MATRIX DOUBLE PRECISION S(4,4) C OUTPUT MATRIX DOUBLE PRECISION SINV(4,4) C OUTPUT ERROR IER C IER=1 DETERMINANT OF MATRIX S TOO SMALL c------------------------------------------------------------------ IMPLICIT NONE c DOUBLE PRECISION EPS6 PARAMETER (EPS6=1.D-6) c------------------------------------------------------------------ c Define the dummy arguments DOUBLE PRECISION S(4,4) DOUBLE PRECISION SINV(4,4) integer*4 IER c------------------------------------------------------------------ integer*4 m(3) integer*4 n(3) c integer*4 id, je integer*4 mi, mj, ni, nj double precision det c------------------------------------------------------------------ DATA M/2,3,1/ DATA N/3,1,2/ IER = 0 DO ID = 1,4 DO JE = 1,4 SINV(ID,JE) = 0.D0 ENDDO ENDDO DET = 0.D0 DO ID = 1,3 DO JE = 1,3 MI = M(ID) MJ = M(JE) NI = N(ID) NJ = N(JE) SINV(JE,ID) = S(MI,MJ)*S(NI,NJ)-S(NI,MJ)*S(MI,NJ) ENDDO DET = DET+S(ID,1)*SINV(1,ID) ENDDO IF(ABS(DET).LT.EPS6) THEN IER = 1 ELSE DO ID = 1,3 DO JE = 1,3 SINV(ID,JE) = SINV(ID,JE)/DET ENDDO ENDDO ENDIF SINV(4,4) = 1.D0 RETURN END c SUBROUTINE ERRESO(MODULE,IER) ! Part of [MAD.SRC]T_RLP.FOR c ============================= c c------------------------------------------------------------------ C SUBROUTINE TO TREAT ERRORS FROM RESOLUTION CALCULATIONS C MODULE = 1 -> SETRLP C MODULE = 2 -> RL2SPV C MODULE = 3 -> EX_CASE c------------------------------------------------------------------ IMPLICIT NONE c integer*4 MXER PARAMETER (MXER = 4) integer*4 MXMOD PARAMETER (MXMOD = 3) c C include 'iolsddef.inc' c------------------------------------------------------------------ c Define the dummy arguments integer*4 module integer*4 ier c------------------------------------------------------------------ integer*4 lmodule, lier CHARACTER*64 MESER(MXER,MXMOD) DATA MESER/ + ' Check lattice spacings (AS,BS,CS)', + ' Check cell angles (AA,BB,CC)', + ' Check scattering plane (AX....BZ)', + ' Check lattice and scattering plane', C + ' Check Lattice and Scattering Plane', + ' Q not in scattering plane', + ' Q modulus too small', + ' KI,KF,Q triangle cannot close', C + ' Error in KI or KF. Check D-spacings & units', + ' KI or KF cannot be obtained', + ' KI or KF too small', + ' KI,KF,Q triangle will not close'/ C--------------------------------------------------------------------------- LIER = MIN(MAX(IER,1),MXER) LMODULE = MIN(MAX(MODULE,1),MXMOD) C WRITE(iolun,501) MESER(LIER,LMODULE) 501 FORMAT(A) RETURN END