curvature Subroutine

public subroutine curvature(tab, nx, ny, dx, dy, S_param_grad, S_param_curv, gcurvt)

Function to calculate the gaussian curvature of a 2D array, its mean quadratic value and the gradient mean quadratic value

Arguments

Type IntentOptional Attributes Name
real(kind=R8), intent(in), dimension(1:nx, 1:ny) :: tab

input 2D array
integer(kind=I4), intent(in) :: nx

number of pixels along x
integer(kind=I4), intent(in) :: ny

number of pixels along x
real(kind=R8), intent(in) :: dx

x lag
real(kind=R8), intent(in) :: dy

y lag
real(kind=R8), intent(out) :: S_param_grad

mean quadratic gradient value
real(kind=R8), intent(out) :: S_param_curv

mean quadratic curvature value
real(kind=R8), intent(out), dimension(1:nx, 1:ny) :: gcurvt

gaussian curvature 2D array

Calls

proc~~curvature~~CallsGraph proc~curvature curvature proc~gauss_curv gauss_curv proc~curvature->proc~gauss_curv proc~gradient gradient proc~curvature->proc~gradient proc~gauss_curv->proc~gradient

Called by

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

Source Code

   subroutine curvature(tab, nx, ny, dx, dy, S_param_grad, S_param_curv, gcurvt)
   !================================================================================================
   !! Function to calculate the gaussian curvature of a 2D array,
   !!        its mean quadratic value and the gradient mean quadratic value
   implicit none
   integer(kind=I4), intent(in ) :: nx                              !! *number of pixels along x*
   integer(kind=I4), intent(in ) :: ny                              !! *number of pixels along x*
   real   (kind=R8), intent(in ) :: dx                              !! *x lag*
   real   (kind=R8), intent(in ) :: dy                              !! *y lag*
   real   (kind=R8), intent(out) :: S_param_grad                    !! *mean quadratic gradient value*
   real   (kind=R8), intent(out) :: S_param_curv                    !! *mean quadratic curvature value*
   real   (kind=R8), intent(in ), dimension(1:nx, 1:ny) :: tab      !! *input 2D array*
   real   (kind=R8), intent(out), dimension(1:nx, 1:ny) :: gcurvt   !! *gaussian curvature  2D array*

      real(kind=R8), allocatable, dimension(:,:) :: gradx    !! *derivative along x 2D array*
      real(kind=R8), allocatable, dimension(:,:) :: grady    !! *derivative along y 2D array*
      real(kind=R8), allocatable, dimension(:,:) :: gradxx   !! *double derivative along x, x 2D array*
      real(kind=R8), allocatable, dimension(:,:) :: gradyy   !! *double derivative along y, y 2D array*
      real(kind=R8), allocatable, dimension(:,:) :: gradxy   !! *double derivative along x, y 2D array*

      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) )

      call gradient(tab   = tab(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

!~       S_param_grad = sum( sqrt( gradx(1:nx, 1:ny)**2 + grady(1:nx, 1:ny)**2 ) ) / (nx * ny)
      S_param_grad = sqrt( sum( gradx(1:nx, 1:ny)**2 + grady(1:nx, 1:ny)**2 ) ) / (nx * ny)

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

      gcurvt(1:nx, 1:ny) = ( gradxx(1:nx, 1:ny) * gradyy(1:nx, 1:ny) - gradxy(1:nx, 1:ny)**2 ) / ( UN + gradx(1:nx, 1:ny)**2 + grady(1:nx, 1:ny)**2 )**2

!~       S_param_curv = sum( sqrt( gradxx(1:nx, 1:ny)**2 + gradyy(1:nx, 1:ny)**2 ) ) / (nx * ny)
      S_param_curv = sqrt( sum( gradxx(1:nx, 1:ny)**2 + gradyy(1:nx, 1:ny)**2 ) ) / (nx * ny)

      deallocate( gradx, grady, gradxy, gradxx, gradyy )

   return
   endsubroutine curvature