apod2 Subroutine

private subroutine apod2(tab_in, tab_out, long, larg, tau1, tau2, ang)

Note

Function that returns an apodized array.
To prevent gaps from appearing after FFT (because of non periodic waves), the surface must be transformed, but not too much. Here a modified Tukey window is determined. The starting surface is not modified below the “correlation lengths”. Above the correlation lengths, a smooth decrease is forced with a cosine squared.

Arguments

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

input acf
real(kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_out

apodized acf
integer(kind=I4), intent(in) :: long

surface acf length
integer(kind=I4), intent(in) :: larg

surface acf width
real(kind=R8), intent(in) :: tau1

surface first correlation length
real(kind=R8), intent(in) :: tau2

surface second correlation length
real(kind=R8), intent(in) :: ang

ellipsis angle

Called by

proc~~apod2~~CalledByGraph proc~apod2 apod2 proc~calc_imp_acf calc_imp_acf proc~calc_imp_acf->proc~apod2 proc~acf_theo acf_theo proc~acf_theo->proc~calc_imp_acf proc~sub_surf sub_surf proc~sub_surf->proc~calc_imp_acf proc~read_job read_job proc~read_job->proc~acf_theo 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 apod2(tab_in, tab_out, long, larg, tau1, tau2, ang)
   !================================================================================================
   !<@note Function that returns an apodized array.<br/>
   !< To prevent gaps from appearing after FFT (because of non periodic waves), the surface must
   !< be transformed, but not too much. Here a modified Tukey window is determined. The starting
   !< surface is not modified below the "correlation lengths". Above the correlation lengths, a
   !< smooth decrease is forced with a cosine squared.
   !<@endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in )                            :: long     !! *surface acf length*
   integer(kind=I4), intent(in )                            :: larg     !! *surface acf width*
   real   (kind=R8), intent(in )                            :: tau1     !! *surface first correlation length*
   real   (kind=R8), intent(in )                            :: tau2     !! *surface second correlation length*
   real   (kind=R8), intent(in )                            :: ang      !! *ellipsis angle*
   real   (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab_in   !! *input acf*
   real   (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_out  !! *apodized acf*

      real   (kind=R8) :: r2, c0, s0, rd, rr, theta, theta_diag, x, y, t, a_min, sum_inn, sum_tab, sum_int
      integer(kind=I4) :: i, j, k, long2, larg2, npt_out, n

      real(kind=R8), allocatable, dimension(:,:) :: tab_tmp

      ! first bissector angle (rad)
      theta_diag = atan2( PARAM%surf_height, PARAM%surf_width )

      ! sine and cosine of the ellipsis angle
      c0 = cos(ang) ; s0 = sin(ang)

      long2 = long / 2
      larg2 = larg / 2

      if ( long == 2 * (long/2) ) long2 = long/2 + 1
      if ( larg == 2 * (larg/2) ) larg2 = larg/2 + 1

      tab_out(1:long, 1:larg) = tab_in(1:long, 1:larg)

      rr = 1.e6_R8

      do j = 1, larg
      do i = 1, long

         x = ( i - long2 ) * PARAM%surf_dx
         y = ( j - larg2 ) * PARAM%surf_dy

         ! Distance as expressed in the theoretical acf. When r2 is 1., the acf is 0.2.
         r2 = ( (+ c0 * x + s0 * y) / tau1 ) ** 2 +  &  !
              ( (- s0 * x + c0 * y) / tau2 ) ** 2

         ! Below the correlation length, no transformation
         if ( r2 <= 1._R8 ) cycle

         ! The correlation length is exceeded.
         ! The angle corresponding to the position (i,j) is calculated
         theta = atan2( y, x ) ; if ( theta < 0. ) theta = theta + 2 * Pi_R8

         t = tan(theta)

         ! According the location of (i,j) (right, top, left or bottom) the line that begin at the surface center, passing by (i,j),
         ! ends on one of the four borders.
         if ( theta > 2 * Pi_R8 - theta_diag .or.  theta <=           + theta_diag ) then ; x = ( long - long2 ) * PARAM%surf_dx ; y = x * t ; endif

         if ( theta >           + theta_diag .and. theta <=     PI_R8 - theta_diag ) then ; y = ( larg - larg2 ) * PARAM%surf_dy ; x = y / t ; endif

         if ( theta >     PI_R8 - theta_diag .and. theta <=     PI_R8 + theta_diag ) then ; x = (    1 - long2 ) * PARAM%surf_dx ; y = x * t ; endif

         if ( theta >     PI_R8 + theta_diag .and. theta <= 2 * Pi_R8 - theta_diag ) then ; y = (    1 - larg2 ) * PARAM%surf_dy ; x = y / t ; endif

         ! The same distance as above is calculated, from the center to the surface edge, then ...
         rd = ( (+ c0 * x + s0 * y) / tau1 ) ** 2 +  &  !
              ( (- s0 * x + c0 * y) / tau2 ) ** 2

         rr = min( rr, sqrt(rd) )

      enddo
      enddo

      rr = 0.99 * rr

      sum_inn = 0
      sum_tab = 0
      sum_int = 0
      npt_out = 0

      n = 2

      do j = 1, larg
      do i = 1, long

         x = ( i - long2 ) * PARAM%surf_dx
         y = ( j - larg2 ) * PARAM%surf_dy

         ! Distance as expressed in the theoretical acf. When r2 is 1., the acf is 0.2.
         r2 = ( (+ c0 * x + s0 * y) / tau1 ) ** 2 +  &  !
              ( (- s0 * x + c0 * y) / tau2 ) ** 2

         ! Below the correlation length, no transformation
         if ( r2 <= 1._R8 ) then

            sum_inn = sum_inn + tab_in(i, j)

            cycle

         endif

         r2 = sqrt(r2)

         if ( r2 >= rr ) then

            npt_out = npt_out + 1

            cycle

         endif

         ! ... the modified Tuckey window can be determined.

         sum_tab = sum_tab + ( cos( 0.5_R8 * PI_R8 * ( r2 - 1. ) / ( rr - 1. ) ) ** n ) * tab_in(i, j)
         sum_int = sum_int + ( r2 - 1. ) / ( rr - 1. )

      enddo
      enddo

      a_min = -( sum_inn + sum_tab ) / ( npt_out + sum_int )

      do j = 1, larg
      do i = 1, long

         x = ( i - long2 ) * PARAM%surf_dx
         y = ( j - larg2 ) * PARAM%surf_dy

         ! Distance as expressed in the theoretical acf. When r2 is 1., the acf is 0.2.
         r2 = ( (+ c0 * x + s0 * y) / tau1 ) ** 2 +  &  !
              ( (- s0 * x + c0 * y) / tau2 ) ** 2

         ! Below the correlation length, no transformation
         if ( r2 <= 1._R8 ) cycle

         r2 = sqrt(r2)

         if ( r2 >= rr ) then

            tab_out(i, j) = a_min

            cycle

         endif

         ! ... the modified Tuckey window can be determined.
         tab_out(i, j) = ( cos( 0.5_R8 * PI_R8 * ( r2 - 1. ) / ( rr - 1. ) ) ** n ) * tab_in(i, j) + a_min * ( r2 - 1. ) / ( rr - 1. )

      enddo
      enddo

      allocate( tab_tmp(1:long, 1:larg) )

      do k = 1, 10

         tab_tmp(1:long, 1:larg) = tab_out(1:long, 1:larg)

         do j = 1 + 1, larg - 1
         do i = 1 + 1, long - 1

            x = ( i - long2 ) * PARAM%surf_dx
            y = ( j - larg2 ) * PARAM%surf_dy

            ! Distance as expressed in the theoretical acf. When r2 is 1., the acf is 0.2.
            r2 = ( (+ c0 * x + s0 * y) / tau1 ) ** 2 +  &  !
                 ( (- s0 * x + c0 * y) / tau2 ) ** 2

            r2 = sqrt(r2)

            if ( r2 <= 0.98_R8 .or. ( r2 >= 1.02_R8 .and. r2 <= 0.98_R8 * rr ) .or. r2 >= 1.02_R8 * rr ) cycle

            tab_out(i, j) = ( 2*tab_tmp(i, j) +tab_tmp(i +1, j   ) +tab_tmp(i -1, j   ) + &  !
                                               tab_tmp(i   , j +1) +tab_tmp(i   , j -1) +( tab_tmp(i +1, j -1) +tab_tmp(i -1, j -1) + &  !
                                                                                           tab_tmp(i -1, j +1) +tab_tmp(i +1, j +1) )/sqrt(2._R8) )/( 6. +4./sqrt(2._R8) )

         enddo
         enddo

      enddo

      deallocate( tab_tmp )

   return
   endsubroutine apod2