linse/sics
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

PSI sics-cvs-psi_pre-ansto

Browse Source
This commit is contained in:
koennecke
2003-06-13 00:00:00 +00:00
committed by Douglas Clowes
parent 2e3ddfb6c6
commit 3ffd0d8af4
1099 changed files with 318432 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

547
difrac/lsormt.f Normal file
View File

@@ -0,0 +1,547 @@
C-----------------------------------------------------------------------
C Linear Least Squares Derivation of Orientation Matrix
C
C The routine is called from the terminal with the MM command, or
C internally from re-orientation during data-collection with OZ.
C
C The data is obtained from ORIENT.DA beginning at record 250.
C There are 10 reflections per record
C h k l 2theta omega chi phi in the arrays
C IHK,NREFB,ILA,BCOUNT,BBGR1,BBGR2,BTIME,BPSI
C The 81st variable NBL is the number of reflections in the record.
C If NBL = 10 the block is full, if not it is the last record.
C
C-----------------------------------------------------------------------
SUBROUTINE LSORMT
INCLUDE 'COMDIF'
DIMENSION RHX(3,3),RHH(3,3),IHI(3),XOBS(3),XCNEW(3),XCOLD(3),
$ RHHI(3,3),DEL(3),RNEW(3,3),
$ IBH(10),IBK(10),IBL(10),BTHETA(10),BOMEGA(10),BCHI(10),
$ BPHI(10)
CHARACTER ANS0*1
MARK = 0
AVEANG = 0.0
DO 100 I = 1,3
DO 100 J = 1,3
ROLD(I,J) = R(I,J)/WAVE
100 CONTINUE
ANS0 = 'N'
DO 110 I = 1,3
DO 110 J = 1,3
RHX(I,J) = 0.0
RHH(I,J) = 0.0
110 CONTINUE
IF (IORNT .EQ. 1) THEN
IHSV = IH
IKSV = IK
ILSV = IL
NBSV = NB
GO TO 200
ENDIF
IF (KI .EQ. 'OP') GO TO 200
WRITE (COUT,10000)
CALL GWRITE (ITP,' ')
WRITE (COUT,13000) WAVE
CALL FREEFM (ITR)
IF (RFREE(1) .NE. 0.0) WAVE = RFREE(1)
C-----------------------------------------------------------------------
C Read from the terminal or from the Idata file
C-----------------------------------------------------------------------
WRITE (COUT,15000)
CALL YESNO ('N',ANS)
IF (ANS .EQ. 'Y') THEN
CALL TERMRD
ENDIF
C-----------------------------------------------------------------------
C Print the current MM data if wanted
C-----------------------------------------------------------------------
115 WRITE (COUT,16000)
CALL FREEFM (ITR)
IOPSHN = IFREE(1)
IF (IOPSHN .EQ. 0) IOPSHN = 3
IF (IOPSHN .EQ. 4) THEN
KI = ' '
RETURN
ENDIF
IF (IOPSHN .EQ. 1) LLIST = LPT
IF (IOPSHN .EQ. 2) THEN
LLIST = IOUNIT(10)
CALL IBMFIL ('MMDATA.DA',LLIST,80,'SU',IERR)
WRITE (LLIST,16100) WAVE
ENDIF
IF (IOPSHN .EQ. 1 .OR. IOPSHN .EQ. 2) THEN
NBLOKO = 250
120 READ (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,BOMEGA,BCHI,BPHI,
$ BPSI,NBL
NBLOKO = NBLOKO + 1
IF (NBL .NE. 0) THEN
DO 130 I = 1,NBL
IDR = BPSI(I) + 0.1
WRITE (LLIST,17000) IBH(I),IBK(I),IBL(I),BTHETA(I),
$ BOMEGA(I),BCHI(I),BPHI(I),IDR
130 CONTINUE
GO TO 120
ENDIF
IF (IOPSHN .EQ. 2) CALL IBMFIL ('MMDATA.DA',-LLIST,80,'SU',IERR)
GO TO 115
ENDIF
C-----------------------------------------------------------------------
C Routine to Delete(1) or Restore(0) reflections from LS
C-----------------------------------------------------------------------
WRITE (COUT,19000)
CALL GWRITE (ITP,' ')
140 WRITE (COUT,29000)
CALL FREEFM (ITR)
JHD = IFREE(1)
KD = IFREE(2)
LD = IFREE(3)
IDR = IFREE(4)
C-----------------------------------------------------------------------
C Find the reflection to be changed
C-----------------------------------------------------------------------
IF (JHD .NE. 0 .OR. KD .NE. 0 .OR. LD .NE. 0) THEN
NBLOKO = 250
150 READ (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,BOMEGA,BCHI,
$ BPHI,BPSI,NBL
NBLOKO = NBLOKO + 1
IF (NBL .EQ. 0) GO TO 140
C-----------------------------------------------------------------------
C Find the reflection, change its status, write back and get the next
C-----------------------------------------------------------------------
DO 160 NB = 1,NBL
IF (IBH(NB) .EQ. JHD .AND. IBK(NB) .EQ. KD .AND.
$ IBL(NB) .EQ. LD) THEN
BPSI(NB) = IDR
NBLOKO = NBLOKO-1
WRITE (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,BOMEGA,BCHI,
$ BPHI,BPSI,NBL
NBLOKO = NBLOKO + 1
GO TO 140
ENDIF
160 CONTINUE
GO TO 150
ENDIF
C-----------------------------------------------------------------------
C Insert new reflections
C-----------------------------------------------------------------------
WRITE (COUT,23000)
CALL YESNO ('N',ANS)
IF (ANS .EQ. 'Y') THEN
NBLOKO = 250
170 READ (ISD,REC=NBLOKO) (JUNK,I = 1,80),NBL
NBLOKO = NBLOKO + 1
IF (NBL .NE. 0) GO TO 170
NBLOKO = NBLOKO - 1
WRITE (COUT,24000) NBLOKO
CALL GWRITE (ITP,' ')
CALL TERMRD
ENDIF
C-----------------------------------------------------------------------
C Get the zero values of omega and chi and possibly use them
C-----------------------------------------------------------------------
CALL WCZERO (ZOMEGA,ZCHI)
WRITE (COUT,24100)
CALL YESNO ('Y',ANS0)
IF (ANS0 .EQ. 'Y') THEN
DOMEGA = DOMEGA + ZOMEGA
DCHI = DCHI + ZCHI
WRITE (COUT,24200) DOMEGA,DCHI
CALL GWRITE (ITP,' ')
ENDIF
C-----------------------------------------------------------------------
C Start of the Least Squares Procedure
C Data is read from the file twice
C MARK = 0 when making and solving the normal equations
C MARK = -1 when forming the deltas for the e.s.d.'s
C-----------------------------------------------------------------------
200 NRF = 0
NBLOKO = 250
IF (MARK .EQ. -1) WRITE (LPT,25000)
210 READ (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,BOMEGA,BCHI,
$ BPHI,BPSI,NBL
NBLOKO = NBLOKO + 1
IF (NBL .NE. 0) THEN
DO 250 NB = 1,NBL
IF (BPSI(NB) .EQ. 0.0) THEN
NRF = NRF + 1
IHI(1) = IBH(NB)
IHI(2) = IBK(NB)
IHI(3) = IBL(NB)
TH = 2.0*SIN(0.5*BTHETA(NB)/DEG)/WAVE
IF (ANS0 .EQ. 'Y') THEN
BOMEGA(NB) = BOMEGA(NB) - ZOMEGA
BCHI(NB) = BCHI(NB) - ZCHI
ENDIF
CCH = COS(BCHI(NB)/DEG)
SCH = SIN(BCHI(NB)/DEG)
COM = COS(BOMEGA(NB)/DEG)
SOM = SIN(BOMEGA(NB)/DEG)
CPH = COS((BPHI(NB))/DEG)
SPH = SIN((BPHI(NB))/DEG)
XOBS(1) = TH*(CCH*CPH*COM - SOM*SPH)
XOBS(2) = TH*(CCH*SPH*COM + SOM*CPH)
XOBS(3) = TH*SCH*COM
C-----------------------------------------------------------------------
C MARK = 0 Form the RHH and RHX elements
C MARK = -1 For IORNT = 1, i.e. re-orientation,
C form the calcd values of X with the old and new matrices;
C then form the differences and hence the angular deviation.
C For IORNT = 0, i.e. normal MM,
C form the calcd value of X at the observed phi,
C and then the differences between the obs and calcd X and
C then the angular deviation.
C-----------------------------------------------------------------------
IF (MARK .EQ. 0) THEN
DO 220 I = 1,3
DO 220 J = 1,3
RHH(I,J) = RHH(I,J) + IHI(I)*IHI(J)
RHX(I,J) = RHX(I,J) + IHI(I)*XOBS(J)
220 CONTINUE
ELSE
IH = IHI(1)
IK = IHI(2)
IL = IHI(3)
DO 230 I = 1,3
XCOLD(I) = 0.0
XCNEW(I) = 0.0
DO 230 J = 1,3
XCNEW(I) = XCNEW(I) + RNEW(I,J)*IHI(J)
XCOLD(I) = XCOLD(I) + ROLD(I,J)*IHI(J)
230 CONTINUE
PHI = BPHI(NB)
CALL ANGPHI (XCNEW)
C-----------------------------------------------------------------------
C Work out the sums for the average angular deviation
C Keep or re-orientation, keep the new R matrix if AVEANG .gt. REOTOL
C-----------------------------------------------------------------------
TOP = 0.0
BOT = 0.0
DO 240 I = 1,3
IF (IORNT .EQ. 1) THEN
DELTA = XCOLD(I) - XCNEW(I)
ELSE
DELTA = XOBS(I) - XCNEW(I)
ENDIF
DELTA = DELTA*DELTA
DEL(I) = DEL(I) + DELTA
TOP = TOP + DELTA
BOT = BOT + XCNEW(I)*XCNEW(I)
240 CONTINUE
ANGLE = DEG*RATAN2(SQRT(TOP),SQRT(BOT))
AVEANG = AVEANG + ANGLE
C-----------------------------------------------------------------------
C Write the results for this reflection
C-----------------------------------------------------------------------
OMG = OMEGA
IF (OMG .GT. 359.995) OMG = 0.0
WRITE (LPT,26000)
$ IHI,BTHETA(NB),BOMEGA(NB),BCHI(NB),BPHI(NB),
$ THETA,OMG,CHI,PHI,ANGLE
ENDIF
ENDIF
250 CONTINUE
IF (ANS0 .EQ. 'Y' .AND. MARK .EQ. -1) THEN
NBLOKO = NBLOKO - 1
WRITE (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,BOMEGA,BCHI,
$ BPHI,BPSI,NBL
NBLOKO = NBLOKO + 1
ENDIF
GO TO 210
ENDIF
C-----------------------------------------------------------------------
C Solve for R matrix elements
C-----------------------------------------------------------------------
IF (MARK .EQ. 0) THEN
CALL MATRIX (RHH,RHHI,RHHI,RHHI,'INVERT')
DO 270 I = 1,3
DO 260 J = 1,3
R(I,J) = 0.0
DO 260 KK = 1,3
R(I,J) = R(I,J) + RHHI(J,KK)*RHX(KK,I)
260 CONTINUE
270 CONTINUE
DET = R(1,1)*(R(2,2)*R(3,3) - R(2,3)*R(3,2)) -
$ R(1,2)*(R(2,1)*R(3,3) - R(2,3)*R(3,1)) +
$ R(1,3)*(R(2,1)*R(3,2) - R(2,2)*R(3,1))
IF (NRC*DET .GT. 0) THEN
WRITE (LPT,27000) NRF
ELSE
WRITE (LPT,28000) NRF
ENDIF
C-----------------------------------------------------------------------
C Extract the real and reciprocal cell parameters
C-----------------------------------------------------------------------
CALL GETPAR
C-----------------------------------------------------------------------
C Clear the esds array
C-----------------------------------------------------------------------
DO 290 J = 1,3
DEL(J) = 0.0
290 CONTINUE
C-----------------------------------------------------------------------
C Store new R matrix times Wavelength
C-----------------------------------------------------------------------
DO 300 I = 1,3
DO 300 J = 1,3
RNEW(I,J) = R(I,J)
R(I,J) = R(I,J)*WAVE
300 CONTINUE
MARK = MARK - 1
GO TO 200
ENDIF
C-----------------------------------------------------------------------
C S.D.'s of the R matrix elements from the diagonal elements of the
C inverted (RHHI) matrix.
C Print the matrix and its e.s.ds
C-----------------------------------------------------------------------
DO 310 I = 1,3
DO 310 J = 1,3
SR(I,J) = SQRT(RHHI(J,J)*DEL(I)/(NRF - 3))
310 CONTINUE
WRITE (LPT,28100) ((RNEW(I,J),J=1,3),(SR(I,J),J=1,3),I = 1,3)
C-----------------------------------------------------------------------
C Call CELLSD to find the s.d.'s of the cell parameters
C-----------------------------------------------------------------------
CALL CELLSD
IF (IORNT .EQ. 1) THEN
AVEANG = AVEANG/NRF
IH = IHSV
IK = IKSV
IL = ILSV
NB = NBSV
IF (AVEANG .LT. REOTOL) THEN
WRITE (COUT,30000) AVEANG
CALL GWRITE (ITP,' ')
WRITE (LPT,30000) AVEANG
DO 320 I = 1,3
DO 320 J = 1,3
R(I,J) = ROLD(I,J)*WAVE
320 CONTINUE
RETURN
ELSE
WRITE (COUT,31000) AVEANG
CALL GWRITE (ITP,' ')
WRITE (LPT,31000) AVEANG
ENDIF
ENDIF
C-----------------------------------------------------------------------
C Calculate the SINABS matrix
C-----------------------------------------------------------------------
IF (KI .NE. 'MM' .AND. KI .NE. 'OP') THEN
ISYS = 1
CALL SINMAT
ENDIF
RETURN
10000 FORMAT (20X,' Least Squares Orientation Matrix')
13000 FORMAT (' Reflection data can be on file or from the terminal.'/
$ ' Wavelength (',F7.5,') ',$)
14000 FORMAT (/10X,' Reflections in the Alignment List')
15000 FORMAT (' Read the data from the terminal (N) ? ',$)
16000 FORMAT (' The following options are available :--'/
$ ' 1. List the MM data on the printer for editting, or'/
$ ' 2. Write the MM data to the ASCII file MMDATA.DA, or'/
$ ' 3. Proceed to the next step, or'/
$ ' 4. Exit from MM.'/
$ ' Which do you want (3) ? ',$)
16100 FORMAT (' MM data from DIFRAC'/F10.6)
17000 FORMAT (3I4,4F8.3,I3)
19000 FORMAT (' Reflections may be deleted or restored to the list',
$ ' by typing :--',/,
$ ' h,k,l,1 for Delete or h,k,l,0 for Restore (End)')
23000 FORMAT (' Do you wish to insert reflections (N) ? ',$)
24100 FORMAT (' Do you want to include these zero values (Y) ? ',$)
24200 FORMAT (' The new zeroes for Omega and Chi are',2F7.3)
24000 FORMAT (' First non-written record: ',I4)
25000 FORMAT (/,22X,'Observed',22X,'Calculated',10X,'Angular'/
$ ' h k l 2Theta Omega Chi Phi ',
$ ' 2Theta Omega Chi Phi ',
$ 'Deviation')
26000 FORMAT (3I4,4F7.2,2X,4F7.2,F8.3)
26100 FORMAT (I4,2X,3I4,4F8.3)
27000 FORMAT (/' Right-handed Orientation Matrix from ',I4,
$ ' Reflections')
28000 FORMAT (/' Left-handed Orientation Matrix from ',I4,
$ ' Reflections')
28100 FORMAT (/9X,'Orientation Matrix',30X,'E.S.Ds'/(3F12.8,6X,3F12.8))
29000 FORMAT (' > ',$)
30000 FORMAT (' The angular deviation is',F6.3,'. The old matrix will',
$ ' be retained.')
31000 FORMAT (' The angular deviation is',F6.3,'. The new matrix will',
$ ' be used.')
END
C-----------------------------------------------------------------------
C Routine to find the zeroes of omega and chi from the alignment data
C ZOMEGA is the average value of omega;
C ZCHI is half the average value of chi for pairs of +++/--- reflns
C-----------------------------------------------------------------------
SUBROUTINE WCZERO (ZOMEGA,ZCHI)
INCLUDE 'COMDIF'
DIMENSION IBH(10),IBK(10),IBL(10),BTHETA(10),BOMEGA(10),BCHI(10),
$ BPHI(10)
C EQUIVALENCE (IHK(1),IBH(1)),(NREFB(1),IBK(1)),(ILA(1),IBL(1)),
C $ (BCOUNT(1),BTHETA(1)),(BBGR1(1),BOMEGA(1)),
C $ (BBGR2(1),BCHI(1)),(BTIME(1),BPHI(1))
SUMOME = 0.0
SUMCHI = 0.0
NOMEGA = 0
NCHI = 0
NBLOKO = 250
IH1 = 0
IK1 = 0
IL1 = 0
CHI1 = 999.0
100 READ (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,BOMEGA,BCHI,BPHI,
$ BPSI,NBL
IF (NBL .NE. 0) THEN
DO 110 NB = 1,NBL
IF (BPSI(NB) .EQ. 0.0) THEN
WOMEGA = BOMEGA(NB)
IF (WOMEGA .GT. 180.0) WOMEGA = WOMEGA - 360.0
SUMOME = SUMOME + WOMEGA
NOMEGA = NOMEGA + 1
IH2 = IBH(NB)
IK2 = IBK(NB)
IL2 = IBL(NB)
CHI2 = BCHI(NB)
IF (IH1 .EQ. -IH2 .AND.
$ IK1 .EQ. -IK2 .AND.
$ IL1 .EQ. -IL2) THEN
CHI12 = CHI1 + CHI2
IF (CHI12 .GT. 350.0) CHI12 = CHI12 - 360.0
IF (CHI12 .GT. 350.0) CHI12 = CHI12 - 360.0
SUMCHI = SUMCHI + CHI12
NCHI = NCHI + 1
ENDIF
IH1 = IH2
IK1 = IK2
IL1 = IL2
CHI1 = CHI2
ENDIF
110 CONTINUE
NBLOKO = NBLOKO + 1
GO TO 100
ENDIF
ZOMEGA = SUMOME/NOMEGA
ZCHI = 0.0
IF (NCHI .NE. 0) ZCHI = 0.5*SUMCHI/NCHI
WRITE (COUT,10000) ZOMEGA,NOMEGA,ZCHI,NCHI
CALL GWRITE (ITP,' ')
WRITE (LPT,10000) ZOMEGA,NOMEGA,ZCHI,NCHI
RETURN
10000 FORMAT (' Omega(0) is',F7.3,' from',I4,' reflections.'/
$ ' Chi(0) is',F7.3,' from',I4,' +/- pairs.')
END
C-----------------------------------------------------------------------
C Get reflection angle input from the terminal
C-----------------------------------------------------------------------
SUBROUTINE TERMRD
INCLUDE 'COMDIF'
DIMENSION IBH(10),IBK(10),IBL(10),BTHETA(10),BOMEGA(10),BCHI(10),
$ BPHI(10)
C EQUIVALENCE (IHK(1),IBH(1)),(NREFB(1),IBK(1)),(ILA(1),IBL(1)),
C $ (BCOUNT(1),BTHETA(1)),(BBGR1(1),BOMEGA(1)),
C $ (BBGR2(1),BCHI(1)),(BTIME(1),BPHI(1))
NBLOKO = 250
WRITE (COUT,10000)
CALL GWRITE (ITP,' ')
100 NBL = 0
DO 110 I = 1,10
WRITE (COUT,11000)
CALL FREEFM (ITR)
IH = IFREE(1)
IK = IFREE(2)
IL = IFREE(3)
IF (IH .EQ. 0 .AND. IK .EQ. 0 .AND. IL .EQ. 0) GO TO 120
IBH(I) = IH
IBK(I) = IK
IBL(I) = IL
BTHETA(I) = RFREE(4)
BOMEGA(I) = RFREE(5)
BCHI(I) = RFREE(6)
BPHI(I) = RFREE(7)
BPSI(I) = 0.
NBL = NBL + 1
110 CONTINUE
120 WRITE (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,BOMEGA,BCHI,
$ BPHI,BPSI,NBL
NBLOKO = NBLOKO + 1
IF (NBL .EQ. 10) GO TO 100
NBL = 0
WRITE (ISD,REC=NBLOKO) IBH,IBK,IBL,BTHETA,BOMEGA,BCHI,
$ BPHI,BPSI,NBL
RETURN
10000 FORMAT (' Type h,k,l,2theta,omega,chi,phi for each refln. (End)')
11000 FORMAT (' > ',$)
END
C-----------------------------------------------------------------------
C Calculate chi and omega for a given phi value
C-----------------------------------------------------------------------
SUBROUTINE ANGPHI (XC)
INCLUDE 'COMDIF'
DIMENSION XC(3)
CP = COS(PHI/DEG)
SP = SIN(PHI/DEG)
TOPC = XC(3)
BOTC = CP*XC(1) + SP*XC(2)
CHI = DEG*RATAN2(TOPC,BOTC)
TOPO = CP*XC(2) - SP*XC(1)
IF (CHI .EQ. 0.0) THEN
OMEGA = 0.0
ELSE
BOTO = XC(3)/SIN(CHI/DEG)
OMEGA = DEG*RATAN2(TOPO,BOTO)
ENDIF
TH = 0.5*WAVE*SQRT(XC(1)*XC(1) + XC(2)*XC(2) + XC(3)*XC(3))
THETA = 2.0*DEG*ATAN(TH/SQRT(1.0 - TH*TH))
TH = 2.0*TH/WAVE
CC = COS(CHI/DEG)
SC = SIN(CHI/DEG)
CO = COS(OMEGA/DEG)
SO = SIN(OMEGA/DEG)
XC(1) = TH*(CC*CP*CO - SO*SP)
XC(2) = TH*(CC*SP*CO + SO*CP)
XC(3) = TH*SC*CO
CALL MOD360 (CHI)
CALL MOD360 (OMEGA)
RETURN
END
C-----------------------------------------------------------------------
C Get around stupid Microsoft compiler problems
C-----------------------------------------------------------------------
FUNCTION RATAN2 (TOP,BOT)
RA = 57.2958
IF (BOT .EQ. 0) THEN
X = 90.0/RA
IF (TOP .LT. 0.0) X = 270.0/RA
ELSE
X = ATAN2(TOP,BOT)
ENDIF
RATAN2 = X
RETURN
END
C-----------------------------------------------------------------------
C Extract the real and reciprocal cell parameters
C-----------------------------------------------------------------------
SUBROUTINE GETPAR
INCLUDE 'COMDIF'
DIMENSION RT(3,3)
DO 100 I = 1,3
DO 100 J = 1,3
RT(I,J) = R(J,I)
100 CONTINUE
CALL MATRIX (RT,R,GI,GI,'MATMUL')
C-----------------------------------------------------------------------
C Use CANGS array for reciprocal angles
C-----------------------------------------------------------------------
CALL PARAMS (GI,APS,CANGS)
C-----------------------------------------------------------------------
C Use RT array for metric tensor G
C-----------------------------------------------------------------------
CALL MATRIX (GI,RT,RT,RT,'INVERT')
C-----------------------------------------------------------------------
C Use CANG array for real angles
C-----------------------------------------------------------------------
CALL PARAMS (RT,AP,CANG)
RETURN
END