Note
Function that returns the acf of an array.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=R8), | intent(in), | dimension(1:long, 1:larg) | :: | tab_in |
input array |
|
| real(kind=R8), | intent(out), | dimension(1:long, 1:larg) | :: | tab_out |
acf of the input array |
|
| integer(kind=I4), | intent(in) | :: | long |
2D array length |
||
| integer(kind=I4), | intent(in) | :: | larg |
2D array width |
||
| logical(kind=I4), | intent(in) | :: | sub_samp |
sampling? |
subroutine acv(tab_in, tab_out, long, larg, sub_samp) !================================================================================================ !< @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 ) :: long !! *2D array length* integer(kind=I4), intent(in ) :: larg !! *2D array width* logical(kind=I4), intent(in ) :: sub_samp !! *sampling?* real (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab_in !! *input array* real (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_out !! *acf of the input array* integer(kind=I4) :: nx2, ny2, iex, iey, ibx, iby, i, j, lo2, la2 real (kind=R8) :: tmp integer(kind=I4), dimension(1:2) :: loc_max complex(kind=R8), dimension(:,:), allocatable :: cmpl real (kind=R8), dimension(:,:), allocatable :: tab_ext1, tab_ext2 type(SCALE_SURF) :: scal_surf type(MOMENT_STAT) :: m_res ! 0-padding if ( PAD_FFT_ANI < 0 ) then nx2 = long ny2 = larg ibx = 1 ; iex = long iby = 1 ; iey = larg else nx2 = 2 * ( nint(PAD_FFT_ANI * long) / 2 ) ny2 = 2 * ( nint(PAD_FFT_ANI * larg) / 2 ) ibx = max( ceiling( (nx2 - long)/2. ), 1 ) ; iex = ibx + long - 1 iby = max( ceiling( (ny2 - larg)/2. ), 1 ) ; iey = iby + larg - 1 endif allocate( tab_ext1(1:nx2, 1:ny2), & ! tab_ext2(1:nx2, 1:ny2) ) ! allocate( cmpl(1:nx2, 1:ny2) ) ! tab_ext1(1:nx2, 1:ny2) = 0 call calc_moments( tab = tab_in(1:long, 1:larg), & ! mx = m_res, & ! nb_mom = 2 ) ! tab_ext1(ibx:iex, iby:iey) = ( tab_in(1:long, 1:larg) - m_res%mu ) / m_res%si if ( APO_FFT_ANI /= "no_apo" ) then call apod( tab_in = tab_ext1(1:nx2, 1:ny2), & ! tab_out = tab_ext2(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2, & ! type_apo = trim(APO_FFT_ANI), & ! param = 0.1_R8 ) ! else tab_ext2(1:nx2, 1:ny2) = tab_ext1(1:nx2, 1:ny2) endif !---------------- if (sub_samp) then call tab_calc_fftw3_real_fwd( tab_in = tab_ext2(1:nx2, 1:ny2), & ! tab_ou = cmpl(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2) ! else call calc_fftw3_real_fwd( tab_in = tab_ext2(1:nx2, 1:ny2), & ! tab_ou = cmpl(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2) ! endif cmpl(1:nx2, 1:ny2) = cmplx( abs( cmpl(1:nx2, 1:ny2) )**2, 0, kind = R8 ) ! théorème de wiener if (sub_samp) then call tab_calc_fftw3_real_bwd( tab_in = cmpl(1:nx2, 1:ny2), & ! tab_ou = tab_ext1(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2) ! else call calc_fftw3_real_bwd( tab_in = cmpl(1:nx2, 1:ny2), & ! tab_ou = tab_ext1(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2) ! endif call trans_corner2center( tab_in = tab_ext1(1:nx2, 1:ny2), & ! tab_out = tab_ext2(1:nx2, 1:ny2), & ! long = nx2, & ! larg = ny2 ) ! tab_out(1:long, 1:larg) = tab_ext2(ibx:iex, iby:iey) ! normalisation loc_max(1:2) = maxloc( tab_out(1:long, 1:larg) ) lo2 = loc_max(1) la2 = loc_max(2) tmp = tab_out(lo2, la2) tab_out(1:long, 1:larg) = tab_out(1:long, 1:larg) / tmp if (.false.) then call init_scal( scal = scal_surf, & ! nx = long, & ! ny = larg, & ! lx = 1._R8, & ! ly = 1._R8, & ! unit_z = 'm ') ! call write_surf( nom_fic = "test_acv.sur", & ! tab_s = tab_out(1:long, 1:larg), & ! scal = scal_surf ) ! endif deallocate(cmpl) deallocate(tab_ext1, tab_ext2) return endsubroutine acv