program main

  use sep
  use wei_myfft_mod
  implicit none
  integer  :: n1, n2, nk1, nk2, esize, i1, i2
  real     :: o1, o2, d1, d2, ok1, ok2, dk1, dk2  
  real, allocatable    :: data(:,:)
  complex, allocatable :: cdata(:,:)  
  real, parameter      :: PI = 3.14159265359

  call initpar()

  call from_history("n1", n1); call from_history("d1", d1); call from_history("o1", o1)
  call from_history("n2", n2); call from_history("d2", d2); call from_history("o2", o2)

  esize = 8

  nk1 = n1
  nk2 = n2
  dk1 = 2*PI/(n1*d1)
  dk2 = 2*PI/(n2*d2)
  ok1 = -PI/d1
  ok2 = -PI/d2

  allocate(data(n1, n2), cdata(n1, n2))

  call sep_read(data, "in", esize=4)

  cdata = data

  call FFT2DX(cdata, 1)

  do i1=1, n1
     call fftshift(cdata(i1,:), n2)
  end do
  do i2=1, n2
     call fftshift(cdata(:,i2), n1)
  end do

  call to_history("n1", n1); call to_history("d1", dk1); call to_history("o1", ok1)
  call to_history("n2", n2); call to_history("d2", dk2); call to_history("o2", ok2)
  call to_history("esize", esize)
  call sep_write(cdata, "out", esize=8)

end program


subroutine fftshift(a,n)
  complex a(n),tmp
  n0=n/2
  do i=1,n0
     tmp=a(i)
     a(i)=a(n0+i)
     a(n0+i)=tmp
  enddo
  return
end subroutine