      PROGRAM DSPACE
c	program written by kenn gardner to calculate the reflections
c	ie d-spacing, for any unit cell
      REAL		A			! UNIT CELL PARAM.
	REAL		AL			! UNIT CELL PARAM. - ALPHA
	REAL		B			! UNIT CELL PARAM.
	REAL		BE			! UNIT CELL PARAM. - BETA
	REAL		C			! UNIT CELL PARAM.
	REAL		D(10000)	! DSPACING
	REAL		GA			! UNIT CELL PARAM. - GAMMA
	INTEGER 	H(10000) 	! MILLER INDEX
      INTEGER		HH(10000)   ! MILLER INDEX 
	REAL		LAMBDA		! WAVELENGTH
      INTEGER		K(10000)    ! MILLER INDEX 
      INTEGER		KK(10000)   ! MILLER INDEX 
      INTEGER		L(10000)    ! MILLER INDEX 
      INTEGER		LL(10000)   ! MILLER INDEX 
	REAL		RO(10000)
	CHARACTER *60	TITLE
	REAL		TT(10000)	! TWO THETA
	REAL		VOL			! UNIT CELL VOLUME
C
C
C
      PI=3.1415926536
      OPEN(UNIT=21,FILE='DSPACE.OUT',TYPE='NEW')
      WRITE(6,35)
35    FORMAT(' TITLE - ',$)
      READ(5,42) TITLE
42	FORMAT (A)
      WRITE(21,37) TITLE
37    FORMAT(1X,A)
	WRITE(6,1)     
1	FORMAT('      1 - TRICLINIC',/,'      2 - MONOCLINIC',/,
     1	   '      3 - HEXAGONAL ',/,'      4 - ORTHORHOMBIC ',/,
     2       '      5 - TETRAGONAL ',/,'      6 - CUBIC ',/,
     3	   ' CRYSTAL CLASS - ',$) 
	READ (5,*) IGRP
100	WRITE(6,2)
2	FORMAT (/,' ASSUME C-AXIS IS UNIQUE ',/,'      A - ',$)
	READ(5,*) A
C	SET DEFAULTS
	B=A
	C=A
	AL=90
	BE=90
	GA=90
	IF (IGRP.EQ.1) THEN	!TRICLINIC
		WRITE(6,3)
3		FORMAT('      B - ',$)
		READ(5,*) B
		WRITE(6,4)
4		FORMAT('      C - ',$)
		READ(5,*) C
		WRITE(6,5)
5		FORMAT('      ALPHA - ',$)
		READ(5,*) AL
		WRITE(6,6)
6		FORMAT('      BETA - ',$)
		READ(5,*) BE
		WRITE(6,7)
7		FORMAT('      GAMMA - ',$)
		READ(5,*) GA
	ENDIF
	IF (IGRP.EQ.2) THEN		!MONOCLINIC
		WRITE(6,3)
		READ(5,*) B
		WRITE(6,4)
		READ(5,*) C
		WRITE(6,7)
		READ(5,*) GA
	ENDIF	
	IF (IGRP.EQ.3) THEN		! HEXAGONAL
		WRITE(6,4)
		READ(5,*) C
		GA=120.0
	ENDIF	
	IF (IGRP.EQ.4) THEN		! ORTHORHOMBIC
		WRITE(6,3)
		READ(5,*) B
		WRITE(6,4)
		READ(5,*) C
	ENDIF	
	IF (IGRP.EQ.5) THEN		!TETRAGONAL
		WRITE(6,4)
		READ(5,*) C
	ENDIF	

      WRITE(6,40) A,B,C,AL,BE,GA
40    FORMAT(/' UNIT CELL PARAMETERS ',/,3F7.3,3F7.2)
      WRITE(6,57)
57    FORMAT(' OK? ',$)
      CALL YESNO(KKEY)
      IF(KKEY.NE.1)GO TO 100 
      WRITE(21,40) A,B,C,AL,BE,GA
      AL=AL*PI/180.
      BE=BE*PI/180.
      GA=GA*PI/180.
      ABC=(AL+BE+GA)/2.
      VOL=2.*A*B*C*SQRT(SIN(ABC)*SIN(ABC-AL)*SIN(ABC-BE)*
     1      SIN(ABC-GA))
      CA=COS(AL)
      CB=COS(BE)
      CG=COS(GA)
      SA=SIN(AL)
      SB=SIN(BE)
      SG=SIN(GA)
      AS=B*C*SA/VOL
      BS=C*A*SB/VOL
      CS=A*B*SG/VOL
      CAS=(CB*CG-CA)/(SB*SG)
      CBS=(CG*CA-CB)/(SG*SA)
      CGS=(CA*CB-CG)/(SA*SB)
      SGS=SQRT(1.-CGS**2)
C	DEFINE WAVELENGTH
	WRITE(6,47)
47    FORMAT(' 1 - Cu ',/,' 2 - Mo ',/,' 3 - Other ',/,
     1' XRAY RADIATION - ',$)
      read(5,*) IRAD
      IF (IRAD.EQ.1) LAMBDA=1.5418
	IF (IRAD.EQ.2) LAMBDA=0.7107
	IF (IRAD.EQ.3) THEN
		WRITE(6,891)
891		FORMAT (' WAVELENGTH (ANGSTROMS) - ',$)
		READ (5,*) LAMBDA
	ENDIF
	WRITE(6,49) LAMBDA
      WRITE(21,49) LAMBDA
49    FORMAT(/' XRAY RADIATION - ',F10.5)
C
	IPOW=0
	WRITE(6,89)
89	FORMAT(/'  ARRANGE DATA AS 1 - POWDER PATTERN ',/,
     1       '                  0 - LAYER LINES ? ',$)
	READ(5,*) IPOW
9     WRITE(6,50)
50    FORMAT(/' DEFINE CUTOFF: 1 - TWO THETA CUTOFF '/
     1       '                2 - DSPACING CUTOFF ? ',$)
      READ (5,*) KKEY
	IF (KKEY.EQ.1) THEN
		WRITE(6,51)
51		FORMAT(/' TWO THETA CUTOFF (DEFAULT=50DEG) - ',$)
		read(5,39) TTHETA
39		FORMAT(F)
		IF(TTHETA.EQ.0)TTHETA=50.
		WRITE(21,52) TTHETA
52		FORMAT(/' TWO THETA CUTOFF - ',F10.4)
		TTHETA=TTHETA*PI/180.
		DCUTOF=.5*LAMBDA/SIN(TTHETA/2.0)
      ELSE
		WRITE(6,54)
54		FORMAT(/' D-SPACING CUTOFF (DEFAULT=1.5A) - '$)
		read(5,39) DCUTOF
		IF(DCUTOF.LE.0)DCUTOF=1.5
		IF(DCUTOF.LT.LAMBDA/2.)DCUTOF=LAMBDA/2.
		WRITE(21,55) DCUTOF
55		FORMAT(/' D-SPACING CUTOFF - ',F10.4)
		TTHETA=2.*ASIN(LAMBDA/(2*DCUTOF))
	END IF
      RHOMAX=1./DCUTOF
C     NEXT DETERMINE THE NUMBER OF REFLECTIONS THAT WILL BE
C     CALCULATED
      SPVOL=4.*PI*(RHOMAX**3)/3.
      NOREFL=SPVOL*VOL
      IF(IGRP.EQ.1)NOREFL=NOREFL/2.
      IF(IGRP.EQ.2.AND.IGRP.EQ.3)NOREFL=NOREFL/4
      IF(IGRP.GE.4)NOREFL=NOREFL/8
C     DETERMINE RANGE OF MILLER INDICES
      DO 12 I1=1,100
	    ROL=SQRT(((I1-1)**2)*CS**2)
	    IF(ROL.GT.RHOMAX)GO TO 13
12    CONTINUE
13    IL=I1-1
      DO 14 I2=1,100
	      ROK=SQRT(((I2-1)**2)*BS**2)
	      IF(ROK.GT.RHOMAX)GO TO 15
14    CONTINUE
15    IK=I2-1
      DO 16 I3=1,100
	      ROH=SQRT(((I3-1)**2)*AS**2)
	      IF(ROH.GT.RHOMAX)GO TO 17
16    CONTINUE
17    IH=I3-1
C      DETERMINE STARTING INDICES
      IHEND=IH
      IKEND=IK
      ILEND=IL
      IHBGN=0
      IKBGN=0
      ILBGN=0
      IF(IGRP.EQ.1) THEN
      IHBGN=-IH
	IKBGN=-IK
	ENDIF
      IF(IGRP.EQ.2) THEN
	IKBGN=-IK
	ENDIF
      IF(IGRP.EQ.3) THEN
	IKBGN=-IK
	ENDIF

      DTHET2=TTHETA*180./PI
      WRITE(6,56) DTHET2,DCUTOF,NOREFL,IHBGN,IHEND,IKBGN,
     1      IKEND,ILBGN,ILEND
56    FORMAT(2F10.4,I5/(2I5/)) 
      WRITE(21,56) DTHET2,DCUTOF,NOREFL,IHBGN,IHEND,IKBGN,
     1      IKEND,ILBGN,ILEND
      IHFIN=(IHEND-IHBGN)+1
      IKFIN=(IKEND-IKBGN)+1
      ILFIN=(ILEND-ILBGN)+1
      IHCOR=IHBGN-1
      IKCOR=IKBGN-1
      ILCOR=ILBGN-1
c
c
c
c	
	IF (IPOW.EQ.1) GO TO 1001
      DO 34 IL=1,ILFIN
	      LL(IL)=ILCOR+IL
	      WRITE(21,64) LL(IL)
64    FORMAT(//,' CALCULATED VALUES FOR LAYER',I4,/
     1      '   H   K   L',5X,'D-SPACING',5X,'TWO THETA')
     		INT=0
		DO 29 IK=1,IKFIN
			KK(IK)=IKCOR+IK
			DO 28 IH=1,IHFIN
				HH(IH)=IHCOR+IH
				INT=INT+1
				H(INT)=HH(IH)
				K(INT)=KK(IK)
				L(INT)=LL(IL)
      ROSAVE=(H(INT)*AS)**2+(K(INT)*BS)**2+(L(INT)*CS)**2+
     1     2.*K(INT)*L(INT)*BS*CS*CAS+2.*L(INT)*H(INT)*
     2      CS*AS*CBS+2.*H(INT)*K(INT)*AS*BS*CGS
			    RO(INT)=SQRT(ROSAVE)
			    IF(RO(INT).LE.RHOMAX)GO TO 28
			    INT=INT-1
28    		    CONTINUE
29	    CONTINUE
C     ARRANGE RHO VALUES IN INCREASING ORDER OF MAGNITUDES
C
		INT1=INT-1
		INT2=INT
		DO 31 J=1,INT1
			INT2=INT2-1
		    DO 30 I=1,INT2
				IF(RO(I+1).GE.RO(I))GO TO 30
			    Z=RO(I+1)
			    RO(I+1)=RO(I)
			    RO(I)=Z
			    IZ=K(I+1)
			    K(I+1)=K(I)
			    K(I)=IZ
			    IZ=H(I+1)
			    H(I+1)=H(I)
			    H(I)=IZ
30		    CONTINUE
31	    CONTINUE
C      PRINT OUT RESULTS
C
1003	CONTINUE
	    DO 33 I=1,INT
			IF(RO(I).EQ.0)GO TO 33
c
c	fix dated nov 85 to remove redundant reflections from hex 
c       and tet cases
			idum = ifix(10000./ro(i))
		
			idum1 = ifix(10000./ro(i+1))
		      IF(idum.NE.idum1)GO TO 32
      IF((ABS(H(I))+ABS(K(I))+ABS(L(I))).EQ.(ABS(H(I+1))+ABS(K(
     1 I+1))+ABS(L(I+1))))GO TO 33
32			CONTINUE
			D(I)=1./RO(I)
			TT(I)=2.*180.*ASIN(LAMBDA/(2.*D(I)))/PI
			WRITE(21,65) H(I),K(I),L(I),D(I),TT(I)
65			FORMAT(3I4,2F15.4)
33	    CONTINUE
34    CONTINUE
      STOP
1001	CONTINUE
	INT=0
	DO 91 IL=1,ILFIN
		LL(IL)=ILCOR+IL
		DO 92 IK=1,IKFIN
			KK(IK)=IKCOR+IK
			DO 93 IH=1,IHFIN
				HH(IH)=IHCOR+IH
				INT=INT+1
				IF (INT.GE.10000) GO TO 1002
	
				H(INT)=HH(IH)
				K(INT)=KK(IK)
				L(INT)=LL(IL)
	ROSAVE=(H(INT)*AS)**2+(K(INT)*BS)**2+(L(INT)*CS)**2+
	1	2.*K(INT)*L(INT)*BS*CS*CAS+2.*L(INT)*H(INT)*
	2	CS*AS*CBS+2.*H(INT)*K(INT)*AS*BS*CGS
				RO(INT)=SQRT(ROSAVE)
				IF (RO(INT).LE.RHOMAX) GO TO 93
				INT=INT-1
93			CONTINUE
92		CONTINUE
91	CONTINUE
1002	CONTINUE
	INT1=INT-1
	INT2=INT
	DO 94 J=1,INT1
		INT2=INT2-1
		DO 95 I=1,INT2
			IF (RO(I+1).GE.RO(I)) GO TO 95
			Z=RO(I+1)
			RO(I+1)=RO(I)
			RO(I)=Z
			IZ=K(I+1)
			K(I+1)=K(I)
			K(I)=IZ
			IZ=H(I+1)
			H(I+1)=H(I)
			H(I)=IZ
			IZ=L(I+1)
			L(I+1)=L(I)
			L(I)=IZ
95		CONTINUE
94	CONTINUE
	GO TO 1003


      END
	SUBROUTINE YESNO(IKEY)
C	IF SYM IS 'Y' THAN IKEY IS 1 OTHERWISE IKEY=0
	CHARACTER *1 SYM
	IKEY=0
	READ(5,1) SYM
1	FORMAT(A)
	IF (SYM.EQ.'Y'.OR.SYM.EQ.'y') iKEY=1
	RETURN
	END
