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