test_label_surf_summits Subroutine

public subroutine test_label_surf_summits()

Function to test the capicity in detecting peaks, pits and saddles in a simple double sinus surface.

Arguments

None

Calls

proc~~test_label_surf_summits~~CallsGraph proc~test_label_surf_summits test_label_surf_summits init_scal init_scal proc~test_label_surf_summits->init_scal proc~label_surf_summits label_surf_summits proc~test_label_surf_summits->proc~label_surf_summits write_surf write_surf proc~test_label_surf_summits->write_surf 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_label_surf_summits~~CalledByGraph proc~test_label_surf_summits test_label_surf_summits program~test_grad_curv test_grad_curv program~test_grad_curv->proc~test_label_surf_summits

Source Code

   subroutine test_label_surf_summits()
   !================================================================================================
   !! Function to test the capicity in detecting peaks, pits and saddles in a simple double sinus
   !!        surface.
   implicit none

      integer(kind=I4), allocatable, dimension(:,:) :: topo
      integer(kind=I4), allocatable, dimension(:,:) :: vall, peak, sadd

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

      integer(kind=I4), dimension(1:3) :: nb_null_derivatives

      integer(kind=I4), parameter :: nx = 750, ny = 1000

      integer(kind=I4) :: i, j, i1, j1

      type(SCALE_SURF) :: scal

      ! create a "digital surf" object
      call init_scal(scal = scal,          &      ! out; creates a surface type, containing ...
                       nx = nx,            &      !  in; ... the number of points along x ...
                       ny = ny,            &      !  in; ... the number of points along y ...
                       lx = nx * 1.e-6_R8, &      !  in; ... the length (default unit : m) ...
                       ly = ny * 1.e-6_R8, &      !  in; ... the width ...
                   unit_z = 'm'       )           !  in; ... and the unit along z.

      allocate( heights(1:nx, 1:ny) )
      allocate( topo   (1:nx, 1:ny) )

      do i = 1, nx
      do j = 1, ny

         heights(i, j) = sinsin(i, j, nx, ny)

      enddo
      enddo

      call write_surf( nom_fic = "out/test_sinus.sur",     &  !
                       tab_s   = heights(1:nx, 1:ny),      &  !
                       scal    = scal )                       !

      call label_surf_summits(tab        = heights(1:nx, 1:ny),      &  !
                              nx         = nx,                       &  !
                              ny         = ny,                       &  !
                              valleys    = vall,                     &  !
                              peaks      = peak,                     &  !
                              saddles    = sadd,                     &  !
                              nb_summits = nb_null_derivatives(1:3) )   !

      write(*,*) nb_null_derivatives(1:3)

      ! csv to compare computed and theoretical values
      open(10, file = 'out/extrema.csv')

         do i = 1, nb_null_derivatives(1)
            write(10,*) 'comput,', vall(i, 1), ',', vall(i, 2), ',', '1,VALLEY'
         enddo

         do i = 1, nb_null_derivatives(2)
            write(10,*) 'comput,', sadd(i, 1), ',', sadd(i, 2), ',', '2,SADDLE'
         enddo

         do i = 1, nb_null_derivatives(3)
            write(10,*) 'comput,', peak(i, 1), ',', peak(i, 2), ',', '3,PEAK'
         enddo

         do i = 0, 2
         i1 = int( nx * (0.5 + i) / 3. )

            do j = 0, 5
            j1 = int( ny * (0.5 + j) / 6.)

               if ( (j == 2 * (j/2)     .and. i == 2 * (i/2)    ) .or.     &  !
                    (j == 2 * (j/2) + 1 .and. i == 2 * (i/2) + 1) ) then      !

                  write(10,*) 'theory,', i1, ',', j1, ',', '3,PEAK'

               else

                  write(10,*) 'theory,', i1, ',', j1, ',', '1,VALLEY'

               endif

         enddo

      enddo

      do i = 0, 2
      i1 = int( nx * (1.0 + i) / 3. )

         do j = 0, 5
         j1 = int( ny * (1.0 + j) / 6.)

            if ( j1 /= ny .and. i1 /= nx ) write(10,*) 'theory,', i1, ',', j1, ',', '2,SADDLE'

         enddo

      enddo

      close(10)

      deallocate( heights, topo, vall, peak, sadd)

   contains

      real(kind=R8) function sinsin(i, j, nx, ny)
      implicit none
      integer(kind=I4), intent(in) :: i, j, nx, ny

         sinsin = sin( 3 * PI_R8 * i / nx + 1.e-8_R8 ) * sin( 6 * PI_R8 * j / ny + 2.e-8_R8 )

      return
      endfunction sinsin

   endsubroutine test_label_surf_summits