make_tex Subroutine

private subroutine make_tex()

Note

Function that creates a periodic macro-texture: knowing the FFT of an analytical texture

Arguments

None

Calls

proc~~make_tex~~CallsGraph proc~make_tex make_tex calc_fftw3 calc_fftw3 proc~make_tex->calc_fftw3 calc_moments calc_moments proc~make_tex->calc_moments

Called by

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

Source Code

   subroutine make_tex()
   !================================================================================================
   !< @note Function that creates a periodic macro-texture: knowing the FFT of an analytical texture
   !        a surface is created and then transformed by a FFT. As a result it is periodic.
   !  @endnote
   !------------------------------------------------------------------------------------------------
   implicit none

      real(kind=R8), allocatable, dimension(:)   :: xloc, yloc, haut

      real   (kind=R8) :: ca, sa, a, b, c, x0, y0, h0
      integer(kind=I4) :: kk, w2, h2

      complex  (kind=R8) :: eix, eiy, ftmp
      real     (kind=R8) :: x1, y1, ax, ay, sx, sy, coeff
      integer  (kind=I4) :: ksi, eta

      complex(kind=R8), allocatable, dimension(:,:) :: cmpl_in, cmpl_ou, ft

      real(kind=R8) :: alpha                 ! texture orientation
      real(kind=R8) :: sigmax, sigmay        ! texture stdv

      integer(kind=I4) :: texture            ! number of texture elements

      character(len=016) :: dimples          ! texture elementar shapes: "circle" or "square"

      real(kind=R8) :: dx, dy

      integer(kind=I4) :: w, h, l

      type(MOMENT_STAT) :: m_res

      w = PARAM%width
      h = PARAM%height
      l = PARAM%npts

      dx = PARAM%surf_dx
      dy = PARAM%surf_dy

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

      sigmax = PARAM%curr_surf%cl1 / dx / w * 2
      sigmay = PARAM%curr_surf%cl2 / dx / h * 2
      alpha  = PARAM%curr_surf%ang

      w2 = w / 2
      h2 = h / 2
      if ( w == 2 * (w / 2) + 1 ) w2 = w / 2 +1
      if ( h == 2 * (h / 2) + 1 ) h2 = h / 2 +1

      allocate( cmpl_in(1:w, 1:h),     &  !
                cmpl_ou(1:w, 1:h) )       !

      allocate( ft(-w/2:w2 -1,         &  !
                   -h/2:h2 -1) )          !

      allocate( xloc(1:TEXTURE),       &  !
                yloc(1:TEXTURE),       &  !
                haut(1:TEXTURE) )         !


      ca = cos(alpha * PI_R8 / 180)
      sa = sin(alpha * PI_R8 / 180)

      ! Gaussian coefficients
      a = 0.5*( (ca / sigmax) * (ca / sigmax) + (sa / sigmay) * (sa / sigmay) )
      c = 0.5*( (sa / sigmax) * (sa / sigmax) + (ca / sigmay) * (ca / sigmay) )
      b = 1.0*( (sa / sigmax) * (ca / sigmax) - (sa / sigmay) * (ca / sigmay) )

      call random_number( xloc(1:TEXTURE ) ) ! localisation in [0,1]x[0,1]
      call random_number( yloc(1:TEXTURE ) )

      call random_number( haut(1:TEXTURE ) )

      haut(1:TEXTURE ) = 2 * haut(1:TEXTURE ) - 1._R8

      coeff = 2 * pi_R8

      do ksi = -w / 2, w2 - 1
      do eta = -h / 2, h2 - 1

         x1 = (+ksi * ca + eta * sa) * coeff
         y1 = (-ksi * sa + eta * ca) * coeff

         ftmp = 0
         sx = sigmax
         sy = sigmay

         do kk = 1, TEXTURE ! nombre de pics sur le niveau niv

            x0 = xloc(kk)
            y0 = yloc(kk)
            h0 = haut(kk)   ! d'un niveau à l'autre, on atténue les pics d'un facteur "per_tex"

            ax = -x1 * x0
            ay = -y1 * y0

            eix = cmplx(cos(ax), sin(ax), kind=R8)
            eiy = cmplx(cos(ay), sin(ay), kind=R8)

            select case( DIMPLES(1:6) )

               case("square")
                  ftmp = ftmp + 2 * h0 * pi_R8 * sx * sy * sinc( sx * x1 ) * eix &  !
                                                         * sinc( sy * y1 ) * eiy
               case("circle")
                  ftmp = ftmp + 2 * h0 * pi_R8 * sx * sy * exp( -0.5 * (sx**2) * (x1**2) ) * eix &  !
                                                         * exp( -0.5 * (sy**2) * (y1**2) ) * eiy
               case default
                  stop "make_terra, no dimple type selected"

            endselect

         enddo

         ft(ksi, eta) = ftmp

      enddo
      enddo

      cmpl_in(1:w2,      1:h2) = ft(0:w2 - 1,      0:h2 - 1)   ! SO = NE
      cmpl_in(1:w2, h2 + 1:h ) = ft(0:w2 - 1, -h / 2:   -1 )   ! SE = NO

      cmpl_in(w2 + 1:w,     1:h2) = ft(-w/2:-1,    0:h2 -1)    ! NO = SE
      cmpl_in(w2 + 1:w, h2 +1:h ) = ft(-w/2:-1, -h/2:   -1)    ! NE = SO

      FFT_DIM = 0 ! forces fftw desalloc and reinit

      call calc_fftw3(sens = BACKWARD,                &  !
                    tab_in = cmpl_in(1:w, 1:h),       &  !
                    tab_ou = cmpl_ou(1:w, 1:h),       &  !
                      long = w,                       &  !
                      larg = h)                          !

      FFT_DIM = 0 ! forces fftw desalloc and reinit

      PARAM%surf(1:w, 1:h) = real(cmpl_ou(1:w, 1:h), kind=R8) * sqrt(real(w * h, kind=R8))

      call calc_moments(    tab = PARAM%surf(1:w, 1:h),      &  ! IN
                             mx = m_res,                     &  ! OUT
                         nb_mom = 3 )                           ! IN

      if ( m_res%sk > 0. ) PARAM%surf(1:w, 1:h) = - PARAM%surf(1:w, 1:h)

      deallocate( xloc, yloc, haut )

      deallocate( cmpl_in, cmpl_ou )

   contains
      !-----------------------------------------
      function sinc(x)
      implicit none
      real(kind=R8) :: sinc
      real(kind=R8), intent(in) :: x
         if ( abs(x)>1.0e-8) then
            sinc = sin(x)/x
         else
            sinc = UN
         endif
      return
      endfunction sinc
      !-----------------------------------------
   endsubroutine make_tex