module sod2adcig_mod

  use fft_data_mod
  implicit none

  type axis
     integer 					:: n
     real    					:: o, d, v ! v=(i-1)*d+o
     real, allocatable				:: all(:)
  end type axis

  public :: axis
  public :: kh,kx,kz, gama
  private
  public :: fftvar_init
  public :: so2adcig_init, so2adcig_exe, adcig_fft, adcig_ifft, do_adcig
  
  real,    parameter 			:: pi=3.1415926535898
  complex, parameter 			:: i = cmplx(0.,1.)
  integer            			:: nx, nh, nz
  integer	     			:: npadx, npadz, npadh
  real               			:: dx, dh, dz
  real               			:: ox, oh, oz
  type(fft_type)     			:: fft_1d_idx, fft_1d_idz, fft_1d_idh, fft_2d_idxh, fft_2d_idzh,fft_3d_idzhx
  type(axis)         			:: kh, kx, kz, gama

  ! Bug danger! The common variables should be modified ONLY by init & clean. Or at worst, by a single subroutine only.

contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  function fftvar_init(gama1, x, z, h) result(stat)

    type(axis), intent(in)     :: x, z, h, gama1
    integer                :: ikh, ikx, ig, ikz, stat

    nx = x%n;    dx = x%d
    nz = z%n;    dz = z%d
    nh = h%n;    dh = h%d

    stat=1
!!!
    gama%n=gama1%n;gama%o=gama1%o;gama%d=gama1%d
    allocate(gama%all(gama%n)); gama%all=0.
    do ig = 1, gama%n
      gama%all(ig) = gama%o+(ig-1)*gama%d
    end do
!!!
    kh%n = nh+npadh
    kh%d = 2*pi/(kh%n*dh)
    kh%o = -(kh%n/2)*kh%d
    allocate(kh%all(kh%n)); kh%all=0.
    do ikh = 1, kh%n
      kh%all(ikh) = kh%o+(ikh-1)*kh%d
    end do
!!!
    kx%n = nx+npadx
    kx%d = 2*pi/(kx%n*dx)
    kx%o = -(kx%n/2)*kx%d
    allocate(kx%all(kx%n)); kx%all=0.
    do ikx = 1, kx%n
      kx%all(ikx) = kx%o+(ikx-1)*kx%d
    end do
!!!
    kz%n = nz+npadz
    kz%d = 2*pi/(kz%n*dz)
    kz%o = -(kz%n/2)*kz%d
    allocate(kz%all(kz%n)); kz%all=0.
    do ikz = 1, kz%n
      kz%all(ikz) = kz%o+(ikz-1)*kz%d
    end do
!!!
    fft_1d_idx = init_fft_1d(kx%n)
    fft_1d_idz = init_fft_1d(kz%n)
    fft_2d_idxh = init_fft_2d(kh%n,kx%n)
    fft_2d_idzh = init_fft_2d(kz%n,kh%n)
    fft_3d_idzhx = init_fft_3d(kz%n,kh%n,kx%n)

!    write(0,*) kz%n,kz%o,kz%d
!    write(0,*) kh%n,kh%o,kh%d
!    write(0,*) kx%n,kx%o,kx%d

!    write(0,*) kz%all
!    write(0,*) kh%all
!    write(0,*) kx%all

stat = 0

  end function fftvar_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine so2adcig_init(npadz_in,npadx_in,npadh_in)

 integer			:: npadz_in,npadx_in,npadh_in

 npadz=npadz_in; npadx=npadx_in; npadh=npadh_in

 end subroutine so2adcig_init 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function so2adcig_exe(Ccig) result(cig_out)

 complex, intent(in)			:: Ccig(:,:,:)
 real, allocatable			:: cig_out(:,:,:)
 complex,allocatable			:: Rcig(:,:,:)
 integer				:: ix, iz, ih

 allocate(Rcig(kz%n,gama%n,nx),cig_out(nz,gama%n,nx))

 Rcig=0.;cig_out=0.
 write(0,*) "passei exec00"
 Rcig = do_adcig(Ccig)
 write(0,*) "passei exec01"
 cig_out = adcig_ifft(Rcig)
 write(0,*) "passei exec02"

 deallocate (Rcig)
 end function so2adcig_exe
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function adcig_fft(cig) result(Ccig)

 complex, allocatable			:: Ccig(:,:,:)
 real 					:: cig(:,:,:)
 integer				:: ix, iz, ih

 allocate(Ccig(kz%n,kh%n,nx))

 do ix=1,nx
    do ih=1,nh
       do iz=1,nz
          cig(iz,ih,ix)=cig(iz,ih,ix)*(-1)**(ih+iz)
       end do
    end do
 end do

 Ccig=0.
 Ccig(1:nz,1:nh,1:nx)=cmplx(cig,0)

 do ix=1,nx
    call fft_2d_data(fft_2d_idzh,  1,  kz%n,kh%n,Ccig(:,:,ix))
 end do

 end function adcig_fft

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function adcig_ifft(Ccig) result(cig) 

 complex				:: Ccig(:,:,:)
 real, allocatable			:: cig(:,:,:)
 complex, allocatable			:: temp(:,:,:)
 integer				:: ix, iz, ih, ig

 allocate(cig(nz,gama%n,nx))
 write(0,*) "passei ifft00"
 allocate(temp(kz%n,gama%n,nx))
 write(0,*) "passei ifft00.5"
! temp=Ccig
 cig=0.
 write(0,*) "passei ifft01"
! do iz=1,kz%n
!    if (iz.lt.int(kz%n/2)+1) then
!       temp(iz+int(kz%n/2),:,:)=Ccig(iz,:,:)
!    else if (iz.ge.int(kz%n/2)+1) then
!       temp(iz-int(kz%n/2),:,:)=Ccig(iz,:,:)
!    end if
! end do
 write(0,*) "passei ifft02"
! Ccig=-temp

 deallocate(temp)

 do ix=1,nx
    do ig=1,gama%n
       call fft_1d_data(fft_1d_idz,  -1,  kz%n,Ccig(:,ig,ix))
   end do
 end do
 write(0,*) "passei ifft03"

 cig=real(Ccig(1:nz,1:gama%n,1:nx))

 end function adcig_ifft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function do_adcig(Ccig) result(Acig)

 complex				:: Ccig(:,:,:)
 complex, allocatable			:: Acig(:,:,:)
 integer				:: ig, ix, iz, ikh,k
 real					:: gw,fw, kh1

 allocate(Acig(kz%n,gama%n,nx));Acig=0.

 do ix=1,kx%n
    do ig=1,gama%n
       do iz=1,kz%n

	  kh1=-kz%all(iz)*tand(gama%all(ig))
	  ikh=int((kh1-kh%o)/kh%d)+1

	  fw=kh1-ikh-1
	  gw=1.-fw

  	  if (ikh.ge.1 .and. ikh.lt.kh%n) then
	     k=k+1
  	     Acig(iz,ig,ix)=Acig(iz,ig,ix)+gw*Ccig(iz,ikh,ix)+fw*Ccig(iz,ikh,ix)
	  end if
       end do
    end do
 end do
 write(0,*) k

 end function do_adcig
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module sod2adcig_mod

