c************
c Description: Create a data waveform in a Midas file
c
c @author Jeff Schoen
c @version $Id: waveform.for,v 1.4 2002/07/19 17:51:03 schoenj Exp $
c  
c************
	program WAVEFORM
c
	include 'primitive.inc'
	record /HEADER/ HCB
c
	parameter (L_bbuf=32768)
	include 'buffers.inc'
	record /BUFFER/ R
c
	real*8 amp,freq,poff,phase,fnyq,chirp
	real*4 famp,fphasel,fpchirp,rx,delay
	complex*8 cx,cxd,cxc
	logical*4 pulse
	integer*4 throttle
	character*4 format
	character*3 type
	character*46 shapes 
     &	  /'SHAPES;Sinusoid,Square,Triangle,Sawtooth,Pulse'/
	equivalence (cx,rx)
c
	call M$CHECKIN
c
c find the format of waveform
	format = M$APICK(2)
	format = M$TEST_FORMAT(format,'S#,C#,-SP,-CP')
c
c find the type of waveform
	type = M$APICK(4)
	if (type.eq.'SIN') then	
	  itype = 1
	elseif (type.eq.'SQU') then	
	  itype = 2
	elseif (type.eq.'TRI') then	
	  itype = 3
	elseif (type.eq.'SAW') then	
	  itype = 4
	elseif (type.eq.'PUL') then	
	  itype = 5
	else
	  call M$ERROR ('Unsupported waveform type: '//type)
	endif
c
c fill out header and open ouput file
	call M$INITIALIZE (HCB)
	HCB.file_name =	M$APICK(1)
	HCB.type   = 1000
	HCB.format = format
	HCB.xstart = M$DPICK(8)
	HCB.xdelta = M$DPICK(9)
	HCB.xunits = 1
	HCB.size   = M$DPICK(3)
	call M$ALLOCATE (HCB)
c
c get waveform parameters
	amp  = M$DPICK(5)
	freq = M$DPICK(6)
	poff = M$DPICK(7) / 360.  ! internally in cycles
	chirp = 0.d0
	phase = 0.d0
	throttle = M$GET_SWITCH_DEF('RT',2)
c
c setup real time parameter controls
	fnyq = .5/HCB.xdelta
	idt = M$MWINIT ('WAVE SHAPE', itype, shapes)
	ida = M$DWINIT ('AMPLITUDE',  amp, 1.d0, -1.d0, 0.01*amp)
	idf = M$DWINIT ('FREQUENCY', freq, 1.d0, -1.d0, 0.001*fnyq)
	idp = M$DWINIT ('PHASE OFF', poff, 1.d0, -1.d0, 0.001d0)
	idc = M$DWINIT ('CHIRP', chirp, 1.d0, -1.d0, 0.001*fnyq)
	idr = M$MWINIT ('THROTTLE', throttle, 'RealTime,Full,Block')
c
c find buffer size
	LPER = M$XFER_LENGTH (8,L_bbuf)
	HCB.buf_type = 'F'
	pulse = poff .eq. 0.d0
c
c signal start of processing block
	call M$SYNC()
c
c loop on buffers depending on data type
	call M$TODO (HCB.size)
	do while ( M$DO(LPER,NN) )

	  idum = M$MWGET(idt,itype)
	  idum = M$DWGET(ida,amp)
	  idum = M$DWGET(idf,freq)
	  idum = M$DWGET(idp,poff)
	  idum = M$DWGET(idc,chirp)
	  idum = M$MWGET(idr,throttle)

	  do while (throttle.eq.3 .and. .not.Mc.break)
	    call M$PAUSE (0.1)
	    idum = M$MWGET(idr,throttle)
	  enddo

	  famp = amp
	  fphase = mod( phase+poff, 1.d0)
	  fpdelta = freq * HCB.xdelta
	  fpchirp = chirp * HCB.xdelta * HCB.xdelta

	  if (itype.eq.1) then				! sinusoid
	    fphase = fphase*twopi
	    fpdelta = fpdelta*twopi
	    fpchirp = fpchirp*twopi
	    cx  = amp * CMPLX( cos(fphase), sin(fphase) )
	    cxd = CMPLX( cos(fpdelta), sin(fpdelta) )
	    cxc = CMPLX( cos(fpchirp), sin(fpchirp) )
	    if (format(1:1).eq.'S') then
	      do i = 1,NN
		R.fbuf(i) = rx
		cx = cx * cxd
		cxd = cxd * cxc
	      enddo
	    else
	      do i = 1,NN
		R.cfbuf(i) = cx
		cx = cx * cxd
		cxd = cxd * cxc
	      enddo
	    endif

	  elseif (itype.eq.2) then			! square
	    do i = 1,NN
	      if (fphase.lt.0.5) then
		R.fbuf(i) = famp
	      else
		R.fbuf(i) = -famp
	      endif
	      fpdelta = fpdelta + fpchirp
	      fphase = fphase + fpdelta
	      if (fphase.ge.1.0) fphase = fphase - 1.0
	    enddo

	  elseif (itype.eq.3) then			! triangle
	    famp = famp * 4
	    do i = 1,NN
	      if (fphase.lt.0.25) then
		R.fbuf(i) = famp*fphase
	      elseif (fphase.lt.0.75) then
		R.fbuf(i) = famp*(.5-fphase)
	      else
		R.fbuf(i) = famp*(fphase-1.)
	      endif
	      fpdelta = fpdelta + fpchirp
	      fphase = fphase + fpdelta
	      if (fphase.ge.1.0) fphase = fphase - 1.0
	    enddo

	  elseif (itype.eq.4) then			! sawtooth
	    famp = famp * 2
	    do i = 1,NN
	      R.fbuf(i) = famp*(fphase-.5)
	      fpdelta = fpdelta + fpchirp
	      fphase = fphase + fpdelta
	      if (fphase.ge.1.0) fphase = fphase - 1.0
	    enddo

	  elseif (itype.eq.5) then			! pulse
	    do i = 1,NN
	      if (fphasel.gt.fphase) then
		R.fbuf(i) = famp
	      else
		R.fbuf(i) = 0.
	      endif
	      fphasel = fphase
	      fpdelta = fpdelta + fpchirp
	      fphase = fphase + fpdelta
	      if (fphase.ge.1.0) fphase = fphase - 1.0
	    enddo

	  endif

c	  complex of all but sinusoids
	  if (format(1:1).eq.'C' .and. itype.ne.1) then	
	    do i = NN,1,-1
	      R.cfbuf(i) = CMPLX( R.fbuf(i), R.fbuf(i) )
	    enddo
	  endif

	  call M$FILAD (HCB,R,NN)

	  if (throttle.ne.2) then
	    delay = HCB.xdelta*NN
	    call M$PAUSE(delay)
	  endif

	  phase = DMOD( phase+freq*HCB.xdelta*NN, 1.d0 )

	enddo
c
c close output file
	call M$CLOSE (HCB)
c
	call M$RETURN
	end

