	SUBROUTINE SIMULATE(files,numf,status)
	INTEGER numf,files(*),status
C   SIMULATES DIFFRACTOMETER TRACE FROM JCPDS CARD DATA
C   
C     WRITTEN MAY,1982  M.F.GARBAUSKAS
C     MODIFIED FOR INVERTED CU 2-THETA STRUCTURE  AUGUST, 1982  MFG
C	IMPROVED AND MODIFIED TO WRITE DUPONT DATA FILES 8/29/86 RLH
C	Major adaptations 12/86-4/87 BHT
C	updated to use new JCPDS format 12/87 BHT
C
C   SIMULATION ADAPTED FROM POWD7 'CCHART' SUBROUINE
C
	include 'cmpr.common'

	INTEGER CLI$GET_VALUE,CLI$PRESENT
	EXTERNAL CLI$_PRESENT,CLI$_DEFAULTED,CLI$_ABSENT,CLI$_NEGATED
	LOGICAL GET_NUMBER,flag,GET_FNUMBER

	CHARACTER*1	IANS
        real CMULT(10)

	LOGICAL DTOTT,TTTOD
	CHARACTER*80	newext
	LOGICAL newflag
	character target*2

	LOGICAL*1 WAVEF,NOISEF,peakf 
C
	STATUS = .false.
	
	IF (numf .eq. 1 .and. NPOINT(files(1)) .le. 0) THEN
	  newfil = files(1)
	  newflag = .false.
	ELSEIF (MFILE .lt. M$FIL) THEN
	  newfil = MFILE + 1
C Clear arrays of peak information
	  CALL ZEROPK(newfil,1,m$pks)
	  newflag = .true.
	ELSE
	  write(lo,'(A)') ' Unable to create a new file'
	  return
	ENDIF
C
C  Clear out the arrays
C
	DO I=1,m$pts
	  BESD(I,newfil)=0.0
	  B(I,newfil)=0.0
	  BX(I,newfil) = 0.0
	ENDDO

	NCURVE = 0
	IWDEF = 4
	peakf = .false.
	NOISEF = .false. 
	SCMAX = 5000.
	W = 0.0
	NBKG = 0
C......................................................................
C............. Parse the qualifiers from the command string............

C WAVELENGTH qualifier
	wavnew = 0.0
	flag = .true.
	IF (CLI$PRESENT('WAVELENGTH')) 
	1	FLAG = GET_FNUMBER('WAVELENGTH',wav1new)
	IF (.not. FLAG) THEN 
	  WAVEF = .true.
	  iw = 1
	  IRAD(newfil) = 'SY'
	  wavnew=wav1new
	  wav2new=wav1new
	ELSEIF (CLI$PRESENT('TARGET') .or. CLI$PRESENT('WTYPE')) THEN 
	  IRAD(newfil)  = 'CU'
C TARGET qualifier
	  IF (CLI$PRESENT('TARGET')) THEN 
	    flag = CLI$GET_VALUE('TARGET',target,ln)
	    read(target(1:ln),'(a2)') IRAD(newfil)
	    WAVEF = .true.
	  ENDIF
C WTYPE qualifier
	  IF (CLI$PRESENT('WTYPE')) THEN 
	    FLAG = GET_NUMBER('WTYPE',IW)
	    IF (IW .le. 0 .or. IW .gt. 4) IW = 0
	  ENDIF
	ELSEIF (DEFAULT .and. wave0 .ne. 0) THEN
	  wavnew = wave0
	  wav1new = wave0
	  wav2new = wave0
	  IW = 1
	ELSEIF (DEFAULT) THEN
34	  IW = 3
	  IRAD(newfil) = 'CU'
	ENDIF

C  GET K ALPHA1 AND K ALPHA2 WAVELENGTHS FOR TARGET
	IF (wavnew .eq. 0.0 .and. IRAD(newfil) .ne. '  ') THEN
	  CALL WAVESET(IRAD(newfil),wavnew,wav1new,wav2new,*35)
	ELSEIF (wavnew .eq. 0.0) THEN
35	  CALL WAVEINQ(LO,LI,IRAD(newfil),wavnew,wav1new,wav2new,*34)
	ENDIF

	IF (wav1new .EQ. wav2new) then 
	  IW = 1
	ELSEIF (IW .eq. 0) then 
273	  WRITE (LO,9319) IRAD(newfil),wav1new,
	1	IRAD(newfil),wav2new,IRAD(newfil),wavnew,
	2	IRAD(newfil),wav1new,wav2new,IWDEF
9319	  FORMAT ('$Choose wavelength(s) for simulated pattern',/,
	1	' 	1 for ',A2,' K alpha 1 (w=',F9.6,'A)',/,
	2	' 	2 for ',A2,' K alpha 2 (w=',F9.6,'A)',/,
	3	' 	3 for ',A2,' K alpha bar (weighted avg of 1+2, w=',
	4	F9.6,'A)',/,' 	4 for both ',A2,
	5	' K alpha 1 and K alpha 2 (w=',F9.6,'A, ',F9.6,'A)',/,
	6	'$Enter digit for choice [Default=',i1,']: ')
	  READ (LI,9351,err=273,end=999) IW
	  IF (IW .le. 0 .or. IW .gt. 4) IW = IWDEF
	ENDIF
	IF (IW .EQ. 1) THEN 
	  wavnew = wav1new
	  wav2new = wav1new
	ELSEIF (IW .EQ. 2) THEN 
	  wavnew = wav2new
	  wav1new = wav2new
	ELSEIF (IW .EQ. 3) THEN 
	  wav1new = wavnew
	  wav2new = wavnew
	ENDIF

C RANGE qualifier
	FLAG = .TRUE.
	IF (CLI$PRESENT('RANGE')) THEN 
	  FLAG = GET_FNUMBER('RANGE',rmin)
	  IF (.not. FLAG) FLAG = GET_FNUMBER('RANGE',RMAX)
	ENDIF
	
	IF (FLAG) THEN
C get the max and min d-space from the peak list(s)
	  dmin = dspace(1,files(1))
	  dmax = dspace(1,files(1))
	  DO n=1,numf
	    nfil = files(n)
	    DO j=1,NPKS(nfil)
	      IF (PKHT(j,nfil) .gt. 0) THEN
	        dmin = MIN(dmin,dspace(j,nfil))
	        dmax = max(dmax,dspace(j,nfil))
	      ENDIF
	    ENDDO
	  ENDDO
	  FLAG = DTOTT(DMIN,RMAX,wav1new)
	  FLAG = DTOTT(DMAX,RMIN,wav2new)
C expand to the nearest degree allow at least half a degree on either side
	  rmin = NINT(rmin - 1.)
	  rmax = NINT(rmax + 1.)
	  IF (.not. DEFAULT) THEN
321	    WRITE(LO,9310) 'Minimum',rmin
9310	    FORMAT('$Enter ',a,' 2-theta values [Default=',f8.3,
	1	' degrees]: ')
	    READ(LI,'(q,F10.0)',END=999,err=321) N,x1
	    IF (n .gt. 0) rmin = x1
	    WRITE(LO,9310) 'Maximum',rmax
	    READ(LI,'(q,F10.0)',END=999,err=321) N,x1
	    IF (n .gt. 0) rmax = x1
	  ENDIF
	ENDIF

	FLAG = TTTOD(DMIN,RMAX,wav1new)
	FLAG = TTTOD(DMAX,RMIN,wav1new)

C STEP qualifier
	delr = 0.01
	IF (CLI$PRESENT('STEP')) THEN 
	  FLAG = GET_FNUMBER('STEP',DELR)
	ELSEIF (.not. DEFAULT) THEN
310	  WRITE(LO,9312) delr
9312	  FORMAT('$Enter stepsize [Default=',f7.3,' degrees]: ')
	  READ(LI,'(q,F10.0)',END=999,err=321) N,x1
	  IF (n .gt. 0) delr = x1
	  IF (delr .EQ. 0.0) delr = DELTA(files(1))
	  IF (delr .EQ. 0.0) goto 310
	ENDIF
	npts = (rmax-rmin)/delr+1
	IF (npts.GT.M$PTS) THEN 
	  WRITE(LO,9318) rmin,rmax,delr,npts,(rmax-rmin)/(M$PTS-1.) 
9318	  FORMAT(' *** Too many points: From',f6.1,' to',f6.1,
	1	' degrees with step of',f6.3,/,
	2	' *** requires',i8,
	3	' points. Minimum allowed step size is',f8.4)
	  DO WHILE (2*npts .GT. M$PTS) 
	    delr = delr * 2.
	    npts = (rmax-rmin)/delr+1
	  ENDDO
	  IF (.not. DEFAULT) GOTO 310
	ENDIF

C get the weights
	IF (CLI$PRESENT('WEIGHT')) THEN 
	  DO i=1,numf
	     FLAG = GET_FNUMBER('WEIGHT',CMULT(i))
	     IF (FLAG) CMULT(i) = 100.
	  ENDDO
	ELSEIF (DEFAULT .or. numf .eq. 1) THEN
	  DO i=1,numf
	     CMULT(i) = 100.
	  ENDDO
	ELSE
	  DO i=1,numf
210	     write (lo,'(A,i3,A)') '$Enter the weight for file',
	1	files(i),'[Default=100.]: '
	     READ (lo,'(F10.0)',err=210,end=998) CMULT(i)
	     IF (CMULT(i) .le. 0) CMULT(i) = 100.0
	  ENDDO
	ENDIF

	npkt = 0
	FLAG = .FALSE.
	DO n=1,numf
	  nfil = files(n)
	  npk = 0
	  DO j=1,NPKS(nfil)
	    IF (dspace(j,nfil).gt.dmin .and. dspace(j,nfil).le.dmax
	1	.and. PKHT(j,nfil) .gt. 0 .and. CMULT(n) .gt. 0)  
	1	npk = npk + 1
	  ENDDO
	  IF (cmult(N) .le. 0.) THEN
	    WRITE(lo,'(A,i3)') ' Multiplier for file',nfil,
	1	' is zero -- phase skipped'
	    FLAG = .true.
	  ELSEIF (npk .eq. 0) THEN
	    WRITE(lo,'(A,i3)') ' No peaks in range for file',nfil
	    FLAG = .true.
	  ENDIF
	  npkt = npkt + npk
	ENDDO
	IF (Npkt .EQ. 0) THEN
	  WRITE (LO,412) rmin,rmax
412	  FORMAT (' There are no peaks within the specified '
	1  '2theta range:'/,1x,f8.2,' to',f8.2,
	1  ', Reenter range or exit (^Z).')
	  IF (DEFAULT) GOTO 998
	  GOTO 321
	ELSEIF (FLAG .and. .not. DEFAULT) THEN
	  WRITE (lo,'(A)') '$Proceed with simulation? ([Y]/N): '
	  READ (li,'(A)',end=998) ians
	  IF (ians .eq. 'N' .or. ians .eq. 'n') goto 321
	ENDIF

C MAXIMUM qualifier
	IF (CLI$PRESENT('MAXIMUM') .eq. %LOC(CLI$_PRESENT)) THEN 
	  FLAG = GET_NUMBER('MAXIMUM',I)
	  IF (I .gt. 0 .and. .not. FLAG) SCMAX = I
	ELSEIF (.not. DEFAULT) THEN
274	  WRITE(LO,9320)
9320	  FORMAT('$ENTER MAXIMUM INTENSITY FOR PLOT [Default=5000.]: ')
	  READ(LI,'(F10.0)',END=999,err=274) SCMAX
	  IF (SCMAX.EQ.0.) SCMAX=5000.
	ENDIF

C CAUCHY qualifier
	IF (CLI$PRESENT('CAUCHY')) THEN 
	  NCURVE = 0
	  W = 0.
	  FLAG = GET_FNUMBER('CAUCHY',TTHETA)
	  IF (.not. FLAG) FLAG = GET_FNUMBER('CAUCHY',W)
C GAUSSIAN qualifier
	ELSEIF (CLI$PRESENT('GAUSSIAN')) THEN 
	  NCURVE = 1
	  W = 0.
	  FLAG = GET_FNUMBER('GAUSSIAN',TTHETA)
	  IF (.not. FLAG) FLAG = GET_FNUMBER('GAUSSIAN',W)
	ELSEIF (.not. DEFAULT) THEN
271	  WRITE(LO,9350)
9350	  FORMAT(' ENTER PROFILE: 0 for CAUCHY,  1 for GAUSSIAN',
	1	/,'$[Default=0 (CAUCHY)]: ')
	  READ(LI,9351,END=999,err=271) NCURVE
9351	  FORMAT(I)
272	  WRITE(LO,9355)
9355	  FORMAT(' Enter a peak half-width and the 2-theta value ',
	1	'for this width',/,'$[Default=0 (use NBS widths)]: ')
	  READ(LI,9356,END=999,err=272) W,TTHETA
9356	  FORMAT(2F8.0)
	  IF (TTHETA .eq. 0) W = 0.0
	ENDIF

C NOISE qualifier
	NOISEF = .false. 
	IF (CLI$PRESENT('NOISE')) THEN 
	  NOISEF = .true. 
	  SIG = 2.
	  FLAG = GET_FNUMBER('NOISE',SIG)
	ELSEIF (.not. DEFAULT) THEN
	  WRITE(LO,9500)
9500	  FORMAT('$ADD NOISE? (Y/[N]): ')
	  READ(LI,9510)IANS
9510	  FORMAT(A1)
	  IF (IANS.EQ.'Y' .or. IANS.EQ.'y') THEN 
	    NOISEF = .true. 
277	    WRITE(LO,9520)
9520	    FORMAT('$  NOISE STANDARD DEVIATION LEVEL [Default=2.0]: ')
	    READ(LI,'(F10.0)',END=999,err=277) SIG
	    IF (SIG.EQ.0.0) SIG = 2.
	  ENDIF
	ENDIF

C AMORPHOUS qualifier
	IF (CLI$PRESENT('AMORPHOUS') .eq. %LOC(CLI$_PRESENT)) THEN 
	  NBKG = 2
	  PEAK = 25.
	  PWID = 30.
	  PINT = 2.
	  FLAG = GET_FNUMBER('AMORPHOUS',PEAK)
	  IF (.not. FLAG) FLAG = GET_FNUMBER('AMORPHOUS',PWID)
	  IF (.not. FLAG) FLAG = GET_FNUMBER('AMORPHOUS',PINT)
C OFFSET qualifier
	ELSEIF (CLI$PRESENT('OFFSET') .eq. %LOC(CLI$_PRESENT)) THEN 
	  NBKG = 1
	  CONST = 10.
	  FLAG = GET_NUMBER('OFFSET',I)
	  IF (I .gt. 0 .and. .not. FLAG) CONST=I
	ELSEIF (.not. DEFAULT) THEN 
	  WRITE(LO,9400)
9400	  FORMAT('$ADD BACKGROUND? (Y/[N]): ')
	  READ(LI,9510,END=999)IANS
	  IF (IANS.EQ.'Y' .or. IANS.EQ.'y') THEN
	    WRITE(LO,9410)
9410	    FORMAT('$  USE AMORPHOUS PEAK? (Y/[N]): ')
	    READ(LI,9510,END=999) IANS
	    IF (IANS.EQ.'Y' .or. IANS.EQ.'y') THEN 
	      NBKG = 2
275	      WRITE(LO,9430)
9430	      FORMAT('   ENTER PARAMETERS FOR AMORPHOUS COMPONENT:')
	      WRITE(LO,9432)
9432	      FORMAT('$     PEAK POSITION (Degrees)  [Default=25.0]: ')
	      READ(LI,'(F10.0)',END=999,err=275) PEAK
	      IF(PEAK.EQ.0.)PEAK=25.
	      WRITE(LO,9434)
9434	      FORMAT('$     PEAK WIDTH (Degrees)  [Default=30.0]: ')
	      READ(LI,'(F10.0)',END=999,err=275) PWID
	      IF (PWID.EQ.0.0)PWID=30.
	      WRITE(LO,9436)
9436	      FORMAT('$     FRACTION OF INTENSITY  [Default=2.0]: ')
	      READ(LI,'(F10.0)',END=999,err=275) PINT
	      IF(PINT.EQ.0.0)PINT=2.
	    ELSE
276	      NBKG = 1
	      WRITE(LO,9420)
9420		FORMAT('$  NUMBER OF BACKGROUND COUNTS DESIRED [Default=10]: ')
	      READ(LI,'(F10.0)',END=999,err=276) CONST
	      IF(CONST.EQ.0.0)CONST=10.
	    ENDIF
	  ENDIF
	ENDIF
C----------------------------------------------------------------------
C----------------------------------------------------------------------
	OPEN(UNIT=out,TYPE='SCRATCH',FORM='UNFORMATTED',disp='delete')

	npkt = 0
	DO n=1,numf
	  nfil = files(n)
	  npk = 0
	  DO j=1,NPKS(nfil)
	    IF (dspace(j,nfil).gt.dmin .and. dspace(j,nfil).le.dmax
	1	.and. PKHT(j,nfil) .gt. 0 .and. CMULT(n) .gt. 0) THEN
	    WRITE (out) dspace(j,nfil),CMULT(n) * PKHT(j,nfil)
	    npk = npk + 1
	    ENDIF
	  ENDDO
	  npkt = npkt + npk
	ENDDO
	IF (Npkt .EQ. 0) THEN
	  WRITE (LO,412) rmin,rmax
	  close(UNIT=out,disp='delete')
	  GOTO 998
	ENDIF
C----------------------------------------------------------------------------
C----------------------------------------------------------------------------
C  Convert the list of peaks to a diffractogram
C
	REWIND out
	CALL SIMUL8(NCURVE,W,TTHETA,wav1new,wav2new,out,npkt,
	1	npts,rmin,rmax,delr,SCMAX,B(1,newfil))
C add background
      IF (NBKG.NE.0) CALL BACKGROUND(NBKG,CONST,PEAK,PWID,PINT,
	1	B(1,newfil),npts,rmin,delr)
C add noise
      IF (NOISEF) CALL NOISE(npts,B(1,newfil),SIG)
C rescale to SCMAX
      CALL SCALE(rmax,rmin,delr,SCMAX,B(1,newfil),npts)
C simulate ESD's
	DO I=1,npts
	  BX(I,newfil) = (I-1)*DELR + RMIN
	  BESD(I,newfil)=SQRT(B(I,newfil))
	ENDDO
C----------------------------------------------------------------------------
	REWIND out
	NPEAKS = 0
	J = M$PKS / 2
	IF (wav1new .eq. wav2new) J = M$PKS
	DO I = 1,J
	  READ (out,end=501) dsp,xint
	  NPEAKS = NPEAKS + 1
	  FLAG = DTOTT(dsp,TT(NPEAKS,newfil),wav1new)
	  RELINT(NPEAKS,newfil)=XINT/100.
	  PKHT(NPEAKS,newfil) = XINT*(SCMAX-CONST)/10000.
	  BKG(NPEAKS,newfil) = CONST
	  DSPACE(NPEAKS,newfil)=dsp
	  IF (wav1new .eq. wav2new) THEN
	    IWAVE(NPEAKS,newfil) = 0
	  ELSE
	    IWAVE(NPEAKS,newfil) = 1
 	    NPEAKS = NPEAKS + 1
	    FLAG = DTOTT(dsp,TT(NPEAKS,newfil),wav2new)
	    RELINT(NPEAKS,newfil)=XINT/200.
	    PKHT(NPEAKS,newfil) = XINT * (SCMAX-CONST) / 20000.
	    BKG(NPEAKS,newfil) = CONST
	    DSPACE(NPEAKS,newfil)=dsp
	    IWAVE(NPEAKS,newfil) = 2
	  ENDIF
	ENDDO
	READ (out,end=501) dsp,xint
	WRITE (LO,'(2A)') 'Warning: simulated pattern has too many ',
	1	'peaks for the peak list'
501	CLOSE(UNIT=out,disp='delete')
C----------------------------------------------------------------------------
	MONO(newfil)= .false.
	wave(newfil)=wavnew
	WKA1(newfil)=wav1new
	WKA2(newfil)=wav2new
	IF (IW .EQ. 1) THEN 
	    MONO(newfil)= .true.
	    wave(newfil)=wav1new
	    WKA2(newfil)=wav1new
	ELSEIF (IW .EQ. 2) THEN 
	    MONO(newfil)= .true.
	    wave(newfil)=wav2new
	    WKA1(newfil)=wav2new
	ELSEIF (IW .EQ. 3) THEN 
	    MONO(newfil)= .true.
	    WKA1(newfil)=wavnew
	    WKA2(newfil)=wavnew
	ENDIF
	CALL IDATE(DATE(1,newfil),DATE(2,newfil),DATE(3,newfil))
	NPKS(newfil) = NPEAKS
	NPOINT(newfil) = NPTS
	INST(newfil) = 'CR'
	SOURCE(newfil)='JC'
	START(newfil) = RMIN
	FINISH(newfil)= RMAX
	DELTA(newfil) = DELR
	TIME(newfil)= 0.0
	TEMP1(newfil)= 0
	TEMP2(newfil)= 0.
	REFINT(newfil)=1.0
	NBKGP(newfil) = 0
	FXTYP(newfil) = 1
	FYTYP(newfil) = 1
	FEXISTS(newfil)=.false.

C create file name
	acfnam = newext(CFILNAM(files(1)), '.SIM')
	IQ = 0
	IF (.not. DEFAULT) THEN
	  ln = lench(acfnam)
	  WRITE(LO,9000) acfnam(:ln)
9000	  FORMAT('$Enter name for the new file [DEFAULT=',a,']: ')
	  READ(LI,'(Q,A)',END=999) IQ,CFILNAM(newfil)
	ENDIF
	IF (IQ .eq. 0) CFILNAM(newfil) = acfnam

C create title
	ATITLE(newfil) = ' '
	DO n=1,numf
	  ln = lench(ATITLE(newfil))
	  IF (n .ne. 1) ATITLE(newfil) = ATITLE(newfil)(:ln) // ' + ' 
	  nfil = files(n)
	  ln = lench(ATITLE(newfil))
	  ln0 = MAX(1,lench(CFILNAM(nfil))-(80/numf))
	  write (ATITLE(newfil)(ln+1:),'(i5,2A)') NINT(cmult(nfil)),
	1		' * ',CFILNAM(nfil)(ln0:lench(CFILNAM(nfil)))
	ENDDO
	MFILE = MAX(MFILE,newfil)
	STATUS = .true.
	RETURN
C-----------------------------------------------------------------------------
998	Type *,'*** Unable to prompt for required input ***'
999	Type *,'*** Ignoring SIMULATE command ***'
3	CLOSE(UNIT=out,disp='delete')
	STATUS = .FALSE.
	RETURN 
	END
