simple_anisotropy Subroutine

   subroutine simple_anisotropy(tabin, long, larg, vec_len, vec_pks, vec_slp, scale_xy, multi_fft)
   !================================================================================================
   !< @note Function that returns some anisotropy parameters calculated on a polar representation,
   !< for each angle from 0 to 179
   !<
   !< + the excess of length
   !< + the RMS slope
   !< + the mean peak width
   !<
   !  @endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   integer(kind=I4), intent(in )                             :: long       !! *surface length*
   integer(kind=I4), intent(in )                             :: larg       !! *surface width*
   logical(kind=I4), intent(in )                             :: multi_fft  !! *parallel ffts?*
   real   (kind=R8), intent(in ), dimension(1:long, 1:larg)  :: tabin      !! *surface acf array*
   real   (kind=R8), intent(in ), dimension(1:2)             :: scale_xy   !! *lag along x, y*
   real   (kind=R8), intent(out), dimension(0:179)           :: vec_len    !! *vector containing path length ratios*
   real   (kind=R8), intent(out), dimension(0:179)           :: vec_pks    !! *vector containing peak mean width*
   real   (kind=R8), intent(out), dimension(0:179)           :: vec_slp    !! *vector containing RMS slopes*

      integer(kind=I4) :: lo, la, ll, nx, ny, ird
      integer(kind=I4) :: p, q, k, nb_cross, nb_peak
      real   (kind=R8) :: scx, scy, dr, s, c
      real   (kind=R8) :: mx, my, reg_a, reg_b
      real   (kind=R8) :: h1, h2, h3, h4, hh
      real   (kind=R8) :: inc_a, theta, x, y, xb, yb, xm, ym, xp, yp, tmp

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

      real   (kind=R8), allocatable, dimension(:)   :: vec_tmp, slope

      integer(kind=I4), allocatable, dimension(:)   :: cross, peaks

      ! lags in micron
      scx = scale_xy(1)  ! x lag
      scy = scale_xy(2)  ! y lag

      ! define the surface center location
      if ( long == 2 * (long / 2) ) then ; lo = long/2 + 1 ; else ; lo = long/2 ; endif
      if ( larg == 2 * (larg / 2) ) then ; la = larg/2 + 1 ; else ; la = larg/2 ; endif

      ! define the square length embedded in the surface
      ll = min(lo, la) - 1

      allocate( height_disc(0:ll, 0:359) )

      allocate( vec_tmp(0:2*ll) )
      allocate( cross(1:2*ll) )
      allocate( slope(1:2*ll) )
      allocate( peaks(1:2*ll) )

      ! angle increment
      inc_a = 2 * PI_R8 / 360

      ! determination of heights on a diameter of the rosette, obtained by linear interpolation:
      ! a point "falls" into a rectangular element [h1,h2,h3,h4], its height is determined by linear interpolation.
      !* Calculate heights along the radius for each angular position.

      do p = 0, ll                ! point number on a radius
         ird = p                  ! point radius

         do q = 0, 359            ! point angular location (°)
            theta = q * inc_a     ! angular location (rad)

            ! projection on x and y of the point determined by its radius and angle
            ! by taking the floor, we get the number of the lower row and the left column of the rectangle
            ! the remainder (x-nx) thus represents the abscissa of the point in the rectangle with sides=1

            x = lo + ird * cos(theta)

            if ( abs(x-nint(x)) < 1.e-3_R8) then

               nx = nint(x)

            else

               nx = floor(x)

            endif
            xb = x - nx

            y = la + ird * sin(theta)

            if ( abs(y-nint(y)) < 1.e-3_R8) then

               ny = nint(y)

            else

               ny = floor(y)

            endif
            yb = y - ny

            xm = UN - xb ; xp = xb
            ym = UN - yb ; yp = yb

            if ( nx+1 <= long .and. ny+1 <= larg .and.      &  !
                 nx   >= 1    .and. ny   >= 1) then
               ! note ird may be greater than lo or la
               h1 = tabin(nx    , ny    )
               h2 = tabin(nx + 1, ny    )
               h3 = tabin(nx + 1, ny + 1)
               h4 = tabin(nx    , ny + 1)

               hh = h1 * xm * ym +  &  !
                    h2 * xp * ym +  &  !
                    h3 * xp * yp +  &  !
                    h4 * xm * yp       !

               height_disc(p, q) = hh

            endif

         enddo ! q = 0, 359

      enddo ! p = 0, ll

      ! centering
      tmp = sum( height_disc(0:ll, 0:359) ) / ( (ll + 1) * 360 )

      height_disc(0:ll, 0:359) = height_disc(0:ll, 0:359) - tmp

      ! for each degree, determine the mean peak width and the path length ratio
      do q = 0, 359 - 180

         do p = ll, 1, -1 ; vec_tmp(ll - p) = height_disc(p, q + 180) ; enddo ! left part of the height vector
         do p = 0, ll, +1 ; vec_tmp(ll + p) = height_disc(p, q      ) ; enddo ! right part

         theta = q * inc_a ; s = sin(theta) ; c = cos(theta)

         ! element length, usually scx=scy
         dr = sqrt( (s * scx)**2 + (c * scy)**2 )

         ! path relative length  -----------------------------------------------------------------------------------------------
         !
         tmp = 0.
         do p = 1, 2*ll

            tmp = tmp + sqrt( dr**2 + (vec_tmp(p) - vec_tmp(p - 1))**2 )

         enddo

         vec_len(q) = tmp / (2*ll*dr) - 1.

         ! path absolute slope  ------------------------------------------------------------------------------------------------
         !
         do p = 1, 2*ll

            slope(p) = ( vec_tmp(p) - vec_tmp(p - 1) ) / dr

         enddo
         tmp = sum( slope(1:2*ll) ) / (2*ll)
         slope(1:2*ll) = slope(1:2*ll) - tmp

         vec_slp(q) = sqrt( sum( slope(1:2*ll)**2 ) ) / (2*ll)

         ! find each abscissa point where the height crosses z=0 -----------------------------------------------------------------
         !

         !------- subtract mean profile (least square)
         mx = (ll + 1) * dr                       ! x mean x = 0, dr, 2*dr, ..., 2*ll*dr
         my = sum( vec_tmp(0:2*ll) ) / (2*ll + 1) ! y mean

         reg_a = 0
         tmp   = 0
         do p = 0, 2*ll

            reg_a = reg_a + (p * dr - mx) * (vec_tmp(p) - my)
            tmp   = tmp   + (p * dr - mx)**2

         enddo
         reg_a = reg_a / tmp        ! first regressor

         reg_b = my - reg_a * mx    ! second regressor

         do p = 0, 2*ll

            vec_tmp(p) = vec_tmp(p) - (reg_a * p * dr + reg_b)

         enddo

         !------- find peaks
         cross(1:2*ll) = 0
         peaks(1:2*ll) = 0

         k = 1
         cross(k) = 0
         do p = 1, 2*ll

            if ( vec_tmp(p - 1) * vec_tmp(p) < 0. ) then ! the height crosses z=0

               k = k + 1
               cross(k) = p

            endif

         enddo
         k = k + 1
         cross(k) = 2*ll
         nb_cross = k

         ! determine the peak width
         !
         k = 0
         do p = 1, nb_cross - 1

            if ( vec_tmp( (cross(p + 1) + cross(p))/2 ) > 0 ) then ! it must be a peak

               k = k + 1

               peaks(k) = cross(k + 1) - cross(k)

            endif

         enddo
         nb_peak = k

         ! mean peak width
         !
         if ( nb_peak > 0 ) then

            vec_pks(q) = dr * sum( peaks(1:nb_peak) ) / nb_peak

         else

            vec_pks(q) = dr * 2*ll

         endif

      enddo

      deallocate( height_disc, vec_tmp, cross, peaks, slope )

   return
   endsubroutine simple_anisotropy