linse/sicspsi
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
07cd8a0215003421c6a6432083f28c99d2aef42e
sicspsi/t_rlp.f
T
cvs cd13637987 - Updated Makefiles 
- Moved TAS code to psi
- Updated programmers documentation
2003-06-30 11:51:39 +00:00

464 lines
13 KiB
FortranFixed
Executable File
Raw Blame History

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
Reference in New Issue View Git Blame Copy Permalink