module sogrot_mod3e

  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
  private
  public :: fftvar_init
  public :: sogrot_exe, sog_fft, sog_ifft, rotsou_geodip, rotrec_geodip
  
  real,    parameter 			:: pi=3.1415926535898
  complex, parameter 			:: i = cmplx(0.,1.)
  complex,pointer  			:: tfield(:,:,:)
  integer            			:: nx, nh, nz,nzs(3)
  real               			:: dx, dh, dz
  real               			:: ox, oh, oz
  type(axis)         			:: kh, kx, kz
  real,pointer                 		:: kxv(:),khv(:),kzv(:)

  integer, parameter   			::  FFTW_R2HC=0,FTW_HC2R=1,FFTW_DHT=2,FFTW_REDFT00=3,FFTW_REDFT01=4
  integer, parameter   			::  FFTW_REDFT10=5,FFTW_REDFT11=6,FFTW_RODFT00=7,FFTW_RODFT01=8,FFTW_RODFT10=9
  integer, parameter   			::  FFTW_RODFT11=10,FFTW_FORWARD=-1,FFTW_BACKWARD=+1,FFTW_DESTROY_INPUT=1
  integer, parameter   			::  FFTW_UNALIGNED=2,FFTW_CONSERVE_MEMORY=4,FFTW_EXHAUSTIVE=8,FTW_PRESERVE_INPUT=16
  integer, parameter   			::  FFTW_PATIENT=32,FFTW_ESTIMATE=64,FFTW_ESTIMATE_PATIENT=128
  integer, parameter   			::  FFTW_BELIEVE_PCOST=256,FFTW_DFT_R2HC_ICKY=512,FFTW_NONTHREADED_ICKY=1024
  integer, parameter   			::  FFTW_NO_BUFFERING=2048,FFTW_NO_INDIRECT_OP=4096,FFTW_ALLOW_LARGE_GENERIC=8192
  integer, parameter   			::  FFTW_NO_RANK_SPLITS=16384,FFTW_NO_VRANK_SPLITS=32768,FFTW_NO_VRECURSE=65536
 
!  integer*8, pointer,private   		::  planF,planI
  integer*8, private	   		::  planF,planI
!  integer*8, target,private   		::  tgt_plan
  

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

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

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

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

    stat=1
!!!
    kh%n = nh
    kh%d = 2*pi/(kh%n*dh)
!!!
    kx%n = nx
    kx%d = 2*pi/(kx%n*dx)
!!!
    kz%n = nz
    kz%d = 2*pi/(kz%n*dz)
!!!
    nzs=(/nz,nh,nx/)
    allocate(tfield(nzs(1),nzs(2),nzs(3)))
 write(0,*) "pre sfftw_plan"
    call sfftw_plan_dft(planF,3,nzs,tfield,tfield,fftW_FORWARD,fftW_PATIENT)
 write(0,*) "passei sfftw_plan"
    call sfftw_plan_dft(planI,3,nzs,tfield,tfield,fftw_BACKWARD,fftW_PATIENT)
 write(0,*) "passei sfftw_plan"
    allocate(kxv(nx),kzv(nz),khv(nh))
    call klinemap(kxv,x%n,x%d,0)
    call klinemap(kzv,z%n,z%d,0)
    call klinemap(khv,h%n,h%d,0)
 write(0,*) "passei klinemap"
 write(0,*) kzv,khv,kxv
!!!

stat = 0

  end function fftvar_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function sogrot_exe(Ccig,souwav) result(cig_out)

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

 allocate(Rcig(kz%n,kh%n,kx%n),cig_out(nz,nh,nx))

 Rcig=0.;cig_out=0.

 if (souwav) then
    Rcig = rotsou_geodip(Ccig)
 else
    Rcig = rotrec_geodip(Ccig)
 end if
 cig_out = sog_ifft(Rcig)

 deallocate (Rcig)
 end function sogrot_exe
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function sog_fft(cig) result(Ccig)

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

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

 tfield=cmplx(cig,0.)
 call sfftw_execute(planF)

 Ccig=tfield*sqrt(1./(1.*nzs(1)*nzs(2)*nzs(3)));

 end function sog_fft

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function sog_ifft(Ccig) result(cig) 

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

 allocate(cig(nz,nh,nx))

 tfield=Ccig;
 call sfftw_execute(planI)
 cig=tfield*sqrt(1./(1.*nzs(1)*nzs(2)*nzs(3)));

 end function sog_ifft
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function rotsou_geodip(Ccig) result(Scig)

 complex				:: Ccig(:,:,:)
 complex, allocatable			:: Scig(:,:,:)
 integer				:: ix, iz, ih, stat, ikhrot,k,ib,ie
 real					:: tgama,talfa,tteta_s,gwrot,fwrot,khrot

 allocate(Scig(kz%n,kh%n,kx%n));Scig=0.

k=0
 do ix=1,kx%n
    do iz=1,kz%n
       if(kzv(iz).eq.0.) then
          if(kxv(ix).lt.0.) then
             talfa=1.E+37
          else if(kxv(ix).gt.0.) then
             talfa=-1.E+37
          else
             talfa=-kx%d/kz%d
          end if
       else
          talfa=-kxv(ix)/kzv(iz)
       end if

       do ih=1,kh%n
          if(kzv(iz).eq.0.) then
!             if(khv(ih).lt.0.) then
!                tteta_s=-1.E+37
!             else if(khv(ih).gt.0.) then
!                tteta_s=-1.E+37
!	     else
!	        tteta_s=-kh%d/kz%d
!             end if
                 khrot=sign(khv(ih),-kxv(ix))
          else
             tteta_s=-khv(ih)/kzv(iz)
  	     if (abs(1.+talfa*tteta_s).lt.(1.E-37)) then
                 khrot=sign(khv(ih),-kxv(ix))
!	        if(tteta_s-talfa.lt.0.) then
!	          tgama=-1.E+37
!	        else if(tteta_s-talfa.gt.0.) then
!	       tgama=1.E+37
!	     else
!	       tgama=0.
!	     end if
	  else
             tgama=(tteta_s-talfa)/(1.+talfa*tteta_s)
	     khrot=-kzv(iz)*tgama
	  end if

          end if


	  ib=nint(khrot/kh%d)+1
	  ib=int(khrot/kh%d)+1
	  ie=ib+1
	  if(ib<0) ib=ib+kh%n
	  if(ie<0) ie=ie+kh%n

!          if(ib<0) ib=ib+kh%n
!          ie=ib+1

  	  if (ib.ge.1.and.ib.le.kh%n.and.ie.ge.1.and.ie.le.kh%n) then
!             if (khrot < 0) then
                fwrot=1-(khrot-khv(ib))/kh%d
!             else
!                fwrot=khrot/kh%d-ib
!	     end if
             gwrot=1.-fwrot
	     k=k+1
  	     Scig(iz,ih,ix)=gwrot*Ccig(iz,ib,ix)+fwrot*Ccig(iz,ie,ix)+Scig(iz,ih,ix)
	  end if

       end do
    end do
 end do
 write(0,*) k

 end function rotsou_geodip

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function rotrec_geodip(Ccig) result(Gcig)

 complex				:: Ccig(:,:,:)
 complex, allocatable			:: Gcig(:,:,:)
 integer				:: ix, iz, ih, stat, ikhrot,k,ib,ie
 real					:: tgama,talfa,tteta_g,gwrot,fwrot,khrot

 allocate(Gcig(kz%n,kh%n,kx%n));Gcig=0.

k=0
 do ix=1,kx%n
    do iz=1,kz%n
       if(kzv(iz).eq.0.) then
          if(kxv(ix).lt.0.) then
             talfa=1.E+37
          else if(kxv(ix).gt.0.) then
             talfa=-1.E+37
          else
             talfa=-kx%d/kz%d
          end if
       else
          talfa=-kxv(ix)/kzv(iz)
       end if

       do ih=1,kh%n
          if(kzv(iz).eq.0.) then
             if(khv(ih).lt.0.) then
                tteta_g=-1.E+37
             else if(khv(ih).gt.0.) then
                tteta_g=1.E+37
             else
		tteta_g=kh%d/kz%d
             end if
          else
             tteta_g=khv(ih)/kzv(iz)
          end if

	  if (abs(1.+talfa*tteta_g).lt.(1.E-37)) then
	     if((talfa-tteta_g).lt.0.) then
	       tgama=-1.E+37
	     else if((talfa-tteta_g).gt.0.) then
	       tgama=1.E+37
	     else
	       tgama=0.
	     end if
	  else
             tgama=(talfa-tteta_g)/(1.+talfa*tteta_g)
	  end if

	  khrot=-kzv(iz)*tgama
	  ib=int(khrot/kh%d)+1
	  ie=ib+1
	  if(ib<0) ib=ib+kh%n
	  if(ie<0) ie=ie+kh%n

  	  if (ib.ge.1.and.ib.le.kh%n.and.ie.ge.1.and.ie.le.kh%n) then
	     k=k+1
             fwrot=1.-(khrot-khv(ib))/kh%d
             gwrot=1.-fwrot
  	     Gcig(iz,ih,ix)=gwrot*Ccig(iz,ib,ix)+fwrot*Ccig(iz,ie,ix)+Gcig(iz,ih,ix)
	  end if

       end do
    end do
 end do
 write(0,*) k

 end function rotrec_geodip
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
   subroutine klinemap(k,n,d,shift)
    real, dimension(:), pointer :: k
    integer, intent(in) :: n,shift
    real,    intent(in) :: d
    integer             :: i,m
#ifdef DEBUG
    write(0,*) "IN KLINEMAP"
#endif
    k=0.
    if(n/=1) then
       m = n/2
       do i=1,m
          k(i+shift  ) =    i-1
          k(i+m+shift) = -m+i-1
       end do

 
!       do i=1,n
!          if(i <= m) then
!            k(i+shift) = i -1
!          else
!            k(i+shift) = i - n -1
!          end if
!       end do


    else
       k=0. !! make sure k=0 if n=1
    end if

    k = k / (n*d)      !! 1/length
    k = 2 * atan(1.)*4 * k !! radians
#ifdef DEBUG
    write(0,*) "OUT KLINEMAP"
#endif
  end subroutine klinemap
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module sogrot_mod3e
