
	SUBROUTINE DISP(DT,XTHE,AVDISP,IDN,DTR)
	IMPLICIT REAL*8 (A-H,O-Z)
	CHARACTER*1 IDN
	IF(IDN.EQ.'1')GO TO 222
C
C  HANDLE A REFLECTION SAMPLE DISPLACEMENT FOR DIFFRACTOMETER RADIUS OF R mm
C
	R = 175.0
	IF(IDN.NE.'2')R = 200.0
	DT = (2.0 * AVDISP * COS(XTHE/2.*DTR)) / (R * DTR)
	GO TO 223
C
C  HANDLE A DISPLACEMENT FOR THE MOVING SAMPLE TRANSMISSION GUINIER
C***********NOTE: A SAMPLE ANGLE OF 45 DEG TO THE BEAM IS ASSUMED (IE ALPHA)
C
222	ALPHA = 45.0
	R = 280. * COS( ALPHA * DTR)
	DT = 2.0 * AVDISP / (DTR * R *
     1  (SIN((XTHE+ALPHA)*DTR) - ((COS((XTHE+ALPHA)*DTR) + 1.0)
     1  / TAN(XTHE*DTR))))
223	RETURN
	END

	SUBROUTINE POLY(TRUE,THET,MCO,YINT,COEF)
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION COEF(5)
	TRUE = YINT
	DO 103 J=1,MCO
103	TRUE = TRUE + COEF(J) * THET**J
	RETURN
	END

	SUBROUTINE SI(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR SILICON PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(12)
C
	DATA DSI/ 28.44271,  47.30334,  56.12309,  69.13093,  
     1 76.37729, 88.03148,  94.95399, 106.71030, 114.09426, 
     1 127.54728, 136.89663, 158.63824/
C
	NUM=12
C
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE AG(DINT,NUM)
C
C ** THIS ROUTINE RETURNS SILVER TWO-THETA VALUES IN DINT
C **	NUM IS THE NUMBER OF PEAKS
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DAG(10)
	DATA DAG/38.112,44.295,64.437,77.390,81.533,97.875,
     +	         110.499,114.914,134.871,156.737/
	NUM=10
	DO 10 I=1,NUM
10	DINT(I)=DAG(I)
	RETURN
	END

	SUBROUTINE BI(DINT,NUM)
C
C ** THIS ROUTINE RETURNS BI12GEO20 TWO-THETA VALUES IN DINT
C **	NUM IS THE NUMBER OF PEAKS
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DAG(80)
	DATA DAG/
     1 12.3281, 17.4685, 21.4365, 24.8018, 27.7848, 30.4981, 33.0088,
     1 35.3604, 37.5831, 39.6991, 41.7247, 43.6730, 45.5543, 49.1474,
     1 50.8719, 52.5552, 54.2014, 55.8142, 57.3966, 58.9517, 60.4817,
     1 61.9889, 63.4753, 64.9427, 66.3927, 67.8268, 69.2465, 70.6530,
     1 73.4309, 74.8046, 76.1693, 77.5261, 78.8759, 80.2194, 81.5575,
     1 82.8910, 84.2206, 85.5471, 86.8712, 88.1936, 89.5151, 90.8363,
     1 92.1580, 94.8055, 96.1327, 97.4633, 98.7979,100.1374,101.4824,
     1102.8340,104.1928,105.5598,108.3223,109.7199,111.1297,112.5531,
     1113.9913,116.9181,118.4099,119.9229,121.4594,123.0214,124.6117,
     1126.2330,127.8887,129.5825,131.3188,133.1026,134.9400,136.8382,
     1138.8060,140.8543,145.2514,147.6419,150.2015,152.9787,156.0495,
     1159.5464,163.7428,169.4431/
	NUM=80
	DO 10 I=1,NUM
10	DINT(I)=DAG(I)
	RETURN
	END


	SUBROUTINE BS(DINT,NUM)
C
C ** THIS ROUTINE RETURNS BI12GEO20 + SI TWO-THETA VALUES IN DINT
C **	NUM IS THE NUMBER OF PEAKS
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DAG(92)
	DATA DAG/
     1 12.3281, 17.4685, 21.4365, 24.8018, 27.7848, 30.4981, 33.0088,
     1 35.3604, 37.5831, 39.6991, 41.7247, 43.6730, 45.5543, 49.1474,
     1 50.8719, 52.5552, 54.2014, 55.8142, 57.3966, 58.9517, 60.4817,
     1 61.9889, 63.4753, 64.9427, 66.3927, 67.8268, 69.2465, 70.6530,
     1 73.4309, 74.8046, 76.1693, 77.5261, 78.8759, 80.2194, 81.5575,
     1 82.8910, 84.2206, 85.5471, 86.8712, 88.1936, 89.5151, 90.8363,
     1 92.1580, 94.8055, 96.1327, 97.4633, 98.7979,100.1374,101.4824,
     1102.8340,104.1928,105.5598,108.3223,109.7199,111.1297,112.5531,
     1113.9913,116.9181,118.4099,119.9229,121.4594,123.0214,124.6117,
     1126.2330,127.8887,129.5825,131.3188,133.1026,134.9400,136.8382,
     1138.8060,140.8543,145.2514,147.6419,150.2015,152.9787,156.0495,
     1159.5464,163.7428,169.4431,
     D 28.44271,  47.30334,  56.12309,  69.13093,  
     E 76.37729, 88.03148,  94.95399, 106.71030, 114.09426, 
     F 127.54728, 136.89663, 158.63824/
	NUM=92
C
C  SORT THEM TOGETHER
C
	DO 5 I=1,NUM-1
	DO 5 J=I,NUM
	IF(DAG(I).LE.DAG(J))GO TO 5
	T = DAG(I)
	DAG(I) = DAG(J)
	DAG(J) = T
5	CONTINUE
	DO 10 I=1,NUM
10	DINT(I)=DAG(I)
	RETURN
	END

	SUBROUTINE NI(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR NICKEL PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(8)
	DATA DSI/
     1 44.4970, 51.8505, 76.3829, 92.9395, 98.4446,121.9471,144.6714,
     1155.7024/
C
	NUM=8
C
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE NO(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR NICKEl OXIDE PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(10)
	DATA DSI/
     1 37.2561, 43.2880, 62.8815, 75.4177, 79.4118, 95.0678, 
     2 107.0016, 111.1259, 129.2350, 146.7792/
C
	NUM=10
C
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE CU(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR COPPER PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(8)
C
	DATA DSI/ 43.31686, 50.44927, 74.12602, 89.93736, 95.14729,
     2 116.93344, 136.50079, 144.70555/
C
	NUM=8
C
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE FL(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR FLUOROPHLOGOPITE PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(12)
	DATA DSI/
     1  8.8525, 17.7585, 26.7742, 35.9624, 45.3969, 55.1691, 65.3991,
     1 76.2545, 87.9885,101.0248,116.1928,135.6741/
	NUM = 12
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE FC(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR FLUOROPHLOGOPITE AND CU PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(20)
	DATA DSI/
     1  8.8525, 17.7585, 26.7742, 35.9624, 43.3169, 45.3969, 50.44927, 
     2 55.1691, 65.3991, 74.1260, 76.2545, 87.9885, 89.9374, 95.14729,
     3101.0248,116.1928,116.9334, 135.674,136.5008, 144.70555/
	NUM = 20
	DO 10 I=1,NUM
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

	SUBROUTINE FS(DINT,NUM)
C
C **	THIS ROUTINE RETURNS TWO-THETA VALUES FOR SILICON AND
C **	FLUOROPHLOGOPITE (NBS SRM) PEAKS
C **	NUM IS THE NUMBER OF PEAKS
C
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION DINT(1),DSI(24)
C
	DATA DSI/ 28.44271,  47.30334,  56.12309,  69.13093,  
     1 76.37729, 88.03148,  94.95399, 106.71030, 114.09426, 
     1 127.54728, 136.89663, 158.63824,
     1  8.8525, 17.7585, 26.7742, 35.9624, 45.3969, 55.1691, 65.3991,
     1 76.2545, 87.9885,101.0248,116.1928,135.6741/
C
	NUM = 24
C
C  SORT THEM TOGETHER
C
	DO 5 I=1,NUM-1
	DO 5 J=I,NUM
	IF(DSI(I).LE.DSI(J))GO TO 5
	T = DSI(I)
	DSI(I) = DSI(J)
	DSI(J) = T
5	CONTINUE
C
	DO 10 I=1,12
	DINT(I)=DSI(I)
10	CONTINUE
	RETURN
	END

      FUNCTION REGR (M, X, Y, N, A, B)
C
C     THIS PROGRAM ORIGINALLY CAME FROM THE BOOK:
C
C     B. CARNAHAN, H. LUTHER, AND J. WILKES, APPLIED NUMERICAL
C     METHODS,  JOHN WILEY AND SONS, NEW YORK, 1969, P580.
C
C
C            FUNCTION REGR, WHICH CARRIES OUT AN NTH-ORDER
C            POLYNOMIAL REGRESSION ON M DATA POINTS CONTAINED
C            IN THE X AND Y ARRAYS.  THE FUNCTION NORMALLY
C            RETURNS THE STANDARD DEVIATION S OF THE POINTS
C            ABOUT THE REGRESSION CURVE.  THE REGRESSION
C            COEFFICIENTS ARE PLACED IN A AND B(1)...B(N).
C            HOWEVER, IF THE SIMULTANEOUS EQUATION SOLVING
C            ROUTINE SIMUL ENCOUNTERS A NEAR-SINGULAR MATRIX,
C            THE FUNCTION RETURNS THE VALUE 0.0 AS REGR.
C
C
C     PARAMETER DEFINITIONS:
C
C     M - THE NUMBER OF DATA POINTS IN X AND Y.
C
C     X,Y - THE ARRAYS CONTAINING THE ABSCISSA AND ORDINATE
C           DATA VALUES RESPECTIVELY
C
C     N - THE ORDER OF THE POLYNOMIAL DESIRED FOR THE FIT.
C
C     A - Y INTERCEPT OF THE POLYNOMIAL.
C
C     B - ARRAY CONTAINING THE COEFFICIENTS OF THE
C         POLYNOMIAL.
C
C     REGR - RETURNS AS THE VALUE OF THE STANDARD
C            DEVIATION OF THE FIT TO THE DATA. IF
C            A NEAR SINGULAR OR IF A SINGULAR
C            MATRIX IS ENCOUNTERED THE VALUE OF
C            0 (ZERO) IS RETURNED AS REGR.
C
C
C     *****  NOTE  *****
C
C     THE POLYNOMIAL HAS THE FORM OF:
C
C     ORDINATE = A + B(I) * (ABSCISSA ** I)
C        FOR I=1 TO N
C
C
      IMPLICIT REAL*8 (A - H, O - Z)
      REAL*8 A, B, REGR, X, Y
      DIMENSION C(11,11), SX(20), SYX(10), CYX(10), X(100),Y(100),
     1B(10)
      DATA EPS/ 1.0D-20/
C
C     .....  COMPUTE SUMS OF POWERS AND PRODUCTS  .....
      IO6 = 6
      NTWO = 2*N
      NP1 = N + 1
      SY = 0.0
      SYY = 0.0
      DO 1 I=1,N
      NPI = N + I
      SX(I) = 0.0
      SX(NPI) = 0.0
      SYX(I) = 0.0
1     CONTINUE
      DO 3 I=1,M
      SY = SY + Y(I)
      SYY = SYY + Y(I)**2
      DUM=1.0
      DO 2 J=1,N
      DUM = DUM*X(I)
      SX(J) = SX(J) + DUM
      SYX(J) = SYX(J) + Y(I)*DUM
2     CONTINUE
      DO 3 J=NP1,NTWO
      DUM = DUM*X(I)
      SX(J) = SX(J) + DUM
3     CONTINUE
C
C     .....  COMPUTE COEFFICIENTS C(I,J)  .....
C
      FM = M
      CYY = SYY - SY*SY/FM
      DO 4 I=1,N
      CYX(I) = SYX(I) - SY*SX(I)/FM
      C(I,NP1) = CYX(I)
      DO 4 J=1,N
      IPJ = I + J
      C(I,J) = SX(IPJ) - SX(I)*SX(J)/FM
4     CONTINUE
C
C     .....  CALL ON SIMUL TO SOLVE SIMULTANEOUS EQUATIONS
C
      DET = SIMUL (N, C, B, EPS, 1, 11)
C
C     I INSERTED THE COMMENTS INFRONT OF THE WRITE
C     STATEMENT BECAUSE I DIDN'T WANT TO PRINT IT OUT.
C
C     WRITE(IO6,200) DET
C200  FORMAT(13H0       DET = ,E14.8)
      IF(DET.NE.0.0) GO TO 6
      REGR = 0.0
      RETURN
C
C     .....  COMPUTE INTERCEPT A AND STANDARD DEVIATION S  .....
C
6     DUM = SY
      TEMP = CYY
      DO 7 I=1,N
      DUM = DUM - B(I)*SX(I)
      TEMP = TEMP - B(I)*CYX(I)
7     CONTINUE
      A = DUM/FM
      DENOM = M - N - 1
C
C     *****  NOTE  *****
C
C     I INSERTED THE FOLLOWING IF STATEMENT TO ELIMINATE
C     TEMP BEING NEAR ZERO BUT SLIGHTLY NEGATIVE.
C     AN EXAMPLE IS IF TEMP EQUALS -.12345E-15.
C
C     THIS STATEMENT IS NOT SUPPLIED IN
C     THE BOOK VERSION, BUT RATHER IS A HOMEBREW PATCH.
C
C     ******************
      IF(TEMP.LE.EPS) TEMP=EPS
      S = DSQRT(TEMP/DENOM)
      REGR = S
      RETURN
      END

      FUNCTION SIMUL( N, A, X, EPS, INDIC, NRC)
C
C     THIS PROGRAM ORIGINALLY CAME FROM THE BOOK:
C
C     B. CARNAHAN, H. LUTHER, AND J. WILKES, APPLIED NUMERICAL
C     METHODS,  JOHN WILEY AND SONS, NEW YORK, 1969, P290
C
C
C        WHEN INDIC IS NEGATIVE, SIMUL COMPUTES THE INVERSE
C        OF THE N BY N MATRIX A IN PLACE.  WHEN INDIC IS ZERO,
C        SIMUL COMPUTES THE N SOLUTIONS X(1)...X(N) CORRESPONDING
C        TO THE SET OF LINEAR EQUATIONS WITH AUGMENTED MATRIX
C        OF COEFFICIENTS IN THE N BY N+1 ARRAY A AND IN
C        ADDITION COMPUTES TH INVERSE OF THE COEFFICIENT
C        MATRIX IN PLACE AS ABOVE.  IF INDIC IS POSITIVE,
C        THE SET OF LINEAR EQUATIONS IS SOLVED BUT THE
C        INVERSE IS NOT COMPUTED IN PLACE.  THE GAUSS-
C        JORDAN COMPLETE ELIMINATION METHOD IS EMPLOYED
C        WITH THE MAXIMUM PIVOT STRATEGY.  ROW AND
C        COLUMN SUBSCRIPTS OF SUCCESIVE PIVOT ELEMENTS ARE
C        SAVED IN ORDER IN THE IROW AND JCOL ARRAYS
C        RESPECTIVELY.  K IS THE PIVOT COUNTER, PIVOT IS
C        THE ALGEBRAIC VALUE OF THE PIVOT ELEMENT. MAX
C        THE NUMBER OF COLUMNS IN A AND DETER THE DETERMINANT
C        OF THE COEFFICIENT MATRIX.  THE SOLUTIONS ARE
C        DETERMINED IN THE N+1TH COLUMN OF A AND THEN
C        UNSCRAMBLED AND PUT IN ORDER IN X(1)...X(N)
C        USING THE PIVOT SUBSCRIPT INFORMATION AVAILABLE
C        IN THE IROW AND JCOL ARRAYS.  THE SIGN OF THE
C        DETERMINANT IS ADJUSTED, IF NECESSARY, BY DETERMINING
C        IF AN EVEN OR ODD NUMBER OF PAIRWISE INTERCHANGES IS
C        REQUIRED TO PUT THE ELEMENTS OF THE JORD ARRAY IN
C        ASCENDING SEQUENCE WHERE JORD(IROW(I))  =  JCOL(I).
C        IF THE INVERSE IS REQUIRED, IT IS UNSCRAMBLED IN PLACE
C        USING Y(1)...Y(N) AS TEMPORARY STORAGE.  THE VALUE OF
C        THE DETERMINANT IS GIVEN AS THE VALUE OF THE FUNCTION.
C        SHOULD THE POTENTIAL PIVOT OF LARGEST MAGNITUDE BE
C        SMALLER IN MAGNITUDE THAN EPS, THE MATRIX IS
C        CONSIDERED TO BE SINGULAR AND A TRUE ZERO IS RETURNED
C        AS THE VALUE OF THE FUNCTION.
C
      IMPLICIT REAL*8(A-H, O-Z)
      REAL*8 A, X, EPS, SIMUL
      DIMENSION IROW(50),JCOL(50),JORD(50),Y(50),A(NRC,NRC),X(N)
C
      MAX = N
      IO6 = 6
      IF(INDIC.GE.0) MAX = N+1
C
C     .....  IS N LARGER THAN 50  .....
      IF(N.LE.50) GO TO 5
      WRITE(IO6,200)
      SIMUL = 0.0
      RETURN
C
C     .....  BEGIN ELIMINATION PROCEDURE  .....
5     DETER = 1.
      DO 18 K = 1,N
      KM1 = K - 1
C
C     .....  SEARCH FOR THE PIVOT ELEMENT  .....
      PIVOT = 0.0
      DO 11 I = 1,N
      DO 11 J = 1,N
C     .....  SCAN IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS  .....
      IF(K.EQ.1) GO TO 9
      DO 8 ISCAN = 1,KM1
      DO 8 JSCAN = 1,KM1
      IF(I.EQ.IROW(ISCAN)) GO TO 11
      IF(J.EQ.JCOL(JSCAN)) GO TO 11
8     CONTINUE
9     IF(DABS(A(I,J)).LE.DABS(PIVOT)) GO TO 11
      PIVOT = A(I,J)
      IROW(K) = I
      JCOL(K) = J
11    CONTINUE
C
C     .....  INSURE THAT SELECTED PIVOT IS LARGER THAN EPS  .....
      IF(DABS(PIVOT).GT.EPS) GO TO 13
      SIMUL = 0.0
      RETURN
C
C     .....  UPDATE THE DETERMINANT VALUE  .....
13    IROWK = IROW(K)
      JCOLK = JCOL(K)
      DETER = DETER*PIVOT
C
C     .....  NORMALIZE PIVOT ROW ELEMENTS  .....
      DO 14 J = 1,MAX
      A(IROWK,J) = A(IROWK,J)/PIVOT
14    CONTINUE
C
C     .....  CARRY OUT ELIMINATION AND DEVELOP INVERSE .....
      A(IROWK,JCOLK) = 1./PIVOT
      DO 18 I = 1,N
      AIJCK = A(I,JCOLK)
      IF(I.EQ.IROWK) GO TO 18
      A(I,JCOLK) = - AIJCK/PIVOT
      DO 17 J = 1,MAX
      IF(J.NE.JCOLK) A(I,J) = A(I,J) - AIJCK*A(IROWK,J)
17    CONTINUE
18    CONTINUE
C
C     .....  ORDER SOLUTION VALUES (IF ANY) AND CREATE JORD ARRAY  .....
      DO 20 I = 1,N
      IROWI = IROW(I)
      JCOLI = JCOL(I)
      JORD(IROWI) = JCOLI
      IF(INDIC.GE.0) X(JCOLI) = A(IROWI,MAX)
20    CONTINUE
C
C     .....  ADJUST SIGN OF DETERMINANT  .....
      INTCH = 0
      NM1 = N - 1
      DO 22 I = 1,NM1
      IP1 = I + 1
      DO 22 J = IP1,N
      IF(JORD(J).GE.JORD(I)) GO TO 22
      JTEMP = JORD(J)
      JORD(J) = JORD(I)
      JORD(I) = JTEMP
      INTCH = INTCH + 1
22    CONTINUE
      IF(INTCH.NE.INTCH/2*2) DETER = - DETER
C
C     .....  IF INDIC IS POSITIVE RETURN WITH RESULTS  .....
24    IF(INDIC.LE.0) GO TO 26
      SIMUL = DETER
      RETURN
C
C     .....  IF INDIC IS NEGATIVE OR ZERO, UNSCRAMBLE THE INVERSE
C            FIRST BY ROWS  .....
26    DO 28 J = 1,N
      DO 27 I = 1,N
      IROWI = IROW(I)
      JCOLI = JCOL(I)
      Y(JCOLI) = A(IROWI,J)
27    CONTINUE
      DO 28 I = 1,N
      A(I,J) = Y(I)
28    CONTINUE
C     .....  THEN BY COLUMNS  .....
      DO 30 I = 1,N
      DO 29 J = 1,N
      IROWJ = IROW(J)
      JCOLJ = JCOL(J)
      Y(IROWJ) = A(I,JCOLJ)
29    CONTINUE
      DO 30 J = 1,N
      A(I,J) = Y(J)
30    CONTINUE
C
C     .....  RETURN FOR INDIC NEGATIVE OR ZERO  .....
      SIMUL = DETER
      RETURN
C
C     .....  FORMAT FOR OUTPUT STATEMENT  .....
200   FORMAT(10H0N TOO BIG )
      END

	SUBROUTINE HKL(I1,I2,I3,T1,D1,N)
	IMPLICIT REAL*8(A-H,O-Z)
C
C  2-THETA AND HKL FINDER
C
	INTEGER H
	DIMENSION T1(200),D1(200),I1(200),I2(200),I3(200)
	DATA WAVE / 1.5405981/, IREAD/5/, IWRIT/6/
	DTR = 3.141592654 / 180.
110	WRITE(IWRIT,1000)
1000	FORMAT(/' HKL AND 2 THETA COMPUTATION'//
     1' INPUT ONLY THE NEEDED LATTICE PARAMETERS FOR YOUR SYMMETRY'/
     1'$a, b, c, alpha, beta, gamma: ')
	READ(IREAD,1010)A,B,C,AL,BE,GA
1010	FORMAT(6G10.5)
	WRITE(IWRIT,1005)
1005	FORMAT(/' INPUT THE WAVELENGTH (if other than Cu Ka1)')
	READ(IREAD,1010)W
	IF(W.NE.0.0)WAVE = W
	WRITE(IWRIT,1011)
1011	FORMAT(/' INPUT LATTICE TYPE (P,I OR F)')
	IF(B.EQ.0.0)B = A
	IF(C.EQ.0.0)C = A
	IF(AL.EQ.0.0)AL = 90.
	IF(BE.EQ.0.0)BE = 90.
	IF(GA.EQ.0.0)GA = 90.
	READ(IREAD,1012)LAT
1012	FORMAT(A1)
	ILAT = 1
	IF(LAT.EQ.1HI)ILAT = 2
	IF(LAT.EQ.1HF)ILAT = 3
	WRITE(IWRIT,82)A,B,C,AL,BE,GA
   82 FORMAT(//' REAL CELL DIMENSIONS'/
     1' A = ',F10.5,' B = ',F10.5,' C = ',F10.5/
     2' ALPHA = ',F10.3,' BETA = ',F10.3,' GAMMA = ',F10.3)
	AL = AL * DTR
	BE = BE * DTR
	GA = GA * DTR
  802 S = AL+BE+GA
      S=S/2.0
      V=2.0*A*B*C*SQRT(SIN(S)*SIN(S-AL)*SIN(S-
     1BE)*SIN(S-GA))
      WRITE(IWRIT,83) V, LAT
   83 FORMAT(' VOLUME = ',F10.5, ' CUBIC ANGSTROMS   LATTICE = ',A1)
      AS = B*C*SIN(AL)/V
      BS = A*C*SIN(BE)/V
      CS = A * B* SIN(GA)/V
      V = 1./V
	AL = 180.- AL/DTR
	BE = 180.- BE/DTR
	GA = 180.- GA/DTR
      WRITE(IWRIT,72)AS,BS,CS,AL,BE,GA,V
   72 FORMAT(/' RECIPROCAL CELL DIMENSIONS'/
     1' A* = ',F10.6,' B* = ',F10.6,' C* = ',F10.6/
     2' ALPHA* = ',F10.3,' BETA* = ',F10.3,'GAMMA* = ',F10.6/
     3' RECIPROCAL VOLUME = ',G20.8//)
	AL = AL * DTR
	BE = BE * DTR
	GA = GA * DTR
C      WRITE(IWRIT,92)
C   92 FORMAT(//'  H  K  L     D*        D      2*THETA')
C
C  SET LIMITS FOR DIFFERENT CRYSTAL SYSTEMS
C
C  FIRST TRICLINIC
C
	KMIN = 0
	KMAX = 8
	IF(ABS(AL-90.) .LT. .001)KMIN=-8
C
C  NOW MONOCLINIC
C
	LMIN = 0
	LMAX = 8
	IF(ABS(BE-90.) .LT. .001 .OR. ABS(GA-90.)
     1 .LT. .001)LMIN = -8
	N = 0
      DO 900 J=0,8
      DO 900 K=KMIN,KMAX
      DO900 L=LMIN,LMAX
      IF(J .EQ.0.AND.K.EQ.0.AND.L.EQ.0)GO TO 900
	GO TO(850,830,840),ILAT
C
C  TEST FOR BODY CENTERING
C
830	H = JMOD(J+K+L,2)
	IF(H.EQ.0)GO TO 850
	GO TO 900
C
C  TEST FOR F CENTERING CONDITION
C
840	H = JMOD(J,2) + JMOD(K,2) + JMOD(L,2)
	IF(H.EQ.0.OR.H.EQ.3)GO TO 850
	GO TO 900
850      DS = J*J*AS*AS + K*K*BS*BS + L*L*CS*CS +
     1 2*J*K*COS(GA)*AS*BS +
     2 2*J*L*AS*CS*COS(BE) +
     3 2*K*L*BS*CS*COS(AL)
      DS = SQRT(DS)
      D = 1.0/DS
	T2 = WAVE / (2.0 * D)
	IF(T2.GT.1.0)GO TO 900
C
C  HAS THIS D ALREADY BEEN FOUND?
C
	IF(N.EQ.0)GO TO 160
	DO 150 I=1,N
	IF(D1(I).EQ.D)GO TO 900
150	CONTINUE
160	N = N + 1
	IF(N.GT.200)GO TO 250
	IF(IPRIR.EQ.'T')GO TO 250
	D1(N) = D
	I1(N) = J
	I2(N) = K
	I3(N) = L
	T2 = ASIN( T2 ) / DTR * 2.0
	T1(N) = T2
C      WRITE(IWRIT,93)J,K,L,DS,D,T1(N)
C   93 FORMAT(3I3,3F10.5)
  900 CONTINUE
200	CONTINUE
C
C  SORT ON DECREASING d
250	DO 220 I=1,N-1
	DO 220 J=I+1,N
	IF(D1(I).GT.D1(J))GO TO 220
	T = D1(I)
	D1(I) = D1(J)
	D1(J) = T
	T = T1(I)
	T1(I) = T1(J)
	T1(J) = T
	IT = I1(I)
	I1(I) = I1(J)
	I1(J) = IT
	IT = I2(I)
	I2(I) = I2(J)
	I2(J) = IT
	IT = I3(I)
	I3(I) = I3(J)
	I3(J) = IT
220	CONTINUE
	RETURN
	END

