c************
c  Description:  Create a Midas file containing a sine wave
c
c  Author:	JGS   4//90
c
c  $Revision: 1.1 $
c  $Modtime:   29 Nov 1993 18:28:58  $
c
c  HISTORY:
c     DR 16125 JSLee with HP:
c        Fixes errors when output file format is of type Double.
c        The output format is set using the file qualifier (ft=d).
c        More specifics:
c          Before this fix, the following command cause the "cos"
c          output file to be SD. The format of the file "cexp" is CF.
c          The "cos" output file's waveform is incorrect.
c
c             SINCOS,sin(ft=d),cos,100,1000,0,4096,0,1/8000,cexp
c
c          Before this fix, the following command causes Signal 11 crash.
c             SINCOS,sin(ft=d),cos(ft=d),100,1000,0,4096,0,1/8000,cexp(ft=d)
c
c************
	program SINCOSINE
c
	include 'primitive.inc'
	record /HEADER/ HCBS,HCBC,HCBX
c
	parameter (LSEC=128,L_fbuf=4096)
c  ORIGINAL CODE: DR 16125
c 	real*4    bo(LSEC*2),rx(L_fbuf),ix(L_fbuf),rphase
c  NEW CODE: DR 16125
	real*4    bo(LSEC*2),rx(L_fbuf*2),ix(L_fbuf*2),rphase
c
	real*8    freq,chirp,delta,phase,dphase,fnyq,flow,freql,ephase
c  COMMENT: DR 16125
c       In case complex output is 'D', made cx array size allocation.
c       This fixes the Signal 11 crash if the complex ouput is
c       specified as cmpOutFile(FT=D).
c  ORIGINAL CODE: DR 16125
c	complex*8 bw(LSEC),cx(L_fbuf),cz
c  NEW CODE: DR 16125
	complex*8 bw(LSEC),cz
	complex*8 cx(L_fbuf*2)
c
	complex*16 dcx,dcz
	logical*4  exact
c
	call M$CHECKIN
c
c fill out header and open output file
	call M$INITIALIZE (HCBS)
	HCBS.type   = 1000
	HCBS.format = 'SF'
	HCBS.size   = DNINT( M$DPICK(6) )
	HCBS.xstart = M$DPICK(7)
	HCBS.xdelta = M$DPICK(8)
	HCBS.xunits = 1
c
c output files ?
	HCBS.file_name = M$APICK(1)
	HCBS.switch = 'OPTIONAL'
	call M$ALLOCATE (HCBS)

	call M$PROPAGATE (HCBS,HCBC)
	HCBC.file_name = M$APICK(2)
c  COMMENT: DR 16125
c       The change fixes the following problem:
c         If sinfile(ft=d), cos file data format also became,
c         even if "ft=d" not used, scalar double.
c         The trouble was that the cosine wave did not look right.
c
c  NEW CODE: DR 16125
	HCBC.format = 'SF'          !jslee added.
c
	HCBC.switch = 'OPTIONAL'
	call M$ALLOCATE (HCBC)


	call M$PROPAGATE (HCBS,HCBX)
	HCBX.file_name = M$APICK(9)
	HCBX.format = 'CF'
	HCBX.switch = 'OPTIONAL'
	call M$ALLOCATE (HCBX)
c
c check for degenerate output case
	if (.not.(HCBS.open .or. HCBC.open .or. HCBX.open))
     &	  call M$WARNING ('No output files specified')
c
c check for EXACT case for testing purposes
	exact = M$GET_PSWITCH('EXACT')
c
c get freq, amplitude and phase values
	amp = M$FPICK(3)
	freq = M$DPICK(4)
	phase = M$DPICK(5) / 360.d0
	delta = HCBS.xdelta
c
c set real-time controls and signal start of processing block
	fnyq = .5/delta
	flow = 0
	if (HCBX.open) flow = -fnyq
	idf = M$DWINIT ('FREQUENCY', freq, flow, fnyq, 0.001d0*fnyq)
	ida = M$FWINIT ('AMPLITUDE',  amp, 1.0, 
     &                   -1.0, 0.01*amp)
c
c get buffer length and sync
	LPER = M$XFER_LENGTH (4,L_fbuf*4)  ! always float data
	freql = freq + 1.0d0  ! force 1st time generation of waveform
	call M$SYNC()
c
c loop through buffers calculating complex phasor
	call M$TODO (HCBS.size)
	do while ( M$DO(LPER,N) )

	  idum = M$FWGET(ida,amp) 
	  idum = M$DWGET(idf,freq) 
c
c create a new waveform buffer if necessary
	  dphase = freq*delta
	  if (freq .ne. freql .and. .not.exact) then
	    freql = freq
	    dcz = dCMPLX( dcos(dtwopi*dphase), dsin(dtwopi*dphase) )
	    dcx = dCMPLX( 1.d0, 0.d0)
	    do i = 1,LSEC
	      bw(i) = dcx
	      dcx = dcx * dcz  ! rotate phasor
	    enddo
	  endif
c
c loop through waveform buffer shifting and scaling
	  i = 0
	  do while (i.lt.N)

	    ntodo = min( LSEC,N-i )

	    if (exact) then
	      ephase = DMOD( phase+dphase*i, 1.d0 )
	      do ii = 1,ntodo*2,2
		bo(ii  ) = amp*dcos(ephase*dtwopi)
		bo(ii+1) = amp*dsin(ephase*dtwopi)
		ephase = ephase + dphase
	      enddo
	    else
	      rphase = dtwopi * DMOD( phase+dphase*i, 1.d0 )
	      cz = amp * CMPLX( cos(rphase), sin(rphase) )
	      call M$CVSMUL (bw,cz,bo,ntodo)
	    endif

	    if (HCBS.open) call M$VMOV ( bo(2), 2, ix(i+1), 1, ntodo)
	    if (HCBC.open) call M$VMOV ( bo(1), 2, rx(i+1), 1, ntodo)
	    if (HCBX.open) call M$VMOV ( bo(1), 1, cx(i+1), 1, ntodo*2)

	    i=i+ntodo
	  enddo
c
c write output buffers
	  call M$FILAD (HCBS,ix,N)
	  call M$FILAD (HCBC,rx,N)
	  call M$FILAD (HCBX,cx,N)
c
	  phase = DMOD( phase+dphase*N, 1.d0 )

	enddo
c
c close output files
	call M$CLOSE (HCBS)
	call M$CLOSE (HCBC)
	call M$CLOSE (HCBX)
c
	call M$RETURN
	end

