From 03b0619991c37a4789c6cbc3d22e34244ceb3369 Mon Sep 17 00:00:00 2001
From: JHendrikx <j.hendrikx@student.tue.nl>
Date: Mon, 26 May 2025 10:31:54 +0200
Subject: [PATCH] Add recombination and ionization to model

---
 vlaplex.f90 | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 64 insertions(+), 2 deletions(-)

diff --git a/vlaplex.f90 b/vlaplex.f90
index 5031648..04d3163 100644
--- a/vlaplex.f90
+++ b/vlaplex.f90
@@ -20,14 +20,36 @@
 
 end module eos
 
+ module Temperature
+   use constantParameters, only: dp
+   implicit none
+
+   private
+   public:: T_function
+
+   contains
+     pure function T_function(T0,n0, ne, gamma_e) result(T)
+       use constantParameters, only: eV_to_K
+       implicit none
+
+       real(dp), intent(in):: T0, n0, ne, gamma_e
+       real(dp):: T
+
+      T = T0 * (ne/n0)**(gamma_e -1.0_dp)
+
+     end function T_function
+
+end module Temperature
 program VlaPlEx
   use constantParameters, only: dp, kb, qe, eps_0, ev_to_K, cm3_to_m3, PI
   use output
   use referenceValues
   use eos,                only: T_to_Z
+  use Temperature,         only: T_function
   use moduleTableBC
   use moduleTableTtoZ
   use moduleTableTtoZne
+  use moduleTableRateCoeff
   use omp_lib
   implicit none
 
@@ -60,6 +82,10 @@ program VlaPlEx
   character(:), allocatable:: TtoZ_file
   type(tableTtoZne):: TtoZne
   character(:), allocatable:: TtoZne_file
+  type(tableRateCoeff):: recombination_rateCoeff
+  character(:), allocatable:: recombination_rate_file
+  type(tableRateCoeff):: ionization_rateCoeff
+  character(:), allocatable:: ionization_rate_file
 
   real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old
   real(dp), allocatable, dimension(:):: f0   ! Boundary at r = x_0
@@ -81,6 +107,8 @@ program VlaPlEx
   real(dp):: phiConv
   real(dp):: phi0
   real(dp):: T_e
+  real(dp):: recombination_rate
+  real(dp):: ionization_rate
   ! real(dp):: phiF
   integer:: k
 
@@ -107,6 +135,10 @@ program VlaPlEx
   call TtoZ%init(TtoZ_file)
   TtoZne_file = 'Zave_values_per_Ne.csv'
   call TtoZne%init(TtoZne_file)
+  recombination_rate_file = 'dr.csv'
+  call recombination_rateCoeff%init(recombination_rate_file)
+  ionization_rate_file = 'ci.csv'
+  call ionization_rateCoeff%init(ionization_rate_file)
 
   ! Set domain boundaries (non-dimensional units)
   r0 =  10.0e-6_dp / L_ref
@@ -184,6 +216,8 @@ program VlaPlEx
   T_i     = 0.0_dp
   n_e     = 0.0_dp
   T_e     = 0.0_dp
+  recombination_rate = 0.0_dp
+  ionization_rate = 0.0_dp
   Zave    = 0.0_dp
   Z_list  = (/ 6.0, 12.0 /)
   Zave_bc_old = 0.0_dp
@@ -277,9 +311,37 @@ program VlaPlEx
     f_i_old(:,:,1)  = 0.0_dp
     f_i_old(:,:,nv) = 0.0_dp
     sum_ni = 0.0_dp
+
+    ! --- Recombination Step --- !
+    do iz = 2, nz
+      do i = 2, nr
+        call recombination_rateCoeff%get(T_function(Temp_bc, n_e(1), n_e(i), gamma_e), Z_list(iz), recombination_rate)
+        ! Recombination of current charge state
+        f_i(iz,i,:) = f_i_old(iz,i,:) - f_i_old(iz,i,:) * dt * recombination_rate * n_e(i)
+
+        ! Transfer recombination density to lower charge state
+        f_i(iz - 1,i,:) = f_i_old(iz - 1,i,:) + f_i_old(iz,i,:) * dt * recombination_rate * n_e(i)
+      end do
+      f_i_old = f_i
+    end do
+
+
+    ! --- Ionization Step --- !
+    do iz = 1, nz - 1
+      do i = 2, nr
+        call ionization_rateCoeff%get(T_function(Temp_bc, n_e(1), n_e(i), gamma_e), Z_list(iz), ionization_rate)
+        ! Ionization of current charge state
+        f_i(iz,i,:) = f_i_old(iz,i,:) - f_i_old(iz,i,:) * dt * ionization_rate * n_e(i)
+
+        ! Transfer ionization density to higher charge state
+        f_i(iz + 1,i,:) = f_i_old(iz + 1,i,:) + f_i_old(iz,i,:) * dt * ionization_rate * n_e(i)
+      end do
+      f_i_old = f_i
+    end do
+
     ! Advect in the r direction
     do iz = 1, nz
-      if (all(n_i(iz,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
+      if (all(f_i(iz,:,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
         cycle
       end if
       !$omp parallel do
@@ -386,7 +448,7 @@ program VlaPlEx
     f_i_old = f_i
 
     do iz = 1, nz
-      if (all(n_i(iz,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
+      if (all(f_i(iz,:,:) < 1.0e-16_dp) .and. iz .NE. z_inj) then
         cycle
       end if
       ! Advect in the v direction