!=======================================================================
!     FUNCTION   --- w ---                                           ===
!=======================================================================
!     Function for "Kernel Function"                                 ===
!                                                                    ===
!         In this function program,                                  ===
!         spline function with compact support is used.              ===
!                        ref. Monaghan & Lattanzio (1985)            ===
!                                                                    ===
!         Kernel function is given as                                ===
!                             1-3/2*(r/h)**2+3/4*(r/h)**3   0<r/h<1  ===
!         W(r,h)=1/pai/h**3   1/4*[2-(r/h)]**3              1<r/h<2  ===
!                             0                               r/h>2  ===
!                                                                    ===
!                                                         2005/07/31 ===
!                                                         2011/06/25 ===
!                                                         2011/08/22 ===
!=======================================================================

real(8) function w(r2,h)

  use mod_parameter

  implicit none

!==== Arguments ========================================================
  real(8), intent(IN) :: r2    ! Distance Squared
  real(8), intent(IN) :: h     ! Smoothing Length

!==== Variables Used in This Function ==================================
  real(8) :: h2      ! h**2
  real(8) :: h2_rev  ! 1/h**2
!=======================================================================

  h2=h*h

  if (r2 .lt. h2) then
     h2_rev=1.0d0/h2
     w=1.0d0-1.5d0*r2*h2_rev+0.75d0*dsqrt(r2*h2_rev)**3
     w=w*PAI_REV*h2_rev/h
  else if (r2 .lt. 4.0d0*h2) then
     h2_rev=1.0d0/h2
     w=0.25d0*(2.0d0-dsqrt(r2*h2_rev))**3
     w=w*PAI_REV*h2_rev/h
  else
     w=0.0d0
  end if


  return
end function w


!=======================================================================
!     FUNCTION   --- dw ---                                          ===
!=======================================================================
!     Function for differentiation of Kernel Function                ===
!                                                         2005/07/31 ===
!                                                         2011/06/25 ===
!                                                         2011/08/22 ===
!=======================================================================

real(8) function dw(r2,h)

  use mod_parameter

  implicit none

!==== Arguments ========================================================
  real(8), intent(IN) :: r2    ! Distance Squared
  real(8), intent(IN) :: h     ! Smoothing Length

!==== Variables Used in This Function ==================================
  real(8) :: h2    ! h**2
  real(8) :: r     ! sqrt(r2)
  real(8) :: h_rev ! 1/h
!=======================================================================

  h2=h*h

  if (r2 .lt. h2) then
     h_rev=1.0d0/h
     dw=3.0d0-2.25d0*dsqrt(r2)*h_rev
     dw=-dw*PAI_REV*h_rev**5
  else if (r2 .lt. 4.0d0*h2) then
     r=dsqrt(r2)
     dw=0.75d0*(2.0d0-r/h)**2/r
     dw=-dw*PAI_REV/(h2**2)
  else
     dw=0.0d0
  end if

  return
end function dw


!=======================================================================
!     FUNCTION   --- dwdh ---                                        ===
!=======================================================================
!     Function for "dw/dh"                                           ===
!                                                                    ===
!                                                         2011/08/30 ===
!=======================================================================

real(8) function dwdh(r,h)

  use mod_parameter

  implicit none

!==== Arguments ========================================================
  real(8), intent(IN) :: r     ! Distance
  real(8), intent(IN) :: h     ! Smoothing Length

!==== Variables Used in This Function ==================================
  real(8) :: q       ! r/h
  real(8) :: q2      ! (r/h)**2
  real(8) :: q3      ! (r/h)**3
!=======================================================================

  if (r .lt. h) then
     q=r/h
     q2=q**2
     q3=q*q2
     dwdh=-3.0d0+7.5d0*q2-4.5d0*q3
     dwdh=dwdh*PAI_REV/h**4
  else if (r .lt. 2.0d0*h) then
     q=r/h
     q2=q**2
     q3=q*q2
     dwdh=-6.0d0+12.0d0*q-7.5d0*q2+1.5d0*q3
     dwdh=dwdh*PAI_REV/h**4
  else
     dwdh=0.0d0
  end if

  return
end function dwdh



!=======================================================================
!     FUNCTION   --- dphidh ---                                      ===
!=======================================================================
!     Function for                                                   ===
!        "the derivative of the potential with respect to h"         ===
!                                                                    ===
!                                                         2011/08/30 ===
!=======================================================================

real(8) function dphidh(r,h)

  use mod_parameter

  implicit none

!==== Arguments ========================================================
  real(8), intent(IN) :: r     ! Distance
  real(8), intent(IN) :: h     ! Smoothing Length

!==== Variables Used in This Function ==================================
  real(8) :: q       ! r/h
  real(8) :: q2      ! (r/h)**2
  real(8) :: q3      ! (r/h)**3
  real(8) :: q4      ! (r/h)**4
  real(8) :: q5      ! (r/h)**5
!=======================================================================

  if (r .lt. h) then
     q=r/h
     q2=q**2
     q4=q2**2
     q5=q4*q
     dphidh=-2.0d0*q2+1.5d0*q4-0.6d0*q5+1.4d0
     dphidh=dphidh/h**2
  else if (r .lt. 2.0d0*h) then
     q=r/h
     q2=q**2
     q3=q2*q
     q4=q2**2
     q5=q4*q
     dphidh=-4.0d0*q2+4.0d0*q3-1.5d0*q4+0.2d0*q5+1.6d0
     dphidh=dphidh/h**2
  else
     dphidh=0.0d0
  end if

  return
end function dphidh