C
C  PROGRAM PREF   TO DETERMINE THE DIRECTIONS OF PREFERRED
C  ORIENTATION IN A POWDER PATTERN
C
C  PROGRAMMED BY ROBERT L. SNYDER
C                 NEW YORK STATE COLLEGE OF CERAMICS
C                  ALFRED UNIVERSITY
C                   ALFRED N.Y. 14802  USA
C  PLEASE CONTACT ME WITH ANY COMMENTS, PROBLEMS OR SUGGESTIONS
C
C  This is based on a method described in a paper by Snyder and Carr
Cin a book entitled "Surfaces and Interfaces in Ceramics" Published
C  by Plenum press New York 1972.
C
C  INPUT SET UP  (There is an interactive SUPREF program to do this)
C
C
C  ITEM 1 COL 1-80  ANY TITLE INFORMATION
C       NOTE: THE FIRST 10 COLS WILL BE PLOTTED ON ANY PLOTS
C       THE COMPOUND FORMULA OR NAME SHOULD BE CENTERED
C       AROUND COL 5. 17-80 WILL ONLY APPEAR ON PRINTOUT
C  ITEM 2 CELL PARAMETERS A,B,C,ALPHA,BETA,GAMMA (6F10.4)
C  ITEM 3 CONTROL PARAMETER CARD
C     COL 1 ICON(1) = 0  CALCULATED PATTERN TO BE INPUT
C                   = 1  NO CALC PATTERN TO BE INPUT
C     COL 2 ICON(2)       NOT PRESENTLY USED
C     COL 3 ICON(3)       NOT PRESENTLY USED
C     COL 4 ICON(4) = 0  FIND Q FACTORS FOR ALL GENERATED ZONES AND
C                        PLOT AGAINST ANGLE FROM EACH PRINCIPAL
C                        ZONE
C                   = 1  SKIP PLOT
C                   = 2  SKIP THIS ROUTINE
C     COL 5 ICON(5) = 0  LAST DATA SET
C                   = 1  ANOTHER DATA SET TO FOLLOW ITEM 8
C  COL 6-8 ONLY APPLY IF COL 4 = 0  IF SO ANY OR ALL MAY BE 1
C  7 OF EACH KIND OF PLOT ARE PRODUCED ONE FOR EACH PRINCIPAL
C  ZONE
C
C       COL 6 = 1 PLOT Q-OBS VS ANGLE
C             = 0 NO Q-OBS PLOTS
C       COL 7 = 1 PLOT Q-CALC VS ANGLE
C             = 0 NO Q-CALC PLOTS
C       COL 8 = 1 PLOT DELTA Q VS ANGLE
C             = 0 NO DELTA Q PLOTS
C       COL 9 = 0 NORMAL OPERATION
C             = 1 WRITE QO, QC AND DELTA Q'S ON UNIT 20
C  ITEM 4   OBSERVED PATTERN SETTINGS
C     COL 1-5  MAXPKN =  MAXIMUM NUMBER OF UNMATCHED
C                        OBSERVED PEAKS TOLERABLE FOR MEANINGFUL
C                        EXAMINATION BLANK DEFAULTS TO 200.
C     COL  6-15 SINT  =  FULL SCALE SETTING OF RECORDER EXPRESSED
C                        IN COUNTS PER SECOND
C                        BLANK DEFAULTS TO 1000
C     COL 16-25 BCKGRD = AVERAGE BACKGROUND VALUE FROM STRIP CHART
C                        BLANK DEFAULTS TO 1.0
C     COL 26-35 FACTR- A MULTIPLYER FOR THE SIZE OF ANY PLOTS
C               IF 1.0 PLOTS WILL BE 10X10 INCHES
C               IF 0.5 PLOTS WILL BE 5X5 INCHES
C     COL 36-45 X2 = MAXIMUM ANGLE OUTPUT FOR TABLE LISTING OF Q VALUES
C                    BLANK DEFAULTS TO 180. (I.E. ALL DATA)
C     COL 46-55 X3 = MAXIMUM ANGLE OUTPUT ON PLOTS OF Q VALUES
C                    BLANK DEFAULTS TO 90.
C     COL 56-65 WAVE = WAVELENGTH OF RADIATION USED
C                BLANK DEFAULTS TO CU KA1  1.54051
C  ITEM 5  OBSERVED PATTERN
C     COL 1-10 2-THETA
C     COL 11-20 UNSCALED INTENSITY
C  ITEM 6 BLANK CARD ENDS OBSERVED DATA SET
C  ITEM 7  CALCULATED PATTERN, IF ENTERED, USES FORMAT TO READ
C          PUNCHED CARDS WHICH ARE OUTPUT FROM POWDER CALC PROGRAM
C          OF DEANE SMITH  I.E.
C     COL  1-17 2-THETA
C     COL 26-32 D-CALCULATED
C     COL 34-36 H    RIGHT ADJUSTED
C     COL 37-39 K    RIGHT ADJUSTED
C     COL 40-42 L    RIGHT ADJUSTED
C     COL 51-54 I-CALCULATED  RIGHT ADJUSTED
C  ITEM 8  BLANK CARD ENDS CALC DATA SET
C
C
C         INSTRUCTIONS FOR IMPLEMENTING ON OTHER COMPUTERS
C
C  THIS PROGRAM IS WRITTEN IN ANSI STANDARD FORTRAN THE
C  ONLY MACHINE DEPENDANT FEATURES ARE THE FOLLOWING TWO
C
C  1. THE INPUT UNIT (IR) IS CURRENTLY SET ONCE ONLY AS
C      THE FIRST EXECUTEABLE STATEMENT TO 105
C     THE OUTPUT UNIT (IW) IS CURRENTLY SET ONCE ONLY AS
C      THE SECOND EXECUTABLE STATMENT TO 108
C  THESE TWO VARIABLES MUST BE SET TO CONFORM WITH LOCAL
C   DEFINITIONS
C
C  2. PLOTTER ROUTINES: ALL ROUTINES USED BY SUBROUTINE
C     GPLOT ARE STANDARD CALCOMP SUBROUTINES.  BUT THE MAIN
C     PROGRAM AT THE PLACES MARKED WITH ****** USES THREE
C     LOCAL ROUTINES WHICH USUALLY VARY BETWEEN COMPUTERS
C  START(X,Y) IS A PLOTTER INITIALIZATION ROUTINE THAT
C              STARTS THE PLOTTER AT THE COORDINATES X AND Y
C              REPLACE WITH PLOTS OR ANY SIMILAR ROUTINE
C  END(X) IS A PLOTTER TERMINATION ROUTINE WHICH ROLLS OFF
C      X INCHES OF PAPER.  ANY SIMILAR TERMINATION ROUTINE
C         MAY BE SUBSTITUTED
C  FACTOR(X) IS A ROUTINE WHICH ALLOWS ALL PLOTS TO BE SCALLED
C          UP OR DOWN BY THE FACTOR (X) SPECIFIED AS THE
C          ARGUMENT.  THIS MAY BE SIMPLY DUMMIED OUT OR
C          THE CALL MAY BE REMOVED
C
C
C  OVERLAY SCHEME
C
C  IF DESIRED THE SUBROUTINES MAY BE OVERLAYED AS FOLLOWS
C
C     ******            ********
C     *    *  <------>  *BRAGG *
C     *    *            ********
C     *    *
C     *    *
C     *    *
C     *    *            ********
C     *    *  <------>  * SCAL *
C     *    *            ********
C     *    *
C     *    *
C     *    *            **********************
C     *    *  <------>  * START, END, FACTOR *
C     *    *            **********************
C     *    *
C     *    *
C     *    *
C     *    *            *******          *******
C     *    *            *     *  <---->  * MAX *
C     *MAIN*            *     *          *     *
C     *    *            *     *
C     *    *            *     *
C     *    *            *     *          ********
C     *    *            *     *  <---->  * SORT *
C     *    *            *     *          ********
C     *    *            *     *
C     *    *  <------>  *ZONE *
C     *    *            *     *
C     *    *            *     *          ********          *******
C     *    *            *     *          *      *  <---->  * LLS *
C     *    *            *     *          *      *          *******
C     *    *            *     *          *      *
C     *    *            *     *  <---->  *QPLOT *
C     *    *            *     *          *      *          **********
C     *    *            *     *          *      *          *CALCOMP *
C     *    *            *     *          *      *  <---->  *ROUTINES*
C     *    *            *     *          *      *          **********
C     ******            *******          ********
C
C
C
C  GLOSSARY OF SYMBOLS
C
C
C  OBSERVED VALUES OV(N,1) = 2-THETA
C                  OV(N,2) = DO
C                  OV(N,3) = IO
C                  OV(N,4) = 0 IF UNMATCHED IN CALC LIST(CV)
C                          = M WHERE M = # OF MATCHED PEAK IN CV LIST
C
C  CALCULATED VALUES HAVE SAME MEANING IN CV ARRAY EXCEPT:
C    CV(N,4) = 0 IF UNMATCHED
C            = 1 IF MATCHED
C
C   IO = TOTAL NUMBER OF OBSERVED PEAKS READ
C   IB = TOTAL NUMBER OF MATCHES
C
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      COMMON /BLOK/ DELTA(300)
      DIMENSION IOVCV(3438),PCT(2,7)
      EQUIVALENCE (OV(1),IOVCV(1))
C
C  FORTRAN SOURCE PROGRAM - SNYDER - 8/4/73.
C
C  DEFINE I/O UNITS  IR= READER   IW= PRINTER
      IR=5
      IW=6
      AIW=0.02
      LV=200
      LV1=3438
C
C  INITIALIZE ARRAYS
C
5     DO 10 I=1,LV1
10    IOVCV(I)=0
C
      DTR=3.141592/180.
      RTD=1.0/DTR
C
C     INPUT DATA
C
      READ(IR,200) ITITLE
      READ(IR,201) A,B,C,AL,BE,GA
      READ(IR,202)ICON
      READ(IR,204)MAXPKN,SINT,BCKGRD,FACTR,X2,X3,WAVE
      IF(FACTR.EQ.0.0.OR.FACTR.GT.1.0) FACTR = 1.0
      IF(X2.EQ.0.0) X2 = 180.
      IF(X3.EQ.0.0) X3 = 90.
      IF(MAXPKN.EQ.0) MAXPKN = 200
      IF(SINT.EQ.0.0) SINT = 1000.
      IF(BCKGRD.EQ.0.0) BCKGRD = 1.0
      IF(WAVE.EQ.0.0) WAVE = 1.54051
      WRITE(IW,300)ITITLE,A,B,C,AL,BE,GA,ICON,X2,X3,
     1 MAXPKN,SINT,BCKGRD,FACTR,WAVE
  300 FORMAT(1X,20A4//16H CELL CONSTANTS ,6F10.4/
     130H CONTROL INTEGERS (ICON(1-10) ,10I2/
     456H MAX ANGLE OUTPUT FOR Q-FACTOR LISTING (ICON(4).NE.2) IS,F10.3,
     5 28H  MAX ANGLE USED IN PLOTS IS,F10.3//
     642H  MAX NUMBER OF UNMATCHED PEAKS ALLOWED = ,I4//
     735H  FULL SCALE SETTING OF RECORDER IS,F10.2,10H CPS  WITH,
     816H A BACKGROUND OF,F6.2,12H CHART UNITS/
     920H PLOT SIZE FACTOR = ,F10.2/
     A14H WAVELENGTH = ,F10.5)
C
C  READ OBSERVED DATA IN SUBROUTINE BRAGG
C
      CALL BRAGG(IO,IR,LV,WAVE)
C
C     IF OPTION IS NO CALC INPUT GO TO Q FACTOR ROUTINE ONLY
C
      IF(ICON(1).EQ.1)  GO TO 33
C
C     READ CALCULATED PATTERN
C
      DO 15 IC=1,LV
      READ(IR,205)(CV(IC,J),J=1,2),(IHKL(IC,J),J=1,3),IOY
      CV(IC,3)=IOY
      IF(CV(IC,1).EQ.0.0) GO TO 17
15    CONTINUE
17    IF(CV(LV,1).EQ.0.) IC=IC-1
C
C  MATCH OBSERVED TO CALCULATED PEAKS ACCORDING TO THE SMALLEST
C  DIFFERENCE IN D-SPACING - SPACING MUST BE LESS THAN THE
C  ERROR WINDOW = AIW =.02
C
      DO 30 M=1,IO
      AMIN=100.
      DO 20 N=1,IC
      IF(ABS(OV(M,2)-CV(N,2)).GE.AMIN) GO TO 20
      AMIN=ABS(OV(M,2)-CV(N,2))
      J=N
20    CONTINUE
      IF(AMIN.GT.AIW) GO TO 30
      AJ=J
C
C  TEST OV LIST FOR PREVIOUS CV(J,2) MATCH AND DELETE
C
      DO 25 I1=1,IO
25    IF(OV(I1,4).EQ.AJ) OV(I1,4)=0.
      OV(M,4)=J
30    CONTINUE
C
C  SCALE OV AND CV INTENSITY PEAKS TO THEIR HIGHEST PEAKS
C
  33  CALL SCAL(IO,IC,SCALEO,SCALEC)
      IF(ICON(1).EQ.1) GO TO 48
C
C  SCALEO = HIGHEST OBSERVED PEAK
C  SCALEC = HIGHEST CALCULATED PEAK
C
C  SCALE BACKGROUND AND FIND BACKGROUND IN COUNTS/SECOND
C
      BCKSCO = BCKGRD*100./SCALEO
      BCKCPS=BCKGRD*SINT/100.
      BCKCPS=2.5*SQRT(BCKCPS)+BCKCPS
C
C  PRINT MATCHED AND UNMATCHED PEAK LISTINGS
C
      WRITE(IW,206) ITITLE
      MATCH=0
      MAXU = 0
C
C  MATCH=NUMBER OF MATCHES
C  MAXU=NUMBER OF OBSERVED PEAKS UNMATCHED
C
      DO 38 M=1,IO
      N=OV(M,4)
      IF(N.EQ.0) GO TO 38
      MATCH=MATCH+1
      CV(N,4)=M
C
      DD=OV(M,2)-CV(N,2)
      DI = OV(M,3) - CV(N,3)
      WRITE(IW,207)M,N,(OV(M,J),CV(N,J),J=1,2),DD,OV(M,3),
     1CV(N,3),DI,(IHKL(N,J),J=1,3)
38    CONTINUE
C
C  UNMATCHED LISTS
C
C  UNMATCHED OBSERVED LIST
C
      WRITE(IW,208)
C
C
      DO 40 M=1,IO
      IF(OV(M,4).GT.0.0) GO TO40
      MAXU=MAXU+1
      WRITE(IW,209) M,(OV(M,J),J=1,3)
40    CONTINUE
      IF(MAXU.EQ.0) WRITE(IW,210)
C
C  UNMATCHED CALC LIST
C
      WRITE(IW,211)
C
      DO 45 N=1,IC
      IF(CV(N,4).GT.0.0) GO TO 45
      WRITE(IW,212) N,(CV(N,J),J=1,3),(IHKL(N,J),J=1,3)
45    CONTINUE
C
C  LIST OF OBSERVED INTENSITIES LESS THAN 2.5SIGMA OF BACKGROUND
C
      WRITE(IW,213)
      KK=0
      DO 47 M=1,IO
      CPSIO=((OV(M,3)*SCALEO/100.)+BCKGRD)*SINT/100.
      IF(CPSIO.GT.BCKCPS) GO TO 47
      KK=1
      WRITE(IW,209)M,(OV(M,J),J=1,3)
47    CONTINUE
C
      IF(KK.EQ.0) WRITE(IW,210)
C
      IF(MAXU.LE.MAXPKN) GO TO 48
      WRITE(IW,215)
      GO TO 100
C
C  FIND SUM IO > 0. ; SUM IOB > BCKGRD ; SUM IC
C
48    SUMIO=0.
      SUMIOB = 0.
      SUMIC = 0.
      DO 50 M=1,IO
      IF(OV(M,4).EQ.0.) GO TO50
      SUMIO=SUMIO+OV(M,3)
      CPSIO=((OV(M,3)*SCALEO/100.)+BCKGRD)*SINT/100.
      IF(CPSIO.GT.BCKCPS) SUMIOB=SUMIOB+OV(M,3)
50    CONTINUE
C
C
      DO 55 N=1,IC
55    SUMIC=SUMIC+CV(N,3)
C
      WRITE(IW,216)SCALEO,SCALEC,BCKGRD,BCKSCO,BCKCPS,
     1SUMIO,SUMIOB,SUMIC
C
C
C  OPTIONS
C
C
   60 IF(ICON(3).NE.2.OR.ICON(4).NE.0)GO TO 65
C
C**********************************************************
C  START IS A LOCAL PLOTTER INITIATION ROUTINE AND SHOULD
C  BE REPLACED BY YOUR  INSTALLATIONS ROUTINE
C*********************************************************
C  FACTOR IS A LOCAL ROUTINE WHICH REDUCES THE SIZE OF THE
C  PLOT BY ITS ARGUMENT--- IT MAY SAFELY BE REMOVED
C*********************************************************
C
      CALL START(2.0,5.5)
      CALL FACTOR(FACTR)
   65 IF(ICON(1).EQ.1)  GO TO 75
C
C  Q FACTORS
C
  75  IF(ICON(4).EQ.2)  GO TO 80
      CALL ZONE(IW,IO,X2,X3,SUMIO,SUMIC,IC,PCT)
C
C************************************************************
C  END IS A LOCAL PLOTTER TERMINATION ROUTINE AND SHOULD
C  BE REPLACED BY YOUR VERSION
C*************************************************************
C
80    IF(ICON(4).EQ.0)                   CALL END(12.)
      IF(ICON(1).EQ.1)  GO TO 100
C
C  EVALUATION OF PATTERN
C
      WRITE(IW,217)
C
      IP=0
      DO 81 M=1,IO
      N=OV(M,4)
      IF(N.EQ.0) GO TO 81
      CPSIO=((OV(M,3)*SCALEO/100.)+BCKGRD)*SINT/100.
      CPSIC=((CV(N,3)*SCALEC/100.)+BCKGRD)*SINT/100.
      IF(ABS(CPSIC-CPSIO).LE.(2.5*SQRT(CPSIO))) GO TO 81
      IP=1
81    CONTINUE
C
C  WRITE PATTERN IS RANDOM IF IP=0.0
C
      IF(IP.NE.0) GO TO 815
      WRITE(IW,218)
      GO TO 100
C
  815 DO 82 M=1,IO
      IF(OV(M,3).EQ.100.) GO TO 83
   82 CONTINUE
83    DO 84 N=1,IC
      IF(CV(N,3).EQ.100.) GO TO 85
   84 CONTINUE
85    IF(N.NE.M) GO TO 100
C
C
      IF(SUMIC.GT.SUMIO) WRITE(IW,219)
      IF(SUMIO.GT.SUMIC) WRITE(IW,220)
      IF(SUMIC.GT.SUMIOB) WRITE(IW,221)
      IF(SUMIOB.GT.SUMIC) WRITE(IW,222)
C
      WRITE(IW,223)(IHKL(N,J),J=1,3)
C
C  CALCULATE PHYSICAL PROPERTY ENHANCEMENT COEFFICIENTS
C
  100 IP = 0
      IT = 0
      WRITE(IW,225)
  225 FORMAT(1H0,42HA VECTOR LIKE PHYSICAL PROPERTY DEPENDANT
     1,41HON ANY OF THE FOLLOWING 7 PRINCIPAL ZONES/
     257H WILL BE INCREASED OR DECREASED BY THE FOLLOWING PERCENTS
     3    //16H  ZONE   PERCENT)
      DO120 J=IP,1
      DO120 K=IP,1
      DO120 L=IP,1
      IF(J+K+L.EQ.0)  GO TO 120
      IT = IT + 1
      IF(PCT(1,IT).LT.0.0) GO TO 110
      T = PCT(1,IT)*100./(1.0-PCT(2,IT))
      GO TO 115
  110 T = PCT(1,IT)*100./PCT(2,IT)
  115 WRITE(IW,226)J,K,L,T
  226 FORMAT(1X,3I2,F8.2)
  120 CONTINUE
      IF(ICON(5).EQ.1) GO TO 5
      IF(IDATA.EQ.1) GO TO 5
C
C  FORMATS
200   FORMAT(20A4)
201   FORMAT(6F10.4)
202   FORMAT(10I1)
204   FORMAT(I5,6F10.4)
  205 FORMAT(F17.2,8X,F7.4,1X,3I3,8X,I4)
206   FORMAT(1H1,20A4/23H OBSERVED PEAKS MATCHED,
     120H TO CALCULATED PEAKS//
     226H NO   NC  2THETAO  2THETAC,5X,2HDO,7X,2HDC,
     34X,6HDO-DC ,5X,2HIO,7X,2HIC,4X,6HIO-IC ,9X,1HH,
     42X,1HK,2X,1HL/)
207   FORMAT(I3,2X,I3,8(2X,F7.3),5X,3I3)
208   FORMAT(1H1,24H UNMATCHED OBSERVED LIST//
     11X,2HNO,4X,6H2THETA,8X,2HDO,10X,2HIO/)
209   FORMAT(I3,3(2X,F10.4))
210   FORMAT(/8H**NONE**/)
211   FORMAT(1H1 , 31H UNMATCHED CALCULATED PEAK LIST /
     12X,2HNC,4X,6H2THETA,10X,2HDC,10X,2HIC,6X,1HH,2X,1HK,2X,1HL/)
212   FORMAT(1X,I3,3(2X,F10.4),2X,3I3)
213   FORMAT(1H1 , 48H OBSERVED PEAKS LESS THAN 2.5SIGMA OF BACKGROUND /
     1/1X,2HNO,6X,6H2THETA,8X,2HDO,10X,2HIO/)
215   FORMAT(1H1 , 38H NUMBER OF UNMATCHED OBSERVED PEAKS IS /
     152H TOO LARGE FOR TOLERANCE....PROCESSING WILL CONTINUE/
     217H TO NEXT DATA SET//)
216   FORMAT(1H1,   34H THE FOLLOWING VALUES ARE USED IN /
     140H CALCULATING SCALE FACTORS AND Q-FACTORS/
     2    / 24H HIGHEST OBSERVED PEAK =, F10.4 /
     326H HIGHEST CALCULATED PEAK =,F10.4/
     428H BACKGROUND VALUE REPORTED =,F10.4/
     536H BACKGROUND VALUE SCALED TO SCALEO =,F10.4/
     651H BACKGROUND + 2.5 STANDARD DEVIATION IN COUNTS PER ,
     79HSECONDS =,F10.4/
     854H TOTAL SUM OF INTENSITY OBSERVED (MATCHED PEAKS ONLY) ,
     917HWHICH IS .GE. 0 =,F10.4/
     C54H TOTAL SUM OF INTENSITY OBSERVED (MATHCED PEAKS ONLY) ,
     C26HWHICH IS .GE. BACKGROUND =,F10.4/
     C45H TOTAL SUM OF INTENSITY OF CALCULATED PEAKS =,F10.4/)
C
217   FORMAT(1H1 ,  30HEVALUATION OF OBSERVED PATTERN//)
218   FORMAT( 27H0OBSERVED PATTERN IS RANDOM/)
219   FORMAT(/ 41H INTENSITY .GE. 0    DIRECTION IS NONPREF,
     117HERRED ORIENTATION//)
220   FORMAT(/ 37H INTENSITY .GE. 0   DIRECTION IS PREF,
     118HERRED  ORIENTATION//)
221   FORMAT(/ 43H INTENSITY .GE. BCKGRD DIRECTION IS NONPREF,
     117HERRED ORIENTATION//)
222   FORMAT(/ 44HINTENSITY .GE. BCKGRD DIRECTION IS PREFERRED,
     112H ORIENTATION/)
223   FORMAT(/9HZONE AXIS,3I3)
C
      STOP
      END
      SUBROUTINE SORT(K1,LH)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
C
C  SORT ROUTINE     LH=1 SORT GEN(K,1) FROM LOW TO HIGH
C                   LH=2 SORT GEN(K,1) FROM HIGH TO LOW
C  ARRAY HKL IS SORTED ALONG WITH GEN FOR USE IN SUBROUTINE
C  ZON
C
1     K=1
2     IF(LH.EQ.2) GO TO 3
      IF(GEN(K+1,1).GE.GEN(K,1)) GO TO 5
      GO TO 4
3     IF(GEN(K,1).GE.GEN(K+1,1)) GO TO 5
    4 DO 6 I = 1,3
      AA = GEN(K,I)
      GEN(K,I) = GEN(K+1,I)
      GEN(K+1,I) = AA
    6 CONTINUE
      AA = HKL(K)
      HKL(K) = HKL(K+1)
      HKL(K+1) =AA
      GO TO 1
5     K=K+1
      IF(K.LT.K1) GO TO 2
      RETURN
      END
      SUBROUTINE DSP(I1,I2,I3,J1,J2,J3,ANG)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL1,BE1,GA1,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      AL=AL1*DTR
      BE=BE1*DTR
      GA=GA1*DTR
C
C  FLOAT 2 SETS OF HKL'S
C
      ZH = I1
      ZK = I2
      ZL = I3
      CH = J1
      CK = J2
      CL = J3
C
C     CALCULATION OF ANGLE BETWEEN TWO HKL'S.
C
      X=1.-(COS(AL))**2-(COS(BE))**2-(COS(GA))**2
      X=X+2.*(COS(AL))*(COS(BE))*(COS(GA))
      V=A*B*C*SQRT(ABS(X))
      S11=(B*C*SIN(AL))**2
      S22=(A*C*SIN(BE))**2
      S33=(A*B*SIN(GA))**2
      S12=A*B*C**2
      S12=S12*((COS(AL)*COS(BE))-(COS(GA)))
      S23=B*C*A**2
      S23=S23*(COS(BE)*COS(GA)-COS(AL))
      S13=A*C*B**2
      S13=S13*(COS(GA)*COS(AL)-COS(BE))
      ZD=(1.0/V**2)*(S11*ZH**2+S22*ZK**2+S33*ZL**2+2.*S12*
     1ZH*ZK+2.*S23*ZK*ZL+2.*S13*ZH*ZL)
      ZD=SQRT(ABS(1.0/ZD))
      CD=(1.0/V**2)*(S11*CH**2+S22*CK**2+S33*CL**2+2.*S12*
     1CH*CK+2.*S23*CK*CL+2.*S13*CH*CL)
      CD=SQRT(ABS(1.0/CD))
      COSA=(ZD*CD/V**2)*(S11*ZH*CH+S22*ZK*CK+S33*ZL*CL+
     1S23*(ZK*CL+CK*ZL)+S13*(ZL*CH+CL*ZH)+S12*(ZH*CK+
     2CH*ZK))
      IF(COSA.GT.1.0) COSA = 1.0
      ANG = ACOS(COSA)*RTD
      RETURN
      END
      SUBROUTINE QPLOT(K,N,X,I1,I2,I3,IW)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      COMMON /BLOK/ DELTA(300)
      DIMENSION XC(2),YC(2)
      A1=I1
      A2 = I2
      A3 = I3
      AMAX=0.
      BMAX = 0.0
      CMAX = 0.0
      LOOP = 0
      IP = 0
      IPL = ICON(6) + ICON(7) + ICON(8)
C
C  FIND MAX.Y-VALUE
C
      DO 5 J=1,N
      IF(GEN(J,1).GT.X) GO TO 6
      IF(ABS(DELTA(J)).GT.CMAX) CMAX = ABS(DELTA(J))
      IF(ABS(GEN(J,3)).GT.BMAX) BMAX = ABS(GEN(J,3))
5     IF(ABS(GEN(J,2)).GT.AMAX)AMAX=ABS(GEN(J,2))
C
C  SCALE X AND Y COORDINATES
C
    6 E=X/10.
      D=AMAX/5.
      F = BMAX / 5.0
      G = CMAX / 5.0
      AM = -1.0 * AMAX
      DO 10 J=1,N
      GEN(J,1)=GEN(J,1)/E
      DELTA(J) = DELTA(J) / G
      GEN(J,3) = GEN(J,3)/F
10    GEN(J,2)=GEN(J,2)/ D
C
C
C  DRAW HEADING, LABELS, AND AXES
C
   12 CALL NUMBER(1.0,-3.08,.77,A1,90.,-1)
      CALL NUMBER(999.,999.,.77,A2,90.,-1)
      CALL NUMBER(999.,999.,.77,A3,90.,-1)
      CALL SYMBOL(999.,999.,.77,5H ZONE,90.,+5)
      CALL PLOT(2.5,0.0,-3)
      CALL SYMBOL(.75,-5.3,.77,ITITLE,0.0,10)
   15 IF(ICON(6).EQ.1.AND.LOOP.EQ.0) GO TO 16
      IF(ICON(7).EQ.1.AND.LOOP.NE.2) GO TO 17
      GO TO 18
   16 CALL AXIS(0.0,-5.0,1H ,1,10.,90.,AM,D)
      CALL SYMBOL(-0.4,-3.85,.77,10HQ OBSERVED,90.,10)
      LOOP = 1
      LO = 2
      WRITE(IW,160)
  160 FORMAT(//42H LEAST SQUARES RESULTS FOR OBSERVED Q PLOT)
      GO TO 20
   17 AM =  -1.0 * BMAX
      CALL AXIS(0.0,-5.0,1H ,1,10.,90.,AM,F)
      CALL SYMBOL(-0.4,-4.62,.77,12HQ CALCULATED,90.,12)
      LOOP = 2
      LO = 3
      WRITE(IW,170)
  170 FORMAT(//44H LEAST SQUARES RESULTS FOR CALCULATED Q PLOT)
      GO TO 20
   18 LOOP = 3
      LO = 2
      AM = -1.0 * CMAX
      WRITE(IW,180)
  180 FORMAT(//39H LEAST SQUARES RESULTS FOR DELTA Q PLOT)
      CALL SYMBOL(-0.4,-5.005,.77,13HQ OBS - Q CAL,90.,13)
      CALL AXIS(0.0,-5.0,1H ,1,10.,90.,AM,G)
      DO 185 J = 1,N
      T = GEN(J,2)
      GEN(J,2) = DELTA(J)
      DELTA(J) = T
  185 CONTINUE
C
   20 CALL AXIS(0.0,0.0,1H ,1,10.,0.0,0.0,E)
      CALL SYMBOL(0.1,-.80,.63,15HANGLE FROM ZONE,0.0,15)
      N1=0
      DO 25 J=1,N
      IF(GEN(J,1).GT.10.) GO TO 25
      N1=N1+1
      CALL SYMBOL(GEN(J,1),GEN(J,LO),.07,1,0.0,-1)
25    CONTINUE
C
      CALL LLS(IW,N1,XC,YC,SLP,YINT,LO)
C
C  THE FOLLOWING STATMENTS WILL PUT THE LEAST SQUARES
C  EQUATION ON THE PLOTS--WHOLE LEAST SQUARES NOT VERY
C  USEFUL SO I HAVE COMMENTED THIS OUT  R. SNYDER
C
CC
C     CALL PLOT(XC(1),YC(1),3)
C     CALL PLOT(XC(2),YC(2),2)
CC
CC
C     CALL SYMBOL(-2.5,-5.0,.14,'SLOPE=',0.0,6)
C     CALL NUMBER(999.,999.,.14,SLP,0.0,3)
CC
C     CALL SYMBOL(-2.5,-5.45,.14,'Y-CEPT =',0.0,8)
C     CALL NUMBER(999.,999.,.14,YINT,0.0,3)
C
      CALL PLOT(12.0,0.0,-3)
      IP = IP + 1
      GO TO(26,26,27),LOOP
   26 IF(IP.GT.IPL) GO TO 29
      GO TO 12
  27  DO 28 J = 1,N
   28 GEN(J,2) = DELTA(J)
   29 DO 30 J=1,N
      GEN(J,1)=GEN(J,1)*E
      GEN(J,2)=GEN(J,2)*D
      GEN(J,3) = GEN(J,3)*F
30    CONTINUE
C
      RETURN
      END
      SUBROUTINE LLS(IW,NUM,XC,YC,SLOPE,YINCP,LO)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      DIMENSION XC(2), YC(2)
C
C  GENERAL LINEAR REGRESSION CIRCA 1965  BY R. SNYDER
C
      AN=NUM
      XB=0.0
      YB=0.0
      SUMX=0.0
      SUMY=0.0
      SUMXY=0.0
      SUMYY=0.0
      SUMXX=0.0
C     XB MEANS X-BAR
C     YB MEANS Y-BAR
      DO 10 J = 1,NUM
      SUMXX=SUMXX+GEN(J,1)*GEN(J,1)
      SUMYY=SUMYY+GEN(J,LO)*GEN(J,LO)
      SUMX=SUMX+GEN(J,1)
      SUMY=SUMY+GEN(J,LO)
   10 SUMXY=SUMXY+GEN(J,1)*GEN(J,LO)
      XB = SUMX / AN
      YB = SUMY / AN
      SLOPE=(AN*SUMXY-SUMX*SUMY)/(AN*SUMXX-SUMX*SUMX)
C     YINCP MEANS Y-INTERCEPT
      YINCP = YB - SLOPE * XB
C     SQRS MEANS S SQUARED OR THE VARIANCE OF A SINGLE MEASURMENT
      SQRS=SUMYY-SUMY*SUMY/AN
      SQRS=SQRS-((SUMXY-SUMX*SUMY/AN)**2)/(SUMXX-SUMX*SUMX/AN)
      SQRS=SQRS/(AN-2.)
      S=SQRT(SQRS)
C     VARS MEANS THE VARIANCE FOR THE SLOPE
      VARS=SQRS/(SUMXX-AN*XB*XB)
      SS=SQRT(VARS)
C     SINCP MEANS THE STANDARD DEVIATIONOF THE Y-INTERCEPT
      SINCP=SQRT(SQRS*SUMXX/(AN*SUMXX-SUMX*SUMX))
      T=YINCP/SINCP
      NF = NUM - 2
      WRITE(IW,100)  YINCP,SLOPE
  100 FORMAT(//31H THE LEAST SQUARES EQUATION IS  /6H  YC =,G14.8,3H + ,
     1G14.8,2H X//)
      WRITE(IW,110)  SINCP,SS
  110 FORMAT(
     147H THE STANDARD DEVIATION OF THE Y INTERCEPT IS =,G14.8/
     241H THE STANDARD DEVIATION OF THE SLOPE IS =,G14.8)
      TOTAL = 0.0
      R = 0.0
      XC(1) = .1E30
      XC(2) = -.1E30
      DO 20 J=1,NUM
      YCAL = YINCP + SLOPE * GEN(J,1)
      IF ( GEN(J,1) .GE. XC(1)) GO TO 12
      XC(1) = GEN(J,1)
      YC(1) = YCAL
   12 IF ( GEN(J,1) .LE. XC(2)) GO TO 14
      XC(2) = GEN(J,1)
      YC(2) = YCAL
   14 DIFF = YCAL - GEN(J,2)
      R = R + DIFF*DIFF
      TOTAL = TOTAL + ABS(DIFF)
   20 CONTINUE
   40 TOTAL = TOTAL / AN
      WRITE(IW,140)  TOTAL,R
  140 FORMAT(//48H THE AVERAGE OF THE ABSOLUTE VALUES OF YC-Y IS =,G14.8
     1 /44H THE SUM OF THE SQUARES OF THE RESIDUALS R =,G14.8)
      RETURN
      END
      SUBROUTINE BRAGG(IO,IR,LV,WAVE)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      RAD=DTR/2.0
C
C  INPUT 2THETA AND INTENSITY
C
      DO 3 IO=1,LV
      READ(IR,2)(OV(IO,J),J=1,3,2)
2     FORMAT(2F10.4)
      IF(OV(IO,1).EQ.0.0) GO TO 4
C
C  BRAGG'S EQUATION FOR FINDING D-SPACING GIVEN RADIATION AND
C  2THETA OF X-RAY BEAM
C
      OV(IO,2)=WAVE/(2.0*SIN(OV(IO,1)*RAD))
3     CONTINUE
4     IF(OV(LV,1).EQ.0.0) IO=IO-1
      RETURN
      END
      SUBROUTINE SCAL(IO,IC,AMAX1,AMAX2)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      AMAX1=0.
      AMAX2 = 0.0
C
C  FIND MAX VALUE INTENSITY PEAKS OF BOTH OV AND CV LIST
C
      DO 1 I=1,IO
      IF(AMAX1.LT.OV(I,3).AND.OV(I,4).NE.0.0) AMAX1=OV(I,3)
    1 CONTINUE
      DO 2 I=1,IO
2     OV(I,3)=OV(I,3)*100./AMAX1
C
      IF(ICON(1).EQ.1) GO TO 5
C
      DO 3 I=1,IC
      IF(AMAX2.LT.CV(I,3))AMAX2=CV(I,3)
    3 CONTINUE
      DO 4 I=1,IC
4     CV(I,3)=CV(I,3)*100./AMAX2
C
5     RETURN
      END
      SUBROUTINE ZONE(IW,IO,X2,X3,SUMIO,SUMIC,IC,PCT)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      COMMON /BLOK/DELTA(300)
      DIMENSION PCT(2,7)
      INTEGER HZ,HPZ
      LH=1
      K=2
      I=0
      IZ = 0
      IPCT = 0
C
C  FIND LIMITS ON H, K, AND L WHICH WILL PRODUCE REFLECTIONS
C  LESS THAN 10. DEGREES FROM ALL ZONES
C
      CALL MAX(JMAX,KMAX,LMAX)
      LIM = (JMAX +1) * (KMAX +1)  * (LMAX + 1) -1
      WRITE(IW,103)JMAX,KMAX,LMAX,LIM
  103 FORMAT(//35H MAXIMUM VALUES FOR H, K AND L ARE ,3I3,
     140H  THE NUMBER OF REFLECTIONS GENERATED IS,I5)
      IF(LIM.GT.300)WRITE(IW,104)
  104 FORMAT(/46H ONLY THE FIRST 300 WILL BE USED FOR Q-FACTORS,
     131H (HIGH H VALUES WILL BE OMITTED)
C
C  FIND Q FACTORS AND TABULATE AND PLOT AGAINST ANGLE
C  AWAY FROM PRINCIPAL ZONES
C
C  GENERATE FIRST ALL HKL'S AND DETERMINE THEIR Q FACTORS
C  WHICH WILL BE STORED IN GEN(I,2)
C  THE CORRESPONDING HKL WILL BE STORED IN ARRAY HKL TO
C  FACILITATE SORTING
C
      DO 5 HZ=IZ,JMAX
      DO 5 KZ=IZ,KMAX
      DO 5 LZ=IZ,LMAX
      IF(HZ.EQ.0.AND.KZ.EQ.0.AND.LZ.EQ.0) GO TO 5
C
      I=I+1
      IF(I.LE.300)  GO TO 1000
      I = 300
      GO TO 8
 1000 SUM=0.
      SUV = 0.
C
C  FIND Q
C
      DO 1 M=1,IC
C  FIND ANGLE BETWEEN THIS VECTOR AND EACH OBSERVED VECTOR
      CALL DSP(HZ,KZ,LZ,IHKL(M,1),IHKL(M,2),IHKL(M,3),ANGLE)
      ANGL = COS(ANGLE*DTR)
      SUV = SUV + CV(M,3) * ANGL
      N = CV(M,4)
      IF(N.EQ.0) GO TO 1
      SUM = SUM + OV(N,3) * ANGL
1     CONTINUE
C
      GEN(I,2)=SUM/SUMIO
      GEN(I,3) = SUV / SUMIC
C
C  PACK HKL INTO ARRAY
C
      HKL(I) = HZ*10000 + KZ*100 + LZ
5     CONTINUE
C
C  GENERATE PRINCIPAL ZONES
C
    8 DO 20  HPZ=IZ,1
      DO 20  KPZ=IZ,1
      DO 20  LPZ=IZ,1
      IF(HPZ.EQ.0.AND.KPZ.EQ.0.AND.LPZ.EQ.0) GO TO 20
      WRITE(IW,100) ITITLE,HPZ,KPZ,LPZ
      IPCT = IPCT + 1
C
C  FIND ANGLE OF Q FROM ZONE HPZ,KPZ,LPZ
C
      DO 10 IG=1,I
C
C  UNPACK HKL FROM ARRAY HKL
C
      CALL UNPACK(IG,HZ,KZ,LZ)
C
      CALL DSP(HPZ,KPZ,LPZ,HZ,KZ,LZ,ANGLE)
      IF(ANGLE.LT.1.0) ANGLE = 0.0
      GEN(IG,1)=ANGLE
C
10    CONTINUE
C
      CALL SORT(I,LH)
C
C  WRITE TABLE OF HKL ANGLE Q FACTOR
C
C THE FOLLOWING STATMENT AND NUMBER 11 WITH COMMENTS
C  IN FRONT WILL CAUSE THE DELTA Q PLOTS TO START
C  AT THE ORIGIN---CURRENT VERSION WILLL PRODUCE
C  ABSOLUTE VALUES OF DELTA Q
C     SCL = GEN(1,3) / GEN(1,2)
      DO 11 IG = 1,I
C  11 DELTA(IG) = GEN(IG,2) * SCL - GEN(IG,3)
   11 DELTA(IG) = GEN(IG,2)       - GEN(IG,3)
      PCT(1,IPCT) = DELTA(1)
      PCT(2,IPCT) = GEN(1,3)
      NN = I / 3
      DO 15 IG=1,NN
      N = IG + NN
      M = IG + NN * 2
      CALL UNPACK(IG,J1,K1,L1)
      CALL UNPACK(N,J2,K2,L2)
      CALL UNPACK(M,J3,K3,L3)
      IF(ICON(9).GT.0)WRITE(20,101)J1,K1,L1,(GEN(IG,LL),LL=1,3),
     1 DELTA(IG),
     1 J2,K2,L2,(GEN(N,LL),LL=1,3),DELTA(N),
     2 J3,K3,L3,(GEN(M,LL),LL=1,3),DELTA(M)
   15 WRITE(IW,101)J1,K1,L1,(GEN(IG,LL),LL=1,3),DELTA(IG),
     1 J2,K2,L2,(GEN(N,LL),LL=1,3),DELTA(N),
     2 J3,K3,L3,(GEN(M,LL),LL=1,3),DELTA(M)
      MM = I - NN*3
      DO16 IG = 1,MM
      LL = NN*3 + IG
      CALL UNPACK(LL,J1,K1,L1)
      IF(ICON(9).GT.0)WRITE(20,102)J1,K1,L1,(GEN(LL,MN),MN=1,3),
     1 DELTA(LL)
  16  WRITE(IW,102)J1,K1,L1,(GEN(LL,MN),MN=1,3),DELTA(LL)
C
      IF(ICON(4).EQ.0)   CALL QPLOT(K,I,X3,HPZ,KPZ,LPZ,IW)
   20 CONTINUE
C
C  FORMATS
C
100   FORMAT(1H1,20A4/10HZONE AXIS   ,3I3/
     13(3X,1HH,2X,1HK,2X,1HL,2X,5HANGLE,3X,5HQ-OBS,2X,6HQ-CALC,
     22X,7HDELTA Q)/)
101   FORMAT(3(1X,3I3,1X,F6.2,2X,F6.4,2X,F6.4,2X,F6.4,1X))
102   FORMAT(84X,1X,3I3,1X,F6.2,1X,3(1X,F6.4,1X))
C
      RETURN
      END
      SUBROUTINE UNPACK(M,J,K,L)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      IF = HKL(M)
      J = IF / 10000
      K = IF / 100 - J * 100
      L = IF - J * 10000 - K * 100
      RETURN
      END
      SUBROUTINE MAX(JMAX,KMAX,LMAX)
      COMMON OV(200,4),CV(200,4),IHKL(200,3),ITITLE(20),
     1A,B,C,AL,BE,GA,GEN(300,3),HKL(300),ICON(10),DTR,RTD
      JMAX = 0
      KMAX = 0
      LMAX = 0
      K1 = 0
      K2 = 0
      K3 = 0
      DO 100 J = 1,50
      I1 = 0
      I2 = 0
      I3 = 0
      J1 = 0
      J2 = 0
      J3 = 0
      IF(A-B)10,10,20
   10 IF(B-C)30,30,40
   30 I3 = 1
      J3 = J
      LMAX = J
   61 J2 = 1
      GO TO 70
   20 IF(A-C)50,50,60
   50 I3 = 1
      J3 = J
      LMAX = J
   63 J1 = 1
      GO TO 70
   60 I1 = 1
      J1 = J
      JMAX = J
      IF(B-C)61,61,62
   40 I2 = 1
      J2 = J
      KMAX = J
      IF(A-C)62,62,63
   62 J3 = 1
   70 CALL   DSP(I1,I2,I3,J1,J2,J3,ANG)
      IF(ANG.LT.10.)GO TO 120
  100 CONTINUE
  120 IF(I1.NE.0)K1 = 1
      IF(I2.NE.0)K2=1
      IF(I3.NE.0)K3=1
      DO 200 J = 1,50
      I1 = 0
      I2 = 0
      I3 = 0
      J1 = 0
      J2 = 0
      J3 = 0
      IF(K1)140,140,130
  130 IF(B-C)132,132,133
  132 I3 = 1
      LMAX = J
      J3 = J
  131 J1 = 1
      GO TO 195
  133 I2 = 1
      KMAX = J
      J2 = J
      GO TO 131
  140 IF(K2)160,160,150
  150 IF(A-C)152,152,154
  152 I3 = 1
      LMAX = J
      J3 = J
  153 J2 = 1
      GO TO 195
  154 I1 = 1
      JMAX = J
      J1 = J
      GO TO 153
  160 IF(A-B)180,180,190
 180  I2 = 1
      J2 = J
      KMAX = J
  185 J3 = 1
      GO TO 195
  190 I1 = 1
      J1 = J
      JMAX = J
      GO TO 185
  195 CALL DSP(I1,I2,I3,J1,J2,J3,ANG)
      IF(ANG.LT.10.)GO TO 210
  200 CONTINUE
  210 IF(I1.NE.0)K1 = 2
      IF(I2.NE.0)K2 = 2
      IF(I3.NE.0)K3 = 2
      DO 300 J = 1,50
      I1 = 0
      I2 = 0
      I3 = 0
      J1 = 0
      J2 = 0
      J3 = 0
      IF(K1.NE.0) GO TO 220
      I1 = 1
      J1 = J
      JMAX = J
      IF(K2.EQ.1)J2 = 1
  215 IF(K3.EQ.1)J3 = 1
      GO TO 250
  220 IF(K2.NE.0) GO TO 230
      I2 = 1
      J2 = J
      KMAX = J
      IF(K1.EQ.1) J1 = 1
      GO TO 215
  230 I3 = 1
      J3 = J
      LMAX = J
      IF(K1.EQ.1) J1 = 1
      IF(K2.EQ.1) J2 = 1
  250 CALL DSP(I1,I2,I3,J1,J2,J3,ANG)
      IF(ANG.LT.10.) GO TO 310
  300 CONTINUE
  310 RETURN
      END
