int_g1_roma Function

    pure function int_g1_roma(r) result(val)
      !> Integral of Roma and Peskin's filter from 0 to r.
      implicit none
      real(wp), intent(in) :: r                                                !! Distance from filter center
      real(wp)             :: val                                              !! Result
      ! Work variables
      real(wp), parameter :: pi=4.0_wp*atan(1.0_wp)
      real(wp) :: part1,part2,part3

      if (abs(r).le.1.0_wp/3.0_wp) then
        part1 =  1.0_wp/8.0_wp*r*(sqrt(4.0_wp-27.0_wp*r**2)+4.0_wp)
        part2 =  sqrt(3.0_wp)/18.0_wp*asin(1.5_wp*sqrt(3.0_wp)*r)
        part3 =  0.0_wp
      else if ((abs(r).ge.1.0_wp/3.0_wp).and.(abs(r).le.1.0_wp)) then
        part1 =  1.25_wp*r-9.0_wp/16.0_wp*r*abs(r)-sqrt(3.0_wp)/36.0_wp*asin(sqrt(3.0_wp)*(1.5_wp*r-r/abs(r)))
        part2 = -1.0_wp/24.0_wp*sqrt(-27.0_wp*r**2+36.0_wp*abs(r)-8.0_wp)*(1.5_wp*r-r/abs(r))
        part3 =  r/abs(r)*(sqrt(3.0_wp)/108.0_wp*pi-1.0_wp/6.0_wp)
      else
        part1 = 0.5_wp*r/abs(r)
        part2 = 0.0_wp
        part3 = 0.0_wp
      end if
      val = part1 + part2 + part3

      return
    end function int_g1_roma