acf_theo Subroutine

private subroutine acf_theo()

Note

Function that returns the theoretical acf PARAM%imp_acf.

If the surface to generate is non periodic, the starting surface is extended. The final surface will be a part of it. Indeed the extended surface will be periodic, because the use of FFTs.

If a roughness orientation is chosen, in addition with long correlation lengths, a windowing should be applied to the acf to prevent from artifacts (vertical and horizontal lines)

Arguments

None

Calls

proc~~acf_theo~~CallsGraph proc~acf_theo acf_theo ellipse_acf ellipse_acf proc~acf_theo->ellipse_acf proc~calc_imp_acf calc_imp_acf proc~acf_theo->proc~calc_imp_acf proc~apod2 apod2 proc~calc_imp_acf->proc~apod2 proc~autocov_impo autocov_impo proc~calc_imp_acf->proc~autocov_impo

Called by

proc~~acf_theo~~CalledByGraph proc~acf_theo acf_theo proc~read_job read_job proc~read_job->proc~acf_theo proc~prg_surf prg_surf proc~prg_surf->proc~read_job program~main main program~main->proc~prg_surf

Source Code

   subroutine acf_theo()
   !================================================================================================
   !<@note Function that returns the theoretical acf PARAM%imp_acf.
   !<
   !< If the surface to generate is non periodic, the starting surface is extended. The final surface
   !< will be a part of it. Indeed the extended surface will be periodic, because the use of FFTs.
   !<
   !< If a roughness orientation is chosen, in addition with long correlation lengths, a windowing
   !< should be applied to the acf to prevent from artifacts (vertical and horizontal lines)
   !<
   !<@endnote
   !------------------------------------------------------------------------------------------------
   implicit none

      integer(kind=I4) :: w, h
      logical(kind=I4) :: with_apod

      real(kind=R8) :: ratio, a, b, c, s, lx, ly

      real(kind=R8), dimension(1:8) :: res

      read(JOB,*) PARAM%l_acf1, PARAM%l_acf2, PARAM%acf__z ; LINE_READ = LINE_READ +1 ; write(SPY,*) "line: ", LINE_READ, "l_acf ", PARAM%l_acf1, PARAM%l_acf2, PARAM%acf__z

      if (PARAM%l_acf1 < PARAM%l_acf2) stop "inverser l_acf1, l_acf2"

      read(JOB,*) PARAM%a_acf  ; LINE_READ = LINE_READ +1 ; write(SPY,*) "line: ", LINE_READ, "a_acf ", PARAM%a_acf

      write(SPY,*) 'acf_theo -> determines the theoretical acf with a padded size, so correlation lengths are adjusted accordingly'

      read(JOB,*) with_apod    ; LINE_READ = LINE_READ +1 ; write(SPY,*) "line: ", LINE_READ, "apod ", with_apod

      PARAM%apod = with_apod

      w = PARAM%width
      h = PARAM%height

      if ( .not.PARAM%periodic ) then

         PARAM%sub_width  = w
         PARAM%sub_height = h
         PARAM%sub_npts   = w * h

         PARAM%sub_surf_width  = PARAM%surf_width
         PARAM%sub_surf_height = PARAM%surf_height

         ! the surface must be extended

         a = PARAM%l_acf1
         b = PARAM%l_acf2

         c = cos( PARAM%a_acf * PI_R8 / 180 )
         s = sin( PARAM%a_acf * PI_R8 / 180 )

         lx = sqrt( (a * c)**2 + (b * s)**2 )
         ly = sqrt( (a * s)**2 + (b * c)**2 )

         w = w + nint( lx / PARAM%surf_dx )
         h = h + nint( ly / PARAM%surf_dy )

         ! update sizes
         PARAM%width  = w
         PARAM%height = h
         PARAM%npts   = w * h

         PARAM%surf_width  = w * PARAM%surf_dx
         PARAM%surf_height = h * PARAM%surf_dy

      endif

      allocate( PARAM%surf(1:w, 1:h) )

      allocate( PARAM%fhi(1:w, 1:h) )

      allocate( PARAM%imp_acf(1:w, 1:h) )

      allocate( PARAM%order(1:w * h) )

      allocate( PARAM%vect_h(1:w * h) )

      allocate( PARAM%acf_surf(1:w, 1:h) )

      ratio = PARAM%l_acf2 / PARAM%l_acf1

      ! acf_theo is the theoretical acf for the normal surface
      call calc_imp_acf( long = w,                                &  ! IN
                         larg = h,                                &  ! IN
                         apod = with_apod,                        &  ! IN
                         tau1 = PARAM%l_acf1,                     &  ! IN
                         tau2 = PARAM%l_acf2,                     &  ! IN
                        alpha = log(PARAM%acf__z),                &  ! IN
                          ang = PARAM%a_acf * PI_R8/180,          &  ! IN
                      tab_acf = PARAM%imp_acf(1:w, 1:h))             ! OUT

      call ellipse_acf( tabin = PARAM%imp_acf(1:w, 1:h),          &  ! IN
                         long = w,                                &  ! IN
                         larg = h,                                &  ! IN
                        p_acv = res(1:8),                         &  ! OUT
                          cut = PARAM%acf__z,                     &  ! IN
                     scale_xy = [PARAM%surf_dx, PARAM%surf_dy],   &  ! IN
                          omp = .true. )                             ! IN

   return
   endsubroutine acf_theo