c************
c  Description:	Host version of MFFT routine
c
c  Author:	JGS   4/18/90
c
c  $Revision:  6 Sep 2001 DSL Added support for type 2000 input files
c              22 Oct 2002 DSL Corrected calculation of output file size
c              26 Jun 2020 Added GFC
c  $Modtime:   20 Mar 1996 10:48:10  $
c  
c DR 11859 merged in DSL's file size calculation fix
c************
	program MFFT
c
	implicit none
	include 'primitive.inc'
        integer*4	nfft, nin, nout, nrot, ngot, navg, modn, nconsume
	integer*4	idw, ida, idc, iw, isok, ind, mode, i
	real*8	  	delx, delf, davg, cellf, dcelf, toffset, overlap
	real*4		sum, rnavg, factor
	logical*4	ac, is_complex, t2000
	character*1	coupling
	character*20	wname
	complex*8	csum
c
	real*4    sig(1)
	complex*8 csig(1)
	pointer (psig,sig)
	pointer (pcsig,csig)

        real*4	win(1), asig(1)
	pointer (pwin,win), (pasig,asig)
c
	record /HEADER/ HIN, HOUT, HCELL

c
	call M$CHECKIN
c
c open and test signal input file #1
	call M$INIT (HIN, M$APICK(1), '1000,2000', 'S#,C#', 0)
	call M$OPEN (HIN, HCBF_INPUT)
	is_complex = HIN.format(1:1) .eq. 'C'
	t2000 = int(HIN.type/1000) .eq. 2
	call M$FORCE1000 (HIN)
c
c warn of loss of precision with DP data
	if (HIN.format(2:2) .eq. 'D') call M$WARNING
     &	  ('Double precision input data will be reformatted to float')

c
c get the fftsize/packing and clip
	nfft = M$LPICK(3)
c	override with frame size for type 2000 input
	if(t2000)nfft=hin.subsize
	delx = HIN.xdelta
	delf = 1.d0/(delx*nfft)
	if (is_complex) then
	  nout = nfft
	  factor = 1./(nfft)  ! fft scaling in window
	  nrot = nfft/2-1
	  if(m$get_switch('NOROTATE').gt.0)nrot=0
	else
	  mode = max(1,min(3,M$GET_SWITCH ('FFTP')))	! packing switch
	  nout = nfft/2
	  if (mode.eq.1) nout = nfft/2 + 1
	  factor = 1./(2*nfft)  ! fft scaling in window
	  nrot = 0
	endif
c
c get the window
	wname = M$APICK(4)
	iw = M$WINDOW_ID(wname)
c
c get input zero padding parameter
	nin = M$GET_SWITCH('TL')
	if (nin.le.0) nin = nfft
	if (nin.gt.nfft) call M$ERROR ('Xfer length must be <= FFT')
c
c get overlap factor and test
	overlap = M$DPICK(5)
c	force zero overlap for type 2000 input
	if(t2000)overlap=0
	if (overlap.ge.0.0d0) then
	  nconsume = max( 1, min( nin,NINT((1.-overlap)*nin) ))
	else
	  nconsume = max( nin, min( 4*nin, NINT(nin*(-overlap)) ))   ! corrected in 3.9.3
	endif
c
c get coupling mode
	coupling = M$APICK(6)
	ac = coupling .eq.'A'
c
c get # ffts to average
	navg = M$LPICK(7)
	navg = max(1,navg)
c
c setup and allocate output fft file
	call M$PROPAGATE (HIN,HOUT)
	HOUT.file_name = M$APICK(2)
	HOUT.type = 2000
	HOUT.format = 'CF'
	HOUT.xstart = -delf*nrot
	HOUT.xdelta = delf
	HOUT.xunits = 3
	HOUT.subsize = nout
	HOUT.ystart = HIN.xstart + delx*( (nin+nconsume*(navg-1))/2 )
	if(t2000)then
	 HOUT.ystart = HOUT.ystart + HIN.ystart
	 HOUT.ydelta = HIN.ydelta * navg
	else
	 HOUT.ydelta = delx*nconsume*navg
	 HOUT.yunits = 1
	endif
c+++
c DR 11859
c Old
c	HOUT.size = DINT((HIN.size-nin+nconsume)/(nconsume*navg*nfft))
c New
	HOUT.size = DINT((HIN.size-nin+nconsume)/nconsume/navg)
c---
	call M$OPEN (HOUT, HCBF_OUTPUT+HCBF_OPTIONAL)
	call M$FORCE1000 (HOUT)
c
c setup and allocate cell output file
	call M$PROPAGATE (HIN,HCELL)
	HCELL.file_name = M$APICK(8)
	HCELL.xstart = HIN.xstart + delx*( HIN.offset + nin/2 )
	HCELL.xdelta = delx*nconsume
	HCELL.xunits = HIN.xunits
	HCELL.subsize = 1
	HCELL.type = 1000
	HCELL.format = 'CF'
	HCELL.size = HOUT.size*navg
	call M$OPEN (HCELL, HCBF_OUTPUT+HCBF_OPTIONAL)
c
c check for degenerate output case
	if (.not. (HOUT.open .or. HCELL.open) ) 
     &	  call M$WARNING('No output files specified')
c
c set cell parameter
	cellf = M$DPICK(9)
c
c set xfer/cons lengths for dataflow grabs
	HIN.xfer_len = nin
	HIN.cons_len = nconsume
	HIN.buf_type = 'F'
c
c get dynamic buffers
	call M$MALLOC (nin*4, pwin)
	call M$MALLOC (max(nout*8,nin*HIN.bpe), psig)
	pcsig = psig
	call M$MALLOC (nout*8, pasig)
c
c setup controls and signal start of processing block
	idw = M$MWINIT ('WINDOW', iw, standard_windows)
c+++
c DR 2616
c	ida = M$LWINIT ('AVERAGING', navg, 1, 50, 1)
	ida = M$LWINIT ('AVERAGING', navg, 1, 262144, 1)
c---
	idc = M$DWINIT ('CELL FREQ', cellf, 0.d0, 1.d0/delx, delf)
c
c signal start of processing block
	call M$SYNC()
c
c make the window
	iw = M$WINDOW (win, wname, nin)
#if _MLIB
	call SSCAL (nin, factor, win, 1)
#else
	call M$VSMUL (win,factor,win,nin)
#endif
c
c loop through buffers
	modn = 0
	do while ( M$DO(nin,HIN.xfer_len) )

	  isok = M$LWGET(ida,navg)
	  modn = mod(modn+1,navg)
	  rnavg = 1./navg
c
c new window ?
	  if (M$MWGET(idw,iw) .ne. 0) then
	    call M$WINDOW_BY_ID (win, iw, nin)
#if _MLIB
	    call SSCAL (nin, factor, win, 1)
#else
	    call M$VSMUL (win,factor,win,nin)
#endif
	  endif
c
c FFT signal
	  if (HIN.open) then
	    call M$GRABX (HIN,sig,ngot)
	    if (ngot.ne.nin) goto 99 
	    if (is_complex) then
	      if (ac) call M$CVSUM (sig,csum,ngot)
	      if (ac) call M$CVSADD (sig,-csum/ngot,sig,ngot)
c	      call M$CRVMUL (sig,win,sig,ngot)
c	      call M$VFILL (sig(ngot*2+1),0.0,(nfft-ngot)*2)
	      do i=1,ngot
	       csig(i)=csig(i)*win(i)
	      enddo
	      do i=ngot+1,nfft
	       csig(i)=(0.0,0.0)
	      enddo
	      call M$CFFT (sig,nfft,1)
	    else
	      if (ac) call M$VSUM (sig,sum,ngot)
	      if (ac) call M$VSADD (sig,-sum/ngot,sig,ngot)
	      call M$VMUL (sig,win,sig,ngot)
	      call M$VFILL (sig(ngot+1),0.0,nfft-ngot)
c	      do i=1,ngot
c	       sig(i)=sig(i)*win(i)
c	      enddo
c	      do i=ngot+1,nfft
c	       sig(i)=0.0
c	      enddo
	      call M$RFFT (sig,nfft,mode)
	    endif
	  endif

c
c compute/output thinned cell
	  if (HCELL.open) then
	    isok = M$DWGET(idc,cellf)
	    ind = mod( nfft+2*NINT(cellf/delf), nfft) + 1
	    call M$FILAD (HCELL, sig(ind), 1)
	  endif
c
c compute/output averaged fft
	  if (HOUT.open) then
	    if (navg.eq.1) then
	      if(nrot.gt.0)then
	        call M$FILAD (HOUT,sig((nout-nrot)*2+1),nrot)
	        call M$FILAD (HOUT,sig,nout-nrot)
	      else
	        call M$FILAD (HOUT,sig,nout)
	      endif
	    else
	      if (modn.eq.1) then
		call M$VMOV (sig,1,asig,1,nout*2)
	      else
	        call M$VADD (sig,asig,asig,nout*2)
	      endif	
	      if (modn.eq.0) then
#if _MLIB
		call SSCAL (nout*2, rnavg, asig, 1)
#else
		call M$VSMUL (asig,rnavg,asig,nout*2)
#endif
	        call M$FILAD (HOUT,asig((nout-nrot)*2+1),nrot)
	        call M$FILAD (HOUT,asig,nout-nrot)
	      endif
	    endif
	  endif
c
  99	enddo
c
c close files
	call M$CLOSE (HIN)
	call M$CLOSE (HOUT)
	call M$CLOSE (HCELL)

c
	call M$RETURN
c
	end

