diff --git a/src/beams.f90 b/src/beams.f90
index 3210ee4..efc924d 100644
--- a/src/beams.f90
+++ b/src/beams.f90
@@ -95,7 +95,7 @@ contains
     !       ellipses in the transverse plane at the launch point (deg)
 
     use gray_params,  only : antenna_parameters
-    use splines,      only : spline_simple
+    use splines,      only : spline_1d, linear_1d
     use utils,        only : locate
     use logger,       only : log_error
 
@@ -105,13 +105,11 @@ contains
 
     ! local variables
     integer :: u, nisteer, i, k, ii
-    real(wp_) :: steer, dal
-    real(wp_), dimension(:), allocatable ::                &
-      alphastv, betastv, x00v, y00v,                       &
+    character(256) :: msg
+    real(wp_) :: steer
+    real(wp_), dimension(:), allocatable ::  &
+      alphastv, betastv, x00v, y00v,         &
       z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v
-    type(spline_simple) :: beta, waist1, waist2,   &
-                           rci1, rci2, phi1, phi2, &
-                           x0, y0, z0
 
     open(newunit=u, file=params%filenm, status='old', action='read', iostat=err)
     if (err /= 0) then
@@ -145,47 +143,95 @@ contains
     rci1v   = 10 * rci1v
     rci2v   = 10 * rci2v
 
-    call beta%init(alphastv, betastv, nisteer)
-    call waist1%init(alphastv, waist1v, nisteer)
-    call waist2%init(alphastv, waist2v, nisteer)
-    call rci1%init(alphastv, rci1v, nisteer)
-    call rci2%init(alphastv, rci2v, nisteer)
-    call phi1%init(alphastv, phi1v, nisteer)
-    call phi2%init(alphastv, phi2v, nisteer)
-    call x0%init(alphastv, x00v, nisteer)
-    call y0%init(alphastv, y00v, nisteer)
-    call z0%init(alphastv, z00v, nisteer)
+    call locate(alphastv, params%alpha, k)
 
-    if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then
-      call locate(alphastv, params%alpha, k)
-      dal           = params%alpha - alphastv(k)
-      params%beta   = beta%raw_eval(k, dal)
-      params%pos(1) = x0%raw_eval(k, dal)
-      params%pos(2) = y0%raw_eval(k, dal)
-      params%pos(3) = z0%raw_eval(k, dal)
-      params%w(1)   = waist1%raw_eval(k, dal)
-      params%w(2)   = waist2%raw_eval(k, dal)
-      params%ri(1)  = rci1%raw_eval(k, dal)
-      params%ri(2)  = rci2%raw_eval(k, dal)
-      params%phi(1) = phi1%raw_eval(k, dal)
-      params%phi(2) = phi2%raw_eval(k, dal)
-    else
-      ! params%alpha outside table range
-      if(params%alpha >= alphastv(nisteer)) ii=nisteer
-      if(params%alpha <= alphastv(1)) ii=1
-      params%alpha  = alphastv(ii)
-      params%beta   = betastv(ii)
-      params%pos(1) = x00v(ii)
-      params%pos(2) = y00v(ii)
-      params%pos(3) = z00v(ii)
-      params%w(1)   = waist1v(ii)
-      params%w(2)   = waist2v(ii)
-      params%ri(1)  = rci1v(ii)
-      params%ri(2)  = rci2v(ii)
-      params%phi(1) = phi1v(ii)
-      params%phi(2) = phi2v(ii)
+    ! No interpolation: 0D table or out of range
+    if (nisteer == 1 .or. k == 0 .or. k == nisteer) then
+      if (k == 0) k = k + 1
+      params%alpha  = alphastv(k)
+      params%beta   = betastv(k)
+      params%pos(1) = x00v(k)
+      params%pos(2) = y00v(k)
+      params%pos(3) = z00v(k)
+      params%w(1)   = waist1v(k)
+      params%w(2)   = waist2v(k)
+      params%ri(1)  = rci1v(k)
+      params%ri(2)  = rci2v(k)
+      params%phi(1) = phi1v(k)
+      params%phi(2) = phi2v(k)
+      return
     end if
 
+    if (nisteer < 4) then
+      ! Linear interpolation for small tables
+      block
+        type(linear_1d) :: beta, waist1, waist2,   &
+                           rci1, rci2, phi1, phi2, &
+                           x0, y0, z0
+
+        call beta%init(alphastv, betastv, nisteer)
+        call waist1%init(alphastv, waist1v, nisteer)
+        call waist2%init(alphastv, waist2v, nisteer)
+        call rci1%init(alphastv, rci1v, nisteer)
+        call rci2%init(alphastv, rci2v, nisteer)
+        call phi1%init(alphastv, phi1v, nisteer)
+        call phi2%init(alphastv, phi2v, nisteer)
+        call x0%init(alphastv, x00v, nisteer)
+        call y0%init(alphastv, y00v, nisteer)
+        call z0%init(alphastv, z00v, nisteer)
+
+        params%beta   = beta%raw_eval(k, params%alpha)
+        params%pos(1) = x0%raw_eval(k, params%alpha)
+        params%pos(2) = y0%raw_eval(k, params%alpha)
+        params%pos(3) = z0%raw_eval(k, params%alpha)
+        params%w(1)   = waist1%raw_eval(k, params%alpha)
+        params%w(2)   = waist2%raw_eval(k, params%alpha)
+        params%ri(1)  = rci1%raw_eval(k, params%alpha)
+        params%ri(2)  = rci2%raw_eval(k, params%alpha)
+        params%phi(1) = phi1%raw_eval(k, params%alpha)
+        params%phi(2) = phi2%raw_eval(k, params%alpha)
+      return
+      end block
+    end if
+
+    ! Cubic interpolation otherwise
+    block
+      type(spline_1d) :: beta, waist1, waist2,   &
+                         rci1, rci2, phi1, phi2, &
+                         x0, y0, z0
+
+      call beta%init(alphastv, betastv, nisteer, err=err)
+      call waist1%init(alphastv, waist1v, nisteer, err=err)
+      call waist2%init(alphastv, waist2v, nisteer, err=err)
+      call rci1%init(alphastv, rci1v, nisteer, err=err)
+      call rci2%init(alphastv, rci2v, nisteer, err=err)
+      call phi1%init(alphastv, phi1v, nisteer, err=err)
+      call phi2%init(alphastv, phi2v, nisteer, err=err)
+      call x0%init(alphastv, x00v, nisteer, err=err)
+      call y0%init(alphastv, y00v, nisteer, err=err)
+      call z0%init(alphastv, z00v, nisteer, err=err)
+
+      if (err > 0) then
+        write (msg, '(a, a, g0)') 'interpolation failed', &
+                                  '`curfit` returned ier=', err
+        call log_error(msg, mod='beams', proc='read_beam1')
+        return
+      end if
+      err = 0
+
+      call locate(beta%knots, params%alpha, k)
+      params%beta   = beta%raw_eval(k, params%alpha)
+      params%pos(1) = x0%raw_eval(k, params%alpha)
+      params%pos(2) = y0%raw_eval(k, params%alpha)
+      params%pos(3) = z0%raw_eval(k, params%alpha)
+      params%w(1)   = waist1%raw_eval(k, params%alpha)
+      params%w(2)   = waist2%raw_eval(k, params%alpha)
+      params%ri(1)  = rci1%raw_eval(k, params%alpha)
+      params%ri(2)  = rci2%raw_eval(k, params%alpha)
+      params%phi(1) = phi1%raw_eval(k, params%alpha)
+      params%phi(2) = phi2%raw_eval(k, params%alpha)
+    end block
+
   end subroutine read_beam1
 
 
diff --git a/src/gray_equil.f90 b/src/gray_equil.f90
index 38c413c..8e34295 100644
--- a/src/gray_equil.f90
+++ b/src/gray_equil.f90
@@ -4,7 +4,7 @@
 !
 module gray_equil
   use const_and_precisions, only : wp_, comp_huge
-  use splines,              only : spline_simple, spline_1d, spline_2d, linear_1d
+  use splines,              only : spline_1d, spline_2d, linear_1d
   use types,                only : contour
 
   implicit none
@@ -26,12 +26,12 @@ module gray_equil
     real(wp_) :: z_boundary(2)  ! z range of the plasma boundary
 
     ! Flux average splines (see `flux_average`)
-    type(spline_simple) :: spline_area, spline_volume
-    type(spline_simple) :: spline_dAdrho_t, spline_dVdrho_t
-    type(spline_simple) :: spline_R_J, spline_B_avg, spline_B_min, spline_B_max
-    type(spline_simple) :: spline_f_c, spline_q, spline_I_p, spline_J_phi_avg
-    type(spline_simple) :: spline_astra_tor, spline_jintrac_tor, spline_paral_tor
-    type(spline_2d)     :: spline_cd_eff
+    type(spline_1d) :: spline_area, spline_volume
+    type(spline_1d) :: spline_dAdrho_t, spline_dVdrho_t
+    type(spline_1d) :: spline_R_J, spline_B_avg, spline_B_min, spline_B_max
+    type(spline_1d) :: spline_f_c, spline_q, spline_I_p, spline_J_phi_avg
+    type(spline_1d) :: spline_astra_tor, spline_jintrac_tor, spline_paral_tor
+    type(spline_2d) :: spline_cd_eff
 
     contains
       procedure(pol_flux_sub),     deferred :: pol_flux
@@ -147,11 +147,10 @@ module gray_equil
     private
     real(wp_) :: fpol_a  ! Poloidal current at the edge (r=a), F_a
     ! Splines
-    type(spline_2d)     :: psi_spline
-    type(contour)       :: psi_domain
-    type(spline_1d)     :: fpol_spline
-    type(spline_simple) :: q_spline
-    type(linear_1d)     :: rhop_spline, rhot_spline
+    type(spline_2d) :: psi_spline
+    type(contour)   :: psi_domain
+    type(spline_1d) :: fpol_spline, q_spline
+    type(linear_1d) :: rhop_spline, rhot_spline
     contains
       procedure :: pol_flux      => numeric_pol_flux
       procedure :: pol_curr      => numeric_pol_curr
@@ -773,6 +772,7 @@ contains
     use gray_params,          only : EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
     use gray_params,          only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING
     use utils,                only : vmaxmin
+    use numint,               only : simpson
     use logger,               only : log_debug, log_info
 
     ! subroutine arguments
@@ -1017,13 +1017,11 @@ contains
 
       call self%q_spline%init(data%psi, data%q, npsi)
 
-      ! Toroidal flux as Φ(ψ) = 2π ∫ q(ψ)dψ
-      ! Specifically, Φ(ψ_n) = 2π|ψ_a| ∫₀ψ_n q(ψ_n)dψ_n, with ψ_n ∈ [0,1]
+      ! Toroidal flux as Φ(ψ) = 2π ∫₀^ψ q(ψ) dψ
       phi(1) = 0
-      do k = 1, npsi-1
-        dx = self%q_spline%data(k+1) - self%q_spline%data(k)
-        phi(k+1) = phi(k) + dx*(self%q_spline%coeffs(k,1)   + dx*(self%q_spline%coeffs(k,2)/2 +  &
-                            dx*(self%q_spline%coeffs(k,3)/3 + dx* self%q_spline%coeffs(k,4)/4) ) )
+      do k = 2, npsi
+        call simpson(2, h=data%psi(2)-data%psi(1), fi=data%q(k-1:k), s=phi(k))
+        phi(k) = phi(k) + phi(k-1)
       end do
 
       self%phi_a = phi(npsi)
diff --git a/src/gray_plasma.f90 b/src/gray_plasma.f90
index 30a198e..1392f5a 100644
--- a/src/gray_plasma.f90
+++ b/src/gray_plasma.f90
@@ -8,7 +8,7 @@
 !
 module gray_plasma
   use const_and_precisions, only : wp_, zero
-  use splines,              only : spline_simple, spline_1d
+  use splines,              only : spline_1d
 
   implicit none
 
@@ -79,10 +79,10 @@ module gray_plasma
   ! Numerical plasma description
   type, extends(abstract_plasma) :: numeric_plasma
     private
-    type(spline_1d)      :: dens_spline
-    type(spline_simple)  :: temp_spline
-    type(spline_simple)  :: zeff_spline
-    type(density_tail)   :: tail
+    type(spline_1d)    :: dens_spline
+    type(spline_1d)    :: temp_spline
+    type(spline_1d)    :: zeff_spline
+    type(density_tail) :: tail
   contains
     procedure :: init    => init_splines
     procedure :: density => numeric_density
diff --git a/src/splines.f90 b/src/splines.f90
index dc2d599..14f25a6 100644
--- a/src/splines.f90
+++ b/src/splines.f90
@@ -1,7 +1,6 @@
 ! This module provides a high-level interface for creating and evaluating
 ! several kind of splines:
 !
-!   `spline_simple` is a simple interpolating cubic spline,
 !   `spline_1d` and `spline_2d` are wrapper around the DIERCKX cubic splines,
 !   `linear_1d` is a linear interpolation.
 module splines
@@ -9,27 +8,16 @@ module splines
 
   implicit none
 
-  ! A 1D interpolating cubic spline
-  type spline_simple
-    integer :: ndata                      ! Number of data points
-    real(wp_), allocatable :: data(:)     ! Data points         (ndata)
-    real(wp_), allocatable :: coeffs(:,:) ! Spline coefficients (ndata, 4)
-    contains
-      procedure :: init     => spline_simple_init
-      procedure :: eval     => spline_simple_eval
-      procedure :: raw_eval => spline_simple_raw_eval
-      procedure :: deriv    => spline_simple_deriv
-  end type
-
   ! A 1D smoothing/interpolating cubic spline
   type spline_1d
     integer :: nknots                   ! Number of spline knots
     real(wp_), allocatable :: knots(:)  ! Knots positions
     real(wp_), allocatable :: coeffs(:) ! B-spline coefficients
     contains
-      procedure :: init   => spline_1d_init
-      procedure :: eval   => spline_1d_eval
-      procedure :: deriv  => spline_1d_deriv
+      procedure :: init     => spline_1d_init
+      procedure :: eval     => spline_1d_eval
+      procedure :: raw_eval => spline_1d_raw_eval
+      procedure :: deriv    => spline_1d_deriv
   end type
 
   ! Wrapper to store pointers in an array
@@ -70,178 +58,10 @@ module splines
   end type
 
   private
-  public spline_simple, spline_1d, spline_2d, linear_1d
+  public spline_1d, spline_2d, linear_1d
 
 contains
 
-  subroutine spline_simple_init(self, x, y, n)
-    ! Initialises the spline
-
-    ! subroutine arguments
-    class(spline_simple),    intent(inout) :: self
-    integer,                 intent(in)    :: n
-    real(wp_), dimension(n), intent(in)    :: x, y
-
-    self%ndata = n
-    allocate(self%data(n))
-    allocate(self%coeffs(n, 4))
-
-    self%data = x
-    call spline_simple_coeffs(x, y, n, self%coeffs)
-  end subroutine spline_simple_init
-
-
-  pure function spline_simple_eval(self, x) result(y)
-    ! Evaluates the spline at x
-    use utils, only : locate
-
-    ! subroutine arguments
-    class(spline_simple), intent(in) :: self
-    real(wp_),            intent(in) :: x
-    real(wp_)                        :: y
-
-    ! local variables
-    integer   :: i
-    real(wp_) :: dx
-
-    call locate(self%data, x, i)
-    i = min(max(1, i), self%ndata - 1)
-    dx = x - self%data(i)
-    y = self%raw_eval(i, dx)
-  end function spline_simple_eval
-
-
-  pure function spline_simple_raw_eval(self, i, dx) result(y)
-    ! Evaluates the i-th polynomial of the spline at dx
-
-    ! subroutine arguments
-    class(spline_simple), intent(in) :: self
-    integer,              intent(in) :: i
-    real(wp_),            intent(in) :: dx
-    real(wp_)                        :: y
-
-    y = self%coeffs(i,1) + dx*(self%coeffs(i,2) &
-                         + dx*(self%coeffs(i,3) + dx*self%coeffs(i,4)))
-  end function spline_simple_raw_eval
-
-
-  pure function spline_simple_deriv(self, x) result(y)
-    ! Computes the derivative of the spline at x
-    use utils, only : locate
-
-
-    ! subroutine arguments
-    class(spline_simple), intent(in) :: self
-    real(wp_),            intent(in) :: x
-    real(wp_)                        :: y
-
-    ! local variables
-    integer   :: i
-    real(wp_) :: dx
-
-    call locate(self%data, x, i)
-    i = min(max(1, i), self%ndata - 1)
-    dx = x - self%data(i)
-    y = self%coeffs(i,2) + dx*(2*self%coeffs(i,3) + 3*dx*self%coeffs(i,4))
-  end function spline_simple_deriv
-
-
-  pure subroutine spline_simple_coeffs(x, y, n, c)
-    ! Computes the cubic coefficients of all n polynomials
-
-    ! subroutine arguments
-    integer,   intent(in)    :: n
-    real(wp_), intent(in)    :: x(n), y(n)
-    real(wp_), intent(inout) :: c(n*4)
-
-    ! local variables
-    integer :: jmp, i, ii, j, j1, j2, j3, k
-    real(wp_) :: xb, xc, ya, yb, h, a, r, dya, dyb, dy2
-    jmp = 1
-
-    if (n <= 1) return
-
-    ! initialisation
-    xc = x(1)
-    yb = y(1)
-    h = 0
-    a = 0
-    r = 0
-    dyb = 0
-    dy2 = c(2)
-
-    ! form the system of linear equations
-    ! and eliminate subsequentially
-    j = 1
-    do i = 1, n
-      j2 = n + i
-      j3 = j2 + n
-      a = h*(2 - a)
-      dya = dyb + h*r
-      if (i >= n) then
-        ! set derivative dy2 at last point
-        dyb = dy2
-        h = 0
-        dyb = dya
-        goto 13
-      else
-        j  = j+jmp
-        xb = xc
-        xc = x(j)
-        h  = xc-xb
-
-        ! ii=0 - increasing abscissae
-        ! ii=1 - decreasing abscissae
-        ii = 0
-        if (h == 0) return
-        if (h < 0) ii = 1
-        ya = yb
-        yb = y(j)
-        dyb = (yb - ya)/h
-
-        if (i <= 1) then
-          j1 = ii
-          goto 13
-        end if
-      end if
-
-      if (j1-ii /= 0) return
-      a = 1 / (2*h + a)
-      13 continue
-      r = a*(dyb - dya)
-      c(j3) = r
-      a = h*a
-      c(j2) = a
-      c(i)  = dyb
-    end do
-
-    ! back substitution of the system of linear equations
-    ! and computation of the other coefficients
-    a  = 1
-    j1 = j3+n+ii-ii*n
-    i  = n
-    do k = 1, n
-      xb = x(j)
-      h  = xc - xb
-      xc = xb
-      a  = a+h
-      yb = r
-      r  = c(j3)-r*c(j2)
-      ya = 2*r
-      c(j3) = ya + r
-      c(j2) = c(i) - h*(ya+yb)
-      c(j1) = (yb - r)/a
-      c(i)  = y(j)
-      a = 0
-      j = j-jmp
-      i = i-1
-      j2 = j2-1
-      j3 = j3-1
-      j1 = j3+n+ii
-    end do
-  end subroutine spline_simple_coeffs
-
-
   subroutine spline_1d_init(self, x, y, n, range, weights, tension, err)
     ! Initialises a spline_1d.
     ! Takes:
@@ -260,7 +80,7 @@ contains
     real(wp_),        intent(in)            :: x(n)
     real(wp_),        intent(in)            :: y(n)
     integer,          intent(in)            :: n
-    real(wp_),        intent(in)            :: range(2)
+    real(wp_),        intent(in), optional  :: range(2)
     real(wp_),        intent(in),  optional :: weights(n)
     real(wp_),        intent(in),  optional :: tension
     integer,          intent(out), optional :: err
@@ -273,10 +93,12 @@ contains
 
     ! default values
     integer   :: err_def
-    real(wp_) :: weights_def(n), tension_def
+    real(wp_) :: range_def(2), weights_def(n), tension_def
 
+    range_def = [x(1), x(n)]
     weights_def = 1
     tension_def = 0
+    if (present(range)) range_def = range
     if (present(weights)) weights_def = weights
     if (present(tension)) tension_def = tension
 
@@ -286,9 +108,10 @@ contains
       allocate(self%knots(nknots_est), self%coeffs(nknots_est))
     end if
 
-    call curfit(0, n, x, y, weights_def, range(1), range(2), 3, tension_def, &
-                nknots_est, self%nknots, self%knots, self%coeffs, residuals, &
-                work_real, size(work_real), work_int, err_def)
+    call curfit(0, n, x, y, weights_def, range_def(1), range_def(2), &
+                3, tension_def, nknots_est, self%nknots, self%knots, &
+                self%coeffs, residuals, work_real, size(work_real),  &
+                work_int, err_def)
 
     if (present(err)) err = err_def
   end subroutine spline_1d_init
@@ -312,6 +135,31 @@ contains
   end function spline_1d_eval
 
 
+  pure function spline_1d_raw_eval(self, i, x) result(y)
+    ! Evaluates the spline at: t(i) ≤ x < t(i+1), where t are the knots
+    use dierckx, only : fpbspl
+
+    ! subroutine arguments
+    class(spline_1d), intent(in) :: self
+    integer,          intent(in) :: i
+    real(wp_),        intent(in) :: x
+    real(wp_)                    :: y
+
+    ! local variables
+    integer   :: err, j
+    real(wp_) :: h(6)
+
+    ! Evaluate the basis splines
+    call fpbspl(self%knots, self%nknots, 3, x, i, h)
+
+    ! Take the linear combination with the coefficients
+    y = 0
+    do j = 1, 4
+      y = y + self%coeffs(i-4+j)*h(j)
+    end do
+  end function spline_1d_raw_eval
+
+
   pure function spline_1d_deriv(self, x, order) result(y)
     ! Evaluates the spline n-th order derivative at x
     use dierckx, only : splder