acf_wiener Subroutine

public subroutine acf_wiener(tab_in, tab_out, w, h, multi_fft)

Note

Function that returns the acf of an array.

Arguments

Type IntentOptional Attributes Name
real(kind=R8), intent(in), dimension(1:w, 1:h) :: tab_in

input array
real(kind=R8), intent(out), dimension(1:w, 1:h) :: tab_out

acf of the input array
integer(kind=I4), intent(in) :: w

2D array length
integer(kind=I4), intent(in) :: h

2D array width
logical(kind=I4), intent(in), optional :: multi_fft

run parallel acfs?

Calls

proc~~acf_wiener~~CallsGraph proc~acf_wiener acf_wiener calc_fftw3_real_bwd calc_fftw3_real_bwd proc~acf_wiener->calc_fftw3_real_bwd calc_fftw3_real_fwd calc_fftw3_real_fwd proc~acf_wiener->calc_fftw3_real_fwd tab_calc_fftw3_real_bwd tab_calc_fftw3_real_bwd proc~acf_wiener->tab_calc_fftw3_real_bwd tab_calc_fftw3_real_fwd tab_calc_fftw3_real_fwd proc~acf_wiener->tab_calc_fftw3_real_fwd trans_corner2center trans_corner2center proc~acf_wiener->trans_corner2center

Called by

proc~~acf_wiener~~CalledByGraph proc~acf_wiener acf_wiener proc~calc_acf calc_acf proc~calc_acf->proc~acf_wiener proc~sub_surf sub_surf proc~sub_surf->proc~acf_wiener proc~read_job read_job proc~read_job->proc~calc_acf proc~read_job->proc~sub_surf proc~prg_surf prg_surf proc~prg_surf->proc~read_job program~main main program~main->proc~prg_surf

Source Code

   subroutine acf_wiener(tab_in, tab_out, w, h, multi_fft)
   !================================================================================================
   !<@note Function that returns the *acf* of an array.
   !< \[
   !< \begin{align*}
   !<    acf(i,j) &= (z \ast z)(i,j) = \sum_{k,l}^{n,n} z(k+1-i,l+1-j)z(k,l)  \\
   !<    TF(acf)  &= ACF = Z \cdot Z                                          \\
   !<    acf      &= TF^{-1}(ACF) = TF^{-1}(Z^2)
   !< \end{align*}
   !< \]
   !<
   !<@endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in )                      :: w           !! *2D array length*
   integer(kind=I4), intent(in )                      :: h           !! *2D array width*
   real   (kind=R8), intent(in ), dimension(1:w, 1:h) :: tab_in      !! *input array*
   real   (kind=R8), intent(out), dimension(1:w, 1:h) :: tab_out     !! *acf of the input array*
   logical(kind=I4), intent(in ), optional            :: multi_fft   !! *run parallel acfs?*

      integer(kind=I4) :: lo2, la2

      real   (kind=R8) :: tmp

      logical(kind=I4) :: parallel_fft

      integer(kind=I4), dimension(1:2) :: loc_max

      complex(kind=R8), dimension(:,:), allocatable :: tab_cmpl
      real   (kind=R8), dimension(:,:), allocatable :: tab_real

      allocate( tab_cmpl(1:w, 1:h) )

      allocate( tab_real(1:w, 1:h) )

      ! check for simultaneous fftw calculations
      !.........................................
      parallel_fft = .false.

      if ( present(multi_fft) ) parallel_fft = multi_fft
      !.........................................

      ! DFFT real -> complex
      !.........................................
      if ( parallel_fft ) then

         call tab_calc_fftw3_real_fwd( tab_in = tab_in  (1:w, 1:h),     &  ! IN
                                       tab_ou = tab_cmpl(1:w, 1:h),     &  ! OUT
                                         long = w,                      &  ! IN
                                         larg = h )                        ! IN

      else

         call calc_fftw3_real_fwd(      tab_in = tab_in  (1:w, 1:h),    &  ! IN
                                        tab_ou = tab_cmpl(1:w, 1:h),    &  ! OUT
                                          long = w,                     &  ! IN
                                          larg = h,                     &  ! IN
                                  planner_flag = FFTW_MEASURE )            ! IN

      endif
      !.........................................

      tab_cmpl(1:w, 1:h) = cmplx( abs( tab_cmpl(1:w, 1:h) )**2, 0, kind = R8 )


      ! IFFT complex -> real
      !.........................................
     if ( parallel_fft ) then

         call tab_calc_fftw3_real_bwd( tab_in = tab_cmpl(1:w, 1:h),     &  ! IN
                                       tab_ou = tab_real(1:w, 1:h),     &  ! OUT
                                         long = w,                      &  ! IN
                                         larg = h )                        ! IN
      else

         call calc_fftw3_real_bwd(      tab_in = tab_cmpl(1:w, 1:h),    &  ! IN
                                        tab_ou = tab_real(1:w, 1:h),    &  ! OUT
                                          long = w,                     &  ! IN
                                          larg = h,                     &  ! IN
                                  planner_flag = FFTW_MEASURE )            ! IN

      endif
      !.........................................

      ! the maximum is placed in the array center
      !.........................................
      call trans_corner2center(  tab_in  = tab_real(1:w, 1:h),  &  ! IN
                                 tab_out = tab_out (1:w, 1:h),  &  ! OUT
                                 long    = w,                   &  ! IN
                                 larg    = h  )                    ! IN
      !.........................................


      ! the maximum is 1
      !.........................................
      loc_max(1:2) = maxloc( tab_out(1:w, 1:h) )
      lo2 = loc_max(1)
      la2 = loc_max(2)

      tmp = tab_out(lo2, la2)

      tab_out(1:w, 1:h) = tab_out(1:w, 1:h) / tmp
      !.........................................

      deallocate(tab_cmpl)
      deallocate(tab_real)

   return
   endsubroutine acf_wiener