**==aa0001.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C********************************* C** Program: LCHK ** C********************************* C C C LCHK: Analyse various characteristics of Laue pictures C C AUTHOR : Original PERQ version - Mike Elder 1986 C VAX/VMS version - Ian Clifton Mar 1988 C S E R C Daresbury WA4 4AD UK. C Increased portability etc. - John Campbell 1989, 1991 C UPDATES : C INCLUDE 'implicit_none.finc' C COMMON /AMAT/A,AMATX COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT COMMON /SPOTS/NR COMMON /SYS/ISYS,KSYS COMMON /XYZ/IHKL REAL CF,DMIN,SPNDLE,WVM,WVS,WVX INTEGER ICENT,ILGP,ISYS,MBEAM,MROT,MSYS,NR,NUMA,NUMC REAL A(3,3),AMATX(3,3),CCELL(6),CELL(6),PHIS(3) INTEGER IHKL(20000),KSYS(3) C INCLUDE 'chanel.finc' C C====== FUNCTION TYPES: C INTEGER LUNSTI, LUNSTO C C====== LOCALS: C LOGICAL QUIT REAL DSTMAX INTEGER I,ICHANG,IENT,IFLAG,J CHARACTER*1 CHAR INTEGER IFAIL C C------------------------------------------------------------------------------- C C ITERM = LUNSTI(0) IPRINT = LUNSTO(0) IUNIT7 = 3 CALL TRWSET(ITERM,IPRINT,0,0,0) CALL QGINIT(0,IPRINT,0,ITERM,IPRINT) C C====== Log file (uses IUNIT6, and units 1&2) C IFAIL = 0 CALL CCPDPN(-1,'LCHKLOG','NEW','F',0,IFAIL) CALL CCPDPN(-2,'LCHKJUNK','NEW','F',0,IFAIL) C IENT = 0 20 CALL QSCLRA CALL QSSCRL(20,24) CALL QSREVV CALL QSSETC(1,1) CALL TSTRNG( + 'LCHK - Last update 17/10/91 Authors: M.Elder, I.J. Clifton' + ) WRITE (1,*) + ' LCHK - Last update 17/10/91 Authors: M.Elder, I.J.Clifton' CALL QSNORV CALL SCREEN(IENT,ICHANG) CALL QSCLRA CALL QSSCRL(1,24) C C====== Calculate the master matrix, a, from the cell and phi angles. C C C====== Apply the scale factor to the reciprocal lattice lengths first. C CELL(1) = CELL(1)*WVS CELL(2) = CELL(2)*WVS CELL(3) = CELL(3)*WVS J = 0 DO 30 I = 1 , 3 IF (CELL(I).NE.0.) THEN J = J + 1 CCELL(J) = CELL(I) END IF 30 CONTINUE CALL ACALC CELL(1) = CELL(1)/WVS CELL(2) = CELL(2)/WVS CELL(3) = CELL(3)/WVS C C====== ICHANG denotes a change in a parameter which necessitates C====== recalculating the reciprocal lattice coords. C IF (ICHANG.NE.0) THEN C DSTMAX = 1./DMIN C C====== Put the indices of possible reflections into /xyz/ C WRITE (IPRINT,6010) CALL HKLGEN(DSTMAX,IFLAG) WRITE (IPRINT,6030) NR WRITE (IPRINT,6020) NR END IF C C====== Execute the analysis commands C IF (IENT.EQ.2) THEN C C====== Analyse the effect of stepping phi-y and phi-z C CALL SEARCH ELSE IF (IENT.EQ.3) THEN C C====== Analyse the wavelength variations for symmetry equivs C CALL WAVVAR ELSE IF (IENT.EQ.4) THEN C C====== Analyse the effect of changing crystal-film distance C CALL VARYCF ELSE IF (IENT.EQ.5) THEN C C====== Analyse the fraction of reflections recorded against resolution C CALL RESOLN ELSE IF (IENT.EQ.6) THEN C C====== Analyse the effect of rotating about an offset spindle C CALL MUROT ELSE IF (IENT.EQ.7) THEN C C====== Analyse the effect of rotating about µ and U C CALL TWOROT ELSE C C====== Analyse the effect of rotating about the spindle in this orientation C CALL ROTSPI END IF C CALL TPROMP('Type to continue...') CALL TREAD(CHAR,QUIT) IF (CHAR.EQ.'Q' .OR. CHAR.EQ.'q') STOP GO TO 20 C C====== Format statements C 6010 FORMAT ( + ' Calculate and store the symmetry unique set of reflections' + ) 6020 FORMAT (I6,' Reflections stored') 6030 FORMAT (///I11) END **==acalc.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** ACALC ** C*************************** C C SUBROUTINE ACALC C INCLUDE 'implicit_none.finc' C C Purpose: Calculate the master matrix, a, from the orientation matrix C premultiplied by the phi matrices. Code from IDXREF. C COMMON /AMAT/A,AA COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN REAL CF,DMIN,SPNDLE,WVM,WVS,WVX INTEGER NUMA,NUMC REAL A(3,3),AA(3,3),CCELL(6),CELL(6),PHIS(3) C C====== LOCALS: C REAL CPX,CPY,CPZ,RAPI,SPX,SPY,SPZ REAL AC1(3,3),AC2(3,3),AC3(3,3),AC4(3,3),AC5(3,3),AC6(3,3),C(3,3), + CX(3,3),CY(3,3),CZ(3,3),PC(3,3) C C------------------------------------------------------------------------------- C C====== Pass the packed cell parameters in ccell to the mats routine C====== of IDXREF. The AC arrays hold the (unnecessary) derivative matrices. C CALL CLEAR(AA) CALL MATS(AA,AC1,AC2,AC3,AC4,AC5,AC6,CCELL) C C====== Define all rotation matrices C RAPI = 3.1415926536/180. CPX = COS(PHIS(1)*RAPI) SPX = SIN(PHIS(1)*RAPI) CPY = COS(PHIS(2)*RAPI) SPY = SIN(PHIS(2)*RAPI) CPZ = COS(PHIS(3)*RAPI) SPZ = SIN(PHIS(3)*RAPI) C CALL CLEAR(CX) CX(1,1) = 1.0 CX(2,2) = CPX CX(2,3) = -SPX CX(3,2) = SPX CX(3,3) = CPX C CALL CLEAR(CY) CY(1,1) = CPY CY(1,3) = SPY CY(2,2) = 1.0 CY(3,1) = -SPY CY(3,3) = CPY C CALL CLEAR(CZ) CZ(1,1) = CPZ CZ(1,2) = -SPZ CZ(2,1) = SPZ CZ(2,2) = CPZ CZ(3,3) = 1.0 C C====== Matrix a is the transformation from indices to orthogonal C====== coordinates C CALL MATM33(CY,CX,C) CALL MATM33(CZ,C,PC) CALL MATM33(PC,AA,A) C RETURN END **==celltest.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** CELLTEST ** C*************************** C C SUBROUTINE CELLTEST(MSYS,ILGP,ICENT,CELL,LABEL,IOK) C INCLUDE 'implicit_none.finc' C C Purpose: Check to see if the reciprocal cell parameters are consistent with C the lattice type. If inconsistent (eg length of zero in ortho) then C give warning message and return iok=1 so user can adjust them. If C ok then leave iok and return. To encourage users to understand the C principle, set redundant parameters to zero (eg b* and c* in cubic). C Note that program assumes b unique for monoclinic, c unique for hex C and tetragonal. C INTEGER ICENT,ILGP,IOK,MSYS REAL CELL(6) CHARACTER LABEL(7,5)*6 C INCLUDE 'chanel.finc' C C====== LOCALS: C INTEGER I C C------------------------------------------------------------------------------- C IF (MSYS.EQ.2) THEN C C====== Monoclinic C DO 20 I = 1 , 3 IF (CELL(I).EQ.0.0) THEN WRITE (IPRINT,6010) GO TO 150 END IF 20 CONTINUE IF (CELL(5).EQ.0.0) THEN WRITE (IPRINT,6020) GO TO 150 END IF DO 30 I = 4 , 6 , 2 IF (CELL(I).NE.0.0) THEN IF (CELL(I).NE.90.) WRITE (IPRINT,6030) CELL(I) = 0.0 END IF 30 CONTINUE IF (ILGP.NE.1) THEN WRITE (IPRINT,6040) GO TO 150 END IF ELSE IF (MSYS.EQ.3) THEN C C====== Orthorhombic C DO 40 I = 1 , 3 IF (CELL(I).EQ.0.0) THEN WRITE (IPRINT,6050) GO TO 150 END IF 40 CONTINUE DO 50 I = 4 , 6 IF (CELL(I).NE.0.0) THEN IF (CELL(I).NE.90.) WRITE (IPRINT,6060) CELL(I) = 0.0 END IF 50 CONTINUE IF (ILGP.NE.2) THEN WRITE (IPRINT,6070) GO TO 150 END IF ELSE IF (MSYS.EQ.4) THEN C C======== Tetragonal C DO 60 I = 1 , 3 , 2 IF (CELL(I).EQ.0.0) THEN WRITE (IPRINT,6080) GO TO 150 END IF 60 CONTINUE IF (CELL(2).NE.0.0) THEN IF (CELL(2).NE.CELL(1)) WRITE (IPRINT,6090) CELL(2) = 0.0 END IF DO 70 I = 4 , 6 IF (CELL(I).NE.0.0) THEN IF (CELL(I).NE.90.) WRITE (IPRINT,6100) CELL(I) = 0.0 END IF 70 CONTINUE IF (ILGP.NE.3 .AND. ILGP.NE.4) THEN WRITE (IPRINT,6110) GO TO 150 END IF ELSE IF (MSYS.EQ.5) THEN C C======== Hexagonal C DO 80 I = 1 , 3 , 2 IF (CELL(I).EQ.0.0) THEN WRITE (IPRINT,6120) GO TO 150 END IF 80 CONTINUE IF (CELL(2).NE.0.0) THEN IF (CELL(2).NE.CELL(1)) WRITE (IPRINT,6130) CELL(2) = 0.0 END IF DO 90 I = 4 , 5 IF (CELL(I).NE.0.0) THEN IF (CELL(I).NE.90.) WRITE (IPRINT,6140) CELL(I) = 0.0 END IF 90 CONTINUE IF (CELL(6).NE.0.) THEN IF (CELL(6).NE.60.) WRITE (IPRINT,6130) CELL(6) = 0. END IF IF (ILGP.LT.5 .OR. ILGP.GT.9) THEN WRITE (IPRINT,6150) GO TO 150 END IF ELSE IF (MSYS.EQ.6) THEN C C====== Rhombohedral C IF (CELL(1).EQ.0.0) THEN WRITE (IPRINT,6160) GO TO 150 END IF DO 100 I = 2 , 3 IF (CELL(I).NE.0.0) THEN IF (CELL(I).NE.CELL(1)) WRITE (IPRINT,6170) CELL(I) = 0.0 END IF 100 CONTINUE IF (CELL(4).EQ.0.0) THEN WRITE (IPRINT,6180) GO TO 150 END IF DO 110 I = 5 , 6 IF (CELL(I).NE.0.0) THEN IF (CELL(I).NE.CELL(4)) WRITE (IPRINT,6190) CELL(I) = 0.0 END IF 110 CONTINUE IF (ILGP.LT.7 .OR. ILGP.GT.9) THEN WRITE (IPRINT,6200) GO TO 150 END IF ELSE IF (MSYS.EQ.7) THEN C C====== Cubic C IF (CELL(1).EQ.0.0) THEN WRITE (IPRINT,6210) GO TO 150 END IF DO 120 I = 2 , 3 IF (CELL(I).NE.0.0) THEN IF (CELL(I).NE.CELL(1)) WRITE (IPRINT,6220) CELL(I) = 0.0 END IF 120 CONTINUE DO 130 I = 4 , 6 IF (CELL(I).NE.0.0) THEN IF (CELL(I).NE.90.) WRITE (IPRINT,6230) CELL(I) = 0.0 END IF 130 CONTINUE IF (ILGP.NE.10 .AND. ILGP.NE.11) THEN WRITE (IPRINT,6240) GO TO 150 END IF ELSE C C====== Triclinic C DO 140 I = 1 , 6 IF (CELL(I).EQ.0.0) THEN WRITE (IPRINT,6250) GO TO 150 END IF 140 CONTINUE IF (ILGP.NE.0) THEN WRITE (IPRINT,6260) GO TO 150 END IF END IF C IOK = 0 RETURN C 150 IOK = 1 RETURN C C====== Format statements C 6010 FORMAT (' Zero cell length in monoclinic case') 6020 FORMAT (' Zero beta-star in monoclinic case') 6030 FORMAT (' Monoclinic case must have unique axis b') 6040 FORMAT (' Laue group should be 2/m for monoclinic') 6050 FORMAT (' Zero cell length in orthorhombic case') 6060 FORMAT (' No angles needed for orthorhombic case') 6070 FORMAT (' Laue group should be mmm for orthorhombic') 6080 FORMAT (' Zero cell length in tetragonal case') 6090 FORMAT (' C must be unique axis for tetragonal case') 6100 FORMAT (' No angles needed for tetragonal case') 6110 FORMAT (' Laue group should be 4/m or 4/mmm for tetragonal') 6120 FORMAT (' Zero cell length in hexagonal case') 6130 FORMAT (' C must be unique axis for hexagonal case') 6140 FORMAT (' No angles needed for hexagonal case') 6150 FORMAT (' Laue group should be 6/m or 6/mmm for hexagonal'/ + ' 3-bar, 3m1 or 31m for trigonal on hexagonal axes') 6160 FORMAT (' Zero cell length in rhombohedral case') 6170 FORMAT (' Only one cell length needed for rhombohedral case') 6180 FORMAT (' Zero cell angle in rhombohedral case') 6190 FORMAT (' Only one angle needed for rhombohedral case') 6200 FORMAT ( +' Laue group should be 3-bar, 3m1 or 31m for trigonal'' on rhomboh +edral axes') 6210 FORMAT (' Zero cell length in cubic case') 6220 FORMAT (' Only one cell length needed for cubic case') 6230 FORMAT (' No angles needed for cubic case') 6240 FORMAT (' Laue group should be m3 or m3m for cubic') 6250 FORMAT (' Zero cell parameter in triclinic case') 6260 FORMAT (' Laue group should be 1-bar for triclinic') END **==clear.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** CLEAR ** C*************************** C C SUBROUTINE CLEAR(A) C INCLUDE 'implicit_none.finc' C REAL A(9) INTEGER I C C------------------------------------------------------------------------------- C DO 20 I = 1 , 9 A(I) = 0.0 20 CONTINUE RETURN END **==gensym.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** GENSYM ** C*************************** C C SUBROUTINE GENSYM(NSYM,NREF,IH,IK,IL) C INCLUDE 'implicit_none.finc' C COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT INTEGER IH,IK,IL,NSYM INTEGER NREF(3,48) INTEGER ICENT,ILGP,MBEAM,MROT,MSYS C C====== LOCALS: C INTEGER I,J C C------------------------------------------------------------------------------- C NREF(1,1) = IH NREF(2,1) = IK NREF(3,1) = IL NSYM = 1 IF (ILGP.NE.0) THEN C C====== Deal with symmetries higher than triclinic C IF (ILGP.EQ.3) GO TO 30 IF (ILGP.EQ.4) GO TO 30 IF (ILGP.EQ.5) GO TO 60 IF (ILGP.EQ.6) GO TO 60 IF (ILGP.EQ.7) GO TO 90 IF (ILGP.EQ.8) GO TO 90 IF (ILGP.EQ.9) GO TO 90 IF (ILGP.EQ.10) GO TO 120 IF (ILGP.EQ.11) GO TO 120 C C====== 2/m C IF (IK.NE.0) THEN IF (IL.NE.0 .OR. IH.NE.0) THEN NREF(1,2) = IH NREF(2,2) = -IK NREF(3,2) = IL NSYM = 2 END IF END IF IF (ILGP.NE.1) THEN C C====== mmm C IF (IH.NE.0 .AND. IL.NE.0) THEN DO 20 I = 1 , NSYM NREF(1,I+NSYM) = NREF(1,I) NREF(2,I+NSYM) = NREF(2,I) NREF(3,I+NSYM) = -NREF(3,I) 20 CONTINUE NSYM = 2*NSYM END IF END IF END IF GO TO 160 C C====== 4/m C 30 IF (IH.NE.0 .OR. IK.NE.0) THEN NREF(1,2) = IK NREF(2,2) = -IH NREF(3,2) = IL NSYM = 2 IF (IL.NE.0) THEN DO 40 I = 1 , NSYM NREF(1,I+NSYM) = NREF(1,I) NREF(2,I+NSYM) = NREF(2,I) NREF(3,I+NSYM) = -NREF(3,I) 40 CONTINUE NSYM = NSYM*2 END IF END IF IF (ILGP.NE.3) THEN C C====== 4/mmm C IF (IK.NE.0 .AND. IH.NE.IK) THEN DO 50 I = 1 , NSYM NREF(1,I+NSYM) = NREF(2,I) NREF(2,I+NSYM) = NREF(1,I) NREF(3,I+NSYM) = NREF(3,I) 50 CONTINUE NSYM = NSYM*2 END IF END IF GO TO 160 C C====== 6/m C 60 IF (IH.NE.0 .OR. IK.NE.0) THEN NREF(1,2) = -(IH+IK) NREF(2,2) = IH NREF(3,2) = IL NREF(1,3) = IK NREF(2,3) = -(IH+IK) NREF(3,3) = IL NSYM = 3 IF (IL.NE.0) THEN DO 70 I = 1 , NSYM NREF(1,I+NSYM) = NREF(1,I) NREF(2,I+NSYM) = NREF(2,I) NREF(3,I+NSYM) = -NREF(3,I) 70 CONTINUE NSYM = NSYM*2 END IF END IF IF (ILGP.NE.5) THEN C C====== 6/mmm C IF (IK.NE.0 .AND. IH.NE.IK) THEN DO 80 I = 1 , NSYM NREF(1,I+NSYM) = NREF(2,I) NREF(2,I+NSYM) = NREF(1,I) NREF(3,I+NSYM) = NREF(3,I) 80 CONTINUE NSYM = NSYM*2 END IF END IF GO TO 160 C C====== 3bar (hexagonal axes) C 90 IF (MSYS.NE.5) THEN C C====== 3bar (rhombohedral axes) C IF (MSYS.EQ.6) THEN IF (IH.NE.IK .OR. IH.NE.IL) THEN NREF(1,2) = NREF(2,1) NREF(2,2) = NREF(3,1) NREF(3,2) = NREF(1,1) NREF(1,3) = NREF(3,1) NREF(2,3) = NREF(1,1) NREF(3,3) = NREF(2,1) NSYM = 3 IF (ILGP.NE.7) THEN C C====== 3m1 and 31m (rhombohedral axes) C IF (IH.NE.IK) THEN DO 100 I = 1 , NSYM NREF(1,I+NSYM) = NREF(2,I) NREF(2,I+NSYM) = NREF(1,I) NREF(3,I+NSYM) = NREF(3,I) 100 CONTINUE NSYM = NSYM*2 END IF END IF END IF END IF GO TO 160 ELSE IF (IH.NE.0 .OR. IK.NE.0) THEN NREF(1,2) = -(IH+IK) NREF(2,2) = IH NREF(3,2) = IL NREF(1,3) = IK NREF(2,3) = -(IH+IK) NREF(3,3) = IL NSYM = 3 END IF IF (ILGP.EQ.7) GO TO 160 C C====== 3m1 and 31m C IF (IL.EQ.0) THEN IF (IH.EQ.IK) GO TO 160 IF (IK.EQ.0) GO TO 160 ELSE IF (ILGP.EQ.8) THEN IF (IK.EQ.0) GO TO 160 ELSE IF (ILGP.EQ.9) THEN IF (IH.EQ.IK) GO TO 160 END IF END IF DO 110 I = 1 , NSYM NREF(1,I+NSYM) = NREF(2,I) NREF(2,I+NSYM) = NREF(1,I) IF (ILGP.EQ.8) THEN NREF(3,I+NSYM) = -NREF(3,I) ELSE NREF(3,I+NSYM) = NREF(3,I) END IF 110 CONTINUE NSYM = NSYM*2 GO TO 160 C C====== m3 C 120 IF (IK.NE.0 .OR. IL.NE.0) THEN NREF(1,2) = -IH NREF(2,2) = IK NREF(3,2) = IL NSYM = 2 IF (IK.NE.0) THEN DO 130 I = 1 , NSYM NREF(1,I+NSYM) = NREF(1,I) NREF(2,I+NSYM) = -NREF(2,I) NREF(3,I+NSYM) = NREF(3,I) 130 CONTINUE NSYM = NSYM*2 END IF END IF IF (IH.NE.IK .OR. IH.NE.IL) THEN DO 140 I = 1 , NSYM NREF(1,I+NSYM) = NREF(2,I) NREF(2,I+NSYM) = NREF(3,I) NREF(3,I+NSYM) = NREF(1,I) NREF(1,I+2*NSYM) = NREF(3,I) NREF(2,I+2*NSYM) = NREF(1,I) NREF(3,I+2*NSYM) = NREF(2,I) 140 CONTINUE NSYM = NSYM*3 END IF IF (ILGP.NE.10) THEN C C====== m3m C IF (IH.NE.IK) THEN DO 150 I = 1 , NSYM NREF(1,I+NSYM) = NREF(2,I) NREF(2,I+NSYM) = NREF(1,I) NREF(3,I+NSYM) = NREF(3,I) 150 CONTINUE NSYM = NSYM*2 END IF END IF C 160 DO 180 I = 1 , NSYM DO 170 J = 1 , 3 NREF(J,I+NSYM) = -NREF(J,I) 170 CONTINUE 180 CONTINUE NSYM = 2*NSYM C RETURN END **==hklgen.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** HKLGEN ** C*************************** C C SUBROUTINE HKLGEN(DSTMAX,IFLAG) C INCLUDE 'implicit_none.finc' C C Purpose: Find all reflections within dstar-max in the hemisphere +h, all k,l. C In essence the routine flow is to C (1) increase h from 0 to max C (2) increase k from 0 to max, then C decrease k from -1 to min C (3) increase l from 0 to max, then C decrease l from -1 to min C with h varying least, l most rapidly. The max and min values are C determined by a dstar which exceeds the limit. Complications are C caused by non-orthogonal cells: for given h,k l=0 may exceed dstar- C max, but some other l may not. The logic is that if at least one l C for given h,k is found within dstar-max, then the l loops end when C dstar exceeds dstar-max. If such an l is not found then the loops end C in both directions when the calculated dstar is increasing. Similarly C for the control of k loops, where the test is for at least one k,l C inside for given h. C COMMON /AMAT/A,DUMMY COMMON /SPOTS/NR COMMON /SYS/ISYS,KSYS COMMON /XYZ/IHKL C INCLUDE 'chanel.finc' C REAL DSTMAX INTEGER IFLAG INTEGER ISYS,NR REAL A(3,3),DUMMY(3,3) INTEGER IHKL(20000),KSYS(3) C C====== LOCALS: C REAL DST2,DST2KM,DST2LM,DSTEST,XG,YG,ZG INTEGER H,K,KSTART,L,LIMIT,LSTART LOGICAL ABSNT,KIN,KPOS,LIN,LPOS,UNIQ CHARACTER IYN*1 INTEGER IPACK C C------------------------------------------------------------------------------- C C====== Initialisations C LIMIT = 20000 DSTEST = DSTMAX*DSTMAX NR = 1 KSTART = 0 LSTART = 0 C C====== Start h loop C H = 0 C C====== Start k loop C 20 K = KSTART KPOS = .TRUE. KIN = .FALSE. DST2KM = 1000. C C====== Check for this h,k systematically absent, whatever the l value C 30 IF (ISYS.NE.0) THEN CALL SYSAHK(H,K,ABSNT) IF (ABSNT) THEN IF (KPOS) THEN K = K + 1 ELSE K = K - 1 END IF GO TO 30 END IF END IF C C====== Start l loop - skip 000 C L = LSTART IF (H.EQ.0 .AND. K.EQ.0) THEN L = 1 END IF LPOS = .TRUE. LIN = .FALSE. DST2LM = 1000. C C====== Check for this h,k,l systematically absent C 40 IF (ISYS.NE.0) THEN CALL SYSABS(H,K,L,ABSNT) IF (ABSNT) THEN IF (LPOS) THEN L = L + 1 ELSE L = L - 1 END IF GO TO 40 END IF END IF C C====== Calculate dstar squared. C XG = A(1,1)*H + A(1,2)*K + A(1,3)*L YG = A(2,1)*H + A(2,2)*K + A(2,3)*L ZG = A(3,1)*H + A(3,2)*K + A(3,3)*L DST2 = XG**2 + YG**2 + ZG**2 C C====== Test for within resolution C IF (DST2.LE.DSTEST) THEN LIN = .TRUE. KIN = .TRUE. C C====== Check for this h,k,l in the unique set according to Laue symmetry C CALL UNIQUE(H,K,L,UNIQ) IF (UNIQ) THEN IHKL(NR) = IPACK(H,K,L) NR = NR + 1 IF (NR.GT.LIMIT) THEN C C====== Ran out of array space C WRITE (IPRINT,6010) LIMIT , H , K , L WRITE (IUNIT6,6010) LIMIT , H , K , L NR = LIMIT RETURN END IF END IF C C====== Increment l and continue C IF (LPOS) THEN L = L + 1 ELSE L = L - 1 END IF GO TO 40 C C====== hkl outside resolution ... test if any hits for this l so far C ELSE IF (LIN) THEN C C====== Yes ... reverse l direction C IF (LPOS) THEN LPOS = .FALSE. L = LSTART - 1 GO TO 40 C C====== or increment k if already reversed l C ELSE IF (KPOS) THEN K = K + 1 ELSE K = K - 1 END IF GO TO 30 C C====== No ... check if getting closer ... save and continue if so C ELSE IF (DST2.LT.DST2LM) THEN DST2LM = DST2 IF (LPOS) THEN L = L + 1 ELSE L = L - 1 END IF GO TO 40 C C====== Getting further away ... reverse l and continue C ELSE IF (LPOS) THEN LPOS = .FALSE. L = LSTART - 1 GO TO 40 C C====== No hits for this l so far ... test if any for this k so far C ELSE IF (KIN) THEN C C====== Yes ... reverse k direction C IF (KPOS) THEN KPOS = .FALSE. K = KSTART - 1 GO TO 30 C C====== or increment h if already reversed k C ELSE H = H + 1 END IF GO TO 20 C C====== No ... check if getting closer ... save and continue if so C ELSE IF (DST2.LT.DST2KM) THEN DST2KM = DST2 IF (KPOS) THEN K = K + 1 ELSE K = K - 1 END IF GO TO 30 C C====== Getting further away ... reverse k and continue C ELSE IF (KPOS) THEN KPOS = .FALSE. K = KSTART - 1 GO TO 30 C C====== No hits for this h value so stop, returning iflag=0 to mark end C ELSE IF (H.EQ.0) THEN H = H + 1 GO TO 20 ELSE IFLAG = 0 NR = NR - 1 RETURN END IF C C====== Format statements C 6010 FORMAT (' limit of',I6,' reflections reached at h,k,l=', + 3I4/' will have to truncate calculation here') END **==inputv.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** INPUTV ** C*************************** C C SUBROUTINE INPUTV(PROMPT,VAR,VALUES,NVALS) C INCLUDE 'implicit_none.finc' C C Purpose : Extract a value from the VALUES array returned by MENU, or C if empty prompt user for the value. C CHARACTER*(*) PROMPT REAL VAR, VALUES(*) INTEGER NVALS C INCLUDE 'chanel.finc' C C====== FUNCTION TYPES: C REAL FLPVAL C C====== LOCALS: C LOGICAL QUIT CHARACTER*30 TMPSTR LOGICAL ERR C C------------------------------------------------------------------------------- C IF (NVALS.GT.0) THEN VAR = VALUES(1) ELSE CALL TPROMP(PROMPT) CALL TREAD(TMPSTR,QUIT) VAR = FLPVAL(TMPSTR,ERR) END IF END **==inputw.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** INPUTW ** C*************************** C C SUBROUTINE INPUTW(PROMPT,WORD,WORDS,NWORDS) C INCLUDE 'implicit_none.finc' C C Purpose : Extract a word from the WORDS array returned by MENU, or C if empty prompt user for the word. C CHARACTER*(*) PROMPT CHARACTER*4 WORD CHARACTER*4 WORDS(*) INTEGER NWORDS C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL QUIT C C------------------------------------------------------------------------------- C IF (NWORDS.GT.0) THEN WORD = WORDS(1) ELSE CALL TPROMP(PROMPT) CALL TREAD(WORD,QUIT) END IF CALL CCPUPC(WORD) RETURN END **==ipack.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** IPACK ** C*************************** C C INTEGER FUNCTION IPACK(H,K,L) C INTEGER H,K,L C C------------------------------------------------------------------------------- C IPACK = 40000*(H+100) + 200*(K+100) + (L+100) RETURN END **==item.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** ITEM ** C*************************** C C SUBROUTINE ITEM(IOPT,ARRAY,LEN,VIDEO) INCLUDE 'implicit_none.finc' C C Purpose : Displays character string on VDU at given line and column C CHARACTER * (*) ARRAY INTEGER IOPT,LEN LOGICAL VIDEO C C====== LOCALS: C INTEGER LINE,ICOL C C------------------------------------------------------------------------------- C ICOL = (IOPT+16)/17 LINE = IOPT - 17*(ICOL-1) CALL QSSETC(1,24) ICOL = ICOL*41 - 40 LINE = LINE + 1 CALL QSSETC(ICOL,LINE) IF (VIDEO) CALL QSREVV CALL TSTRNG(ARRAY(1:LEN)) CALL QSNORV RETURN END **==laue.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** LAUE ** C*************************** C C SUBROUTINE LAUE(ION) C INCLUDE 'implicit_none.finc' C C Purpose: Calculate but do not display the Laue photograph. C The array ion returned by this routine shows which of the nr C reflections in ihkl would appear on the film plot (ion(i)=n>0) C or not (ion(i)=0). C (n) is the total number of symmetry equivalents for this reflection. C COMMON /AMAT/A,DUMMY COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /RADC/RADI,RADI2 COMMON /SPOTS/NR COMMON /XYZ/IHKL INTEGER ION(*) REAL CF,DMIN,RADI,RADI2,SPNDLE,WVM,WVS,WVX INTEGER NR,NUMA,NUMC REAL A(3,3),CCELL(6),CELL(6),DUMMY(3,3),PHIS(3) INTEGER IHKL(20000) C C====== LOCALS: C REAL C1,DD1,DD2,RF,RMAX,RMAX2,RMIN,RMIN2,S1,X,XF,XG,Y,YF,YG,Z,ZG INTEGER I,IH,IK,IL,N,NSYM INTEGER NREF(3,48) C C------------------------------------------------------------------------------- C C1 = SPNDLE*3.1415926536/180.0 S1 = SIN(C1) C1 = COS(C1) RMIN = 1./WVX RMAX = 1./WVM RMIN2 = RMIN*RMIN RMAX2 = RMAX*RMAX C C====== Loop through all the reflections C DO 30 N = 1 , NR ION(N) = 0 C C====== Unpack reflection indices C CALL UNPACK(IHKL(N),IH,IK,IL) C C====== Generate symmetry equivalent reflections C CALL GENSYM(NSYM,NREF,IH,IK,IL) C C====== Loop through these symmetry equivalents C DO 20 I = 1 , NSYM IH = NREF(1,I) IK = NREF(2,I) IL = NREF(3,I) XG = A(1,1)*IH + A(1,2)*IK + A(1,3)*IL YG = A(2,1)*IH + A(2,2)*IK + A(2,3)*IL ZG = A(3,1)*IH + A(3,2)*IK + A(3,3)*IL X = C1*XG - S1*YG IF (X.LT.0.00001) THEN Y = S1*XG + C1*YG Z = ZG C C====== Dists squared of rl point from centre min and max spheres C DD1 = ((RMIN+X)**2+Y*Y+Z*Z) IF (DD1.GE.RMIN2) THEN DD2 = ((RMAX+X)**2+Y*Y+Z*Z) IF (DD2.LE.RMAX2) THEN RF = (X*X+Y*Y+Z*Z)/(2*X) C C====== rl point can reflect, centre refl sphere at RF,0,0 C IF (CF*(-RF+X).GT.0.) THEN XF = CF*Y/(-RF+X) YF = CF*Z/(-RF+X) C C====== Check that point is inside the circle requested by the user C IF (ABS(XF).LE.RADI) THEN IF (ABS(YF).LE.RADI) THEN IF ((XF*XF+YF*YF).LE.RADI2) THEN ION(N) = ION(N) + 1 END IF END IF END IF END IF END IF END IF END IF 20 CONTINUE 30 CONTINUE END **==lauer.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** LAUER ** C*************************** C C SUBROUTINE LAUER(ION,WRES) C INCLUDE 'implicit_none.finc' C C Purpose: Calculate but do not display the Laue photograph. C The array ion returned by this routine shows which of the nr C reflections in ihkl would appear on the film plot (n) or not (0). C (n) is the total number of symmetry equivalents for this reflection. C The array wres records the resolution for each. C COMMON /AMAT/A,DUMMY COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /RADC/RADI,RADI2 COMMON /SPOTS/NR COMMON /XYZ/IHKL REAL WRES(*) INTEGER ION(*) REAL CF,DMIN,RADI,RADI2,SPNDLE,WVM,WVS,WVX INTEGER NR,NUMA,NUMC REAL A(3,3),CCELL(6),CELL(6),DUMMY(3,3),PHIS(3) INTEGER IHKL(20000) C C====== LOCALS: C REAL C1,DD1,DD2,RF,RMAX,RMAX2,RMIN,RMIN2,S1,X,XF,XG,Y,YF,YG,Z,ZG INTEGER I,IH,IK,IL,N,NSYM INTEGER NREF(3,48) C C------------------------------------------------------------------------------- C C1 = SPNDLE*3.1415926536/180.0 S1 = SIN(C1) C1 = COS(C1) RMIN = 1./WVX RMAX = 1./WVM RMIN2 = RMIN*RMIN RMAX2 = RMAX*RMAX C C====== Loop through all the reflections C DO 30 N = 1 , NR ION(N) = 0 C C====== Unpack reflection indices C CALL UNPACK(IHKL(N),IH,IK,IL) C C====== Generate symmetry equivalent reflections C CALL GENSYM(NSYM,NREF,IH,IK,IL) C C====== Loop through these symmetry equivalents C DO 20 I = 1 , NSYM IH = NREF(1,I) IK = NREF(2,I) IL = NREF(3,I) XG = A(1,1)*IH + A(1,2)*IK + A(1,3)*IL YG = A(2,1)*IH + A(2,2)*IK + A(2,3)*IL ZG = A(3,1)*IH + A(3,2)*IK + A(3,3)*IL X = C1*XG - S1*YG IF (X.LT.0.00001) THEN Y = S1*XG + C1*YG Z = ZG C C====== Dists squared of rl point from centre min and max spheres C DD1 = ((RMIN+X)**2+Y*Y+Z*Z) IF (DD1.GE.RMIN2) THEN DD2 = ((RMAX+X)**2+Y*Y+Z*Z) IF (DD2.LE.RMAX2) THEN RF = (X*X+Y*Y+Z*Z)/(2*X) C C====== rl point can reflect, centre refl sphere at RF,0,0 C IF (CF*(-RF+X).GT.0.) THEN XF = CF*Y/(-RF+X) YF = CF*Z/(-RF+X) C C====== Check that point is inside the circle requested by the user C IF (ABS(XF).LE.RADI) THEN IF (ABS(YF).LE.RADI) THEN IF ((XF*XF+YF*YF).LE.RADI2) THEN ION(N) = ION(N) + 1 WRES(N) = 1./SQRT(XG**2+YG**2+ZG**2) GO TO 30 END IF END IF END IF END IF END IF END IF END IF 20 CONTINUE 30 CONTINUE END **==lauew.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** LAUEW ** C*************************** C C SUBROUTINE LAUEW(ION,WION) C INCLUDE 'implicit_none.finc' C C Purpose: Calculate but do not display the Laue photograph. C The array ion returned by this routine shows which of the nr C reflections in ihkl would appear on the film plot (n) or not (0). C (n) is the total number of symmetry equivalents for this reflection. C The array wion records the variation in wavelength for each. C COMMON /AMAT/A,DUMMY COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /RADC/RADI,RADI2 COMMON /SPOTS/NR COMMON /XYZ/IHKL REAL WION(*) INTEGER ION(*) REAL CF,DMIN,RADI,RADI2,SPNDLE,WVM,WVS,WVX INTEGER NR,NUMA,NUMC REAL A(3,3),CCELL(6),CELL(6),DUMMY(3,3),PHIS(3) INTEGER IHKL(20000) C C====== LOCALS: C REAL C1,DD1,DD2,DIFF,RF,RMAX,RMAX2,RMIN,RMIN2,S1,SUM,X,XF,XG,Y,YF, + YG,Z,ZG INTEGER I,IH,IK,IL,J,M,N,NSYM REAL WAVE(48) INTEGER NREF(3,48) C C------------------------------------------------------------------------------- C C1 = SPNDLE*3.1415926536/180.0 S1 = SIN(C1) C1 = COS(C1) RMIN = 1./WVX RMAX = 1./WVM RMIN2 = RMIN*RMIN RMAX2 = RMAX*RMAX C C====== Loop through all the reflections C DO 50 N = 1 , NR ION(N) = 0 C C====== Unpack reflection indices C CALL UNPACK(IHKL(N),IH,IK,IL) C C====== Generate symmetry equivalent reflections C CALL GENSYM(NSYM,NREF,IH,IK,IL) C C====== Loop through these symmetry equivalents C M = 0 DO 20 I = 1 , NSYM IH = NREF(1,I) IK = NREF(2,I) IL = NREF(3,I) XG = A(1,1)*IH + A(1,2)*IK + A(1,3)*IL YG = A(2,1)*IH + A(2,2)*IK + A(2,3)*IL ZG = A(3,1)*IH + A(3,2)*IK + A(3,3)*IL X = C1*XG - S1*YG IF (X.LT.0.00001) THEN Y = S1*XG + C1*YG Z = ZG C C====== Dists squared of rl point from centre min and max spheres C DD1 = ((RMIN+X)**2+Y*Y+Z*Z) IF (DD1.GE.RMIN2) THEN DD2 = ((RMAX+X)**2+Y*Y+Z*Z) IF (DD2.LE.RMAX2) THEN RF = (X*X+Y*Y+Z*Z)/(2*X) C C====== rl point can reflect, centre refl sphere at rf,0,0 C IF (CF*(-RF+X).GT.0.) THEN XF = CF*Y/(-RF+X) YF = CF*Z/(-RF+X) C C====== Check that point is inside the circle requested by the user C IF (ABS(XF).LE.RADI) THEN IF (ABS(YF).LE.RADI) THEN IF ((XF*XF+YF*YF).LE.RADI2) THEN ION(N) = ION(N) + 1 M = M + 1 WAVE(M) = -1./RF END IF END IF END IF END IF END IF END IF END IF 20 CONTINUE C C====== Calculate wave length variation C IF (M.GT.1) THEN SUM = 0. DO 30 J = 1 , M SUM = SUM + WAVE(J) 30 CONTINUE SUM = SUM/M DIFF = 0. DO 40 J = 1 , M DIFF = ABS(WAVE(J)-SUM) + DIFF 40 CONTINUE WION(N) = DIFF/M END IF 50 CONTINUE END **==lgrp.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** LGRP ** C*************************** C C SUBROUTINE LGRP(MSYS,ILGP) C INCLUDE 'implicit_none.finc' C C Purpose: Check the crystal system in msys and ask the user for the Laue C group (returned in ilgp) if necessary. C INTEGER ILGP,MSYS C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL QUIT INTEGER ITEMP C C------------------------------------------------------------------------------- C IF (MSYS.EQ.2) THEN C C====== Monoclinic C ILGP = 1 RETURN ELSE IF (MSYS.EQ.3) THEN C C====== Orthorhombic C ILGP = 2 RETURN ELSE IF (MSYS.EQ.4) THEN C C====== Tetragonal - may be 4/m or 4/mmm C 20 CALL TWRITE( + 'Select Laue group for symmetry equivalent analysis' + ) ITEMP = 1 CALL TASKI('1 = 4/m, 2 = 4/mmm',ITEMP,QUIT) IF (ITEMP.NE.1 .AND. ITEMP.NE.2) GO TO 20 ILGP = 2 + ITEMP RETURN ELSE IF (MSYS.EQ.5) THEN C C====== Hexagonal - may be 6/m or 6/mmm C====== or 3bar, 3m1 or 31m for trigonal indexed on hexagonal lattice C 30 CALL TWRITE( + 'Select Laue group for symmetry equivalent analysis' + ) ITEMP = 1 CALL TASKI('1 = 6/m, 2 = 6/mmm, 3 = 3bar, 4 = 3m1, 5 = 31m', + ITEMP,QUIT) IF (ITEMP.LT.1 .AND. ITEMP.GT.5) GO TO 30 ILGP = 4 + ITEMP RETURN ELSE IF (MSYS.EQ.6) THEN C C====== Rhombohedral lattice is trigonal - may be 3bar, 3m1 or 31m C 40 CALL TWRITE( + 'Select Laue group for symmetry equivalent analysis' + ) ITEMP = 1 CALL TASKI('1 = 3bar, 2 = 3m1, 3 = 31m',ITEMP,QUIT) IF (ITEMP.LT.1 .OR. ITEMP.GT.3) GO TO 40 ILGP = 6 + ITEMP RETURN ELSE IF (MSYS.EQ.7) THEN C C====== Cubic C 50 CALL TWRITE( + 'Select Laue group for symmetry equivalent analysis' + ) ITEMP = 1 CALL TASKI('1 = m3, 2 = m3m',ITEMP,QUIT) IF (ITEMP.NE.1 .AND. ITEMP.NE.2) GO TO 50 ILGP = 9 + ITEMP RETURN ELSE C C====== Triclinic C ILGP = 0 RETURN END IF END **==matm33.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** MATM33 ** C*************************** C C SUBROUTINE MATM33(A,B,C) C INCLUDE 'implicit_none.finc' C REAL A(3,3),B(3,3),C(3,3) C C====== LOCALS: C INTEGER I,J,K C C------------------------------------------------------------------------------- C DO 40 I = 1 , 3 DO 30 J = 1 , 3 C(I,J) = 0.0 DO 20 K = 1 , 3 C(I,J) = C(I,J) + A(I,K)*B(K,J) 20 CONTINUE 30 CONTINUE 40 CONTINUE RETURN END **==mats.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** MATS ** C*************************** C C SUBROUTINE MATS(AM,AA,AB,AC,AD,AE,AF,CELL) C INCLUDE 'implicit_none.finc' C C Purpose: Generalised matrix and derivative builder C the reciprocal cell parameters come in cell(1-6) but only the C unique ones ... cell(1) is astar, cell(2) alphastar for rhombohedral C for example. C am returns the A matrix ... row 1 is components of reciprocal cell C along X, row 2 along Y, row 3 along Z. C aa, ab etc are the derivative matrices for A wrt the unique cell C parameters in order ... thus aa is dA/dastar, ab dA/alphastar for C rhombohedral. C COMMON /MATRIX/MSYS,MROT,MBEAM,IJUNK1,IJUNK2 REAL AA(3,3),AB(3,3),AC(3,3),AD(3,3),AE(3,3),AF(3,3),AM(3,3), + CELL(6) INTEGER IJUNK1,IJUNK2,MBEAM,MROT,MSYS C C====== LOCALS: C REAL A,C,COSPSI,COSRHO,FUN,S,SI,SJ,SK INTEGER I,J,K,N1,N2,N3 REAL CC(3),DA(3,3,3),DB(3,3,3),RL(3),SS(3) C C------------------------------------------------------------------------------- C C C====== I refers to mrot, the crystal rotation axis, Z C I = IABS(MROT) SI = 1.0 IF (MROT.LT.0) SI = -1.0 C C====== J refers to mbeam, the xray beam axis, X C J = IABS(MBEAM) SJ = 1.0 IF (MBEAM.LT.0) SJ = -1.0 C C====== K is therefore the other axis, Y C K = 6 - I - J C C====== The sign of k relies on i,j,k and X,Y,Z being right-handed C IF (K.EQ.J+1 .OR. K.EQ.J-2) THEN SK = SI*SJ ELSE SK = -SI*SJ END IF C IF (MSYS.EQ.2) THEN C C====== Monoclinic C RL(1) = CELL(1) RL(2) = CELL(2) RL(3) = CELL(3) C CC(1) = 0.0 CC(2) = COS(CELL(4)) CC(3) = 0.0 C SS(1) = 1.0 SS(2) = SIN(CELL(4)) SS(3) = 1.0 C C====== A matrix C AM(1,J) = SJ*RL(J)*SS(K) AM(1,K) = SJ*RL(K)*CC(I) AM(2,K) = SK*RL(K)*SS(I)*SS(J) AM(3,I) = SI*RL(I) AM(3,J) = SI*RL(J)*CC(K) AM(3,K) = SI*RL(K)*CC(J) C C====== Derivatives wrt reciprocal cell lengths C DO 40 N1 = 1 , 3 DO 30 N2 = 1 , 3 DO 20 N3 = 1 , 3 DA(N1,N2,N3) = 0. 20 CONTINUE 30 CONTINUE 40 CONTINUE C DA(3,I,I) = SI DA(1,J,J) = SJ*SS(K) DA(3,J,J) = SI*CC(K) DA(1,K,K) = SJ*CC(I) DA(2,K,K) = SK*SS(I)*SS(J) DA(3,K,K) = SI*CC(J) C DO 70 N1 = 1 , 3 DO 60 N2 = 1 , 3 DO 50 N3 = 1 , 3 IF (DA(N1,N2,N3).NE.0.0) THEN IF (N3.EQ.1) AA(N1,N2) = DA(N1,N2,N3) IF (N3.EQ.2) AB(N1,N2) = DA(N1,N2,N3) IF (N3.EQ.3) AC(N1,N2) = DA(N1,N2,N3) END IF 50 CONTINUE 60 CONTINUE 70 CONTINUE C C====== Derivatives wrt reciprocal cell angles (only ss(2) and cc(2)) C IF (I.EQ.2) THEN AD(1,K) = -SJ*RL(K)*SS(I) AD(2,K) = SK*RL(K)*CC(I) END IF IF (J.EQ.2) THEN AD(2,K) = SK*RL(K)*CC(J) AD(3,K) = -SI*RL(K)*SS(J) END IF IF (K.EQ.2) THEN AD(1,J) = SJ*RL(J)*CC(K) AD(3,J) = -SI*RL(J)*SS(K) END IF RETURN ELSE IF (MSYS.EQ.3) THEN C C====== Orthorhombic C RL(1) = CELL(1) RL(2) = CELL(2) RL(3) = CELL(3) AM(1,J) = SJ*RL(J) AM(2,K) = SK*RL(K) AM(3,I) = SI*RL(I) C C====== Derivatives C IF (J.EQ.1) AA(1,1) = SJ IF (K.EQ.1) AA(2,1) = SK IF (I.EQ.1) AA(3,1) = SI C IF (J.EQ.2) AB(1,2) = SJ IF (K.EQ.2) AB(2,2) = SK IF (I.EQ.2) AB(3,2) = SI C IF (J.EQ.3) AC(1,3) = SJ IF (K.EQ.3) AC(2,3) = SK IF (I.EQ.3) AC(3,3) = SI RETURN ELSE IF (MSYS.EQ.4) THEN C C====== Tetragonal C RL(1) = CELL(1) RL(2) = CELL(1) RL(3) = CELL(2) AM(1,J) = SJ*RL(J) AM(2,K) = SK*RL(K) AM(3,I) = SI*RL(I) C C====== Derivatives C IF (J.EQ.1) AA(1,1) = SJ IF (K.EQ.1) AA(2,1) = SK IF (I.EQ.1) AA(3,1) = SI C IF (J.EQ.2) AA(1,2) = SJ IF (K.EQ.2) AA(2,2) = SK IF (I.EQ.2) AA(3,2) = SI C IF (J.EQ.3) AB(1,3) = SJ IF (K.EQ.3) AB(2,3) = SK IF (I.EQ.3) AB(3,3) = SI RETURN ELSE IF (MSYS.EQ.5) THEN C C====== Hexagonal C RL(1) = CELL(1) RL(2) = CELL(1) RL(3) = CELL(2) C CC(1) = 0.0 CC(2) = 0.0 CC(3) = 0.5 C SS(1) = 1.0 SS(2) = 1.0 SS(3) = SQRT(3.)/2. C C====== A matrix C AM(1,J) = SJ*RL(J)*SS(K) AM(1,K) = SJ*RL(K)*CC(I) AM(2,K) = SK*RL(K)*SS(I)*SS(J) AM(3,I) = SI*RL(I) AM(3,J) = SI*RL(J)*CC(K) AM(3,K) = SI*RL(K)*CC(J) C C====== Derivatives wrt reciprocal cell lengths C DO 100 N1 = 1 , 3 DO 90 N2 = 1 , 3 DO 80 N3 = 1 , 3 DA(N1,N2,N3) = 0. 80 CONTINUE 90 CONTINUE 100 CONTINUE C DA(3,I,I) = SI DA(1,J,J) = SJ*SS(K) DA(3,J,J) = SI*CC(K) DA(1,K,K) = SJ*CC(I) DA(2,K,K) = SK*SS(I)*SS(J) DA(3,K,K) = SI*CC(J) C DO 130 N1 = 1 , 3 DO 120 N2 = 1 , 3 DO 110 N3 = 1 , 3 IF (DA(N1,N2,N3).NE.0.0) THEN IF (N3.EQ.1 .OR. N3.EQ.2) AA(N1,N2) = DA(N1,N2,N3) IF (N3.EQ.3) AB(N1,N2) = DA(N1,N2,N3) END IF 110 CONTINUE 120 CONTINUE 130 CONTINUE RETURN ELSE IF (MSYS.EQ.6) THEN C C====== Rhombohedral C A = CELL(1) C = COS(CELL(2)) S = SIN(CELL(2)) FUN = SQRT(1.-3.*C**2+2.*C**3) AM(1,J) = SJ*A*S AM(1,K) = SJ*A*C*(1.-C)/S AM(2,K) = SK*A*FUN/S AM(3,I) = SI*A AM(3,J) = SI*A*C AM(3,K) = SI*A*C C AA(1,J) = SJ*S AA(1,K) = AM(1,K)/A AA(2,K) = AM(2,K)/A AA(3,I) = SI AA(3,J) = SI*C AA(3,K) = SI*C C AB(1,J) = SJ*A*C AB(1,K) = SJ*A*(C-(1.-C)/S**2) AB(2,K) = SK*A*(3.*C*(1.-C)/FUN-FUN*C/S**2) AB(3,J) = -SI*A*S AB(3,K) = -SI*A*S RETURN ELSE IF (MSYS.EQ.7) THEN C C====== Cubic C A = CELL(1) AM(1,J) = SJ*A AM(2,K) = SK*A AM(3,I) = SI*A C AA(1,J) = SJ AA(2,K) = SK AA(3,I) = SI RETURN ELSE C C====== Triclinic C RL(1) = CELL(1) RL(2) = CELL(2) RL(3) = CELL(3) C CC(1) = COS(CELL(4)) CC(2) = COS(CELL(5)) CC(3) = COS(CELL(6)) C SS(1) = SIN(CELL(4)) SS(2) = SIN(CELL(5)) SS(3) = SIN(CELL(6)) C C====== See Buerger's Xray Crystallography C====== for the next formulae, p348 C COSPSI = (CC(I)-CC(K)*CC(J))/SS(K) COSRHO = SQRT(1.-CC(I)**2-CC(J)**2-CC(K)**2+2.*CC(I)*CC(J) + *CC(K))/SS(K) C C====== A matrix C AM(1,J) = SJ*RL(J)*SS(K) AM(1,K) = SJ*RL(K)*COSPSI AM(2,K) = SK*RL(K)*COSRHO AM(3,I) = SI*RL(I) AM(3,J) = SI*RL(J)*CC(K) AM(3,K) = SI*RL(K)*CC(J) C C====== Derivatives wrt reciprocal cell lengths C DO 160 N1 = 1 , 3 DO 150 N2 = 1 , 3 DO 140 N3 = 1 , 3 DA(N1,N2,N3) = 0. DB(N1,N2,N3) = 0. 140 CONTINUE 150 CONTINUE 160 CONTINUE C DA(3,I,I) = SI DA(1,J,J) = SJ*SS(K) DA(3,J,J) = SI*CC(K) DA(1,K,K) = SJ*COSPSI DA(2,K,K) = SK*COSRHO DA(3,K,K) = SI*CC(J) C DO 190 N1 = 1 , 3 DO 180 N2 = 1 , 3 DO 170 N3 = 1 , 3 IF (DA(N1,N2,N3).NE.0.0) THEN IF (N3.EQ.1) AA(N1,N2) = DA(N1,N2,N3) IF (N3.EQ.2) AB(N1,N2) = DA(N1,N2,N3) IF (N3.EQ.3) AC(N1,N2) = DA(N1,N2,N3) END IF 170 CONTINUE 180 CONTINUE 190 CONTINUE C C====== Derivatives wrt reciprocal cell angles C DB(1,K,I) = -SJ*RL(K)*SS(I)/SS(K) DB(2,K,I) = SK*RL(K)*SS(I)*(CC(I)-CC(J)*CC(K)) + /(SS(K)**2*COSRHO) DB(1,K,J) = SJ*RL(K)*CC(K)*SS(J)/SS(K) DB(2,K,J) = SK*RL(K)*SS(J)*(CC(J)-CC(I)*CC(K)) + /(SS(K)**2*COSRHO) DB(3,K,J) = -SI*RL(K)*SS(J) DB(1,J,K) = SJ*RL(J)*CC(K) DB(1,K,K) = SJ*RL(K)*(CC(J)-CC(K)*CC(I))/SS(K)**2 DB(2,K,K) = SK*RL(K)*((CC(K)-CC(I)*CC(J))/(COSRHO*SS(K))-CC(K) + *COSRHO/SS(K)) DB(3,J,K) = -SI*RL(J)*SS(K) C DO 220 N1 = 1 , 3 DO 210 N2 = 1 , 3 DO 200 N3 = 1 , 3 IF (DB(N1,N2,N3).NE.0.0) THEN IF (N3.EQ.1) AD(N1,N2) = DB(N1,N2,N3) IF (N3.EQ.2) AE(N1,N2) = DB(N1,N2,N3) IF (N3.EQ.3) AF(N1,N2) = DB(N1,N2,N3) END IF 200 CONTINUE 210 CONTINUE 220 CONTINUE RETURN END IF END **==menu.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** MENU ** C*************************** C C SUBROUTINE MENU(WORDS,NWORDS,VALUES,NVALS,IOPT) C INCLUDE 'implicit_none.finc' C REAL VALUES(*) CHARACTER*4 WORDS(*) INTEGER NWORDS,NVALS,IOPT C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL ABORT CHARACTER*80 STR CHARACTER*4 OPTION(34) INTEGER KOPT,INOPT,NOPT,I C DATA OPTION/'A ','B ','C ','ALPH','BETA','GAMM','WMIN', + 'WMAX','DIST','DMIN','RADI','SPDR','BMDR','PX ','PY ', + 'PZ ','SPIN','SYST','SCAL','DATS','DATG','LAUE','LATT', + 'ROTA','SEAR','WAVV','VARC','RESO','OFFR','TWOR','STOP', + 'HELP','LOG ',' '/ C DATA NOPT/34/ C C------------------------------------------------------------------------------- C INOPT = IOPT 20 CALL QSSETC(1,24) NWORDS = 0 NVALS = 0 CALL QSBELL WRITE (STR,6010) OPTION(INOPT) CALL TPROMP(STR) CALL TREAD(STR,ABORT) CALL CCPUPC(STR) CALL STRTKN(STR,' ,',2,WORDS,NWORDS,VALUES,NVALS) IF (NWORDS.GT.0) THEN KOPT = 0 DO 30 I = 1 , NOPT IF (WORDS(1).EQ.OPTION(I)) THEN KOPT = I GO TO 40 END IF 30 CONTINUE END IF 40 IF (.NOT.(ABORT)) THEN IF (NWORDS.EQ.0) KOPT = INOPT IOPT = KOPT IF (IOPT.LE.0) THEN C C====== Display error message C CALL QSREVV WRITE (STR,6020) WORDS(1) CALL TWRITE(STR) CALL QSNORV GO TO 20 END IF END IF C NWORDS = NWORDS - 1 IF (NWORDS.GE.1) WORDS(1) = WORDS(2) RETURN C C====== Format statements C 6010 FORMAT ('>>> Select option [',A4,']: ') 6020 FORMAT ('!!! Invalid option: ',A4) END **==murot.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** MUROT ** C*************************** C C SUBROUTINE MUROT C INCLUDE 'implicit_none.finc' C REAL DEGRAD PARAMETER (DEGRAD=3.14159265358979/180.) INTEGER MAXONE, MAXALL PARAMETER (MAXONE=20000,MAXALL=100000) C COMMON /DATA/CELL, WVM, WVX, WVS, CF, PHIS, SPNDLE, CCELL, NUMA, + NUMC, DMIN COMMON /SPOTS/NR REAL CF, DMIN, SPNDLE, WVM, WVS, WVX INTEGER NR, NUMA, NUMC REAL CCELL(6), CELL(6), PHIS(3) C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL QUIT INTEGER ION(MAXONE) INTEGER KSAVE(MAXALL) REAL PHISAV(3) REAL FRACT, STEP, PHIU, MU, SINPHIU, SINPHIY, COSPHIY INTEGER I, ICOUNT, IST, ITOT, J, K, NSTEP LOGICAL ISAV REAL C(4), SP(20) C C------------------------------------------------------------------------------- C WRITE (IPRINT,6010) MU = 0.0 CALL TASKF('Enter spindle offset angle µ (degrees)',MU,1,QUIT) MU = MU*DEGRAD NSTEP = 1 CALL TASKI('Enter number of steps (counting a spindle)'// + ' setting of 0 as the first)',NSTEP,QUIT) STEP = 10.0 CALL TASKF('Enter step size in degrees',STEP,1,QUIT) WRITE (IPRINT,6020) MU/DEGRAD WRITE (IUNIT6,6020) MU/DEGRAD WRITE (IPRINT,6030) WRITE (IUNIT6,6030) C DO 20 I = 1 , 3 PHISAV(I) = PHIS(I) 20 CONTINUE ISAV = .TRUE. IST = 0 C C====== Loop through the spindle steps C DO 80 I = 1 , NSTEP PHIU = (I-1)*STEP SP(I) = PHIU PHIU = PHIU*DEGRAD SINPHIU = SIN(PHIU) SINPHIY = -SIN(MU)*COS(MU)*(1-COS(PHIU)) PHIS(2) = ASIN(SINPHIY) COSPHIY = COS(PHIS(2)) PHIS(3) = ASIN(COS(MU)*SINPHIU/COSPHIY) PHIS(1) = ASIN(-SIN(MU)*SINPHIU/COSPHIY) DO 30 J = 1 , 3 PHIS(J) = PHIS(J)/DEGRAD + PHISAV(J) 30 CONTINUE C C====== Recalculate the A matrix since phis have changed C CELL(1) = CELL(1)*WVS CELL(2) = CELL(2)*WVS CELL(3) = CELL(3)*WVS J = 0 DO 40 K = 1 , 3 IF (CELL(K).NE.0.) THEN J = J + 1 CCELL(J) = CELL(K) END IF 40 CONTINUE CALL ACALC CELL(1) = CELL(1)/WVS CELL(2) = CELL(2)/WVS CELL(3) = CELL(3)/WVS CALL LAUE(ION) C C====== Find out how many reflections recorded C ICOUNT = 0 ITOT = 0 DO 50 J = 1 , 4 C(J) = 0. 50 CONTINUE DO 60 J = 1 , NR K = ION(J) C C======= If there is room save the ion array for overlap purposes C IF (ISAV) KSAVE(IST+J) = K IF (K.NE.0) THEN ITOT = ITOT + K IF (K.GT.4) K = 4 C(K) = C(K) + 1 ICOUNT = ICOUNT + 1 END IF 60 CONTINUE FRACT = FLOAT(ICOUNT)/FLOAT(NR) DO 70 J = 1 , 4 C(J) = C(J)/NR 70 CONTINUE WRITE (IPRINT,6040) PHIU/DEGRAD , FRACT , ICOUNT , (C(J),J=1,4) + , ITOT WRITE (IUNIT6,6040) PHIU/DEGRAD , FRACT , ICOUNT , (C(J),J=1,4) + , ITOT IST = IST + NR IF ((IST+NR).GT.MAXALL) ISAV = .FALSE. 80 CONTINUE DO 90 I = 1 , 3 PHIS(I) = PHISAV(I) 90 CONTINUE C C====== Investigate overlaps? C IF (NSTEP.GT.1) CALL OVERLP(KSAVE,NSTEP,SP) C C====== Format statements C 6010 FORMAT (' This command analyses the effect of rotating about '/ + 'the spindle axis at a fixed angle µ.') 6020 FORMAT (///' Effect of rotating about spindle at µ = ',F6.1, + ' on the fraction of symmetry'/ + ' unique reflections recorded') 6030 FORMAT (/// +' Spindle Fraction recorded Distribution of symmetry equivalent +s'/32X,' 1 2 3 4+'/) 6040 FORMAT (F6.1,F13.3,' (',I5,')',4X,4F6.3,' (',I5,')') END **==overcn.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** OVERCN ** C*************************** C C SUBROUTINE OVERCN(ICNT,IREC,KSAVE,NR) C INCLUDE 'implicit_none.finc' C C Purpose: Count in icnt the reflections in ksave not recorded in irec C irec=1 denotes not yet recorded, ksave.ne.0 denotes on film C INTEGER ICNT,NR INTEGER KSAVE(NR),IREC(NR) C C====== LOCALS: C INTEGER I C C------------------------------------------------------------------------------- C ICNT = 0 DO 20 I = 1 , NR IF (IREC(I).EQ.1 .AND. KSAVE(I).NE.0) ICNT = ICNT + 1 20 CONTINUE RETURN END **==overlp.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** OVERLP ** C*************************** C C SUBROUTINE OVERLP(KSAVE,NSTEP,SP) C INCLUDE 'implicit_none.finc' C INTEGER MAXONE,MAXALL PARAMETER (MAXONE=20000,MAXALL=100000) C COMMON /SPOTS/NR INTEGER NSTEP REAL SP(20) INTEGER KSAVE(MAXALL) INTEGER NR C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL QUIT REAL FR,FRLEFT INTEGER I,ICNT,IFILM,IMAX,IMAXNO,IOFF,IOFF2,IST,J,K REAL OV(20) INTEGER KSET(20) INTEGER IREC(MAXALL) C C------------------------------------------------------------------------------- C I = MAXALL/NR IF (I.LT.NSTEP) THEN NSTEP = I WRITE (IPRINT,6010) NSTEP WRITE (IUNIT6,6010) NSTEP END IF 20 CALL TWRITE('Analyse the overlap of the symmetry uniques') IST = 1 CALL TASKI('Enter row in above table to use as base (0=end)',IST, + QUIT) IF (IST.EQ.0) RETURN IF (IST.GT.NSTEP) GO TO 20 WRITE (IPRINT,6020) WRITE (IUNIT6,6020) C C====== Zero the array which says which spindle settings used C====== the value in the array gives the number of the film in the series C====== ifilm counts through the series C DO 30 I = 1 , NSTEP KSET(I) = 0 30 CONTINUE KSET(IST) = 1 IFILM = 1 C C====== Set irec array to say all uniques not recorded at start C DO 40 I = 1 , NR IREC(I) = 1 40 CONTINUE C C====== Calculate the ones recorded at the first spindle setting C IOFF = (IST-1)*NR + 1 CALL OVERCN(ICNT,IREC,KSAVE(IOFF),NR) CALL OVERST(IREC,KSAVE(IOFF),NR) FRLEFT = 1.0 FR = FLOAT(ICNT)/FLOAT(NR) FRLEFT = FRLEFT - FR WRITE (IPRINT,6030) SP(IST) , FR , FRLEFT WRITE (IUNIT6,6030) SP(IST) , FR , FRLEFT C C====== Loop through the remaining spindle settings C====== imax records max overlap with unrecordeds, imaxno the spindle for it C 50 IMAX = -1 IMAXNO = 0 DO 60 K = 1 , NSTEP IF (KSET(K).EQ.0) THEN IOFF = (K-1)*NR + 1 CALL OVERCN(ICNT,IREC,KSAVE(IOFF),NR) IF (ICNT.GT.IMAX) THEN IMAX = ICNT IMAXNO = K END IF END IF 60 CONTINUE C C====== Check for end - no spindle settings remain C IF (IMAXNO.EQ.0) GO TO 20 C C====== Output the best overlap spindle setting and repeat C IOFF = NR*(IMAXNO-1) + 1 CALL OVERST(IREC,KSAVE(IOFF),NR) FR = FLOAT(IMAX)/FLOAT(NR) FRLEFT = FRLEFT - FR C C====== Count the overlaps of this spindle setting with the earlier ones C DO 80 J = 1 , IFILM DO 70 K = 1 , NSTEP IF (KSET(K).EQ.J) THEN IOFF2 = NR*(K-1) + 1 CALL OVERPR(ICNT,KSAVE(IOFF),KSAVE(IOFF2),NR) OV(J) = FLOAT(ICNT)/FLOAT(NR) GO TO 80 END IF 70 CONTINUE 80 CONTINUE WRITE (IPRINT,6030) SP(IMAXNO) , FR , FRLEFT , (OV(J),J=1,IFILM) WRITE (IUNIT6,6030) SP(IMAXNO) , FR , FRLEFT , (OV(J),J=1,IFILM) IFILM = IFILM + 1 KSET(IMAXNO) = IFILM C C====== but only keep going if some unrecordeds remain C====== or if the last spindle setting yielded some overlaps C IF (IMAX.EQ.0) GO TO 20 IF (FRLEFT.LE.0.00001) GO TO 20 GO TO 50 C C====== Format statements C 6010 FORMAT (' *** warning *** only room to save first',I3,' rows') 6020 FORMAT (// + ' Spndl Add. fract. Fract. symm. Fractional overlap' + / + ' symm. uniq. uniques not with previous spindles' + / + ' captured recorded 1 2 3 4 5' + /' 1.000') 6030 FORMAT (F6.1,4X,F6.3,7X,F6.3,2X,5F7.3/31X,5F7.3/31X,5F7.3/31X, + 5F7.3) END **==overpr.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** OVERPR ** C*************************** C C SUBROUTINE OVERPR(ICNT,KSAVE,KSAVE2,NR) C INCLUDE 'implicit_none.finc' C C Purpose: Count in icnt the reflections common to ksave and ksave2 C INTEGER ICNT,NR INTEGER KSAVE(NR),KSAVE2(NR) C C====== LOCALS: C INTEGER I C C------------------------------------------------------------------------------- C ICNT = 0 DO 20 I = 1 , NR IF (KSAVE(I).NE.0 .AND. KSAVE2(I).NE.0) ICNT = ICNT + 1 20 CONTINUE RETURN END **==overst.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** OVERST ** C*************************** C C SUBROUTINE OVERST(IREC,KSAVE,NR) C INCLUDE 'implicit_none.finc' C C Purpose: Set the irec array to show as recorded the reflection in ksave C INTEGER NR INTEGER KSAVE(NR),IREC(NR) C C====== LOCALS: C INTEGER I C C------------------------------------------------------------------------------- C DO 20 I = 1 , NR IF (IREC(I).EQ.1 .AND. KSAVE(I).NE.0) IREC(I) = 0 20 CONTINUE RETURN END **==resoln.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** RESOLN ** C*************************** C C SUBROUTINE RESOLN() C INCLUDE 'implicit_none.finc' C INTEGER MAXONE,MAXALL PARAMETER (MAXONE=20000,MAXALL=100000) C COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /SPOTS/NR REAL CF,DMIN,SPNDLE,WVM,WVS,WVX INTEGER NR,NUMA,NUMC REAL CCELL(6),CELL(6),PHIS(3) C INCLUDE 'chanel.finc' C C====== LOCALS: C INTEGER ION(MAXALL) REAL WRES(MAXALL) REAL DMAX,F1,F2,F3,V INTEGER ICOUNT,J,K,L REAL CON(11) INTEGER IBOX(10) C DATA CON/.0,.10,.20,.35,.55,.80,1.2,1.8,2.8,4.5,100./ C C------------------------------------------------------------------------------- C WRITE (IPRINT,6010) WRITE (IPRINT,6020) WRITE (IUNIT6,6020) C CALL LAUER(ION,WRES) C C====== Find out how many reflections recorded C ICOUNT = 0 DMAX = 0. DO 20 J = 1 , NR K = ION(J) IF (K.NE.0) THEN IF (WRES(J).GT.DMAX) DMAX = WRES(J) ICOUNT = ICOUNT + 1 END IF 20 CONTINUE CON(11) = DMAX - DMIN C WRITE (IPRINT,6030) ICOUNT , FLOAT(ICOUNT)/FLOAT(NR) WRITE (IUNIT6,6030) ICOUNT , FLOAT(ICOUNT)/FLOAT(NR) DO 30 K = 1 , 10 IBOX(K) = 0 30 CONTINUE DO 60 J = 1 , NR IF (ION(J).EQ.0) GO TO 60 V = WRES(J) K = 0 DO 40 L = 1 , 11 IF (V.LT.(DMIN+CON(L))) GO TO 50 K = K + 1 40 CONTINUE K = 10 50 IBOX(K) = IBOX(K) + 1 60 CONTINUE F2 = 0.0 F3 = FLOAT(ICOUNT)/FLOAT(NR) DO 70 K = 1 , 10 IF (K.NE.1) F3 = F3 - F1 F1 = FLOAT(IBOX(K))/FLOAT(NR) F2 = F2 + F1 WRITE (IPRINT,6040) F1 , F2 , F3 , CON(K) + DMIN , CON(K+1) + + DMIN WRITE (IUNIT6,6040) F1 , F2 , F3 , CON(K) + DMIN , CON(K+1) + + DMIN 70 CONTINUE RETURN C C====== Format statements C 6010 FORMAT (' This command calculates the resolution distribution', + ' for'/' the symmetry equivalents in this orientation.') 6020 FORMAT (///' Distribution of resolution of reflections'/) 6030 FORMAT (' Symmetry uniques recorded =',I6,' (',F4.2, + ')'//' Fract. of unique Cumulative Resolution range'/ + ' reflections fraction '/) 6040 FORMAT (F12.3,5X,2F7.3,F9.2,' - ',F6.2) END **==rotspi.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** ROTSPI ** C*************************** C C SUBROUTINE ROTSPI() C INCLUDE 'implicit_none.finc' C INTEGER MAXONE,MAXALL PARAMETER (MAXONE=20000,MAXALL=100000) C COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /SPOTS/NR REAL CF,DMIN,SPNDLE,WVM,WVS,WVX INTEGER NR,NUMA,NUMC REAL CCELL(6),CELL(6),PHIS(3) C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL QUIT INTEGER ION(MAXONE) INTEGER KSAVE(MAXALL) REAL FRACT,SPSAV,STEP INTEGER I,ICOUNT,IST,ITOT,J,K,NSTEP LOGICAL ISAV REAL C(4),SP(20) C C------------------------------------------------------------------------------- C WRITE (IPRINT,6010) NSTEP = 1 CALL TASKI('Enter number of steps (counting the present spindle'// + ' setting as the first)',NSTEP,QUIT) STEP = 10.0 CALL TASKF('Enter step size in degrees',STEP,1,QUIT) WRITE (IPRINT,6020) WRITE (IUNIT6,6020) WRITE (IPRINT,6030) WRITE (IUNIT6,6030) C C====== Loop through the spindle steps C SPSAV = SPNDLE SPNDLE = SPNDLE - STEP ISAV = .TRUE. IST = 0 DO 50 I = 1 , NSTEP SPNDLE = SPNDLE + STEP SP(I) = SPNDLE CALL LAUE(ION) C C====== Find out how many reflections recorded C ICOUNT = 0 ITOT = 0 DO 20 J = 1 , 4 C(J) = 0. 20 CONTINUE DO 30 J = 1 , NR K = ION(J) C C====== If there is room save the ion array for overlap purposes C IF (ISAV) KSAVE(IST+J) = K IF (K.NE.0) THEN ITOT = ITOT + K IF (K.GT.4) K = 4 C(K) = C(K) + 1 ICOUNT = ICOUNT + 1 END IF 30 CONTINUE FRACT = FLOAT(ICOUNT)/FLOAT(NR) DO 40 J = 1 , 4 C(J) = C(J)/NR 40 CONTINUE WRITE (IPRINT,6040) SPNDLE , FRACT , ICOUNT , (C(J),J=1,4) , + ITOT WRITE (IUNIT6,6040) SPNDLE , FRACT , ICOUNT , (C(J),J=1,4) , + ITOT IST = IST + NR IF ((IST+NR).GT.MAXALL) ISAV = .FALSE. 50 CONTINUE SPNDLE = SPSAV C C====== Investigate overlaps? C IF (NSTEP.GT.1) CALL OVERLP(KSAVE,NSTEP,SP) RETURN C C====== Format statements C 6010 FORMAT (' This command analyses the effect of rotating about ', + 'the spindle axis.') 6020 FORMAT (/// + ' Effect of rotating about spindle on the fraction of symmetry' + /' unique reflections recorded') 6030 FORMAT (/// +' Spindle Fraction recorded Distribution of symmetry equivalent +s'/32X,' 1 2 3 4+'/) 6040 FORMAT (F6.1,F13.3,' (',I5,')',4X,4F6.3,' (',I5,')') END **==screen.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** SCREEN ** C*************************** C C SUBROUTINE SCREEN(IENT,ICHANG) C INCLUDE 'implicit_none.finc' C COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT COMMON /RADC/RADI,RADI2 COMMON /SYS/ISYS,KSYS INTEGER ICHANG,IENT REAL CF,DMIN,RADI,RADI2,RSUB,SPNDLE,SSUB,WVM,WVS,WVX,XSUB,YSUB INTEGER ICENT,ILGP,ISYS,MBEAM,MROT,MSYS,NUMA,NUMC REAL CCELL(6),CELL(6),PHIS(3) INTEGER KSYS(3) C INCLUDE 'chanel.finc' C C====== FUNCTION TYPES: C INTEGER LENSTR C C====== LOCALS: C LOGICAL QUIT REAL RAPI,VALUES(5),VAL CHARACTER*4 WORDS(5) INTEGER I,IOK,IOPT,J,NVALS,NWORDS SAVE IOPT LOGICAL LDATEX,LOGG,LSTIEX CHARACTER CHARS*4,DATFIL*60,SAVFIL*60 CHARACTER LABEL(35)*6,NAME(5)*15,PROMPT*50 C C------------------------------------------------------------------------------- C DATA LABEL/'a* ','b* ','c* ','alpha*','beta* ','gamma*', + 'wmin ','wmax ','dist ','dmin ','radi ','spdr ', + 'bmdr ','px ','py ','pz ','spin ','syst ', + 'scal ','dats ','datg ','laue ','latt ','rotate', + 'search','wavvar','varcf ','resoln','offrot','tworot', + 'stop ',' ','log ',' ',' '/ C C====== First entry is ient=0. on subsequent entries, jump to menu. C IF (IENT.GT.0) THEN ICHANG = 0 GO TO 30 END IF ICHANG = 1 C C====== Switch logging off initially C LOGG = .FALSE. IUNIT6 = 2 C C====== Set default menu selection to `DATG' C IOPT = 21 - 1 WORDS(1) = 'DATG' C C====== Initialise some values C C C====== Cell parameters to silly values to show which are changed C DO 20 I = 1 , 6 CELL(I) = 0.0 20 CONTINUE C C====== Matrix to zero. msys=system, mrot=rotation axis, mbeam=beam axis C MSYS = 0 MROT = 3 MBEAM = 1 ISYS = 0 ICENT = 1 KSYS(1) = 0 KSYS(2) = 0 KSYS(3) = 0 ILGP = -1 C C====== Geometry to 0.0 C PHIS(1) = 0.0 PHIS(2) = 0.0 PHIS(3) = 0.0 WVM = 0.3 WVX = 2.6 CF = 50. WVS = 5.0 C C====== Default spindle setting to 0.0, dmin to 2.5, radius to 59.3 mm. C SPNDLE = 0.0 DMIN = 2.5 RADI = 59.3 RADI2 = RADI*RADI C C====== Display menu and get next command C 30 CALL UPDMEN(LABEL,LOGG) IOPT = IOPT + 1 CALL MENU(WORDS,NWORDS,VALUES,NVALS,IOPT) C C====== Reciprocal cell parameters C IF (IOPT.LE.6) THEN IF (IOPT.LE.3) WRITE (PROMPT,6010) LABEL(IOPT) IF (IOPT.GT.3) WRITE (PROMPT,6020) LABEL(IOPT) CALL INPUTV(PROMPT,CELL(IOPT),VALUES,NVALS) ICHANG = 1 GO TO 30 END IF C C====== Various limits C IF (IOPT.EQ.7) THEN CALL INPUTV('Wavelength minimum (angstroms): ',WVM,VALUES, + NVALS) ICHANG = 1 GO TO 30 END IF IF (IOPT.EQ.8) THEN CALL INPUTV('Wavelength maximum (angstroms): ',WVX,VALUES, + NVALS) ICHANG = 1 GO TO 30 END IF IF (IOPT.EQ.9) THEN CALL INPUTV('Crystal-to-film distance in mm: ',CF,VALUES,NVALS) GO TO 30 END IF IF (IOPT.EQ.10) THEN CALL INPUTV('Enter dmin - resolution limit in Angstroms: ', + DMIN,VALUES,NVALS) ICHANG = 1 GO TO 30 END IF IF (IOPT.EQ.11) THEN CALL INPUTV('Enter radius limit for plot - mm: ',RADI,VALUES, + NVALS) ICHANG = 1 RADI2 = RADI*RADI GO TO 30 END IF C C====== Orientation C C C====== Spindle axis C IF (IOPT.EQ.12) THEN CALL TWRITE('Axis 1 = astar, -3 = -cstar etc.') CALL INPUTV('Reciprocal axis lying along crystal rotation'// + ' axis: ',VAL,VALUES,NVALS) MROT = VAL IF (MROT.EQ.0 .OR. MROT.LT.-3 .OR. MROT.GT.3) THEN WRITE (IPRINT,6030) MROT = 3 END IF GO TO 30 END IF C C====== Beam axis C IF (IOPT.EQ.13) THEN CALL TWRITE('Reciprocal axis lying nearest to X-ray beam'// + ' axis in the plane') CALL INPUTV('defined by the beam axis and the rotation axis: ', + VAL,VALUES,NVALS) MBEAM = VAL IF (MBEAM.EQ.0 .OR. MBEAM.LT.-3 .OR. MBEAM.GT.3 .OR. + IABS(MBEAM).EQ.IABS(MROT)) THEN WRITE (IPRINT,6030) MBEAM = 1 END IF GO TO 30 END IF C C====== Phi angles C IF (IOPT.EQ.14) THEN CALL INPUTV('Rotation about beam axis, phi(x) in degrees: ', + PHIS(1),VALUES,NVALS) GO TO 30 END IF IF (IOPT.EQ.15) THEN CALL INPUTV('Rotation about vertical axis, phi(y) in degrees: ' + ,PHIS(2),VALUES,NVALS) GO TO 30 END IF IF (IOPT.EQ.16) THEN CALL INPUTV( + 'Rotation about horizontal axis, phi(z) in degrees: ' + ,PHIS(3),VALUES,NVALS) GO TO 30 END IF C C====== Miscellaneous film parameters C IF (IOPT.EQ.17) THEN CALL INPUTV('Spindle setting: ',SPNDLE,VALUES,NVALS) GO TO 30 END IF C C====== Lattice input C IF (IOPT.EQ.18) THEN CALL INPUTW('Crystal system: ',CHARS,WORDS,NWORDS) IF (CHARS.EQ.'TRIC') THEN MSYS = 1 ELSE IF (CHARS.EQ.'MONO') THEN MSYS = 2 ELSE IF (CHARS.EQ.'ORTH') THEN MSYS = 3 ELSE IF (CHARS.EQ.'TETR') THEN MSYS = 4 ELSE IF (CHARS.EQ.'HEXA') THEN MSYS = 5 ELSE IF (CHARS.EQ.'RHOM') THEN MSYS = 6 ELSE IF (CHARS.EQ.'CUBI') THEN MSYS = 7 ELSE WRITE (IPRINT,6040) MSYS = 0 GO TO 30 END IF ICHANG = 1 C C====== May as well get user to set Laue group now C CALL LGRP(MSYS,ILGP) GO TO 30 END IF IF (IOPT.EQ.19) THEN CALL INPUTV( + 'Scale factor that multiplies recip. cell lengths: ' + ,WVS,VALUES,NVALS) ICHANG = 1 GO TO 30 END IF C C====== Save file for control data C IF (IOPT.EQ.20) THEN CALL TPROMP( + 'Enter name of file in which to save control data: ') CALL TREAD(SAVFIL,QUIT) C C====== Look for first blank, append .dat if "." not seen C I = 0 40 I = I + 1 IF (SAVFIL(I:I).NE.' ' .AND. SAVFIL(I:I).NE.'.') GO TO 40 IF (SAVFIL(I:I).EQ.' ') THEN SAVFIL(I:I+3) = '.dat' END IF OPEN (UNIT=IUNIT7,FILE=SAVFIL,STATUS='UNKNOWN', + FORM='FORMATTED') C C====== Write the records ... for contents see the comments C====== above concerning reading them back. C WRITE (IUNIT7,6050) PHIS , CELL WRITE (IUNIT7,6060) WVM , WVX , WVS , CF , SPNDLE , DMIN , RADI WRITE (IUNIT7,6070) MSYS , MROT , MBEAM , ISYS , KSYS , ILGP , + ICENT WRITE (IUNIT7,6080) C C====== Final mosaic spread record C REWIND (UNIT=IUNIT7) CLOSE (UNIT=IUNIT7) GO TO 30 END IF C C====== File definition C IF (IOPT.EQ.21) THEN CALL TPROMP( + 'Enter name of previously saved control data file: ') CALL TREAD(DATFIL,QUIT) C C====== Look for first blank, append .dat if "." not seen C I = 0 50 I = I + 1 IF (DATFIL(I:I).NE.' ' .AND. DATFIL(I:I).NE.'.') GO TO 50 IF (DATFIL(I:I).EQ.' ') THEN DATFIL(I:I+3) = '.dat' END IF INQUIRE (FILE=DATFIL,EXIST=LDATEX) IF (.NOT.LDATEX) THEN WRITE (IPRINT,6090) DATFIL(1:LENSTR(DATFIL)) GO TO 30 END IF OPEN (UNIT=IUNIT7,FILE=DATFIL,STATUS='OLD') C C====== Record 1: orientation angles, cell parameters C READ (IUNIT7,*) PHIS , CELL C C====== Record 2: wavelengths, scale, crystal-to-film distance, dmin, radius. C READ (IUNIT7,*) WVM , WVX , WVS , CF , SPNDLE , DMIN , RADI RADI2 = RADI*RADI C C====== Record 3: matrix - msys, mrot, mbeam, isys, ksys, ilgp, icent C READ (IUNIT7,*) MSYS , MROT , MBEAM , ISYS , KSYS , ILGP , + ICENT C REWIND (UNIT=IUNIT7) CLOSE (UNIT=IUNIT7) ICHANG = 1 GO TO 30 END IF C C====== Laue group C IF (IOPT.EQ.22) THEN IF (MSYS.EQ.0) THEN WRITE (IPRINT,6100) GO TO 30 END IF CALL LGRP(MSYS,ILGP) ICHANG = 1 GO TO 30 END IF C C====== Lattice centering C IF (IOPT.EQ.23) THEN CALL INPUTW('Lattice centering: P,A,B,C,I,F or R [P]: ',CHARS, + WORDS,NWORDS) C C====== Default lattice centering is p C IF (CHARS.EQ.' ') CHARS = 'P' C C====== R only applies to rhombohedral indexed on hexagonal lattice C IF (CHARS.EQ.'R' .AND. MSYS.NE.5) CHARS = 'P' IF (CHARS.EQ.'P') THEN ISYS = 0 ICENT = 1 KSYS(1) = 0 KSYS(2) = 0 KSYS(3) = 0 ICHANG = 1 GO TO 30 END IF IF (CHARS.EQ.'A') THEN ISYS = 2 ICENT = 2 KSYS(2) = 1 KSYS(3) = 1 ICHANG = 1 GO TO 30 END IF IF (CHARS.EQ.'B') THEN ISYS = 2 ICENT = 3 KSYS(1) = 1 KSYS(3) = 1 ICHANG = 1 GO TO 30 END IF IF (CHARS.EQ.'C') THEN ISYS = 2 ICENT = 4 KSYS(1) = 1 KSYS(2) = 1 ICHANG = 1 GO TO 30 END IF IF (CHARS.EQ.'I') THEN ISYS = 2 ICENT = 5 KSYS(1) = 1 KSYS(2) = 1 KSYS(3) = 1 ICHANG = 1 GO TO 30 END IF IF (CHARS.EQ.'F') THEN ISYS = 4 ICENT = 6 ICHANG = 1 GO TO 30 END IF IF (CHARS.EQ.'R') THEN ISYS = 3 ICENT = 7 KSYS(1) = -1 KSYS(2) = 1 KSYS(3) = 1 ICHANG = 1 GO TO 30 END IF WRITE (IPRINT,6110) ISYS = 0 ICENT = 1 KSYS(1) = 0 KSYS(2) = 0 KSYS(3) = 0 GO TO 30 END IF C C====== Stop signal. C IF (IOPT.EQ.31) THEN CALL QSSCRL(1,24) CALL QSCLRA STOP END IF C C====== Help information. (Sorry - removed now) C IF (IOPT.EQ.32) GO TO 30 C C====== Log file status. C IF (IOPT.EQ.33) THEN CALL INPUTW('Output to log file: ',CHARS,WORDS,NWORDS) IF (CHARS.EQ.'ON' .OR. CHARS(1:1).EQ.'Y') THEN LOGG = .TRUE. IUNIT6 = 1 ELSE IF (CHARS.EQ.'OFF' .OR. CHARS(1:1).EQ.'N') THEN LOGG = .FALSE. IUNIT6 = 2 ELSE WRITE (IPRINT,6120) CHARS END IF IOPT = 20 GO TO 30 END IF C C====== Check things before proceeding C IOK = 0 IF (CF.EQ.0.0) THEN WRITE (IPRINT,6130) IOK = 1 END IF C C====== Check cell parameters compatible with lattice type C IF (MSYS.GT.0) CALL CELLTEST(MSYS,ILGP,ICENT,CELL,LABEL,IOK) C C====== Copy cell parameters to contiguous array C RAPI = 3.14159265/180. NUMC = 0 NUMA = 0 J = 0 DO 60 I = 1 , 6 IF (CELL(I).NE.0.) THEN J = J + 1 CCELL(J) = CELL(I) NUMC = NUMC + 1 IF (I.GT.3) THEN NUMA = NUMA + 1 CCELL(J) = CCELL(J)*RAPI END IF END IF 60 CONTINUE IF (NUMC.EQ.0) THEN WRITE (IPRINT,6140) IOK = 1 END IF C C====== Check for matrix C IF (MSYS.EQ.0) THEN WRITE (IPRINT,6150) IOK = 1 END IF IF (ILGP.EQ.-1) THEN WRITE (IPRINT,6160) IOK = 1 END IF IF (IOK.NE.0) GO TO 30 C C====== Return to main routine to execute the command C IENT = IOPT - 23 RETURN C C====== Format statements C 6010 FORMAT ('Recip. lattice length (dimensionless) ',A6,': ') 6020 FORMAT ('Reciprocal lattice angle (degrees) ',A6,': ') 6030 FORMAT (' !!! Invalid axis') 6040 FORMAT (' !!! Invalid crystal lattice') 6050 FORMAT (3F12.6,3F10.6/3F10.4) 6060 FORMAT (3F10.6,F10.4,2F10.2,F10.2) 6070 FORMAT (9I3) 6080 FORMAT (' 0.2000') 6090 FORMAT (' !!! File ',A,' does not exist') 6100 FORMAT (' Please enter crystal lattice type first!') 6110 FORMAT (' !!! Invalid lattice centering') 6120 FORMAT (' !!! Illegal choice: ',A4) 6130 FORMAT (' !!! No value for crystal to film distance') 6140 FORMAT (' !!! No cell parameters') 6150 FORMAT (' !!! No lattice information') 6160 FORMAT (' !!! No laue group') END **==search.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** SEARCH ** C*************************** C C SUBROUTINE SEARCH() C INCLUDE 'implicit_none.finc' C INTEGER MAXONE,MAXALL PARAMETER (MAXONE=20000,MAXALL=100000) C COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /SPOTS/NR REAL CF,DMIN,SPNDLE,WVM,WVS,WVX INTEGER NR,NUMA,NUMC REAL CCELL(6),CELL(6),PHIS(3) C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL QUIT INTEGER ION(MAXONE) REAL PHIYSV,SPSAV,STEPY,STEPZ INTEGER I,ICOUNT,J,K,M,NSTEPY,NSTEPZ REAL FR(30) C C------------------------------------------------------------------------------- C WRITE (IPRINT,6010) WRITE (IUNIT6,6020) CALL TWRITE(' ') NSTEPY = 1 CALL TWRITE('Enter number of steps about the vertical axis ') CALL TASKI('(counting the present phi-y as the first)',NSTEPY, + QUIT) STEPY = 10.0 CALL TASKF('Enter step size in degrees',STEPY,1,QUIT) NSTEPZ = 1 CALL TWRITE('Enter number of steps about the spindle axis ') CALL TASKI('(counting the present phi-z as the first)',NSTEPZ, + QUIT) STEPZ = 10.0 CALL TASKF('Enter step size in degrees',STEPZ,1,QUIT) C SPSAV = SPNDLE SPNDLE = 0. PHIYSV = PHIS(2) WRITE (IPRINT,6030) WRITE (IUNIT6,6030) WRITE (IPRINT,6040) ((PHIS(3)+(I-1)*STEPZ),I=1,NSTEPZ) WRITE (IUNIT6,6040) ((PHIS(3)+(I-1)*STEPZ),I=1,NSTEPZ) WRITE (IPRINT,6050) WRITE (IUNIT6,6050) C C====== Loop through the vertical axis steps C PHIS(2) = PHIS(2) - STEPY DO 50 M = 1 , NSTEPY PHIS(2) = PHIS(2) + STEPY C C====== Recalculate the a matrix since phis have changed C CELL(1) = CELL(1)*WVS CELL(2) = CELL(2)*WVS CELL(3) = CELL(3)*WVS J = 0 DO 20 I = 1 , 3 IF (CELL(I).NE.0.) THEN J = J + 1 CCELL(J) = CELL(I) END IF 20 CONTINUE CALL ACALC CELL(1) = CELL(1)/WVS CELL(2) = CELL(2)/WVS CELL(3) = CELL(3)/WVS C C====== Loop through the spindle steps C SPNDLE = -STEPZ DO 40 I = 1 , NSTEPZ SPNDLE = SPNDLE + STEPZ CALL LAUE(ION) C C====== Find out how many reflections recorded C ICOUNT = 0 DO 30 J = 1 , NR K = ION(J) IF (K.NE.0) THEN ICOUNT = ICOUNT + 1 END IF 30 CONTINUE FR(I) = FLOAT(ICOUNT)/FLOAT(NR) 40 CONTINUE WRITE (IPRINT,6060) PHIS(2) , (FR(I),I=1,NSTEPZ) WRITE (IUNIT6,6060) PHIS(2) , (FR(I),I=1,NSTEPZ) 50 CONTINUE PHIS(2) = PHIYSV SPNDLE = SPSAV RETURN C C====== Format statements C 6010 FORMAT (' This command searches the vertical axis and spindle', + ' axis (phi-y, phi-z)'/ + ' to explore the fraction of symmetry', + ' unique reflections recorded.'// + ' Whatever its value now, the', + ' spindle setting will be taken as zero.') 6020 FORMAT (' Whatever its value now, the', + ' spindle setting will be taken as zero.') 6030 FORMAT (///' Fraction of symmetry unique reflections recorded') 6040 FORMAT (//' phi-z',9F6.1/12X,9F6.1/12X,9F6.1) 6050 FORMAT (' phi-y') 6060 FORMAT (1X,F6.1,5X,9F6.3/12X,9F6.3/12X,9F6.3) END **==sysabs.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** SYSABS ** C*************************** C C SUBROUTINE SYSABS(IH,IK,IL,ABSNT) C INCLUDE 'implicit_none.finc' C COMMON /SYS/ISYS,KSYS INTEGER IH,IK,IL LOGICAL ABSNT INTEGER ISYS INTEGER KSYS(3) C C====== LOCALS: C INTEGER I,J INTEGER H(3) C C------------------------------------------------------------------------------- C H(1) = IH H(2) = IK H(3) = IL IF (ISYS.EQ.4) THEN ABSNT = .TRUE. DO 20 I = 1 , 2 IF (MOD(H(I)+H(I+1),2).NE.0) RETURN 20 CONTINUE ABSNT = .FALSE. RETURN ELSE C C====== This version for a,b,c,r,i & f lattice absences C ABSNT = .FALSE. J = 0 DO 30 I = 1 , 3 J = J + H(I)*KSYS(I) 30 CONTINUE IF (MOD(J,ISYS).NE.0) ABSNT = .TRUE. RETURN END IF END **==sysahk.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** SYSAHK ** C*************************** C C SUBROUTINE SYSAHK(IH,IK,ABSNT) C INCLUDE 'implicit_none.finc' C C Purpose: This variant of sysabs checks to see if a given h,k is absent, C whatever the l value C INTEGER IH,IK LOGICAL ABSNT INTEGER ISYS INTEGER KSYS(3) COMMON /SYS/ISYS,KSYS C C------------------------------------------------------------------------------- C ABSNT = .FALSE. IF ((ISYS.EQ.2 .OR. ISYS.EQ.4) .AND. KSYS(3).EQ.0) THEN IF (MOD(IH+IK,2).NE.0) ABSNT = .TRUE. END IF RETURN END **==tworot.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** TWOROT ** C*************************** C C SUBROUTINE TWOROT C INCLUDE 'implicit_none.finc' C REAL DEGRAD PARAMETER (DEGRAD=3.14159265358979/180.) INTEGER MAXONE, MAXALL PARAMETER (MAXONE=20000,MAXALL=100000) C COMMON /DATA/CELL, WVM, WVX, WVS, CF, PHIS, SPNDLE, CCELL, NUMA, + NUMC, DMIN COMMON /SPOTS/NR REAL CF, DMIN, SPNDLE, WVM, WVS, WVX INTEGER NR, NUMA, NUMC REAL CCELL(6), CELL(6), PHIS(3) C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL QUIT INTEGER ION(MAXONE) INTEGER KSAVE(MAXALL) REAL PHISAV(3), COSPHIY REAL FRACT, USTEP, PHIU, MU, MUSTEP INTEGER I, IU, IMU, ICOUNT, IST, ITOT, J, K, NU, NMU, NALL, IALL LOGICAL ISAV REAL C(4), SP(20), ANGSAV(20,2) C C------------------------------------------------------------------------------- C WRITE (IPRINT,6010) NU = 1 CALL TASKI('Enter number of phi_U steps (counting a spindle'// + ' setting of 0 as the first)',NU,QUIT) USTEP = 10.0 CALL TASKF('Enter step size in degrees',USTEP,1,QUIT) NMU = 1 CALL TASKI('Enter number of µ steps (counting µ=0 as the first)', + NMU,QUIT) MUSTEP = 10.0 CALL TASKF('Enter step size in degrees',MUSTEP,1,QUIT) NALL = NMU*NU IF (NALL.GT.20) THEN WRITE (IPRINT,6020) NALL WRITE (IUNIT6,6020) NALL ELSE WRITE (IPRINT,6030) WRITE (IUNIT6,6030) WRITE (IPRINT,6040) WRITE (IUNIT6,6040) C DO 20 I = 1 , 3 PHISAV(I) = PHIS(I) 20 CONTINUE ISAV = .TRUE. IST = 0 C C====== Loop through the spindle steps C DO 90 IU = 1 , NU PHIU = (IU-1)*USTEP PHIU = PHIU*DEGRAD DO 80 IMU = 1 , NMU IALL = (IU-1)*NMU + IMU SP(IALL) = FLOAT(IALL) MU = (IMU-1)*MUSTEP MU = MU*DEGRAD ANGSAV(IALL,1) = PHIU ANGSAV(IALL,2) = MU PHIS(2) = ASIN(SIN(MU)*COS(PHIU)) COSPHIY = COS(PHIS(2)) PHIS(3) = ASIN(SIN(PHIU)/COSPHIY) PHIS(1) = ACOS(COS(MU)/COSPHIY) DO 30 J = 1 , 3 PHIS(J) = PHIS(J)/DEGRAD + PHISAV(J) 30 CONTINUE C C====== Recalculate the A matrix since phis have changed C CELL(1) = CELL(1)*WVS CELL(2) = CELL(2)*WVS CELL(3) = CELL(3)*WVS J = 0 DO 40 K = 1 , 3 IF (CELL(K).NE.0.) THEN J = J + 1 CCELL(J) = CELL(K) END IF 40 CONTINUE CALL ACALC CELL(1) = CELL(1)/WVS CELL(2) = CELL(2)/WVS CELL(3) = CELL(3)/WVS CALL LAUE(ION) C C====== Find out how many reflections recorded C ICOUNT = 0 ITOT = 0 DO 50 J = 1 , 4 C(J) = 0. 50 CONTINUE DO 60 J = 1 , NR K = ION(J) C C====== If there is room save the ion array for overlap purposes C IF (ISAV) KSAVE(IST+J) = K IF (K.NE.0) THEN ITOT = ITOT + K IF (K.GT.4) K = 4 C(K) = C(K) + 1 ICOUNT = ICOUNT + 1 END IF 60 CONTINUE FRACT = FLOAT(ICOUNT)/FLOAT(NR) DO 70 J = 1 , 4 C(J) = C(J)/NR 70 CONTINUE WRITE (IPRINT,6050) IALL , FRACT , ICOUNT , (C(J),J=1,4) + , ITOT WRITE (IUNIT6,6050) IALL , FRACT , ICOUNT , (C(J),J=1,4) + , ITOT IST = IST + NR IF ((IST+NR).GT.MAXALL) ISAV = .FALSE. 80 CONTINUE 90 CONTINUE DO 100 I = 1 , 3 PHIS(I) = PHISAV(I) 100 CONTINUE WRITE (IPRINT,6060) WRITE (IUNIT6,6060) DO 110 I = 1 , NALL WRITE (IPRINT,6070) I , ANGSAV(I,1)/DEGRAD , ANGSAV(I,2) + /DEGRAD WRITE (IUNIT6,6070) I , ANGSAV(I,1)/DEGRAD , ANGSAV(I,2) + /DEGRAD 110 CONTINUE C C====== investigate overlaps? C IF (NALL.GT.1) CALL OVERLP(KSAVE,NALL,SP) END IF C C====== Format statements C 6010 FORMAT (' This command analyses the effect of changing '/ + ' both the spindle axis setting phi_U', + ' and the spindle offset angle µ.') 6020 FORMAT (' Too many steps altogether (',I3,'), maximum 20') 6030 FORMAT (///' Effect of changing phi_U and µ', + ' on the fraction of symmetry'/ + ' unique reflections recorded') 6040 FORMAT (/// +' Spindle# Fraction recorded Distribution of symmetry equivalent +s'/32X,' 1 2 3 4+'/) 6050 FORMAT (I6,F13.3,' (',I5,')',4X,4F6.3,' (',I5,')') 6060 FORMAT (//' Table of spindle numbers and the corresponding'/ + ' phi_U, µ settings'//' Spindle# phi_U µ') 6070 FORMAT (1X,I3,9X,2F8.1) END **==unique.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** UNIQUE ** C*************************** C C SUBROUTINE UNIQUE(IH,IK,IL,UNIQ) C INCLUDE 'implicit_none.finc' C C Purpose: This routine returns uniq=true if reflection h,k,l is in the C unique set according to the Laue symmetry, otherwise iflag=false. C C To be compatible with the order of reflection generation in hklgen C this routine must ensure that it does not reject all l for given k C (which causes hklgen to stop incrementing k) if a subsequent k C would be in the unique set. Similarly, do not reject all k for given C h if a subsequent would be in the unique set. C COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT INTEGER IH,IK,IL LOGICAL UNIQ INTEGER ICENT,ILGP,MBEAM,MROT,MSYS INTRINSIC IABS C C------------------------------------------------------------------------------- C UNIQ = .FALSE. IF (ILGP.EQ.0) THEN C C====== 1bar C IF (IH.LT.0) RETURN IF (IH.EQ.0) THEN IF (IK.LT.0) RETURN IF (IK.EQ.0) THEN IF (IL.LT.0) RETURN END IF END IF GO TO 30 END IF C IF (ILGP.EQ.2) THEN C C====== mmm C IF (IH.LT.0) RETURN IF (IK.LT.0) RETURN IF (IL.LT.0) RETURN ELSE IF (ILGP.EQ.3) THEN C C====== 4/m C IF (IH.LT.0) RETURN IF (IK.LT.0) RETURN IF (IH.EQ.0 .AND. IK.NE.0) RETURN IF (IL.LT.0) RETURN ELSE IF (ILGP.EQ.4) THEN C C====== 4/mmm C IF (IH.LT.0) RETURN IF (IK.GT.IH) RETURN IF (IK.LT.0) RETURN IF (IL.LT.0) RETURN ELSE IF (ILGP.EQ.5) THEN C C====== 6/m C IF (IH.LT.0) RETURN IF (IK.LT.0) RETURN IF (IL.LT.0) RETURN IF (IH.EQ.0 .AND. IK.NE.0) RETURN ELSE IF (ILGP.EQ.6) THEN C C====== 6/mmm C IF (IH.LT.0) RETURN IF (IK.LT.0) RETURN IF (IK.GT.IH) RETURN IF (IL.LT.0) RETURN ELSE IF (ILGP.EQ.7) THEN C C====== 3bar ( hexagonal or rhombohedral lattices possible) C IF (MSYS.EQ.5) THEN IF (IH.LT.0) RETURN IF (IK.LT.0) RETURN IF (IH.EQ.0) THEN IF (IK.NE.0) RETURN IF (IL.LT.0) RETURN END IF ELSE IF (MSYS.EQ.6) THEN IF (IH.LE.0) RETURN IF (IABS(IL).GT.0) THEN IF (IABS(IL).EQ.IH .AND. IK.NE.IH) RETURN IF (IABS(IL).LT.IH .AND. + (IABS(IK).LT.IABS(IL) .OR. IABS(IK).GE.IH)) RETURN IF (IABS(IL).GT.IH .AND. + (IABS(IK).GT.IABS(IL) .OR. IABS(IK).LE.IH)) RETURN END IF END IF ELSE IF (ILGP.EQ.8) THEN GO TO 20 ELSE IF (ILGP.EQ.9) THEN GO TO 20 ELSE IF (ILGP.EQ.10) THEN C C====== m3 C IF (IH.LE.0) RETURN IF (IL.LT.0) RETURN IF (IK.LT.0) RETURN IF (IL.GT.0) THEN IF (IL.EQ.IH .AND. IK.NE.IH) RETURN IF (IL.LT.IH .AND. (IK.LT.IL .OR. IK.GE.IH)) RETURN IF (IL.GT.IH .AND. (IK.GT.IL .OR. IK.LE.IH)) RETURN END IF ELSE IF (ILGP.EQ.11) THEN C C====== m3m C IF (IH.LE.0) RETURN IF (IL.LT.0) RETURN IF (IK.LT.0) RETURN IF (IK.GT.IH) RETURN IF (IK.EQ.0 .AND. IL.NE.0) RETURN IF (IK.NE.0 .AND. IK.NE.IH .AND. IL.GE.IK) RETURN ELSE C C====== 2/m C IF (IH.LT.0) RETURN IF (IK.LT.0) RETURN IF (IH.EQ.0 .AND. IL.LT.0) RETURN END IF GO TO 30 C C====== 3m1 and 31m ( hexagonal or rhombohedral lattices possible) C 20 IF (MSYS.EQ.5) THEN IF (IH.LT.0) RETURN IF (IK.LT.0) RETURN IF (IK.GT.IH) RETURN IF (IH.EQ.0 .AND. IL.LT.0) RETURN IF (IK.EQ.0 .AND. ILGP.EQ.9 .AND. IL.LT.0) RETURN IF (IK.EQ.IH .AND. ILGP.EQ.8 .AND. IL.LT.0) RETURN ELSE IF (MSYS.EQ.6) THEN STOP ' Sorry this case not yet implemented' END IF C 30 UNIQ = .TRUE. RETURN END **==unpack.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** UNPACK ** C*************************** C C SUBROUTINE UNPACK(N,H,K,L) C INCLUDE 'implicit_none.finc' C INTEGER H,K,L,N C C------------------------------------------------------------------------------- C L = MOD(N,200) - 100 K = MOD(N/200,200) - 100 H = N/40000 - 100 RETURN END **==updmen.f C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** UPDMEN ** C*************************** C C SUBROUTINE UPDMEN(LABEL,LOGG) C INCLUDE 'implicit_none.finc' C C Purpose: Updates option menu with associated values C CHARACTER *6 LABEL(34) LOGICAL LOGG C C LABEL : NAMES OF OPTIONS C INCLUDE 'chanel.finc' C COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /MATRIX/MSYS,MROT,MBEAM,ILGP,ICENT COMMON /RADC/RADI,RADI2 COMMON /SYS/ISYS,KSYS REAL CF,DMIN,RADI,RADI2,SPNDLE,WVM,WVS,WVX INTEGER ICENT,ILGP,ISYS,MBEAM,MROT,MSYS,NUMA,NUMC REAL CCELL(6),CELL(6),PHIS(3) INTEGER KSYS(3) C C====== LOCALS: C CHARACTER *15 TEMP INTEGER I CHARACTER MAXIS(7)*7,MSYST(8)*7,MCENT(7)*7,MLAUE(13)*6 C DATA MAXIS/' -cs ',' -bs ',' -as ',' ',' +as ', + ' +bs ',' +cs '/ DATA MSYST/' ','tricl. ','monocl.','ortho. ','tetrag.', + 'hexag. ','rhomb. ','cubic '/ DATA MCENT/' P ',' A ',' B ',' C ',' I ', + ' F ',' R '/ DATA MLAUE/' ','1-bar ',' 2/m ',' mmm ',' 4/m ','4/mmm ', + ' 6/m ','6/mmm ','3-bar ',' 3m1 ',' 31m ',' m3 ', + ' m3m '/ C C------------------------------------------------------------------------------- C DO 20 I = 1 , 6 IF (CELL(I).EQ.0.) THEN WRITE (TEMP,6030) LABEL(I) CALL ITEM(I,TEMP,15,.FALSE.) ELSE WRITE (TEMP,6040) LABEL(I) , CELL(I) CALL ITEM(I,TEMP,15,.FALSE.) END IF C 20 CONTINUE C WRITE (TEMP,6040) LABEL(7) , WVM CALL ITEM(7,TEMP,15,.FALSE.) C WRITE (TEMP,6040) LABEL(8) , WVX CALL ITEM(8,TEMP,15,.FALSE.) C WRITE (TEMP,6010) LABEL(9) , CF CALL ITEM(9,TEMP,15,.FALSE.) C WRITE (TEMP,6010) LABEL(10) , DMIN CALL ITEM(10,TEMP,15,.FALSE.) C WRITE (TEMP,6010) LABEL(11) , RADI CALL ITEM(11,TEMP,15,.FALSE.) RADI2 = RADI*RADI C WRITE (TEMP,6020) LABEL(12) , MAXIS(MROT+4) CALL ITEM(12,TEMP,15,.FALSE.) C WRITE (TEMP,6020) LABEL(13) , MAXIS(MBEAM+4) CALL ITEM(13,TEMP,15,.FALSE.) C DO 30 I = 1 , 3 WRITE (TEMP,6010) LABEL(13+I) , PHIS(I) CALL ITEM(13+I,TEMP,15,.FALSE.) 30 CONTINUE C WRITE (TEMP,6010) LABEL(17) , SPNDLE CALL ITEM(17,TEMP,15,.FALSE.) C WRITE (TEMP,6020) LABEL(18) , MSYST(MSYS+1) CALL ITEM(18,TEMP,15,.FALSE.) C WRITE (TEMP,6010) LABEL(19) , WVS CALL ITEM(19,TEMP,15,.FALSE.) C WRITE (TEMP,6030) LABEL(20) CALL ITEM(20,TEMP,15,.TRUE.) C WRITE (TEMP,6030) LABEL(21) CALL ITEM(21,TEMP,15,.TRUE.) C WRITE (TEMP,6020) LABEL(22) , MLAUE(ILGP+2) CALL ITEM(22,TEMP,15,.FALSE.) C WRITE (TEMP,6020) LABEL(23) , MCENT(ICENT) CALL ITEM(23,TEMP,15,.FALSE.) C C====== Write out the commands. C DO 40 I = 24 , 32 WRITE (TEMP,6030) LABEL(I) CALL ITEM(I,TEMP,15,.TRUE.) 40 CONTINUE C IF (LOGG) THEN WRITE (TEMP,6020) LABEL(33) , 'ON' ELSE WRITE (TEMP,6020) LABEL(33) , 'OFF' END IF CALL ITEM(33,TEMP,15,.TRUE.) C C====== Now write out parameters to log file C IF (LOGG) THEN DO 50 I = 1 , 6 IF (CELL(I).EQ.0.) THEN WRITE (IUNIT6,6070) LABEL(I) ELSE WRITE (IUNIT6,6080) LABEL(I) , CELL(I) END IF C 50 CONTINUE C WRITE (IUNIT6,6080) LABEL(7) , WVM C WRITE (IUNIT6,6080) LABEL(8) , WVX C WRITE (IUNIT6,6050) LABEL(9) , CF C WRITE (IUNIT6,6050) LABEL(10) , DMIN C WRITE (IUNIT6,6050) LABEL(11) , RADI C WRITE (IUNIT6,6060) LABEL(12) , MAXIS(MROT+4) C WRITE (IUNIT6,6060) LABEL(13) , MAXIS(MBEAM+4) C DO 60 I = 1 , 3 WRITE (IUNIT6,6050) LABEL(13+I) , PHIS(I) 60 CONTINUE C WRITE (IUNIT6,6050) LABEL(17) , SPNDLE C WRITE (IUNIT6,6060) LABEL(18) , MSYST(MSYS+1) C WRITE (IUNIT6,6050) LABEL(19) , WVS C WRITE (IUNIT6,6060) LABEL(22) , MLAUE(ILGP+2) C WRITE (IUNIT6,6060) LABEL(23) , MCENT(ICENT) END IF RETURN C C====== Format statements C 6010 FORMAT (A6,F8.2,1X) 6020 FORMAT (A6,1X,A7,1X) 6030 FORMAT (A6,9X) 6040 FORMAT (A6,F8.4,1X) 6050 FORMAT (1X,A6,F8.2) 6060 FORMAT (1X,A6,1X,A7) 6070 FORMAT (1X,A6) 6080 FORMAT (1X,A6,F8.4) END **==varycf.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** VARYCF ** C*************************** C C SUBROUTINE VARYCF() C INCLUDE 'implicit_none.finc' C INTEGER MAXONE,MAXALL PARAMETER (MAXONE=20000,MAXALL=100000) C COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /SPOTS/NR REAL CF,DMIN,SPNDLE,WVM,WVS,WVX INTEGER NR,NUMA,NUMC REAL CCELL(6),CELL(6),PHIS(3) C INCLUDE 'chanel.finc' C C====== LOCALS: C LOGICAL QUIT INTEGER ION(MAXONE) REAL CFSAV,FRACT,STEP INTEGER I,ICOUNT,ITOT,J,K,NSTEP C C------------------------------------------------------------------------------- C WRITE (IPRINT,6010) NSTEP = 1 CALL TASKI('Enter number of steps (counting the present c-to-f'// + ' setting as the first)',NSTEP,QUIT) STEP = 5.0 CALL TASKF('Enter step size in mm',STEP,1,QUIT) WRITE (IPRINT,6020) WRITE (IUNIT6,6020) WRITE (IPRINT,6030) WRITE (IUNIT6,6030) C C====== Loop through the c-to-f steps C CFSAV = CF CF = CF - STEP DO 30 I = 1 , NSTEP CF = CF + STEP CALL LAUE(ION) C C====== Find out how many reflections recorded C ICOUNT = 0 ITOT = 0 DO 20 J = 1 , NR K = ION(J) IF (K.NE.0) THEN ITOT = ITOT + K ICOUNT = ICOUNT + 1 END IF 20 CONTINUE FRACT = FLOAT(ICOUNT)/FLOAT(NR) WRITE (IPRINT,6040) CF , FRACT , ICOUNT , ITOT WRITE (IUNIT6,6040) CF , FRACT , ICOUNT , ITOT 30 CONTINUE CF = CFSAV RETURN C C====== Format statements C 6010 FORMAT (' This command analyses the effect of changing the', + ' crystal to film distance.') 6020 FORMAT (/// +' Effect of increasing c-to-f on the fraction of symmetry unique' +/' reflections recorded') 6030 FORMAT (///' cf(mm) Fraction of symmetry Total reflections'/ + ' uniques recorded on film'/) 6040 FORMAT (F6.1,F13.3,' (',I5,')',12X,I5) END **==wavvar.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** WAVVAR ** C*************************** C C SUBROUTINE WAVVAR() C INCLUDE 'implicit_none.finc' C INTEGER MAXONE,MAXALL PARAMETER (MAXONE=20000,MAXALL=100000) C COMMON /DATA/CELL,WVM,WVX,WVS,CF,PHIS,SPNDLE,CCELL,NUMA,NUMC,DMIN COMMON /SPOTS/NR REAL CF,DMIN,SPNDLE,WVM,WVS,WVX INTEGER NR,NUMA,NUMC REAL CCELL(6),CELL(6),PHIS(3) C INCLUDE 'chanel.finc' C C====== LOCALS: C INTEGER ION(MAXALL) REAL WION(MAXALL) REAL DEL,WMAX INTEGER ICOUNT,ITOT,J,K INTEGER IBOX(10) C C------------------------------------------------------------------------------- C WRITE (IPRINT,6010) WRITE (IPRINT,6020) WRITE (IUNIT6,6020) C CALL LAUEW(ION,WION) C C====== Find out how many reflections recorded C ICOUNT = 0 WMAX = 0. ITOT = 0 DO 20 J = 1 , NR K = ION(J) IF (K.GT.1) THEN ITOT = ITOT + K IF (WION(J).GT.WMAX) WMAX = WION(J) ICOUNT = ICOUNT + 1 END IF 20 CONTINUE C WRITE (IPRINT,6030) ICOUNT , ITOT WRITE (IUNIT6,6030) ICOUNT , ITOT DEL = WMAX/10. IF (DEL.LT..01) DEL = .01 DO 30 K = 1 , 10 IBOX(K) = 0 30 CONTINUE DO 40 J = 1 , NR IF (ION(J).GT.1) THEN K = WION(J)/DEL + 1 IF (K.LE.0) K = 1 IF (K.GT.10) K = 10 IBOX(K) = IBOX(K) + 1 END IF 40 CONTINUE DO 50 K = 1 , 10 WRITE (IPRINT,6040) IBOX(K) , DEL*(K-1) , DEL*K WRITE (IUNIT6,6040) IBOX(K) , DEL*(K-1) , DEL*K 50 CONTINUE RETURN C C====== Format statements C 6010 FORMAT (' This command calculates the wavelength distribution', + ' for'/' the symmetry equivalents in this orientation.') 6020 FORMAT (///' Distribution of wavelength variations'/) 6030 FORMAT (' Symmetry uniques with duplicates recorded =', + I6/' Total symmetry duplicates recorded =', + I6//' Number of reflections Wavelength variation range'/ + ) 6040 FORMAT (I10,7X,F19.2,' - ',F4.2) END