module off2ang2D
  !! OFF2ANG2D - 2D offset to angle transformation
  !! based on Paul Sava's OFF2ANG program

  !! real input/output
  !! input  is z-mx-hx
  !! output is z-mx-ax

  use sep
  use util_
  use adcig2D
  use odcig2D

  !! Global Variables
  type(myaxis) :: az,amx,ahx,aax
  integer, external,private   :: omp_get_num_threads,omp_get_thread_num
  real,private                :: amax,hmax,Aeps,pi
  logical ,private            :: verb,push,inverse
  integer,private             :: nzpad,nzold,esize,nth,ith,stat
  complex,allocatable ,private:: cdata(:,:,:)

contains

  subroutine off2ang2D_init(az_in,amx_in,ahx,aax,verb_in,inverse_in)
    type(myaxis)         :: az_in,amx_in,ahx,aax
    logical :: verb_in,inverse_in
    call from_param("Aeps" , Aeps, 1.)
    call from_param("push" , push,.true.)
    call from_param("nzpad", nzpad,az_in%n)
    call from_param("hmax" , hmax, 1.)
    call from_param("amax" , amax, 60.)

    verb = verb_in
    inverse = inverse_in
    pi = acos(-1.)
    nzold = az%n

    az%n = nzpad
    az%o = az_in%o
    az%d = az_in%d

    amx%n=amx_in%n
    amx%o=amx_in%o
    amx%d=amx_in%d

    !$OMP PARALLEL
    nth = omp_get_num_threads()
    !$OMP END PARALLEL

    if(verb) write(0,*) " push,nth,nzpad=",push,nth,nzpad
    if(verb) write(0,*) " Aeps,hmax,amax",Aeps,hmax,amax

    if(inverse) then !! DEFINE MAX OFFSET PARAMETERS		
       ahx%n = aax%n   
       ahx%d = 1.
       if(aax%o ==0.) then
          ahx%o = 0.
          ahx%d =   hmax/ahx%n
       else
          ahx%o =  -hmax
          ahx%d = 2.*hmax/ahx%n
       end if

       if(verb) write(0,*) "requested max hx =",hmax

    else !! CALCULATE MAXIMUM ANGLES
       aax%n = ahx%n
       if(ahx%o ==0.) then
          aax%o = 0.
          aax%d =   amax/aax%n
       else
          aax%o =  -amax
          aax%d = 2.*amax/aax%n
       end if
       if(verb)  write(0,*) "requested max ax =",amax
       if(verb)  write(0,*) "believe  only ax <",180*atan(az%d/ahx%d)/pi
       if(verb)  write(0,*) "num ahx%n ",ahx%n

    end if

    if(inverse) then
       allocate(cdata(nzpad,aax%n,amx%n), stat=stat)
       if(stat /=0) call seperr("cant allocate cdata")
    else
       allocate(cdata(nzpad,ahx%n,amx%n), stat=stat)
       if(stat /=0) call seperr("cant allocate cdata")
    end if

    !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ih,ix,iz)
    do ix=1,amx%n
       do ih=1,ahx%n
          do iz=1,az%n
             cdata(iz,ih,ix) = cmplx(0.,0.)
          end do
       end do
    end do
    !$OMP END PARALLEL DO

  end subroutine off2ang2D_init

  !----------------------------------------------------------------

  subroutine off2ang2D_run(Rimg,ahx,aax)
    type(myaxis)  :: ahx,aax
    real, dimension(:,:,:) :: Rimg
    integer :: ix,iz,ih

    if(inverse) then  !! ANGLE TO OFFSET
       call odcig2D_fft_init(nth,(/az%n,ahx%n/))
       call odcig2D_shift_init(az,ahx)

       !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ih,ix,iz)
       do ih =1,ahx%n
          do ix=1,amx%n
             do iz=1,az%n
                cdata(iz,ih,ix) = cmplx( Rimg(ix,iz,ih) )
             end do
          end do
       end do
       !$OMP END PARALLEL DO

       !$OMP  PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ith,ix)
       do ix=1,amx%n
          ith = omp_get_thread_num()+1
          call odcig2D_F_fft(ith,ix,cdata)
          call odcig2D_F_shift(  ix,cdata)
       end do
       !$OMP END  PARALLEL DO

       call odcig2D_off_init(nth,az,ahx,aax,Aeps,push)

       !$OMP  PARALLEL DO  SCHEDULE(DYNAMIC) PRIVATE(ith,ix)
       do ix=1,amx%n
          ith = omp_get_thread_num()+1
          call odcig2D_off_x(ith,ix,cdata,az,ahx,aax)
       end do
       !$OMP END PARALLEL DO

       call odcig2D_off_close()

       !$OMP  PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ith,ix)
       do ix=1,amx%n
          ith = omp_get_thread_num()+1
          call odcig2D_I_shift(  ix,cdata)
          call odcig2D_I_fft(ith,ix,cdata)
       end do
       !$OMP END PARALLEL DO

       !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ih,ix,iz)
       do ih=1,ahx%n
          do ix=1,amx%n
             do iz=1,az%n
                Rimg(ix,iz,ih) = real (cdata(iz,ih,ix))	
             end do
          end do
       end do
       !$OMP END PARALLEL DO

       call odcig2D_shift_close()
       call odcig2D_fft_close()

    else !! OFFSET TO ANGLE

       call adcig2D_fft_init(nth,(/az%n,aax%n/))
       call adcig2D_shift_init(az,ahx)

       !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ih,ix,iz)
       do ix=1,amx%n
          do ih =1,ahx%n
             do iz=1,az%n
                cdata(iz,ih,ix) = cmplx( Rimg(ix,iz,ih) )
             end do
          end do
       end do
       !$OMP END PARALLEL DO

       do ix=1,amx%n
          call adcig2D_F_fft  (1,ix,cdata)
          call adcig2D_F_shift(  ix,cdata)
       end do
       write(0,*) 'INITIALIZING WITH az',az%n,az%d,az%o
       call adcig2D_ang_init(nth,az,ahx,aax,Aeps,push)
       !$OMP  PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ith,ix)
       do ix=1,amx%n
          ith = omp_get_thread_num()+1
          call adcig2D_ang_x(ith,ix,cdata,az,ahx,aax)
       end do
       !$OMP END PARALLEL DO
       call adcig2D_ang_close()

       do ix=1,amx%n
          ith = omp_get_thread_num()+1
          call adcig2D_I_shift(  ix,cdata)
          call adcig2D_I_fft(ith,ix,cdata)
       end do

       !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(ih,ix,iz)
       do ih=1,ahx%n
          do ix=1,amx%n
             do iz=1,az%n
                Rimg(ix,iz,ih) = real (cdata(iz,ih,ix))	
             end do
          end do
       end do
       !$OMP END PARALLEL DO

       call adcig2D_shift_close()
       call adcig2D_fft_close()
    end if

  end subroutine off2ang2D_run

  !----------------------------------------------------------------

  subroutine off2ang2D_close()
    deallocate(cdata, stat=stat)
    if(stat /=0) call seperr("cant deallocate cdata")
  end subroutine off2ang2D_close

end module off2ang2D


