program neo2

  use, intrinsic :: iso_c_binding
  use fgsl ! Fortran interface of the GSL Library

  implicit none

  real(fgsl_double), parameter :: x_low = 1.93d1
  real(fgsl_double), parameter :: x_up = 2.0d1
  ! absolute error
  real(fgsl_double), parameter :: epsabs = 1.0d-13
  ! relative error
  real(fgsl_double), parameter :: epsrel = 1.0d-13

  integer(fgsl_size_t), parameter :: n = 100_fgsl_size_t

  integer(fgsl_size_t) :: neval ! number of function evaluations
  real(fgsl_double) :: ra, rda ! result and absolute error
  integer(fgsl_int) :: statusval
  type(fgsl_function) :: stdfunc
  type(fgsl_integration_cquad_workspace) :: workspace
!~   type(fgsl_integration_workspace) :: workspace

  neval = 0

  stdfunc = fgsl_function_init(func_to_integrate, c_null_ptr)
  workspace = fgsl_integration_cquad_workspace_alloc(n)
!~   workspace = fgsl_integration_workspace_alloc(1000_fgsl_size_t)

  statusval = fgsl_integration_cquad(stdfunc, x_low, x_up, &
     & epsabs, epsrel, workspace, ra, rda, neval)
!~   statusval = fgsl_integration_qags(stdfunc, x_low, x_up, &
!~      & epsabs, epsrel, 1000_fgsl_size_t, workspace, ra, rda)

  call fgsl_function_free(stdfunc)
  call fgsl_integration_cquad_workspace_free(workspace)
!~   call fgsl_integration_workspace_free(workspace)

  write(*,*) ra, rda, neval

  stop

contains

  function func_to_integrate(x, params) bind(c)
    real(c_double), value :: x
    type(c_ptr), value :: params
    real(c_double) :: func_to_integrate

    func_to_integrate =  exp(-x**2)
  end function func_to_integrate

end program neo2