	SUBROUTINE SIMUL8(NCURVE,W,TTHETA,WL1,WL2,NFILE,NREFL,NPT 
	1	,FIRST,END,D2T,SCMAX,Y)
C
C   PROGRAM TO SIMULATE A DIFFRACTOMETER CHART 
C	 FROM AN ARRAY OF 2-THETA POSITIONS AND INTENSITIES
C	 READ FROM SCRATCH FILE ON LOGICAL UNIT 'NFILE'
C
C  MODIFIED FROM POWD7 'CCHART' SUBROUTINE   MAY,1982  M.F.GARBAUSKAS
C
	LOGICAL DTOTT,FLAG

	DIMENSION Y(npt),TTH(2),RINT(2)
	DIMENSION TABLE1(15),TABLE2(15),TABLE(15),XANGLE(15)
	DIMENSION TABLE3(15) 
C
C     INITIALIZE VARIABLES
C 
	DATA PI,RAD,ZERO/3.141593,114.5916,0.0/,RINT/1.,0.5/
C 
	DATA TABLE1/0.060,0.066,0.075,0.089,0.108,0.130,0.160,0.175,0.205, 
	1	0.225,0.250,0.295,0.360,0.410,0.500/   
C    
	DATA TABLE2/0.070,0.075,0.088,0.105,0.130,0.185,0.260,0.315,0.390, 
	1	0.450,0.535,0.750,1.200,1.500,2.000/
C 
C     NBS MEASURED FWHM TABLE 
	DATA TABLE3/0.12,0.12,0.12,0.125,0.13,0.135,0.155,0.185,0.23, 
	1	0.255,0.285,0.315,0.36,0.41,0.5/

	DATA XANGLE/0.0,20.0,40.0,60.0,80.0,100.0,120.0,130.0,140.0,145.0,
	1	150.0,155.0,160.0,162.5,165.0/
C
C     IF W=0.0 USE NBS MEASURED FWHM TABLE 
      IF (W .LE. 0.01) THEN
	W = 0.12
	TTHETA = 40.0 
	DO I=1,15 
	  TABLE(I) = TABLE3(I) 
	ENDDO
      ELSE
	DO 220  I=1,15    
	  IF (XANGLE(I).LT.TTHETA) GO TO 220  
	  ITEMP=I-1  
	  XPART=TTHETA-XANGLE(ITEMP)
	  YPART=XANGLE(I)-XANGLE(ITEMP)  
	  FRAC=XPART/YPART   
	  HWU=(FRAC*(TABLE2(I)-TABLE2(ITEMP)))+TABLE2(ITEMP) 
	  HWL=(FRAC*(TABLE1(I)-TABLE1(ITEMP)))+TABLE1(ITEMP)
	  DIF=(W-HWL)/(HWU-HWL)
	  GO TO 230    
220	CONTINUE
230	DO J=1,15
	  TABLE(J)=TABLE1(J)+DIF*(TABLE2(J)-TABLE1(J))
	ENDDO
      ENDIF
C
	NWV = 2
	IF (WL2 .EQ. WL1 .OR. WL2 .EQ. 0) NWV = 1

300   DO 440  II=1,NREFL
	K=2
C  READ POSITIONS AND INTENSITIES FROM SCRATCH FILE
	READ(NFILE,END=450) dsp,XINT1 	! dspace, INTENSITY

	FLAG = DTOTT(dsp,TTH(1),WL1)
	IF (.not. FLAG) GOTO 440
	TTH(2) = TTH(1)
	IF (NWV .eq. 2) FLAG = DTOTT(dsp,TTH(2),WL2)
C
C   CHECK FOR REFLECTION RANGE FOR HALF WIDTH DETERMINATION
C
340	IF (XANGLE(K).LE.TTH(1)) THEN 
	  K=K+1  
	  GO TO 340    
	ENDIF
	KTEMP=K-1 
	DO 430 I=1,NWV
	  IF (TTH(I).LE.0) GO TO 430
	  XPART=TTH(I)-XANGLE(KTEMP) 
	  YPART=XANGLE(K)-XANGLE(KTEMP)  
	  WREFL=(XPART/YPART)*(TABLE(K)-TABLE(KTEMP))+TABLE(KTEMP)     
	  WSQ=WREFL*WREFL   
	  IF (NCURVE.EQ.0) THEN 
	    ZREFL = 2.0*RINT(I)*XINT1*WREFL/PI
	  ELSE
	    ZREFL = RINT(I)*XINT1/WREFL
	  ENDIF
C
C  SET UP LIMITS FOR CALCULATION 
C
	
C 5 DEGREES ON EITHER SIDE OF REFLECTION
	  I2T = (TTH(I)-FIRST)/D2T + 1.
	  KSTART=I2T-(5.0/D2T)
	  KSTOP=I2T+(5.0/D2T)
	  IF (KSTART.LT.1) KSTART=1
	  KSTART = MAX(1,KSTART)
	  KSTOP = min(npt,KSTOP)
C
C     CALCULATE PROFILE FOR A REFLECTION    
C
	  DO J=KSTART,KSTOP     
	    X=J-1
	    DEL = TTH(I)-(FIRST+X*D2T) 
	    DELSQ = DEL*DEL 
	    IF (NCURVE.EQ.0) THEN 
C     CAUCHY PROFILE FOR A REFLECTION    
		Y(J) = Y(J) + (ZREFL/(WSQ+4.0*DELSQ))   
	    ELSE
C     GAUSSIAN PROFILE FOR A REFLECTION     
		Y(J) = Y(J)+(ZREFL*EXP(-PI*DELSQ/WSQ))     
	    ENDIF
	  ENDDO
430	CONTINUE
440   CONTINUE
C
C     SCALE THE CALCULATED PATTERN
C
450	XX=Y(1)
	DO J=2,NPT
	  XX = MAX(XX,Y(J))
	ENDDO
	IF (XX .le. 0) return
	SCALE=SCMAX/XX
	DO J=1,NPT
	  Y(J)=Y(J)*SCALE
	ENDDO
      RETURN
      END
