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
|
|
|