init_particle_batch_from_config Subroutine

  subroutine init_particle_batch_from_config(cfg, batch_idx, pcls, state, mesh, photo_emission_dq, mpi_rank, mpi_size, mpi)
    type(app_config), intent(in) :: cfg
    integer(i32), intent(in) :: batch_idx
    type(particles_soa), intent(out) :: pcls
    type(injection_state), intent(inout), optional :: state
    type(mesh_type), intent(in), optional :: mesh
    real(dp), intent(out), optional :: photo_emission_dq(:)
    integer(i32), intent(in), optional :: mpi_rank, mpi_size
    type(mpi_context), intent(in), optional :: mpi

    integer(i32) :: s, i, batch_n, max_rank, out_idx, local_rank, n_ranks, global_count
    integer(i32), allocatable :: counts_max(:), counts_actual(:), species_cursor(:), species_id(:), emit_elem_species(:, :)
    real(dp), allocatable :: vmin_normal(:), barrier_normal(:), effective_density_m3(:), w_effective(:), &
                             effective_temperature_k(:), effective_drift_velocity(:, :), &
                             photo_emit_current_density(:), photo_vmin_normal(:), photo_normal_drift_speed(:)
    logical, allocatable :: apply_barrier_energy_shift(:)
    real(dp), allocatable :: x_species(:, :, :), v_species(:, :, :), w_species(:, :)
    real(dp), allocatable :: x(:, :), v(:, :), q(:), m(:), w(:)
    type(sheath_injection_context) :: sheath_ctx

    if (cfg%sim%batch_count <= 0_i32) error stop 'sim.batch_count must be > 0.'
    if (batch_idx < 1_i32 .or. batch_idx > cfg%sim%batch_count) then
      error stop 'Requested batch index is out of range.'
    end if
    call resolve_parallel_rank_size(local_rank, n_ranks, mpi_rank, mpi_size, mpi, 'init_particle_batch_from_config')
    if (present(state)) then
      if (.not. allocated(state%macro_residual)) error stop 'injection_state is not initialized.'
      if (size(state%macro_residual) < cfg%n_particle_species) error stop 'injection_state size mismatch.'
    end if
    if (present(photo_emission_dq)) then
      if (.not. present(mesh)) error stop 'photo_emission_dq requires mesh in init_particle_batch_from_config.'
      if (size(photo_emission_dq) /= mesh%nelem) error stop 'photo_emission_dq size mismatch.'
      photo_emission_dq = 0.0d0
    end if

    allocate (counts_max(cfg%n_particle_species), counts_actual(cfg%n_particle_species))
    allocate (vmin_normal(cfg%n_particle_species), barrier_normal(cfg%n_particle_species))
    allocate (effective_density_m3(cfg%n_particle_species), w_effective(cfg%n_particle_species))
    allocate (effective_temperature_k(cfg%n_particle_species), effective_drift_velocity(3, cfg%n_particle_species))
    allocate (photo_emit_current_density(cfg%n_particle_species), photo_vmin_normal(cfg%n_particle_species))
    allocate (photo_normal_drift_speed(cfg%n_particle_species), apply_barrier_energy_shift(cfg%n_particle_species))
    counts_max = 0_i32
    counts_actual = 0_i32
    vmin_normal = 0.0d0
    barrier_normal = 0.0d0
    effective_density_m3 = 0.0d0
    w_effective = 0.0d0
    effective_temperature_k = 0.0d0
    effective_drift_velocity = 0.0d0
    photo_emit_current_density = 0.0d0
    photo_vmin_normal = 0.0d0
    photo_normal_drift_speed = 0.0d0
    apply_barrier_energy_shift = .true.
    do s = 1, cfg%n_particle_species
      if (.not. cfg%particle_species(s)%enabled) cycle
      select case (trim(lower_ascii(cfg%particle_species(s)%source_mode)))
      case ('volume_seed')
        w_effective(s) = cfg%particle_species(s)%w_particle
        effective_temperature_k(s) = species_temperature_k(cfg%particle_species(s))
        effective_drift_velocity(:, s) = cfg%particle_species(s)%drift_velocity
      case ('reservoir_face')
        effective_density_m3(s) = species_number_density_m3(cfg%particle_species(s))
        w_effective(s) = cfg%particle_species(s)%w_particle
        effective_temperature_k(s) = species_temperature_k(cfg%particle_species(s))
        effective_drift_velocity(:, s) = cfg%particle_species(s)%drift_velocity
      case ('photo_raycast')
        photo_emit_current_density(s) = cfg%particle_species(s)%emit_current_density_a_m2
        photo_normal_drift_speed(s) = cfg%particle_species(s)%normal_drift_speed
      end select
    end do

    call resolve_sheath_injection_context(cfg, sheath_ctx)
    if (sheath_ctx%enabled) then
      s = int(sheath_ctx%electron_species)
      effective_density_m3(s) = sheath_ctx%electron_number_density_m3
      vmin_normal(s) = max(vmin_normal(s), sheath_ctx%electron_vmin_normal)
      apply_barrier_energy_shift(s) = .false.

      if (cfg%particle_species(s)%has_target_macro_particles_per_batch .and. &
          cfg%particle_species(s)%target_macro_particles_per_batch > 0_i32) then
        call resolve_reservoir_target_weight( &
          cfg%sim, cfg%particle_species(s), effective_density_m3(s), vmin_normal(s), &
          effective_temperature_k(s), effective_drift_velocity(:, s), &
          cfg%particle_species(s)%target_macro_particles_per_batch, w_effective(s) &
          )
      end if

      if (sheath_ctx%has_local_reservoir_profile) then
        s = int(sheath_ctx%ion_species)
        effective_density_m3(s) = sheath_ctx%ion_number_density_m3
        effective_temperature_k(s) = 0.0d0
        call apply_normal_speed_override( &
          cfg%particle_species(s)%drift_velocity, sheath_ctx%reference_inward_normal, &
          sheath_ctx%ion_normal_speed_mps, effective_drift_velocity(:, s) &
          )

        if (cfg%particle_species(s)%has_target_macro_particles_per_batch .and. &
            cfg%particle_species(s)%target_macro_particles_per_batch > 0_i32) then
          call resolve_reservoir_target_weight( &
            cfg%sim, cfg%particle_species(s), effective_density_m3(s), vmin_normal(s), &
            effective_temperature_k(s), effective_drift_velocity(:, s), &
            cfg%particle_species(s)%target_macro_particles_per_batch, w_effective(s) &
            )
        end if
      end if

      if (sheath_ctx%has_photo_species) then
        s = int(sheath_ctx%photo_species)
        photo_emit_current_density(s) = sheath_ctx%photo_emit_current_density_a_m2
        photo_vmin_normal(s) = sheath_ctx%photo_vmin_normal
        photo_normal_drift_speed(s) = 0.0d0
      end if

      do s = 1, cfg%n_particle_species
        if (.not. cfg%particle_species(s)%enabled) cycle
        if (trim(lower_ascii(cfg%particle_species(s)%source_mode)) /= 'reservoir_face') cycle
        if (.not. cfg%particle_species(s)%has_target_macro_particles_per_batch) cycle
        if (cfg%particle_species(s)%target_macro_particles_per_batch == -1_i32) then
          w_effective(s) = w_effective(1)
        end if
      end do
    end if
    do s = 1, cfg%n_particle_species
      if (.not. cfg%particle_species(s)%enabled) cycle
      select case (trim(lower_ascii(cfg%particle_species(s)%source_mode)))
      case ('volume_seed')
        global_count = cfg%particle_species(s)%npcls_per_step
        counts_max(s) = mpi_split_count(global_count, local_rank, n_ranks)
      case ('reservoir_face')
        if (.not. present(state)) then
          error stop 'reservoir_face requires injection_state in init_particle_batch_from_config.'
        end if
        call reservoir_face_velocity_correction( &
          cfg%sim, cfg%particle_species(s), vmin_normal(s), barrier_normal(s), mesh &
          )
        call compute_macro_particles_for_species( &
          cfg%sim, cfg%particle_species(s), state%macro_residual(s), counts_max(s), vmin_normal=vmin_normal(s), &
          batch_duration_scale=1.0d0/real(n_ranks, dp), number_density_override=effective_density_m3(s), &
          temperature_k_override=effective_temperature_k(s), drift_velocity_override=effective_drift_velocity(:, s), &
          w_particle_override=w_effective(s) &
          )
      case ('photo_raycast')
        global_count = cfg%particle_species(s)%rays_per_batch
        counts_max(s) = mpi_split_count(global_count, local_rank, n_ranks)
      case default
        error stop 'Unknown particles.species.source_mode.'
      end select
    end do
    max_rank = max(1_i32, maxval(counts_max))
    allocate (x_species(3, max_rank, cfg%n_particle_species))
    allocate (v_species(3, max_rank, cfg%n_particle_species))
    allocate (w_species(max_rank, cfg%n_particle_species))
    allocate (emit_elem_species(max_rank, cfg%n_particle_species))
    x_species = 0.0d0
    v_species = 0.0d0
    w_species = 0.0d0
    emit_elem_species = -1_i32
    do s = 1, cfg%n_particle_species
      if (counts_max(s) <= 0_i32) cycle
      select case (trim(lower_ascii(cfg%particle_species(s)%source_mode)))
      case ('volume_seed', 'reservoir_face')
        call sample_species_state( &
          cfg%sim, cfg%particle_species(s), counts_max(s), &
          x_species(:, 1:counts_max(s), s), v_species(:, 1:counts_max(s), s), &
          barrier_normal_energy=barrier_normal(s), vmin_normal=vmin_normal(s), &
          apply_barrier_energy_shift=apply_barrier_energy_shift(s), &
          temperature_k_override=effective_temperature_k(s), drift_velocity_override=effective_drift_velocity(:, s) &
          )
        counts_actual(s) = counts_max(s)
        w_species(1:counts_actual(s), s) = w_effective(s)
      case ('photo_raycast')
        if (.not. present(mesh)) then
          error stop 'photo_raycast requires mesh in init_particle_batch_from_config.'
        end if
        if (photo_emit_current_density(s) <= 0.0d0) then
          counts_actual(s) = 0_i32
          cycle
        end if
        call sample_photo_species_state( &
          cfg%sim, cfg%particle_species(s), mesh, counts_max(s), x_species(:, 1:counts_max(s), s), &
          v_species(:, 1:counts_max(s), s), w_species(1:counts_max(s), s), counts_actual(s), &
          emit_elem_idx=emit_elem_species(1:counts_max(s), s), &
          global_rays_per_batch=cfg%particle_species(s)%rays_per_batch, &
          emit_current_density_override=photo_emit_current_density(s), &
          normal_drift_speed_override=photo_normal_drift_speed(s), vmin_normal=photo_vmin_normal(s) &
          )
        if (present(photo_emission_dq) .and. cfg%particle_species(s)%deposit_opposite_charge_on_emit) then
          do i = 1, counts_actual(s)
            if (emit_elem_species(i, s) < 1_i32 .or. emit_elem_species(i, s) > size(photo_emission_dq)) then
              error stop 'photo_raycast emitted invalid elem_idx.'
            end if
            photo_emission_dq(emit_elem_species(i, s)) = photo_emission_dq(emit_elem_species(i, s)) - &
                                                         cfg%particle_species(s)%q_particle*w_species(i, s)
          end do
        end if
      case default
        error stop 'Unknown particles.species.source_mode.'
      end select
    end do

    batch_n = sum(counts_actual)
    allocate (species_id(batch_n))
    out_idx = 0_i32
    do i = 1, max_rank
      do s = 1, cfg%n_particle_species
        if (i > counts_actual(s)) cycle
        out_idx = out_idx + 1_i32
        species_id(out_idx) = s
      end do
    end do

    allocate (x(3, batch_n), v(3, batch_n), q(batch_n), m(batch_n), w(batch_n))
    allocate (species_cursor(cfg%n_particle_species))
    species_cursor = 0_i32
    do i = 1, batch_n
      s = species_id(i)
      species_cursor(s) = species_cursor(s) + 1_i32
      x(:, i) = x_species(:, species_cursor(s), s)
      v(:, i) = v_species(:, species_cursor(s), s)
      q(i) = cfg%particle_species(s)%q_particle
      m(i) = cfg%particle_species(s)%m_particle
      w(i) = w_species(species_cursor(s), s)
    end do

    call init_particles(pcls, x, v, q, m, w)
  end subroutine init_particle_batch_from_config