diff --git a/input/gray.ini b/input/gray.ini
index 1a2dd97..60a0aee 100644
--- a/input/gray.ini
+++ b/input/gray.ini
@@ -145,12 +145,6 @@ irho = 0
 ; Filepath of the equilibrium (relative to this file)
 filenm = "profiles.txt"
 
-; Value of ψ at the plasma boundary  [multipass module]
-; Notes:
-;   1. boundary means zero density;
-;   2. determines when a ray is considered inside the plasma
-psnbnd = 1.007
-
 ; Tension of the density spline
 ; Note: 0 means perfect interpolation
 sspld = 0.1
diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90
index 615cc11..888b132 100644
--- a/src/coreprofiles.f90
+++ b/src/coreprofiles.f90
@@ -1,136 +1,228 @@
+! This modules handles the loading, interpolation and evaluation of the
+! plasma profiles (density, temperature, effective charge)
+!
+! Two kinds of profiles are supported: analytical (suffix `_an` in the
+! subroutine names) or numerical (suffix `_spline`). For the latter, the
+! the data is interpolated using splines.
 module coreprofiles
   use const_and_precisions, only : wp_, zero, one
+
   implicit none
 
-  integer, save :: npp, nsfd
-  real(wp_), save :: psdbnd, psnpp, denpp, ddenpp, d2denpp
-  real(wp_), dimension(:),  allocatable,  save :: tfn, cfn, psrad
-  real(wp_), dimension(:, :),  allocatable,  save :: ct, cz
-  real(wp_), save :: dens0, aln1, aln2, te0, dte0, alt1, alt2, zeffan
+  ! Parameters of the plasma profiles splines
+  type spline_parameters
+    integer :: ndata   ! Number of data points
+    integer :: nknots  ! Number of spline knots
+    ! Density spline (ψ, knots, B-spline coefficients)
+    real(wp_), dimension(:), allocatable :: knots, coeffs
+    ! Temperature and effective charge arrays (ψ, T(ψ), Zeff(ψ))
+    real(wp_), dimension(:),    allocatable :: psi
+    real(wp_), dimension(:, :), allocatable :: temp, zeff
+  end type
+
+  ! Parameters of the C² polynomial tail of the density spline
+  type density_tail
+    real(wp_) :: start   ! ψ₀, start of the tail
+    real(wp_) :: end     ! ψ₁, end of the end
+    real(wp_) :: value   ! s(ψ₀), value at the start
+    real(wp_) :: deriv1  ! s'(ψ₀), first derivative at the start
+    real(wp_) :: deriv2  ! s"(ψ₀), second derivative at the start
+  end type
+
+  ! Parameters of the analytical profiles model
+  type analytic_model
+    real(wp_) :: dens0     ! Density scaling factor
+    real(wp_) :: n1, n2    ! Density exponents
+    real(wp_) :: te0, te1  ! Temperature at ψ=0, ψ=1
+    real(wp_) :: t1, t2    ! Temperature exponents
+    real(wp_) :: zeff      ! Effective charge
+  end type
+
+  ! Global variable storing the state of the module
+  type(spline_parameters), save :: spline
+  type(density_tail),      save :: tail
+  type(analytic_model),    save :: model
+
+  private
+  public read_profiles, read_profiles_an       ! Reading data files
+  public scale_profiles                        ! Applying rescaling
+  public density, temp, fzeff                  ! Accessing interpolated data
+  public set_profiles_spline, set_profiles_an  ! Initialising internal state
+  public unset_profiles_spline                 ! Deinitialising internal state
 
 contains
 
-  subroutine density(psin,dens,ddens)
+  subroutine density(psin, dens, ddens)
+    ! Computes the density its first derivative as a function of
+    ! normalised poloidal flux.
+    !
+    ! Note: density has units of 10¹⁹ m⁻³.
     use gray_params, only : iprof
-    use dierckx,     only : splev,splder
+    use dierckx,     only : splev, splder
     use logger,      only : log_error
 
     implicit none
 
     ! subroutine arguments
-    real(wp_), intent(in) :: psin
-    real(wp_), intent(out) :: dens,ddens
+    real(wp_), intent(in)  :: psin         ! normalised poloidal flux
+    real(wp_), intent(out) :: dens, ddens  ! density and first derivative
 
     ! local variables
-    integer :: ier,nu
-    real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh
-    real(wp_), dimension(1) :: xxs,ffs
-    real(wp_), dimension(npp+4) :: wrkfd
-    character(256) :: msg
+    integer :: ier                      ! dierck error code
+    real(wp_) :: f(1)                   ! dierck output (must be an array)
+    real(wp_) :: wrkfd(spline%ndata+4)  ! dierck working space array
+    character(256) :: msg               ! for log messages formatting
 
-    !
-    ! Computation of density [10¹⁹ m⁻³] and derivative wrt ψ
-    !
-    dens=zero
-    ddens=zero
-    if((psin >= psdbnd) .or. (psin < zero)) return
+    ! Initialise both to zero
+    dens  = zero
+    ddens = zero
 
-    if(iprof == 0) then
-      if(psin > one) return
-      profd=(one-psin**aln1)**aln2
-      dens=dens0*profd
-      dprofd=-aln1*aln2*psin**(aln1-one)           &
-              *(one-psin**aln1)**(aln2-one)
-      ddens=dens0*dprofd
-    else
-      if(psin > psnpp) then
+    ! Outside the tail end both density and its
+    ! derivatives are identically zero
+    if (psin >= tail%end .or. psin < 0) return
 
-      ! Smooth interpolation for psnpp < psi < psdbnd
-      ! dens = fp * fh
-      ! fp: parabola matched at psi=psnpp with given profile density
-      ! fh=(1-t)^3(1+3t+6t^2) is a smoothing function:
-      ! fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1
+    if (iprof == 0) then
+      ! Use the analytical model
       !
-        dpsib=psin-psnpp
-        fp=denpp+dpsib*ddenpp+0.5_wp_*dpsib**2*d2denpp
-        dfp=ddenpp+dpsib*d2denpp
-        tt=dpsib/(psdbnd-psnpp)
-        fh=(one-tt)**3*(one+3.0_wp_*tt+6.0_wp_*tt**2)
-        dfh=-30.0_wp_*(one-tt)**2*tt**2/(psdbnd-psnpp)
-        dens=fp*fh
-        ddens=dfp*fh+fp*dfh
+      !   n(ψ) = dens0⋅(1 - ψ^aln1)^aln2
+      !
+      if (psin > 1) return
+      dens  =  model%dens0 * (1 - psin**model%n1)**model%n2
+      ddens = -model%dens0 * model%n1*model%n2 * psin**(model%n1 - 1)  &
+                           * (1 - psin**model%n1)**(model%n2 - 1)
+    else
+      ! Use the numerical data
+
+      if (psin < tail%start) then
+        ! Use the interpolating spline when in range
+
+        ! Evaluate the spline
+        ier = 0
+        call splev(spline%knots, spline%nknots, spline%coeffs, &
+                   3, [psin], f, 1, ier)
+        dens = f(1)
+
+        ! Evaluate the spline 1st derivative
+        ier = 0
+        call splder(spline%knots, spline%nknots, spline%coeffs, &
+                    3, 1, [psin], f, 1, wrkfd, ier)
+        ddens = f(1)
+        if (abs(dens) < 1.0e-10_wp_) dens = zero
       else
-        xxs(1)=psin
-        ier=0
-        call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
-        dens=ffs(1)
-        nu=1
-        ier=0
-        call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
-        ddens=ffs(1)
-        if(abs(dens) < 1.0e-10_wp_) dens=zero
+        ! Use a C² polynomial extension outside (ψ > ψ₀)
+
+        ! The tail consists of the product p(ψ)⋅h(t), where:
+        !
+        ! - p(ψ) is the 2nd order Taylor polynomial of the spline,
+        !   centered at ψ₀. See set_profiles_spline for details.
+        !
+        ! - h(t) is a "smoothing" polynomial in the variable
+        !   t = (ψ - ψ₀)/(ψ₁ - ψ₀), defined as:
+        !
+        !     h(t) = (1 - t)³⋅(1 + 3t + 6t²)
+        !
+        !   with the following properties:
+        !
+        !     h(0) = 1   h'(0)=0   h"(0)=0
+        !     h(1) = 0   h'(1)=0   h"(1)=0
+        block
+          real(wp_) :: dpsi, t, p, dp, h, dh
+          dpsi = psin - tail%start  ! Δψ = (ψ - ψ₀)
+
+          ! Taylor polynomial p(ψ) and its derivative
+          p  = tail%value + dpsi*tail%deriv1 + dpsi**2*tail%deriv2/2
+          dp = tail%deriv1 + dpsi*tail%deriv2
+
+          ! Smoothing polynomial h(t) and its derivative
+          t  = dpsi/(tail%end - tail%start)
+          h  = (1 - t)**3 * (1 + 3*t + 6*t**2)
+          dh = -30*(1 - t)**2 * t**2 / (tail%end - tail%start)
+
+          dens  = p*h
+          ddens = dp*h + p*dh
+        end block
       end if
-      if(dens < zero) then
+
+      if (dens < 0) then
         write (msg, '("negative density:", 2(x,a,"=",g0.3))') &
           'ne', dens, 'ψ', psin
         call log_error(msg, mod='coreprofiles', proc='density')
       end if
+
     end if
   end subroutine density
 
 
   function temp(psin)
-    use const_and_precisions, only : wp_,zero,one
-    use gray_params, only : iprof
-    use utils, only : locate
-    use simplespline, only :spli
+    ! Computes the temperature as a function of the
+    ! normalised poloidal flux.
+    !
+    ! Note: temperature has units of keV.
+    use gray_params,  only : iprof
+    use utils,        only : locate
+    use simplespline, only : spli
+
     implicit none
-! arguments
+
+    ! subroutine arguments
     real(wp_), intent(in) :: psin
     real(wp_) :: temp
-! local variables
-    integer :: k
-    real(wp_) :: proft,dps
 
-    temp=zero
-    if((psin >= one).or.(psin < zero)) return
-    if(iprof == 0) then
-      proft=(1.0_wp_-psin**alt1)**alt2
-      temp=(te0-dte0)*proft+dte0
+    ! local variables
+    integer :: k
+    real(wp_) :: proft, dps
+
+    temp = zero
+    if (psin >= 1 .or. psin < 0) return
+    if (iprof == 0) then
+      ! Use the analytical model
+      !
+      !   T(ψ) = (te0 - te1)⋅(1 - ψ^t1)^t2 + te1
+      !
+      proft = (1 - psin**model%t1)**model%t2
+      temp = (model%te0 - model%te1)*proft + model%te1
     else
-      call locate(psrad,npp,psin,k)
-      k=max(1,min(k,npp-1))
-      dps=psin-psrad(k)
-      temp=spli(ct,npp,k,dps)
+      ! Use the interpolated numerical data
+      call locate(spline%psi, spline%ndata, psin, k)
+      k = max(1, min(k, spline%ndata - 1))
+      dps = psin - spline%psi(k)
+      temp = spli(spline%temp, spline%ndata, k, dps)
     endif
   end function temp
 
+
   function fzeff(psin)
-    use const_and_precisions, only : wp_,zero,one
-    use gray_params, only : iprof
-    use utils, only : locate
-    use simplespline, only :spli
+    ! Computes the effective charge Zeff as a
+    ! function of the normalised poloidal flux.
+    use gray_params,  only : iprof
+    use utils,        only : locate
+    use simplespline, only : spli
+
     implicit none
-! arguments
+
+    ! subroutine arguments
     real(wp_), intent(in) :: psin
     real(wp_) :: fzeff
-! local variables
+
+    ! local variables
     integer :: k
     real(wp_) :: dps
 
-    fzeff=one
-    if((psin >= one).or.(psin < zero)) return
-    if(iprof == 0) then
-      fzeff=zeffan
+    fzeff = one
+    if (psin >= 1 .or. psin < 0) return
+    if (iprof == 0) then
+      ! Use the analytical model (just a constant)
+      fzeff = model%zeff
     else
-      call locate(psrad,npp,psin,k)
-      k=max(1,min(k,npp-1))
-      dps=psin-psrad(k)
-      fzeff=spli(cz,npp,k,dps)
+      ! Use the interpolated numerical data
+      call locate(spline%psi, spline%ndata, psin, k)
+      k = max(1, min(k, spline%ndata - 1))
+      dps = psin - spline%psi(k)
+      fzeff = spli(spline%zeff, spline%ndata, k, dps)
     endif
   end function fzeff
 
-  
+
   subroutine read_profiles(filenm, data, unit)
     ! Reads the radial plasma profiles from `file` and store them
     ! into `data`. If given, the file is opened in the `unit` number.
@@ -154,10 +246,10 @@ contains
     integer :: err
 
     ! Free the arrays when already allocated
-    if(allocated(data%psrad)) deallocate(data%psrad)
-    if(allocated(data%terad)) deallocate(data%terad)
-    if(allocated(data%derad)) deallocate(data%derad)
-    if(allocated(data%zfc))   deallocate(data%zfc)
+    if (allocated(data%psrad)) deallocate(data%psrad)
+    if (allocated(data%terad)) deallocate(data%terad)
+    if (allocated(data%derad)) deallocate(data%derad)
+    if (allocated(data%zfc))   deallocate(data%zfc)
 
     u = get_free_unit(unit)
 
@@ -190,7 +282,14 @@ contains
     ! from params%filenm.
     ! If given, the file is opened in the `unit` number.
     !
-    ! TODO: add format description
+    ! The file should be formatted as follows:
+    !
+    !   1 dens0 n1 n2
+    !   2 te0 te1 t1 t2
+    !   3 zeff
+    !
+    ! See `density`, `temp`, `fzeff` subroutines for the meaning
+    ! of the parameters (i.e. the formulae for n,T,Zeff).
     use gray_params, only : profiles_data
     use utils,       only : get_free_unit
     use logger,      only : log_error
@@ -205,7 +304,7 @@ contains
     ! local variables
     integer :: u
     integer :: err
-        
+
     u = get_free_unit(unit)
 
     if (allocated(data%terad)) deallocate(data%terad)
@@ -220,9 +319,9 @@ contains
       call exit(1)
     end if
 
-    read (u,*) data%derad(1:3)  ! dens0, aln1, aln2
-    read (u,*) data%terad(1:4)  ! te0, dte0, alt1, alt2
-    read (u,*) data%zfc(1)      ! zeffan
+    read (u,*) data%derad(1:3)  ! dens0, n1, n2
+    read (u,*) data%terad(1:4)  ! te0, te1, t1, t2
+    read (u,*) data%zfc(1)      ! zeff
     close(u)
   end subroutine read_profiles_an
 
@@ -253,13 +352,13 @@ contains
       aan = one
     end if
 
-    if(params%iscal==2) then
+    if (params%iscal==2) then
       ffact = one
     else
       ffact = factb
     end if
 
-    if(params%iprof==0) then
+    if (params%iprof==0) then
       last_te = 2
       last_ne = 1
     else
@@ -272,140 +371,227 @@ contains
   end subroutine scale_profiles
 
 
-  subroutine set_profiles_spline(params, data)
+  subroutine set_profiles_spline(params, data, launch_pos)
     ! Computes splines for the plasma profiles data and stores them
     ! in their respective global variables, see the top of this file.
+    !
+    ! When `launch_pos` (cartesian launch coordinates in cm) is present,
+    ! the subroutine will also check that the wave launcher is strictly
+    ! outside the reconstructed plasma density boundary.
     use simplespline, only : difcs
     use dierckx,      only : curfit, splev, splder
     use gray_params,  only : profiles_parameters, profiles_data
-    use logger,       only : log_info, log_warning
+    use logger,       only : log_debug, log_info, log_warning, log_error
 
     implicit none
 
     ! subroutine arguments
-    type(profiles_parameters), intent(in)    :: params
+    type(profiles_parameters), intent(inout) :: params
     type(profiles_data),       intent(inout) :: data
+    real(wp_), optional,       intent(in)    :: launch_pos(3)
+
+    ! curfit parameters
+    integer, parameter :: iopt = 0  ! smoothing spline mode
+    integer, parameter :: kspl = 3  ! order of spline (cubic)
 
     ! local variables
-    integer, parameter :: iopt=0, kspl=3
-    integer :: n, npest, lwrkf, ier
-    real(wp_) :: xb, xe, fp, xnv, xxp,xxm,delta2,ssplne_loc
-    real(wp_), dimension(:), allocatable :: wf, wrkf
-    integer, dimension(:), allocatable :: iwrkf
-    real(wp_), dimension(1) :: dedge,ddedge,d2dedge
-    character(256) :: msg ! for log messages formatting
+    integer :: n, npest, ier
+    real(wp_) :: xb, xe, fp, ssplne_loc
 
-    n=size(data%psrad)
-    npest=n+4
-    lwrkf=n*4+npest*16
-    allocate(wrkf(lwrkf),iwrkf(npest),wf(n))
+    ! working space arrays for the dierckx functions
+    integer :: lwrkf
+    real(wp_), dimension(:), allocatable :: wf, wrkf
+    integer,   dimension(:), allocatable :: iwrkf
+
+    ! for log messages formatting
+    character(256) :: msg
+
+    n = size(data%psrad)
+    npest = n + 4
+    lwrkf = n*4 + npest*16
+    allocate(wrkf(lwrkf), iwrkf(npest), wf(n))
 
     ssplne_loc=params%sspld
 
-    ! if necessary, reallocate spline arrays
-    if(.not.allocated(psrad)) then
-      allocate(psrad(n),ct(n,4),cz(n,4))
+    ! If necessary, reallocate the spline arrays
+    if (.not. allocated(spline%psi)) then
+      allocate(spline%psi(n), spline%temp(n, 4), spline%zeff(n, 4))
     else
-      if(size(psrad)<n) then
-        deallocate(psrad,ct,cz)
-        allocate(psrad(n),ct(n,4),cz(n,4))
+      if (size(spline%psi) < n) then
+        deallocate(spline%psi, spline%temp, spline%zeff)
+        allocate(spline%psi(n), spline%temp(n, 4), spline%zeff(n, 4))
       end if
     end if
-    if(.not.allocated(cfn)) then
-      allocate(tfn(npest),cfn(npest))
+    if (.not. allocated(spline%coeffs)) then
+      allocate(spline%knots(npest), spline%coeffs(npest))
     else
-      if(size(cfn)<npest) then
-        deallocate(tfn,cfn)
-        allocate(tfn(npest),cfn(npest))
+      if (size(spline%coeffs) < npest) then
+        deallocate(spline%knots, spline%coeffs)
+        allocate(spline%knots(npest), spline%coeffs(npest))
       end if
     end if
 
-    ! spline approximation of  temperature and data%zfc
-    call difcs(data%psrad,data%terad,  n,iopt,ct,ier)
-    call difcs(data%psrad,data%zfc,n,iopt,cz,ier)
-    psrad=data%psrad
-    npp=n
+    ! Spline interpolation of temperature and effective charge
+    call difcs(data%psrad, data%terad, n, iopt, spline%temp, ier)
+    call difcs(data%psrad, data%zfc,   n, iopt, spline%zeff, ier)
+    spline%psi = data%psrad
+    spline%ndata = n
 
-    ! spline approximation of density
-    xb=zero
-    xe=data%psrad(n)
-    wf(:)=one
-    call curfit(iopt,n,data%psrad,data%derad,wf,xb,xe,kspl,ssplne_loc,npest, &
-       nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
+    ! Spline interpolation of density
+    xb = zero
+    xe = data%psrad(n)
+    wf(:) = one
+    call curfit(iopt, n, data%psrad, data%derad, wf, xb, xe, kspl, &
+                ssplne_loc, npest, spline%nknots, spline%knots,    &
+                spline%coeffs, fp,  wrkf, lwrkf, iwrkf, ier)
     ! if ier=-1 data are re-fitted using sspl=0
-    if(ier == -1) then
+    if (ier == -1) then
         call log_warning('curfit failed with error -1: re-fitting with '// &
                          's=0', mod='coreprofiles', proc='density')
-        ssplne_loc=0.0_wp_
-        call curfit(iopt,n,data%psrad,data%derad,wf,xb,xe,kspl,ssplne_loc,npest, &
-           nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
+        ssplne_loc = zero
+        call curfit(iopt, n, data%psrad, data%derad, wf, xb, xe, kspl, &
+                    ssplne_loc, npest, spline%nknots, spline%knots,    &
+                    spline%coeffs, fp, wrkf, lwrkf, iwrkf, ier)
     end if
 
-    ! compute polinomial extrapolation matching the spline boundary up to the
-    ! 2nd order derivative, extending the profile up to psi=psdbnd where
-    ! data%derad=data%derad'=data%derad''=0
-    ! spline value and derivatives at the edge
-    call splev(tfn,nsfd,cfn,kspl,data%psrad(n:n),dedge(1:1),1,ier)
-    call splder(tfn,nsfd,cfn,kspl,1,data%psrad(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier)
-    call splder(tfn,nsfd,cfn,kspl,2,data%psrad(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier)
-    ! determination of the boundary
-    psdbnd=params%psnbnd
+    ! Computation of the polynomial tail parameters
+    !
+    ! Note: The density is the only quantity that needs to be evaluated
+    ! at the edge. The spline thus has to be extended to transition
+    ! smoothly from the last profile point to 0 outside the plasma.
+    block
+      real(wp_), dimension(1) :: s0, s1, s2  ! spline, 1st, 2nd derivative
+      real(wp_), dimension(1) :: delta4      ! discriminant Δ/4 of q(x)
+      real(wp_), dimension(1) :: x0, x1      ! vertex of q(x), solution
 
-    psnpp=data%psrad(n)
-    denpp=dedge(1)
-    ddenpp=ddedge(1)
-    d2denpp=d2dedge(1)
+      ! Compute the coefficients of a 2nd order Taylor polinomial to
+      ! extend the spline beyond the last point:
+      !
+      !   p(ψ) = s(ψ₀) + (ψ - ψ₀)s'(ψ₀) + ½(ψ - ψ₀)²s"(ψ₀)
+      !
+      !  where s(ψ) is the spline and ψ₀ the last point.
+      !
+      call splev(spline%knots, spline%nknots, spline%coeffs, kspl, &
+                 data%psrad(n:n), s0, 1, ier)
+      call splder(spline%knots, spline%nknots, spline%coeffs, kspl, 1, &
+                  data%psrad(n:n), s1, 1, wrkf(1:spline%nknots), ier)
+      call splder(spline%knots, spline%nknots, spline%coeffs, kspl, 2, &
+                  data%psrad(n:n), s2, 1, wrkf(1:spline%nknots), ier)
 
-    delta2=(ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp
-    xnv=psnpp-ddenpp/d2denpp
-    if(delta2 < zero) then
-    ! if(xnv > psnpp) psdbnd=min(psdbnd,xnv)
-    else
-      xxm=xnv-sqrt(delta2)
-      xxp=xnv+sqrt(delta2)
-      if(xxm > psnpp) then
-        psdbnd=min(psdbnd,xxm)
-      else if (xxp > psnpp) then
-        psdbnd=min(psdbnd,xxp)
+      ! Determine where to end the tail (to ensure the density remains
+      ! positive) from the zeros of the Taylor polynomial p(ψ)
+      !
+      ! Define x=(ψ - ψ₀), then p(ψ)=0 is rewritten as
+      !
+      !   q(x) = x² + 2s'/s" x + 2s/s" = 0
+      !
+      ! The discriminant is Δ/4 = (s'/s")² - 2(s/s") and
+      ! the solutions are x = x₀ ± √(Δ/4), with x₀ = -s'/s".
+      !
+      x0 = -s1 / s2                    ! vertex of parabola y=q(x)
+      delta4 = (s1 / s2)**2 - 2*s0/s2  ! Δ/4 of q(x)
+
+      if (delta4(1) > 0) then
+        ! Pick the smallest positive zero (implying >ψ₀)
+        x1 = x0 + sign(sqrt(delta4), sqrt(delta4) - x0)
+      else
+        ! There are no zeros, use the parabola vertex
+        x1 = x0
+        call log_debug('spline extension has no zeros', &
+                       mod='coreprofiles', proc='set_profiles_spline')
       end if
-      write (msg, '(a,g0.3)') 'density boundary: ψ=', psdbnd
-      call log_info(msg, mod="coreprofiles", proc="set_profiles_spline")
+
+      ! Store the tail parameters
+      tail%start  = data%psrad(n)
+      tail%end    = tail%start + x1(1)
+      tail%value  = s0(1)
+      tail%deriv1 = s1(1)
+      tail%deriv2 = s2(1)
+    end block
+
+    ! Make sure the wave launcher does not fall inside the tail
+    ! Note: if it does, the initial wave conditions become
+    ! invalid as they are given assuming a vacuum (N=1)
+    if (present(launch_pos)) then
+      block
+        use equilibrium, only : equinum_psi
+        real(wp_) :: R, Z, psi
+        real(wp_), parameter :: cm = 1.0e-2_wp_
+
+        ! Convert from cartesian to cylindrical coordinates
+        R = hypot(launch_pos(1), launch_pos(2)) * cm  ! R = √(x²+y²)
+        z = launch_pos(3) * cm
+
+        ! Get the poloidal flux at the launcher
+        ! Note: this returns -1 when the data is not available
+        call equinum_psi(R, z, psi)
+
+        if (psi > tail%start .and. psi < tail%end) then
+          ! Fall back to the midpoint of ψ₀ and the launcher ψ
+          tail%end = (tail%start + psi)/2
+            call log_warning('downscaled tail to not reach the wave launcher', &
+                             mod='coreprofiles', proc='set_profiles_spline')
+        end if
+
+        if (psi > 0 .and. psi < tail%start) then
+          ! This must be a user error, stop here
+          write (msg, '(a, a, g0.3, a, g0.3)') &
+            'wave launcher is inside the plasma! ', &
+            'launcher: ψ=', psi, ' boundary: ψ=', tail%end
+          call log_error(msg, mod='coreprofiles', proc='set_profiles_spline')
+          call exit(2)
+        end if
+      end block
     end if
 
-    deallocate(iwrkf,wrkf,wf)
+    ! Set the density boundary ψ
+    ! Note: this is used to detect entrance in the plasma
+    params%psnbnd = tail%end
+
+    write (msg, '(a,g0.4)') 'density boundary: ψ=', tail%end
+    call log_info(msg, mod='coreprofiles', proc='set_profiles_spline')
+
+    deallocate(iwrkf, wrkf, wf)
   end subroutine set_profiles_spline
 
 
   subroutine unset_profiles_spline
+    ! Unsets the splines global variables, see the top of this file.
+
     implicit none
 
-    if (allocated(psrad)) deallocate(psrad)
-    if (allocated(ct))    deallocate(ct)
-    if (allocated(cz))    deallocate(cz)
-    if (allocated(tfn))   deallocate(tfn)
-    if (allocated(cfn))   deallocate(cfn)
+    if (allocated(spline%psi))     deallocate(spline%psi)
+    if (allocated(spline%temp))    deallocate(spline%temp)
+    if (allocated(spline%zeff))    deallocate(spline%zeff)
+    if (allocated(spline%knots))   deallocate(spline%knots)
+    if (allocated(spline%coeffs))  deallocate(spline%coeffs)
   end subroutine unset_profiles_spline
 
 
-  subroutine set_profiles_an(data)
+  subroutine set_profiles_an(params, data)
     ! Stores the analytical profiles data in their respective
     ! global variables, see the top of this file.
-    use gray_params, only : profiles_data
+    use gray_params, only : profiles_parameters, profiles_data
 
     implicit none
 
     ! subroutine arguments
-    type(profiles_data), intent(in) :: data
+    type(profiles_parameters), intent(inout) :: params
+    type(profiles_data),       intent(in)    :: data
+
+    model%te0   = data%terad(1)
+    model%te1   = data%terad(2)
+    model%t1    = data%terad(3)
+    model%t2    = data%terad(4)
+    model%dens0 = data%derad(1)
+    model%n1    = data%derad(2)
+    model%n2    = data%derad(3)
+    model%zeff  = data%zfc(1)
+
+    ! Define the plasma boundary to be exactly ψ=1
+    ! Note: this is used to detect entrance in the plasma
+    params%psnbnd = one
 
-    te0    = data%terad(1)
-    dte0   = data%terad(2)
-    alt1   = data%terad(3)
-    alt2   = data%terad(4)
-    dens0  = data%derad(1)
-    aln1   = data%derad(2)
-    aln2   = data%derad(3)
-    zeffan = data%zfc(1)
-    psdbnd = one
   end subroutine set_profiles_an
-  
+
 end module coreprofiles
diff --git a/src/main.f90 b/src/main.f90
index e46e96d..275b586 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -63,8 +63,9 @@ program main
   ! Read the input data and set the global variables
   ! of the respective module. Note: order matters.
   call init_equilibrium(params, data)
-  call init_profiles(params%profiles, params%equilibrium%factb, data%profiles)
   call init_antenna(params%antenna)
+  call init_profiles(params%profiles, params%equilibrium%factb, &
+                     params%antenna%pos, data%profiles)
   call init_misc(params, data)
 
   ! Change the current directory to output files there
@@ -241,7 +242,7 @@ contains
   end subroutine deinit_equilibrium
 
 
-  subroutine init_profiles(params, factb, data)
+  subroutine init_profiles(params, factb, launch_pos, data)
     ! Reads the plasma kinetic profiles file (containing the elecron
     ! temperature, density and plasma effective charge) and initialises
     ! the respective GRAY data structure.
@@ -253,9 +254,10 @@ contains
     implicit none
 
     ! subroutine arguments
-    type(profiles_parameters), intent(in)  :: params
-    real(wp_),                 intent(in)  :: factb
-    type(profiles_data),       intent(out) :: data
+    type(profiles_parameters), intent(inout) :: params
+    real(wp_),                 intent(in)    :: factb
+    real(wp_),                 intent(in)    :: launch_pos(3)
+    type(profiles_data),       intent(out)   :: data
 
     if (params%iprof == 0) then
       ! Analytical profiles
@@ -280,10 +282,10 @@ contains
     ! Set global variables
     if (params%iprof == 0) then
       ! Analytical profiles
-      call set_profiles_an(data)
+      call set_profiles_an(params, data)
     else
       ! Numerical profiles
-      call set_profiles_spline(params, data)
+      call set_profiles_spline(params, data, launch_pos)
     end if
   end subroutine init_profiles