test_peaks_and_pits_curvatures Subroutine

public subroutine test_peaks_and_pits_curvatures()

Function to test the function “peaks_and_pits_curvatures” on a real rough surface. Outputs surface gradients and curvatures as 2D arrays or single values.

Arguments

None

Calls

proc~~test_peaks_and_pits_curvatures~~CallsGraph proc~test_peaks_and_pits_curvatures test_peaks_and_pits_curvatures end_fftw3 end_fftw3 proc~test_peaks_and_pits_curvatures->end_fftw3 fftw_plan_with_nthreads fftw_plan_with_nthreads proc~test_peaks_and_pits_curvatures->fftw_plan_with_nthreads init_fftw3 init_fftw3 proc~test_peaks_and_pits_curvatures->init_fftw3 omp_get_max_threads omp_get_max_threads proc~test_peaks_and_pits_curvatures->omp_get_max_threads proc~curv2 curv2 proc~test_peaks_and_pits_curvatures->proc~curv2 proc~curvature curvature proc~test_peaks_and_pits_curvatures->proc~curvature proc~fft_filter fft_filter proc~test_peaks_and_pits_curvatures->proc~fft_filter proc~gauss_curv gauss_curv proc~test_peaks_and_pits_curvatures->proc~gauss_curv proc~gradient gradient proc~test_peaks_and_pits_curvatures->proc~gradient proc~peaks_and_pits_curvatures peaks_and_pits_curvatures proc~test_peaks_and_pits_curvatures->proc~peaks_and_pits_curvatures read_surf read_surf proc~test_peaks_and_pits_curvatures->read_surf unit2IUf unit2IUf proc~test_peaks_and_pits_curvatures->unit2IUf write_surf write_surf proc~test_peaks_and_pits_curvatures->write_surf proc~curvature->proc~gauss_curv proc~curvature->proc~gradient calc_fftw3_real_bwd calc_fftw3_real_bwd proc~fft_filter->calc_fftw3_real_bwd calc_fftw3_real_fwd calc_fftw3_real_fwd proc~fft_filter->calc_fftw3_real_fwd extend extend proc~fft_filter->extend proc~gaussian_filter gaussian_filter proc~fft_filter->proc~gaussian_filter tab_calc_fftw3_real_bwd tab_calc_fftw3_real_bwd proc~fft_filter->tab_calc_fftw3_real_bwd tab_calc_fftw3_real_fwd tab_calc_fftw3_real_fwd proc~fft_filter->tab_calc_fftw3_real_fwd proc~gauss_curv->proc~gradient proc~peaks_and_pits_curvatures->proc~curvature proc~label_surf_summits label_surf_summits proc~peaks_and_pits_curvatures->proc~label_surf_summits sort_array2 sort_array2 proc~peaks_and_pits_curvatures->sort_array2 proc~gradient_corner gradient_corner proc~label_surf_summits->proc~gradient_corner proc~labelize_point labelize_point proc~label_surf_summits->proc~labelize_point selectcase selectcase proc~label_surf_summits->selectcase proc~deriv_n deriv_N proc~labelize_point->proc~deriv_n

Called by

proc~~test_peaks_and_pits_curvatures~~CalledByGraph proc~test_peaks_and_pits_curvatures test_peaks_and_pits_curvatures program~test_grad_curv test_grad_curv program~test_grad_curv->proc~test_peaks_and_pits_curvatures

Source Code

   subroutine test_peaks_and_pits_curvatures()
   !================================================================================================
   !! Function to test the function "peaks_and_pits_curvatures" on a real rough surface.
   !!        Outputs surface gradients and curvatures as 2D arrays or single values.
   implicit none

      real(kind=R8), allocatable, dimension(:,:) :: heights, bf_heights
      real(kind=R8), allocatable, dimension(:,:) :: gradx, grady, gradxx, gradxy, gradyy, cvt

      real(kind=R8)  :: dx, dy, dz, pcurv, pgrad, pmin, pmax, fft_cutoff

      integer(kind=I4) :: nx, ny, nx2, ny2, n_th

      type(SCALE_SURF) :: scal_surf

      call read_surf( nom_fic = "sur/test1.sur", &  ! IN
                        tab_s = heights,         &  ! OUT
                         scal = scal_surf )         ! OUT

      nx = scal_surf%xres
      ny = scal_surf%yres

      dx = scal_surf%dx * unit2IUf(scal_surf%dx_unit)
      dy = scal_surf%dy * unit2IUf(scal_surf%dy_unit)
      dz = 0

      fft_cutoff = dx / 5.e-6

      nx2 = 2 * ( nint(PAD_FFT * nx)/2 )
      ny2 = 2 * ( nint(PAD_FFT * ny)/2 )

      write(*,*) 'nx, ny = ', nx, ny
      write(*,*) 'dx, dy = ', dx, dy

      !=================================================================

      allocate( bf_heights(1:nx, 1:ny) )

      allocate( gradx(1:nx, 1:ny) )
      allocate( grady(1:nx, 1:ny) )

      allocate( gradxx(1:nx, 1:ny) )
      allocate( gradxy(1:nx, 1:ny) )
      allocate( gradyy(1:nx, 1:ny) )

      allocate( cvt(1:nx, 1:ny) )

      n_th = omp_get_max_threads()
      call fftw_plan_with_nthreads(nthreads = n_th)

      call init_fftw3(long = nx2,  & !
                      larg = ny2 )   !

      call fft_filter(tab       = heights(1:nx, 1:ny),      & ! in
                      long      = nx,                       & ! in
                      larg      = ny,                       & ! in
                      cutoff    = fft_cutoff,               & ! in
                      bf_tab    = bf_heights(1:nx, 1:ny),   & ! out
                      multi_fft = .FALSE.)                    ! in

      call end_fftw3()

      call write_surf( nom_fic = "out/test_bf_heights.sur", &  ! in
                       tab_s   = bf_heights(1:nx, 1:ny),    &  ! in
                       scal    = scal_surf )                   ! in

      call gradient(tab   = bf_heights(1:nx, 1:ny),   & ! in
                    nx    = nx,                       & ! in
                    ny    = ny,                       & ! in
                    dx    = dx,                       & ! in
                    dy    = dy,                       & ! in
                    gradx = gradx(1:nx, 1:ny),        & ! out
                    grady = grady(1:nx, 1:ny))          ! out

            call write_surf( nom_fic = "out/test_gradq.sur",  &  ! in
                             tab_s   = gradx(1:nx, 1:ny)**2 + &  ! in
                                       grady(1:nx, 1:ny)**2,  &  !
                             scal    = scal_surf )               ! in

            call write_surf( nom_fic = "out/test_gradx.sur",  &  ! in
                             tab_s   = gradx(1:nx, 1:ny),     &  ! in
                             scal    = scal_surf )               ! in

            call write_surf( nom_fic = "out/test_grady.sur",  &  ! in
                             tab_s   = grady(1:nx, 1:ny),     &  ! in
                             scal    = scal_surf )               ! in

      !=================================================================

      call gauss_curv(gradx  = gradx (1:nx, 1:ny),    & !
                      grady  = grady (1:nx, 1:ny),    & !
                      gradxx = gradxx(1:nx, 1:ny),    & !
                      gradxy = gradxy(1:nx, 1:ny),    & !
                      gradyy = gradyy(1:nx, 1:ny),    & !
                      nx     = nx,                    & !
                      ny     = ny,                    & !
                      dx     = dx,                    & !
                      dy     = dy)                      !

            call write_surf( nom_fic = "out/test_gradxx_1.sur",  &  ! in
                             tab_s   = gradxx(1:nx, 1:ny),       &  ! in
                             scal    = scal_surf )                  ! in

            call write_surf( nom_fic = "out/test_gradxy_1.sur",  &  ! in
                             tab_s   = gradxy(1:nx, 1:ny),       &  ! in
                             scal    = scal_surf )                  ! in

            call write_surf( nom_fic = "out/test_gradyy_1.sur",  &  ! in
                             tab_s   = gradyy(1:nx, 1:ny),       &  ! in
                             scal    = scal_surf )                  ! in

      call curv2(tab    = bf_heights(1:nx, 1:ny),    & ! in
                 nx     = nx,                        & ! in
                 ny     = ny,                        & ! in
                 dx     = dx,                        & ! in
                 dy     = dy,                        & ! in
                 gradxx = gradxx    (1:nx, 1:ny),    & ! out
                 gradyy = gradyy    (1:nx, 1:ny))      ! out

            call write_surf( nom_fic = "out/test_gradxx_2.sur",  &  ! in
                             tab_s   = gradxx(1:nx, 1:ny),       &  ! in
                             scal    = scal_surf )                  ! in

            call write_surf( nom_fic = "out/test_gradxy_2.sur",  &  ! in
                             tab_s   = gradxy(1:nx, 1:ny),       &  ! in
                             scal    = scal_surf )                  ! in

            call write_surf( nom_fic = "out/test_gradyy_2.sur",  &  ! in
                             tab_s   = gradyy(1:nx, 1:ny),       &  ! in
                             scal    = scal_surf )                  ! in

      call curvature(tab          = bf_heights(1:nx, 1:ny),   & ! in
                     nx           = nx,                       & ! in
                     ny           = ny,                       & ! in
                     dx           = dx,                       & ! in
                     dy           = dy,                       & ! in
                     S_param_grad = pgrad,                    & ! out
                     S_param_curv = pcurv,                    & ! out
                     gcurvt       = cvt(1:nx, 1:ny))            ! out

      write(*,*) "S_param_grad: ", pgrad

            call write_surf( nom_fic = "out/test_cvt.sur",       &  ! in
                             tab_s   = cvt(1:nx, 1:ny),          &  ! in
                             scal    = scal_surf )                  ! in

      !=================================================================

      call peaks_and_pits_curvatures( heights      = bf_heights(1:nx, 1:ny),   & ! in
                                      nx           = nx,                       & ! in
                                      ny           = ny,                       & ! in
                                      dx           = dx,                       & ! in
                                      dy           = dy,                       & ! in
                                      S_param_grad = pgrad,                    & ! out
                                      S_param_curv = pcurv,                    & ! out
                                      peak_curv    = pmax,                     & ! out
                                      pits_curv    = pmin )                      ! out

      write(*,*) "S_param_curv: ", pcurv
      write(*,*) "peak_curv: ",    pmin
      write(*,*) "pits_curv: ",    pmax

      !=================================================================

      deallocate( heights, bf_heights, gradx, grady, gradxy, gradxx, gradyy, cvt )

   return
   endsubroutine test_peaks_and_pits_curvatures