program FT4d_
  !! 4-D FFT
  !! Paul Sava (paul@sep.stanford.edu)
  !! 12/6/03

  !! complex input/output
  !! input  is  a1- a2- a3- a4-a5
  !! output is ka1-ka2-ka3-ka4-a5

  use sep
  use util_
  use fftwflags

  integer, parameter :: dbl=selected_real_kind(10,50)
  integer(kind=dbl)  :: plan

  logical :: adj, verb
  integer :: stat
  real :: scale, pi, arg, sign

  type(myaxis):: a1,a2,a3,a4,a5
  integer     :: i1,i2,i3,i4,i5

  complex, dimension(:,:,:,:) :: datai,datao
  complex, dimension(:,:)     :: sft12,sft34
  logical :: shift12,shift34

  call sep_init()

  call from_param("adj",  adj,  .false.)
  call from_param("verb", verb, .false.)
  call from_param("shift12",shift12,.true. )
  call from_param("shift34",shift34,.true. )

  call hinaxis(a1,1);call houaxis(a1,1)
  call hinaxis(a2,2);call houaxis(a2,2)
  call hinaxis(a3,3);call houaxis(a3,3)
  call hinaxis(a4,4);call houaxis(a4,4)
  call hinaxis(a5,5);call houaxis(a5,5)
  if(verb) then
     call printaxis(a1)
     call printaxis(a2)
     call printaxis(a3)
     call printaxis(a4)
     call printaxis(a5)
  end if

  call sep_close()

  !----------------------------------------------------------------
  allocate(datai(a1%n,a2%n,a3%n,a4%n), stat=stat)
  if(stat /=0) call seperr("cant allocate datai")
  allocate(datao(a1%n,a2%n,a3%n,a4%n), stat=stat)
  if(stat /=0) call seperr("cant allocate datao")
  !----------------------------------------------------------------
  if(adj) then
     call fftwnd_f77_create_plan(plan,4,(/a1%n,a2%n,a3%n,a4%n/),FFTW_BACKWARD,FFTW_ESTIMATE)
     sign = +1.
  else
     call fftwnd_f77_create_plan(plan,4,(/a1%n,a2%n,a3%n,a4%n/),FFTW_FORWARD, FFTW_ESTIMATE)
     sign = -1.
  end if
  !----------------------------------------------------------------
  scale = 1./sqrt(1.*a1%n*a2%n*a3%n*a4%n)
  !----------------------------------------------------------------
  allocate(sft12(a1%n,a2%n), stat=stat)
  if(stat /=0) call seperr("cant allocate sft12")
  sft12=cmplx(1.,0.)
  allocate(sft34(a3%n,a4%n), stat=stat)
  if(stat /=0) call seperr("cant allocate sft34")
  sft34=cmplx(1.,0.)

  if(shift12) then
     pi = acos(-1.)

     if(a1%o/=0. .or. a2%o/=0.) then
        do i1=1,a1%n
           do i2=1,a2%n
              arg = 2*pi &
              * ( (i1-1.)/a1%n * a1%o/a1%d &
              +   (i2-1.)/a2%n * a2%o/a2%d )

              sft12(i1,i2) = &
              cmplx(cos(arg),(sign)*sin(arg))
           end do
        end do
     end if
  end if
  if(shift34) then
     pi = acos(-1.)

     if(a3%o/=0. .or. a4%o/=0.) then
        do i3=1,a3%n
           do i4=1,a4%n
              arg = 2*pi &
              * ( (i3-1.)/a3%n * a3%o/a3%d &
              +   (i4-1.)/a4%n * a4%o/a4%d )

              sft34(i3,i4) = &
              cmplx(cos(arg),(sign)*sin(arg))
           end do
        end do
     end if
  end if
  !----------------------------------------------------------------
  if(adj) then
     do i5=1,a5%n
        if(verb) write(0,*) i5,'/',a5%n
        call sep_read(datai)

        datai = datai * scale
        if(shift34) then
           do i1=1,a1%n
              do i2=1,a2%n
                 datai(i1,i2,:,:)=datai(i1,i2,:,:) * sft34 !! origin shift
              end do
           end do
        end if
        if(shift12) then
           do i3=1,a3%n
              do i4=1,a4%n
                 datai(:,:,i3,i4)=datai(:,:,i3,i4) * sft12 !! origin shift
              end do
           end do
        end if
        call fftwnd_f77_one(plan,datai,datao)

        call sep_write(datao)
     end do
  else
     do i5=1,a5%n
        if(verb) write(0,*) i5,'/',a5%n
        call sep_read(datai)

        call fftwnd_f77_one(plan,datai,datao)
        if(shift12) then
           do i3=1,a3%n
              do i4=1,a4%n
                 datao(:,:,i3,i4)=datao(:,:,i3,i4) * sft12 !! origin shift
              end do
           end do
        end if
        if(shift34) then
           do i1=1,a1%n
              do i2=1,a2%n
                 datao(i1,i2,:,:)=datao(i1,i2,:,:) * sft34 !! origin shift
              end do
           end do
        end if
        datao = datao * scale

        call sep_write(datao)
     end do
  end if
  !----------------------------------------------------------------
  call exit(0)
end program FT4d_







