From a03d39826d42c1f01f037fb2a61291a61c2632d2 Mon Sep 17 00:00:00 2001
From: JHendrikx <j.hendrikx@student.tue.nl>
Date: Thu, 3 Apr 2025 10:31:37 +0200
Subject: [PATCH] Add TtoZ table data to simulation

---
 vlaplex.f90 | 25 ++++++++++++++++++++++---
 1 file changed, 22 insertions(+), 3 deletions(-)

diff --git a/vlaplex.f90 b/vlaplex.f90
index 9471a63..44550cf 100644
--- a/vlaplex.f90
+++ b/vlaplex.f90
@@ -26,6 +26,8 @@ program VlaPlEx
   use referenceValues
   use eos,                only: T_to_Z
   use moduleTableBC
+  use moduleTableTtoZ
+  use moduleTableTtoZne
   use omp_lib
   implicit none
 
@@ -48,12 +50,16 @@ program VlaPlEx
   integer:: j0 ! First integer of positive velocity
 
   real(dp):: Temp_bc ! Temperature
-  real(dp):: Zave_bc ! Average charge state
+  real(dp):: Zave_bc, Zave_bc_old ! Average charge state
   real(dp):: u_bc    ! Injection velocity
   real(dp):: n_bc    ! Injection density
   real(dp):: c_s     ! Ion sound speed
   type(tableBC):: boundaryConditions
   character(:), allocatable:: bc_file
+  type(tableTtoZ):: TtoZ
+  character(:), allocatable:: TtoZ_file
+  type(tableTtoZne):: TtoZne
+  character(:), allocatable:: TtoZne_file
 
   real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old
   real(dp), allocatable, dimension(:):: f0   ! Boundary at r = x_0
@@ -97,6 +103,10 @@ program VlaPlEx
   c_s      = sqrt(11.0_dp * gamma_i * 1.0_dp)
   bc_file = 'bc.csv'
   call boundaryConditions%init(bc_file)
+  TtoZ_file = 'average_TtoZ_curve.csv'
+  call TtoZ%init(TtoZ_file)
+  TtoZne_file = 'Zave_values_per_Ne.csv'
+  call TtoZne%init(TtoZne_file)
 
   ! Set domain boundaries (non-dimensional units)
   r0 =  10.0e-6_dp / L_ref
@@ -176,6 +186,7 @@ program VlaPlEx
   T_e     = 0.0_dp
   Zave    = 0.0_dp
   Z_list  = (/ 6.0, 12.0 /)
+  Zave_bc_old = 0.0_dp
   phi     = 0.0_dp
   phi_old = 0.0_dp
   E       = 0.0_dp
@@ -231,15 +242,23 @@ program VlaPlEx
   call writeOutputFCum(t, dt, nz, r(rCum_index), nv, v, fCum_i, Z_list)
   call writeOutputPhi(t, dt, nr, r, phi, E, n_e)
   call writeOutputMom(t, dt, nz, nr, r, n_i, u_i, T_i, Z_list)
+  call writeOutputZList(nz, Z_list)
 
   ! Main loop
   do t = 1, nt
     time = t * dt + t0
     call boundaryConditions%get(time, n_bc, u_bc, Temp_bc)
-    z_inj = minloc(abs(Z_list - T_to_Z(Temp_bc)),1)
+    call TtoZ%get(Temp_bc, Zave_bc)
+    z_inj = minloc(abs(Z_list - Zave_bc),1)
     Zave_bc = Z_list(z_inj)
+    do while (Zave_bc - Zave_bc_old > 0.1_dp)
+      Zave_bc_old = Zave_bc
+      call TtoZne%get(Temp_bc, Zave_bc * n_bc, Zave_bc)
+      z_inj = minloc(abs(Z_list - Zave_bc),1)
+      Zave_bc = Z_list(z_inj)
+    end do
+    call writeOutputBoundary(t, dt, Zave_bc*n_bc, u_bc, Temp_bc, Zave_bc, Zave_bc)
     u_bc = sqrt(Zave_bc * Temp_bc)
-    call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, T_to_Z(Temp_bc), Zave_bc)
 
     ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
     f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)