Carpenteri0 · Aug 10, 2026 · Aug 10, 2026 · Aug 20, 2026 · Aug 20, 2026 · Aug 20, 2026
diff --git a/communication/broadcast_phys.f90 b/communication/broadcast_phys.f90
@@ -924,7 +924,6 @@ subroutine broadcast_phys(my_id)
   call MPI_PACK(part_group_configs%homma2020_alpha,            n_part_groups_max,       MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr)
   call MPI_PACK(part_group_configs%ics_group_idx,              n_part_groups_max,     MPI_INTEGER,buffer,bufsize,position,MPI_COMM_WORLD,ierr)
 
-  call MPI_PACK(part_group_configs%num_re,                     n_part_groups_max,     MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr)
   call MPI_PACK(part_group_configs%re_energy,                  n_part_groups_max,     MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr)
   call MPI_PACK(part_group_configs%re_std_energy,              n_part_groups_max,     MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr)
   call MPI_PACK(part_group_configs%re_pitch,                   n_part_groups_max,     MPI_REAL8,buffer,bufsize,position,MPI_COMM_WORLD,ierr)
@@ -1918,7 +1917,6 @@ subroutine broadcast_phys(my_id)
   call MPI_UNPACK(buffer,bufsize,position,part_group_configs%homma2020_alpha,              n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr)
   call MPI_UNPACK(buffer,bufsize,position,part_group_configs%ics_group_idx,              n_part_groups_max,MPI_INTEGER,MPI_COMM_WORLD,ierr)
 
-  call MPI_UNPACK(buffer,bufsize,position,part_group_configs%num_re,                     n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr)
   call MPI_UNPACK(buffer,bufsize,position,part_group_configs%re_energy,                  n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr)
   call MPI_UNPACK(buffer,bufsize,position,part_group_configs%re_std_energy,              n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr)
   call MPI_UNPACK(buffer,bufsize,position,part_group_configs%re_pitch,                   n_part_groups_max,MPI_REAL8,MPI_COMM_WORLD,ierr)

diff --git a/docs/physics/kinetic_models/runaway_electrons/rep_tutorial.md b/docs/physics/kinetic_models/runaway_electrons/rep_tutorial.md
@@ -29,7 +29,7 @@ Apart from the jorek_restart.h5, an part_restart.h5 file is required by any simu
  part_group_configs(1)%type                        = 'particle_kinetic_relativistic'
  part_group_configs(1)%n_particles                 = 1e6    ! number of super particles
  part_group_configs(1)%mass                        = 0.000548579870184805     
- part_group_configs(1)%num_re                      = 1.175403e16      ! number of physical REs, doesn't work if initialized according to current
+ part_group_configs(1)%n_particles_total           = 1.175403e16      ! number of physical REs, doesn't work if initialized according to current
  part_group_configs(1)%use_kin_recombination       = .false.
 
 ! Projecting (particle field -> MHD) smoothing

diff --git a/models/mod_log_params.f90 b/models/mod_log_params.f90
@@ -1148,7 +1148,7 @@ subroutine log_parameters(my_id, short)
 
       ! rep (runaway electrons, only pressure coupling for now) -----
       if (sim%groups(group_num)%coupling_scheme .eq. 'rep') then
-        write(*,REAL_FMT) 'num_re,                 ',part_group_configs(group_num)%num_re
+        write(*,REAL_FMT) 'n_particles_total,      ',part_group_configs(group_num)%n_particles_total
         write(*,REAL_FMT) 're_energy,              ',part_group_configs(group_num)%re_energy
         write(*,REAL_FMT) 're_std_energy,          ',part_group_configs(group_num)%re_std_energy
         write(*,REAL_FMT) 're_pitch,               ',part_group_configs(group_num)%re_pitch

diff --git a/models/phys_module.f90 b/models/phys_module.f90
@@ -1114,7 +1114,6 @@ module phys_module
 
     ! ================ for runaway electrons ('rep' coupling scheme) particles ===============
 
-    real*8              :: num_re                  !< number of runaway electrons in the group
     real*8              :: re_energy               !< energy [eV] of the runaway electrons in the group
     real*8              :: re_std_energy           !< standard deviation of the energy [eV] of the runaway electrons in the group
     real*8              :: re_pitch                !< pitch between RE momentum and magnetic field line (i.e. p_re_par/p_re_tot)
@@ -1125,6 +1124,8 @@ module phys_module
                                                    !< for example n_phi_planes=4, n_particles=1e4 then only 250 particles are initialised
                                                    !< each particle is then copied multiple (3) times around the torus with angle 2pi/n_phi_planes (= pi/2)
                                                    !< if n_phi_planes=int*n_period then projected particle quantities are initialised as 0 for n_tor>1
+
+    ! =============== for rep and epf ===================
     real*8              :: n_particles_total       !< Total number of particles to simulate (ie sum(weights)) !!NOT n_particles - total number of super/numeric-particles
 
 

diff --git a/models/preset_parameters.f90 b/models/preset_parameters.f90
@@ -924,15 +924,16 @@ subroutine preset_parameters
 part_group_configs(:)%homma2020_alpha        = 1.5d0
 part_group_configs(:)%ics_group_idx          = -1
 
-!----- specific to rep 
-part_group_configs(:)%num_re                 = 0.d0
+!----- specific to rep
 part_group_configs(:)%re_energy              = 0.d0
 part_group_configs(:)%re_std_energy          = 0.d0
 part_group_configs(:)%re_pitch               = 0.d0
 
 !----- specific to epf
 part_group_configs(:)%T_maxwell              = 0.d0
 part_group_configs(:)%n_phi_planes           = 0
+
+!----- specific to rep and epf
 part_group_configs(:)%n_particles_total      = 0.d0
 
 

diff --git a/particles/diagnostics/mod_neutral_density.f90 b/particles/diagnostics/mod_neutral_density.f90
@@ -16,6 +16,7 @@ module mod_neutral_density
 
 !> first calculates the neutral density from all particles with ncs, then projects it
 subroutine get_neutral_density(sim, neutral_density_proj)
+  use phys_module, only: n_part_groups
   implicit none
 
   class(particle_sim), target, intent(inout)                :: sim

diff --git a/particles/examples/E_diagnostic.f90 b/particles/examples/E_diagnostic.f90
@@ -12,6 +12,7 @@
 module write_E_time
   use mod_event, only: event, action
   use equil_info, only:find_xpoint
+  use mod_interp
   implicit none
 
   type, extends(action) :: save_E

diff --git a/particles/examples/plasma_volume.f90 b/particles/examples/plasma_volume.f90
@@ -12,7 +12,7 @@ program plasma_volume
 use mod_random_seed, only: random_seed
 use mod_sobseq_rng, only: sobseq_rng
 use constants, only: TWOPI, PI
-use mod_initialise_particles, only: domain_bounding_box
+use initialisers_base, only: domain_bounding_box
 use domains
 use hdf5_io_module
 use mod_sampling, only: transform_uniform_cylindrical

diff --git a/particles/examples/seed_particles.f90 b/particles/examples/seed_particles.f90
@@ -4,9 +4,12 @@
 program seed_particles
 use particle_tracer
 use mod_particle_io
+use mod_rej_f
+use initialisers_base
 implicit none
 
 type(event) :: fieldreader
+type(spatial_pdf) :: space_pdf
 integer :: i, j
 
 ! Start up MPI, jorek
@@ -16,6 +19,10 @@ program seed_particles
 fieldreader = event(read_jorek_fields_interp_linear(basename='jorek', i=-1))
 call with(sim, fieldreader)
 
+! Set up spatial pdf type
+space_pdf%f    => f_psi_inside
+space_pdf%vars =  [1]
+
 ! Set up particles
 sim%groups(:)%Z    = -2 
 sim%groups(:)%mass = 2.d0 !< atomic mass units
@@ -26,16 +33,16 @@ program seed_particles
   call initialise_particles_H_mu_psi(sim%groups(i)%particles, &
       sim%fields, pcg32_rng(), &
       sim%groups(i)%mass, &
-      uniform_space=.true., uniform_space_rej_f=f_psi_inside, &
-      uniform_space_rej_vars=[1], charge=1)
+      uniform_space=.true., space_pdf=space_pdf, charge=1)
   sim%groups(i)%particles(:)%weight = 1.0
 end do
 
 call write_simulation_hdf5(sim, 'part_restart.h5')
 contains
 pure function f_psi_inside(n, P, grad_P) result(f)
-  integer, intent(in) :: n
-  real*8, intent(in) :: P(n), grad_P(3,n)
+  integer,                intent(in) :: n
+  real*8, dimension(n),   intent(in) :: P
+  real*8, dimension(3,n), intent(in) :: grad_P
   real*4 :: f
   f = 5e-3
   if (P(1) .lt. -0.26 .and. P(1) .gt. -0.35) f = 1e0

diff --git a/particles/examples/test_generalised_initialisation_gc_H_mu_psi.f90 b/particles/examples/test_generalised_initialisation_gc_H_mu_psi.f90
@@ -4,6 +4,7 @@ program test_generalised_initialisation_gc_H_mu_psi
 !> initialize_particles_H_mu_psi implemented in mod_initialise_particles
 use constants, only: TWOPI,PI,ATOMIC_MASS_UNIT,EL_CHG
 use particle_tracer
+use mod_interp
 implicit none
 type(pcg32_rng) :: rng_pcg32
 type(event)     :: field_reader,particle_writer

diff --git a/particles/examples/v_ExB_hist.f90 b/particles/examples/v_ExB_hist.f90
@@ -15,7 +15,7 @@ program v_ExB_hist
 use mod_sobseq_rng, only: sobseq_rng
 use mod_random_seed, only: random_seed
 use constants, only: TWOPI
-use mod_initialise_particles, only: domain_bounding_box
+use initialisers_base, only: domain_bounding_box
 use mod_math_operators, only: cross_product
 use domains
 use equil_info, only:find_xpoint

diff --git a/particles/initialisers/initialisers_RE.f90 b/particles/initialisers/initialisers_RE.f90
@@ -5,74 +5,25 @@ module initialisers_RE
   use mod_particle_sim
   use mod_rng
   use initialisers_base
-  use constants, only: EL_CHG, ATOMIC_MASS_UNIT, SPEED_OF_LIGHT, MASS_ELECTRON, TWOPI
+  use constants, only: EL_CHG, ATOMIC_MASS_UNIT, SPEED_OF_LIGHT, MASS_ELECTRON, TWOPI, MU_ZERO
   use mod_pusher_tools, only: get_orthonormals
   use mod_coordinate_transforms, only: vector_cylindrical_to_cartesian
   use mod_model_settings
-  use phys_module, only: CENTRAL_DENSITY, CENTRAL_MASS, ATOMIC_MASS_UNIT, MU_ZERO
+  use phys_module, only: CENTRAL_DENSITY, CENTRAL_MASS
   use equil_info
+  use mod_rej_f
   implicit none
 
   contains
-
-! Quick and rough function to sample markers based on RZ-coordinates
-pure function RZ_pdf(var) result(p)
-  real*8, intent(in)  :: var(2) ! var(1)=j
-  real*8              :: p
-  real*8              :: minor_r
-  real*8              :: R_ax, Z_ax
-
-  R_ax = ES%R_axis
-  Z_ax = ES%Z_axis
-
-  minor_r = sqrt((var(1)-Z_ax)**2 + (var(2)-R_ax)**2)
-
-  p = 1 / (1 + minor_r)**2
-end function RZ_pdf
-
-pure function analytical_pdf(var) result(p)
-  real*8, intent(in)  :: var(2) ! var(1)=j
-  real*8              :: p
-  real*8              :: minor_r
-  real*8              :: nu
-  real*8              :: R_ax, Z_ax
-  real*8              :: LCFS_a
-
-  R_ax   = ES%R_axis
-  Z_ax   = ES%Z_axis
-  LCFS_a = ES%LCFS_a
-  nu     = 2.d0
-
-  minor_r = sqrt((var(1)-Z_ax)**2 + (var(2)-R_ax)**2)
-
-  p = (1.d0 - (minor_r/LCFS_a)**2)**nu
-
-end function analytical_pdf
-
-! Quick and rough function to sample markers proportionally to toroidal current density
-pure function current_pdf(var) result(p)
-  real*8, intent(in)  :: var(2) ! var(2)=j, var(1) = R
-  !real*8, intent(in)  :: var(1) ! var(1)=j
-  real*8              :: p
-  real*8              :: jzmin, jzmax 
-
-  !> temporarily hard coded, but should be able to be obtained from fluid restart file
-  jzmax = 3.0 / 10.0 !1.173 / 10
-  jzmin = 0.0001239 / 11.0 !0.0003166 / 11
-
-  p = (var(2)/var(1)-jzmin)/(jzmax-jzmin)
-
-end function current_pdf
-
-subroutine basic_initialization(sim, group_num, rng, init_pdf, energy, pitch, std_energy)
+subroutine initialise_re_gaussian(sim, group_num, rng, space_pdf, energy, pitch, std_energy)
   use phys_module, only: tstep_particles
   use mod_kinetic_relativistic
   use mod_sampling, only: boxmueller_transform
 
   type(particle_sim),                   intent(inout) :: sim
   integer,                              intent(in)    :: group_num
   class(type_rng),                      intent(in)    :: rng
-  character(len=50),                    intent(in)    :: init_pdf
+  type(spatial_pdf),    optional,       intent(in)    :: space_pdf
   real*8,                               intent(in)    :: energy, pitch ! Kinetic energy in units of eV and pitch
   real*8,               optional,       intent(in)    :: std_energy
   real*8,               allocatable                   :: p_tot(:), p_par(:), p_perp(:)
@@ -83,20 +34,8 @@ subroutine basic_initialization(sim, group_num, rng, init_pdf, energy, pitch, st
   integer                                             :: num_part
   real*8,               allocatable                   :: ran_uniform(:), ran_gaussian(:)
 
-  select case (trim(init_pdf))
-    case ("RZ")
-      call initialise_particles(sim%groups(group_num)%particles, sim%fields%node_list, sim%fields%element_list, rng, variables=[-2,-1], transform=RZ_pdf)
-    case ("current")
-      call initialise_particles(sim%groups(group_num)%particles, sim%fields%node_list, sim%fields%element_list, rng, variables=[-1,var_zj], transform=current_pdf) 
-    case ("analytical")
-      call initialise_particles(sim%groups(group_num)%particles, sim%fields%node_list, sim%fields%element_list, rng, variables=[-2,-1], transform=analytical_pdf)
-    case default
-      if (sim%my_id == 0) then
-        write(*,*) "ERROR: ", trim(init_pdf), " is not a valid pdf/transform function for "
-        write(*,*) "  for group '", sim%groups(group_num)%id, "' when using the 'basic_initialization' function" 
-        endif
-      stop 1
-  end select
+  call initialise_particles(sim%groups(group_num)%particles, &
+    sim%fields%node_list, sim%fields%element_list, rng, space_pdf)
 
   num_part = size(sim%groups(group_num)%particles,1) 
 
@@ -156,6 +95,6 @@ subroutine basic_initialization(sim, group_num, rng, init_pdf, energy, pitch, st
   deallocate(p_par)
   deallocate(p_perp)
 
-end subroutine basic_initialization
+end subroutine initialise_re_gaussian
 
-end module initialisers_RE
+end module initialisers_RE