629 lines
21 KiB
Fortran
629 lines
21 KiB
Fortran
C-----------------------------------------------------------------------
|
||
C
|
||
C Constrained Cell Parameter Least Squares on Theta Data.
|
||
C Adapted from the routine CELLLS of the NRCVAX package.
|
||
C
|
||
C E.J.Gabe Chemistry Division, N.R.C., Ottawa Canada
|
||
C
|
||
C 2theta data is taken from the file ORIENT.DA, which must have been
|
||
C written by the AL command.
|
||
C
|
||
C-----------------------------------------------------------------------
|
||
SUBROUTINE CELLLS
|
||
INCLUDE 'COMDIF'
|
||
DIMENSION IBH(10),IBK(10),IBL(10),BTHETA(10),BOMEGA(10),BCHI(10),
|
||
$ BPHI(10),QOBS(NSIZE)
|
||
EQUIVALENCE (IHK(1),IBH(1)),(NREFB(1),IBK(1)),(ILA(1),IBL(1)),
|
||
$ (BCOUNT(1),BTHETA(1)),(BBGR1(1),BOMEGA(1)),
|
||
$ (BBGR2(1),BCHI(1)),(BTIME(1),BPHI(1)),
|
||
$ (ACOUNT(1),QOBS(1))
|
||
C-----------------------------------------------------------------------
|
||
C File data input. Skip reflections flagged bad in MM (Psi .ne. 0)
|
||
C-----------------------------------------------------------------------
|
||
WRITE (COUT,10000)
|
||
CALL GWRITE (ITP,' ')
|
||
IOUT = -1
|
||
CALL SPACEG (IOUT,0)
|
||
LAUE = LAUENO
|
||
IAXIS = NAXIS
|
||
IF (LAUENO .EQ. 4 .OR. LAUENO .EQ. 5) LAUE = 4
|
||
IF (LAUENO .EQ. 6 .OR. LAUENO .EQ. 7) LAUE = 7
|
||
IF (LAUENO .GE. 8 .AND. LAUENO .LE. 12) LAUE = 6
|
||
IF (LAUENO .EQ. 13 .OR. LAUENO .EQ. 14) LAUE = 5
|
||
NUMD = 0
|
||
NBLOKO = 250
|
||
100 READ (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,(JUNK, I = 41,70),
|
||
$ BPSI,NBL
|
||
NBLOKO = NBLOKO + 1
|
||
IF (NBL .NE. 0) THEN
|
||
DO 110 NB = 1,NBL
|
||
IF (BPSI(NB) .EQ. 0) THEN
|
||
NUMD = NUMD + 1
|
||
S = 2.0*SIN(0.5*BTHETA(NB)/DEG)/WAVE
|
||
QOBS(NUMD) = S*S
|
||
IOH(NUMD) = IBH(NB)
|
||
IOK(NUMD) = IBK(NB)
|
||
IOL(NUMD) = IBL(NB)
|
||
ENDIF
|
||
110 CONTINUE
|
||
GO TO 100
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Do the least squares
|
||
C-----------------------------------------------------------------------
|
||
IF (NUMD .GE. 6) THEN
|
||
WRITE (LPT,11000) WAVE,NUMD
|
||
CALL CLSTSQ
|
||
ELSE
|
||
WRITE (COUT,12000)
|
||
CALL GWRITE (ITP,' ')
|
||
ENDIF
|
||
KI = ' '
|
||
RETURN
|
||
10000 FORMAT (/10X,'Constrained Cell Dimension Least-Squares'/)
|
||
11000 FORMAT (/' Wavelength'F10.6,'; ',I6,' reflections.')
|
||
12000 FORMAT (' Less than 6 reflections. Quit')
|
||
END
|
||
C-----------------------------------------------------------------------
|
||
c General least-squares of lattice parameters
|
||
C-----------------------------------------------------------------------
|
||
SUBROUTINE CLSTSQ
|
||
INCLUDE 'COMDIF'
|
||
DIMENSION AI(6),SIG(7,7),PAR(6),QOBS(NSIZE)
|
||
EQUIVALENCE (ACOUNT(1),QOBS(1))
|
||
EQUIVALENCE (PAR(1),ASO),(PAR(2),BSO),(PAR(3),CSO),
|
||
$ (PAR(4),ALPHA),(PAR(5),BETA),(PAR(6),GAMMA)
|
||
DATA ASIG,BSIG,CSIG,DSIG,ESIG,FSIG/6*0.0/,
|
||
$ AA,AB,AC,ADD,AE,AF/6*0.0/,DETERM/1.0/,AI/6*0.0/
|
||
C-----------------------------------------------------------------------
|
||
C Select the appropriate number of parameters to calculate
|
||
C-----------------------------------------------------------------------
|
||
WC = 1
|
||
N = 2
|
||
IF (LAUE .EQ. 1) N = 6
|
||
IF (LAUE .EQ. 2) N = 4
|
||
IF (LAUE .EQ. 3) N = 3
|
||
IF (LAUE .EQ. 5) N = 1
|
||
L = N
|
||
C-----------------------------------------------------------------------
|
||
C Initialize arrays
|
||
C-----------------------------------------------------------------------
|
||
DO 110 J = 1,7
|
||
DO 100 K = 1,7
|
||
SIG(J,K) = 0.0
|
||
100 CONTINUE
|
||
SIGSQ(J) = 0.0
|
||
SIGMA(J) = 0.0
|
||
110 CONTINUE
|
||
C-----------------------------------------------------------------------
|
||
C Accumulate the sums and make the coeficients of the theta equation
|
||
C-----------------------------------------------------------------------
|
||
DO 140 II = 1,NUMD
|
||
I = II
|
||
IF (IOH(I) .NE. 0 .OR. IOK(I) .NE. 0 .OR. IOL(I) .NE. 0) THEN
|
||
M = L
|
||
CALL ETAI (AI,I)
|
||
N = M
|
||
BI = QOBS(I)
|
||
DO 130 J = 1,N
|
||
DO 120 K = 1,N
|
||
SIG(J,K) = AI(J)*AI(K)*WC + SIG(J,K)
|
||
120 CONTINUE
|
||
SIGMA(J) = SIGMA(J) + WC*BI*AI(J)
|
||
130 CONTINUE
|
||
ENDIF
|
||
140 CONTINUE
|
||
IF (N .EQ. 1) THEN
|
||
SIGMA(1) = SIGMA(1)/SIG(1,1)
|
||
SIG(1,1) = 1.0/SIG(1,1)
|
||
ELSE
|
||
NN = N - 1
|
||
DO 150 J = 1,NN
|
||
JJ = J + 1
|
||
DO 150 K = JJ,N
|
||
SIG(K,J) = SIG(J,K)
|
||
150 CONTINUE
|
||
CALL CMATIN (SIG,N,SIGMA,1,DETERM)
|
||
ENDIF
|
||
IF (DETERM .EQ. 0.0) THEN
|
||
WRITE (COUT,10000)
|
||
CALL GWRITE (ITP,' ')
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Make the sums for the esds
|
||
C-----------------------------------------------------------------------
|
||
SUMWV = 0.0
|
||
SUMW = 0.0
|
||
DO 170 II = 1,NUMD
|
||
I = II
|
||
IF (IOH(I) .NE. 0 .OR. IOK(I) .NE. 0 .OR. IOL(I) .NE. 0) THEN
|
||
T3 = 0.0
|
||
CALL ETAI(AI,I)
|
||
DO 160 K = 1,N
|
||
T3 = T3 + AI(K)*SIGMA(K)
|
||
160 CONTINUE
|
||
VI = T3 - QOBS(I)
|
||
RWGHT = 1
|
||
SUMWV = SUMWV + RWGHT*VI*VI
|
||
SUMW = SUMW + RWGHT
|
||
ENDIF
|
||
170 CONTINUE
|
||
C-----------------------------------------------------------------------
|
||
C Sigma squared
|
||
C-----------------------------------------------------------------------
|
||
DO 180 I = 1,N
|
||
SIGSQ(I) = SUMWV*SIG(I,I)/SUMW
|
||
180 CONTINUE
|
||
C-----------------------------------------------------------------------
|
||
C Calculate a, b, c, alpha, beta, gamma according to the Laue code
|
||
C
|
||
C Triclinic
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 1) THEN
|
||
AF = SIGMA(6)
|
||
AE = SIGMA(5)
|
||
ADD = SIGMA(4)
|
||
FSIG = SIGSQ(6)
|
||
ESIG = SIGSQ(5)
|
||
DSIG = SIGSQ(4)
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Monoclinic - a, b or c unique
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 2) THEN
|
||
IF (IAXIS .EQ. 1) THEN
|
||
AF = SIGMA(4)
|
||
FSIG = SIGSQ(4)
|
||
ENDIF
|
||
IF (IAXIS .EQ. 2) THEN
|
||
AE = SIGMA(4)
|
||
ESIG = SIGSQ(4)
|
||
ENDIF
|
||
IF (IAXIS .EQ. 3) THEN
|
||
ADD = SIGMA(4)
|
||
DSIG = SIGSQ(4)
|
||
ENDIF
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Triclinic, monoclinic or orthorhombic
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 1 .OR. LAUE .EQ. 2 .OR. LAUE .EQ. 3) THEN
|
||
AC = SIGMA(3)
|
||
AB = SIGMA(2)
|
||
AA = SIGMA(1)
|
||
CSIG = SIGSQ(3)
|
||
BSIG = SIGSQ(2)
|
||
ASIG = SIGSQ(1)
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Tetragonal
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 4) THEN
|
||
AA = SIGMA(1)
|
||
AB = AA
|
||
AC = SIGMA(2)
|
||
ASIG = SIGSQ(1)
|
||
BSIG = ASIG
|
||
CSIG = SIGSQ(2)
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Hexagonal and rhombohedral with hexagonal axes
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 6) THEN
|
||
AA = SIGMA(1)
|
||
AB = AA
|
||
ADD = AA/2.0
|
||
AC = SIGMA(2)
|
||
ASIG = SIGSQ(1)
|
||
BSIG = ASIG
|
||
DSIG = ASIG/2.0
|
||
CSIG = SIGSQ(2)
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Rhombohedral with rhombohedral axes
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 7) THEN
|
||
ADD = SIGMA(2)
|
||
AE = ADD
|
||
AF = ADD
|
||
DSIG = SIGSQ(2)
|
||
ESIG = DSIG
|
||
FSIG = DSIG
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Rhombohedral or cubic
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 5 .OR. LAUE .EQ. 7) THEN
|
||
AA = SIGMA(1)
|
||
AB = AA
|
||
AC = AA
|
||
ASIG = SIGSQ(1)
|
||
BSIG = ASIG
|
||
CSIG = ASIG
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Now the actual cell parameters
|
||
C-----------------------------------------------------------------------
|
||
VK = 1.0/SQRT(AA*AB*AC - AA*AF*AF - AB*AE*AE - AC*ADD*ADD +
|
||
$ 2.0*AF*AE*ADD)
|
||
ABC = AB*AC - AF*AF
|
||
AAC = AA*AC - AE*AE
|
||
AAB = AA*AB - ADD*ADD
|
||
ASO = VK*SQRT(ABC)
|
||
BSO = VK*SQRT(AAC)
|
||
CSO = VK*SQRT(AAB)
|
||
ARG1 = AE*ADD - AA*AF
|
||
ARG2 = AAC*AAB
|
||
ARG2 = SQRT(ARG2 - ARG1*ARG1)
|
||
CALL CATAN2 (ARG2,ARG1,ANSWER)
|
||
ALPHA = ANSWER*DEG
|
||
ARG1 = ADD*AF - AB*AE
|
||
ARG2 = AAB*ABC
|
||
ARG2 = SQRT(ARG2 - ARG1*ARG1)
|
||
CALL CATAN2 (ARG2,ARG1,ANSWER)
|
||
BETA = ANSWER*DEG
|
||
ARG1 = AF*AE - AC*ADD
|
||
ARG2 = ABC*AAC
|
||
ARG2 = SQRT(ARG2 - ARG1*ARG1)
|
||
CALL CATAN2 (ARG2,ARG1,ANSWER)
|
||
GAMMA = ANSWER*DEG
|
||
SALPHA = SIN(ALPHA/DEG)
|
||
SBETA = SIN(BETA/DEG)
|
||
SGAMMA = SIN(GAMMA/DEG)
|
||
C-----------------------------------------------------------------------
|
||
C Determine the standard errors using the quantities derived from the
|
||
C least-squares (AA to AF) and their variances
|
||
C
|
||
C Variances of the direct cell parameters a, b and c
|
||
C-----------------------------------------------------------------------
|
||
V2 = AA*AB*AC - AA*AF*AF - AB*AE*AE - AC*ADD*ADD + 2.0*ADD*AE*AF
|
||
V = SQRT(V2)
|
||
V3 = V2*V
|
||
TA2 = AB*AC - AF*AF
|
||
TB2 = AA*AC - AE*AE
|
||
TC2 = AA*AB - ADD*ADD
|
||
TA = SQRT(TA2)
|
||
TB = SQRT(TB2)
|
||
TC = SQRT(TC2)
|
||
C-----------------------------------------------------------------------
|
||
C Variance of a
|
||
C-----------------------------------------------------------------------
|
||
TEM = TA2*TA/(2.0*V3)
|
||
PASO = TEM*TEM*ASIG
|
||
TEM = (V2*AC - TA2*TB2)/(2.0*TA*V3)
|
||
PASO = PASO + TEM*TEM*BSIG
|
||
TEM = (V2*AB - TA2*TC2)/(2.0*TA*V3)
|
||
PASO = PASO + TEM*TEM*CSIG
|
||
TEM = TA*(AE*AF - AC*ADD)/V3
|
||
PASO = PASO + TEM*TEM*DSIG
|
||
TEM = TA*(ADD*AF - AB*AE)/V3
|
||
PASO = PASO + TEM*TEM*ESIG
|
||
TEM = (AF*V2 + TA2*(ADD*AE - AA*AF))/(TA*V3)
|
||
PASO = PASO + TEM*TEM*FSIG
|
||
PASO = SQRT(PASO)
|
||
C-----------------------------------------------------------------------
|
||
C Variance of b
|
||
C-----------------------------------------------------------------------
|
||
TEM = (AC*V2 - TB2*TA2)/(2.0*TB*V3)
|
||
PBSO = TEM*TEM*ASIG
|
||
TEM = TB2*TB/(2.0*V3)
|
||
PBSO = PBSO + TEM*TEM*BSIG
|
||
TEM = (AA*V2 - TB2*TC2)/(2.0*TB*V3)
|
||
PBSO = PBSO + TEM*TEM*CSIG
|
||
TEM = TB*(AE*AF - AC*ADD)/V3
|
||
PBSO = PBSO + TEM*TEM*DSIG
|
||
TEM = (AE*V2 + TB2*(ADD*AF - AB*AE))/(TB*V3)
|
||
PBSO = PBSO + TEM*TEM*ESIG
|
||
TEM = TB*(ADD*AE - AA*AF)/V3
|
||
PBSO = PBSO + TEM*TEM*FSIG
|
||
PBSO = SQRT(PBSO)
|
||
C-----------------------------------------------------------------------
|
||
C Variance of c
|
||
C-----------------------------------------------------------------------
|
||
TEM = (AB*V2 - TC2*TA2)/(2.0*TC*V3)
|
||
PCSO = TEM*TEM*ASIG
|
||
TEM = (AA*V2 - TC2*TB2)/(2.0*TC*V3)
|
||
PCSO = PCSO + TEM*TEM*BSIG
|
||
TEM = TC2*TC/(2.0*V3)
|
||
PCSO = PCSO + TEM*TEM*CSIG
|
||
TEM = (ADD*V2 + TC2*(AE*AF - AC*ADD))/(TC*V3)
|
||
PCSO = PCSO + TEM*TEM*DSIG
|
||
TEM = TC*(ADD*AF - AB*AE)/V3
|
||
PCSO = PCSO + TEM*TEM*ESIG
|
||
TEM = TC*(ADD*AE - AA*AF)/V3
|
||
PCSO = PCSO + TEM*TEM*FSIG
|
||
PCSO = SQRT(PCSO)
|
||
C-----------------------------------------------------------------------
|
||
C Variances of alpha, beta and gamma from their cosines
|
||
C
|
||
C Variance of alpha
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 1 .OR. LAUE .EQ. 2 .OR. LAUE .EQ. 7) THEN
|
||
BOT2 = (AA*AC - AE*AE)*(AA*AB - ADD*ADD)
|
||
BOT = SQRT(BOT2)
|
||
FAC = (AE*ADD - AA*AF)/(2.0*BOT)
|
||
TEM = (AF*BOT + FAC*(2.0*AA*AB*AC-AB*AE*AE-AC*ADD*ADD))/BOT2
|
||
PALPHA = TEM*TEM*ASIG
|
||
TEM = FAC*(AA*AA*AC - AA*AE*AE)/BOT2
|
||
PALPHA = PALPHA + TEM*TEM*BSIG
|
||
TEM = FAC*(AA*AA*AB - AA*ADD*ADD)/BOT2
|
||
PALPHA = PALPHA + TEM*TEM*CSIG
|
||
TEM = (BOT*AE - 2.0*FAC*(ADD*AE*AE - AA*AC*ADD))/BOT2
|
||
PALPHA = PALPHA + TEM*TEM*DSIG
|
||
TEM = (ADD*BOT - FAC*2.0*(ADD*ADD*AE - AA*AB*AE))/BOT2
|
||
PALPHA = PALPHA + TEM*TEM*ESIG
|
||
PALPHA = PALPHA + AA*AA*FSIG/BOT2
|
||
PALPHA = DEG*SQRT(PALPHA/(SALPHA*SALPHA))
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Variance of beta
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 1 .OR. LAUE .EQ. 2) THEN
|
||
BOT2 = (AB*AC - AF*AF)*(AA*AB - ADD*ADD)
|
||
BOT = SQRT(BOT2)
|
||
FAC = (ADD*AF - AB*AE)/(2.0*BOT)
|
||
TEM = FAC*(AB*AB*AC - AB*AF*AF)/BOT2
|
||
PBETA = TEM*TEM*ASIG
|
||
TEM = (BOT*AE + FAC*(2.0*AA*AB*AC-AA*AF*AF-AC*ADD*ADD))/BOT2
|
||
PBETA = PBETA + TEM*TEM*BSIG
|
||
TEM = FAC*(AA*AB*AB - AB*ADD*ADD)/BOT2
|
||
PBETA = PBETA + TEM*TEM*CSIG
|
||
TEM = (BOT*AF - FAC*2.0*(ADD*AF*AF - AB*AC*ADD))/BOT2
|
||
PBETA = PBETA + TEM*TEM*DSIG
|
||
PBETA = PBETA + AB*AB*ESIG/BOT2
|
||
TEM = (BOT*ADD - FAC*2.0*(AF*ADD*ADD - AA*AB*AF))/BOT2
|
||
PBETA = PBETA + TEM*TEM*FSIG
|
||
PBETA = DEG*SQRT(PBETA/(SBETA*SBETA))
|
||
PGAMMA = 0.0
|
||
C-----------------------------------------------------------------------
|
||
C Variance of gamma
|
||
C-----------------------------------------------------------------------
|
||
BOT2 = (AA*AC - AE*AE)*(AB*AC - AF*AF)
|
||
BOT = SQRT(BOT2)
|
||
FAC = (AE*AF - AC*ADD)/(2.0*BOT)
|
||
TEM = FAC*(AB*AC*AC - AC*AF*AF)/BOT2
|
||
PGAMMA = TEM*TEM*ASIG
|
||
TEM = FAC*(AA*AC*AC - AC*AE*AE)/BOT2
|
||
PGAMMA = PGAMMA + TEM*TEM*BSIG
|
||
TEM = (ADD*BOT + FAC*(2.0*AA*AB*AC-AB*AE*AE-AA*AF*AF))/BOT2
|
||
PGAMMA = PGAMMA + TEM*TEM*CSIG
|
||
PGAMMA = PGAMMA + AC*AC*DSIG/BOT2
|
||
TEM = (AF*BOT - FAC*2.0*(AE*AF*AF - AB*AC*AE))/BOT2
|
||
PGAMMA = PGAMMA + TEM*TEM*ESIG
|
||
TEM = (AE*BOT - FAC*2.0*(AE*AE*AF - AA*AC*AF))/BOT2
|
||
PGAMMA = PGAMMA + TEM*TEM*FSIG
|
||
PGAMMA = DEG*SQRT(PGAMMA/(SGAMMA*SGAMMA))
|
||
ENDIF
|
||
CALL DEVLST (PAR)
|
||
WRITE (LPT,11000) ASO, PASO, BSO, PBSO, CSO, PCSO,
|
||
$ ALPHA,PALPHA,BETA,PBETA,GAMMA,PGAMMA
|
||
RETURN
|
||
10000 FORMAT (10X,' Singular Matrix')
|
||
11000 FORMAT (/18X,' Cell Errors '/
|
||
$ 8X,'a ',F12.6,F13.7/
|
||
$ 8X,'b ',F12.6,F13.7/
|
||
$ 8X,'c ',F12.6,F13.7/
|
||
$ 8X,'Alpha ',F9.3,4X,F9.4/
|
||
$ 8X,'Beta ',F9.3,4X,F9.4/
|
||
$ 8X,'Gamma ',F9.3,4X,F9.4/)
|
||
END
|
||
C-----------------------------------------------------------------------
|
||
C Determine the AI values from h, k and l
|
||
C-----------------------------------------------------------------------
|
||
SUBROUTINE ETAI (AI,I)
|
||
INCLUDE 'COMDIF'
|
||
DIMENSION AI(6)
|
||
C-----------------------------------------------------------------------
|
||
C Triclinic
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 1) THEN
|
||
AI(6) = 2*IOK(I)*IOL(I)
|
||
AI(5) = 2*IOH(I)*IOL(I)
|
||
AI(4) = 2*IOH(I)*IOK(I)
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Monoclinic
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 2) THEN
|
||
IF (IAXIS .EQ. 1) AI(4) = 2*IOK(I)*IOL(I)
|
||
IF (IAXIS .EQ. 2) AI(4) = 2*IOH(I)*IOL(I)
|
||
IF (IAXIS .EQ. 3) AI(4) = 2*IOH(I)*IOK(I)
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Triclinic, monoclinic or orthorhombic
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .LE. 3) THEN
|
||
AI(3) = IOL(I)*IOL(I)
|
||
AI(2) = IOK(I)*IOK(I)
|
||
AI(1) = IOH(I)*IOH(I)
|
||
RETURN
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Tetragonal
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 4) THEn
|
||
AI(2) = IOL(I)*IOL(I)
|
||
AI(1) = IOH(I)*IOH(I) + IOK(I)*IOK(I)
|
||
RETURN
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Hexagonal and rhombohedral with hexagonal axes
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 6) THEN
|
||
AI(2) = IOL(I)*IOL(I)
|
||
AI(1) = IOH(I)*IOH(I) + IOK(I)*IOK(I) + IOH(I)*IOK(I)
|
||
RETURN
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Rhombohedral with rhombohedral axes
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 7)
|
||
$ AI(2) = 2*(IOH(I)*IOK(I) + IOH(I)*IOL(I) + IOK(I)*IOL(I))
|
||
C-----------------------------------------------------------------------
|
||
C Rhombohedral or cubic
|
||
C-----------------------------------------------------------------------
|
||
IF (LAUE .EQ. 5 .OR. LAUE .EQ. 7)
|
||
$ AI(1) = IOH(I)*IOH(I) + IOK(I)*IOK(I) + IOL(I)*IOL(I)
|
||
RETURN
|
||
END
|
||
C-----------------------------------------------------------------------
|
||
C List the obs and calc data in the input form
|
||
C-----------------------------------------------------------------------
|
||
SUBROUTINE DEVLST (PAR)
|
||
INCLUDE 'COMDIF'
|
||
DIMENSION PAR(6),REC(6),Q(6),QOBS(NSIZE)
|
||
EQUIVALENCE (ACOUNT(1),QOBS(1))
|
||
C-----------------------------------------------------------------------
|
||
C Make the reciprocal cell, (Int. Tab. Vol. II, p.106.
|
||
C-----------------------------------------------------------------------
|
||
PAR4 = PAR(4)/DEG
|
||
PAR5 = PAR(5)/DEG
|
||
PAR6 = PAR(6)/DEG
|
||
SUM = (PAR4 + PAR5 + PAR6)/2.0
|
||
XPRSS = SIN(SUM)*SIN(SUM - PAR4)*SIN(SUM - PAR5)*SIN(SUM - PAR6)
|
||
VOL = 2.0*PAR(1)*PAR(2)*PAR(3)*SQRT(XPRSS)
|
||
REC(1) = PAR(2)*PAR(3)*SIN(PAR4)/VOL
|
||
REC(2) = PAR(3)*PAR(1)*SIN(PAR5)/VOL
|
||
REC(3) = PAR(1)*PAR(2)*SIN(PAR6)/VOL
|
||
REC(4) = (COS(PAR5)*COS(PAR6) - COS(PAR4))/(SIN(PAR5)*SIN(PAR6))
|
||
REC(5) = (COS(PAR6)*COS(PAR4) - COS(PAR5))/(SIN(PAR6)*SIN(PAR4))
|
||
REC(6) = (COS(PAR4)*COS(PAR5) - COS(PAR6))/(SIN(PAR4)*SIN(PAR5))
|
||
C-----------------------------------------------------------------------
|
||
C Calculate the metric tensor Q
|
||
C-----------------------------------------------------------------------
|
||
Q(1) = REC(1)*REC(1)
|
||
Q(2) = REC(2)*REC(2)
|
||
Q(3) = REC(3)*REC(3)
|
||
Q(4) = REC(2)*REC(3)*REC(4)
|
||
Q(5) = REC(3)*REC(1)*REC(5)
|
||
Q(6) = REC(1)*REC(2)*REC(6)
|
||
C-----------------------------------------------------------------------
|
||
C Derive the Obs and Calc data
|
||
C-----------------------------------------------------------------------
|
||
DO 100 I = 1, NUMD
|
||
QCALC = IOH(I)*IOH(I)*Q(1) + IOK(I)*IOK(I)*Q(2) +
|
||
$ IOL(I)*IOL(I)*Q(3) + 2*IOK(I)*IOL(I)*Q(4) +
|
||
$ 2*IOL(I)*IOH(I)*Q(5) + 2*IOH(I)*IOK(I)*Q(6)
|
||
THOBS = 2.0*DEG*ACOS(SQRT(1.0 - (QOBS(I)*WAVE*WAVE/4.)))
|
||
THCAL = 2.0*DEG*ACOS(SQRT(1.0 - (QCALC *WAVE*WAVE/4.)))
|
||
100 CONTINUE
|
||
RETURN
|
||
END
|
||
C-----------------------------------------------------------------------
|
||
C Find atan(A/B) and put the answer C in the 0 to 180 range
|
||
C-----------------------------------------------------------------------
|
||
SUBROUTINE CATAN2 (A,B,C)
|
||
PI = 3.141592654
|
||
C = PI/2.0
|
||
IF (B .NE. 0) THEN
|
||
C = ATAN(ABS(A/B))
|
||
IF (B .LT. 0) C = PI - C
|
||
ENDIF
|
||
RETURN
|
||
END
|
||
C-----------------------------------------------------------------------
|
||
C Matrix inversion with accompanying solution of linear equations
|
||
C-----------------------------------------------------------------------
|
||
SUBROUTINE CMATIN (A,N,B,M,DETERM)
|
||
DIMENSION IPIVOT(7),A(7,7),B(7,1),INDEX(7,2),PIVOT(7)
|
||
EQUIVALENCE (IROW,JROW),(ICOLUM,JCOLUM),(AMAX,T,SWAP)
|
||
I = 1
|
||
EPS = .0000000001
|
||
DETERM = 1.0
|
||
DO 100 J = 1,N
|
||
IPIVOT(J) = 0
|
||
100 CONTINUE
|
||
C-----------------------------------------------------------------------
|
||
C Search for the pivot element
|
||
C-----------------------------------------------------------------------
|
||
DO 200 I = 1,N
|
||
AMAX = 0.0
|
||
DO 120 J = 1,N
|
||
IF (IPIVOT(J) .NE. 1) THEN
|
||
DO 110 K = 1,N
|
||
IF (IPIVOT(K) .GT. 1) RETURN
|
||
IF (IPIVOT(K) .LT. 1) THEN
|
||
IF (ABS(AMAX) .LT. ABS(A(J,K))) THEN
|
||
IROW = J
|
||
ICOLUM = K
|
||
AMAX = A(J,K)
|
||
ENDIF
|
||
ENDIF
|
||
110 CONTINUE
|
||
ENDIF
|
||
120 CONTINUE
|
||
IPIVOT(ICOLUM) = IPIVOT(ICOLUM) + 1
|
||
C-----------------------------------------------------------------------
|
||
C Interchange rows to put the pivot element on the main diagonal
|
||
C-----------------------------------------------------------------------
|
||
IF (IROW .NE. ICOLUM) THEN
|
||
DETERM = - DETERM
|
||
DO 130 L = 1,N
|
||
SWAP = A(IROW,L)
|
||
A(IROW,L) = A(ICOLUM,L)
|
||
A(ICOLUM,L) = SWAP
|
||
130 CONTINUE
|
||
IF (M .GT. 0) THEN
|
||
DO 140 L = 1,M
|
||
SWAP = B(IROW,L)
|
||
B(IROW,L) = B(ICOLUM,L)
|
||
B(ICOLUM,L) = SWAP
|
||
140 CONTINUE
|
||
ENDIF
|
||
ENDIF
|
||
INDEX(I,1) = IROW
|
||
INDEX(I,2) = ICOLUM
|
||
PIVOT(I) = A(ICOLUM,ICOLUM)
|
||
IF (ABS(PIVOT(I)) .LE. EPS) THEN
|
||
DETERM = 0.0
|
||
RETURN
|
||
ENDIF
|
||
DETERM = DETERM*PIVOT(I)
|
||
C-----------------------------------------------------------------------
|
||
C Divide the pivot row by the pivot element
|
||
C-----------------------------------------------------------------------
|
||
A(ICOLUM,ICOLUM) = 1.0
|
||
DO 150 L = 1,N
|
||
A(ICOLUM,L) = A(ICOLUM,L)/PIVOT(I)
|
||
150 CONTINUE
|
||
IF (M .GT. 0) THEN
|
||
DO 160 L = 1,M
|
||
B(ICOLUM,L) = B(ICOLUM,L)/PIVOT(I)
|
||
160 CONTINUE
|
||
ENDIF
|
||
C-----------------------------------------------------------------------
|
||
C Reduce non-pivot rows
|
||
C-----------------------------------------------------------------------
|
||
DO 200 L1 = 1,N
|
||
IF (L1 .NE. ICOLUM) THEN
|
||
T = A(L1,ICOLUM)
|
||
A(L1,ICOLUM) = 0.0
|
||
DO 170 L = 1,N
|
||
A(L1,L) = A(L1,L) - A(ICOLUM,L)*T
|
||
170 CONTINUE
|
||
IF (M .GT. 0) THEN
|
||
DO 180 L = 1,M
|
||
B(L1,L) = B(L1,L) - B(ICOLUM,L)*T
|
||
180 CONTINUE
|
||
ENDIF
|
||
ENDIF
|
||
200 CONTINUE
|
||
C-----------------------------------------------------------------------
|
||
C Interchange columns
|
||
C-----------------------------------------------------------------------
|
||
DO 220 I = 1,N
|
||
L = N + 1 - I
|
||
IF (INDEX(L,1) .NE. INDEX(L,2)) THEN
|
||
JROW = INDEX(L,1)
|
||
JCOLUM = INDEX(L,2)
|
||
DO 210 K = 1,N
|
||
SWAP = A(K,JROW)
|
||
A(K,JROW) = A(K,JCOLUM)
|
||
A(K,JCOLUM) = SWAP
|
||
210 CONTINUE
|
||
ENDIF
|
||
220 CONTINUE
|
||
RETURN
|
||
END
|
||
|