program MDHspecfac_fast
  ! (MD) Spectral Factorization or (MD) Autocorrelation
  ! using Helix Transform.
  !
  ! Sjoerd de Ridder
  ! Stanford, January 21th, 2009.
  ! Stanford, January 28th, 2009.

  ! Helix is Wrapped and unwrapped as read from disk!
  !
  ! auto   = 0          [1] for MDH AC output
  ! tsym   = 0          [1] for return of positive and negative time-lags
  !                         ( [1]=default when auto=0 )
  use sep

  implicit none

  ! Include for FFTW
  include "fftw3.f" ! That is usually found in /usr/include on ubuntu
                    !    or in /opt/FFTW/include on 64bit machines
  
  integer                          :: n1, n2, n3, nl, nh, nhp
  real                             :: d1, d2, d3, o1, o2, o3
  integer                          :: ii, jj, kk, cent
  integer                          :: auto, tsym, ssym
  real   , allocatable             :: data(:)
  integer                          :: pad2

  ! For FFTW
  real   , allocatable             :: in(:)
  complex, allocatable             :: out(:)
  integer*8                        :: planf,plani
  
  call sep_init()
  
  ! Read from input History
  call from_history('n1',n1,1 ); call from_history('n2',n2,1 ); call from_history('n3',n3,1 )
  call from_history('d1',d1,0.); call from_history('d2',d2,0.); call from_history('d3',d3,0.)
  call from_history('o1',o1,0.); call from_history('o2',o2,0.); call from_history('o3',o3,0.)
  
  ! Read from input Parameters
  call from_par('auto'  , auto  , 0  )
  if ( auto.eq.1 ) then
	call from_par('tsym'  , tsym  , 1  )
  else
	call from_par('tsym'  , tsym  , 0  )
  end if
  call from_par('ssym'  , ssym  , 1  )

  ! Write to output History  
  if (ssym.eq.1) then
	call to_history('o1',-((n1-1)/2)*d1)
     	if (n3.ne.1) then
        	call to_history('o2',-((n2-1)/2)*d2)
	else
        	call to_history('o2',0.)
	end if
  else
	call to_history('o1',o1); call to_history('o2',0.)
  end if
  call to_history('n1',n1); call to_history('n2',n2); call to_history('n3',n3)
  if (tsym.eq.1) then
	if (n3.ne.1) then
        	call to_history('n3',2*n3); call to_history('o3',-((2*n3-1)/2)*d3)
	else if (n2.ne.1) then
        	call to_history('n2',2*n2); call to_history('o2',-((2*n2-1)/2)*d2)
	else
        	call to_history('n1',2*n1); call to_history('o1',-((2*n1-1)/2)*d1)
     	end if
  else
     if (n3.ne.1) then
        call to_history('n3',n3); call to_history('o3',o3)
     else if (n2.ne.1) then
        call to_history('n2',n2); call to_history('o2',o2)
     else
        call to_history('n1',n1); call to_history('o1',o1)
     end if
  end if
  call to_history('d1',d1); call to_history('d2',d2); call to_history('d3',d3)
  
  call sep_close()
  ! Determine Helix lengths
  nl = n1
  if (n3.ne.1) then; nl = nl*n2; end if
  nh   = n1*n2*n3
  nhp  = pad2(2*nh) ! times 2 to Padd with zeros for negative times
  if (n2.ne.1) then
	cent = n1/2-1     ! Find spatial center
  else
    	cent = 0     ! Find spatial center
  end if
  if (n3.ne.1) then; cent = cent + (n2/2-1)*n1; end if

  ! Print Helix length and padd size
  write(0,*) '-------------------------------------------------------------'
  write(0,*) 'Helix length: ',nh,' Padded total length: ',nhp
  
  ! Allocate Memory
  allocate( data(nh),in(nhp),out(nhp/2+1) )
  
  ! Initialize (I)FFTW plans
  call sfftw_plan_dft_r2c_1d(planf,nhp,in,out,FFTW_ESTIMATE)
  call sfftw_plan_dft_c2r_1d(plani,nhp,out,in,FFTW_ESTIMATE)
  
  ! Read data and center on a Helix
  write(0,*) 'Read Data ...'
  call sep_read(data)
  in(            1 : nl  - cent ) = data( cent + 1 : nl   )
  in(  nl-cent + 1 : nh  - cent ) = data(   nl + 1 : nh   )
  in(  nh-cent + 1 : nhp - cent ) = 0.
  in( nhp-cent + 1 : nhp        ) = data(        1 : cent )
  
  ! Perform Fourier Transform
  write(0,*) 'Call FT ...'
  
  in = in/sqrt( 1.*nhp) ! FFT weighting
  call sfftw_execute(planf)
  
  ! Compute Autocorrelation
  write(0,*) 'Compute AC ...'
  out  =  out*conjg(out)
  
  ! Kolmogoroff or IFFT
  if ( auto.eq.0 ) then
    write(0,*) 'Call Kolmogoroff ...'
    call kolmogoroff(planf,plani,nhp,in,out)
  else
    write(0,*) 'Call IFT ...'
    out = out/sqrt( 1.*nhp) ! IFFT weighting
    call sfftw_execute(plani)
  end if
  
  ! Center and write to output
  write(0,*) 'Write Data ...'
  if ( tsym.eq.1 ) then
    data(      1    : cent ) = in( (nhp-nh) + nhp - cent + 1 : nh + nhp        )
    data( 1  + cent : nl   ) = in( (nhp-nh) +              1 : nh +  nl - cent )
    data( nl + 1    : nh   ) = in( (nhp-nh) +  nl - cent + 1 : nh +  nh - cent )
    call sep_write( data )
    data(      1    : cent ) = in( nhp - cent + 1 : nhp        )
    data( 1  + cent : nl   ) = in(              1 :  nl - cent )
    data( nl + 1    : nh   ) = in(  nl - cent + 1 :  nh - cent )
    call sep_write( data )
  else
    data(      1    : cent ) = in( nhp - cent + 1 : nhp        )
    data( 1  + cent : nl   ) = in(              1 :  nl - cent )
    data( nl + 1    : nh   ) = in(  nl - cent + 1 :  nh - cent )
    call sep_write( data )
  end if
  
  write(0,*) '-------------------------------------------------------------'
  
  !! Destroy (I)FFTW plans
  call sfftw_destroy_plan(planf);
  call sfftw_destroy_plan(plani);
  
  ! Deallocate 
  deallocate( data,in,out )
  
  call exit(0)

end program MDHspecfac_fast


subroutine kolmogoroff( planf, plani, n, in,out)  ! Spectral factorization (from pvi)
  integer   n
  real      in (n)
  complex   out(n/2+1)
  integer*8 plani,planf

   out = clog(out)
   out = out/sqrt( 1.*n)  ! IFFTW weighting
   call sfftw_execute(plani)
   in            = in*sqrt(1./n)
   in(1)         = in(1) / 2.
   in(1+n/2)     = in(1+n/2) / 2.
   in(1+n/2+1:n) = 0.
   in = in/sqrt( 1.*n)    ! FFTW weighting
   call sfftw_execute(planf)
   out = out*sqrt(1.*n)
   out = cexp( out)
   out = out/sqrt( 1.*n)  ! IFFTW weighting
   call sfftw_execute(plani)
   in = in*sqrt(1./n)

end subroutine kolmogoroff


integer function pad2(n)
  integer n
  pad2=1
  do
     if (pad2.ge.n) exit
     pad2=pad2*2
  end do
end function pad2
