From d7c23d5577ad0fce90135b00886610e38efc8d1f Mon Sep 17 00:00:00 2001
From: JHendrikx <j.hendrikx@student.tue.nl>
Date: Tue, 4 Feb 2025 12:10:58 +0100
Subject: [PATCH] Change b_i to sum_ni

---
 vlaplex.f90 | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/vlaplex.f90 b/vlaplex.f90
index 8510f57..558957f 100644
--- a/vlaplex.f90
+++ b/vlaplex.f90
@@ -61,7 +61,7 @@ program VlaPlEx
   real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old
   real(dp), allocatable, dimension(:):: f0   ! Boundary at r = x_0
   real(dp), allocatable, dimension(:,:):: n_i
-  real(dp), allocatable, dimension(:):: b_i
+  real(dp), allocatable, dimension(:):: sum_ni
   real(dp), allocatable, dimension(:,:):: u_i
   real(dp), allocatable, dimension(:):: E_i
   real(dp), allocatable, dimension(:,:):: T_i
@@ -161,7 +161,7 @@ program VlaPlEx
   ! Allocate vectors
   allocate(f_i(1:nz,1:nr,1:nv), f_i_old(1:nz,1:nr,1:nv))
   allocate(n_i(1:nz,1:nr))
-  allocate(b_i(1:nr))
+  allocate(sum_ni(1:nr))
   allocate(u_i(1:nz,1:nr), E_i(1:nr), T_i(1:nz,1:nr))
   allocate(Zave(1:nr))
   allocate(n_e(1:nr))
@@ -170,7 +170,7 @@ program VlaPlEx
   f_i     = 0.0_dp
   f_i_old = 0.0_dp
   n_i     = 0.0_dp
-  b_i   = 0.0_dp
+  sum_ni   = 0.0_dp
   u_i     = 0.0_dp
   E_i     = 0.0_dp
   T_i     = 0.0_dp
@@ -266,7 +266,7 @@ program VlaPlEx
     ! set edge velocities to 0
     f_i_old(:,:,1)  = 0.0_dp
     f_i_old(:,:,nv) = 0.0_dp
-    b_i = 0.0_dp
+    sum_ni = 0.0_dp
     ! Advect in the r direction
     !$omp parallel do
     do iz = 1, nz
@@ -296,12 +296,12 @@ program VlaPlEx
         end if
 
       end do
-      b_i = b_i + Zave * n_i(iz,:)
+      sum_ni = sum_ni + Zave * n_i(iz,:)
     end do
     !$omp end parallel do
 
     ! Assume quasi-neutrality to start iterating
-    n_e = 1.0_dp/nz * b_i
+    n_e = 1.0_dp/nz * sum_ni
     db_dphi = 0.0_dp
 
     ! Solve Poission (maximum number of iterations, break if convergence is reached before)
@@ -319,7 +319,7 @@ program VlaPlEx
       diag_low(nr-1) =  2.0_dp / dr**2 - db_dphi(nr) ! Neumann
 
       ! Calculate charge density
-      b = - 1.0_dp/nz * b_i + n_e
+      b = - 1.0_dp/nz * sum_ni + n_e
       ! Apply boundary conditions
       b(1)  = phi0 ! Dirichlet
       ! b(nr) = 0.0_dp ! Dirichlet
@@ -335,11 +335,11 @@ program VlaPlEx
       ! phi0=phi(1) ! Neumann
 
       ! Calculate distribution of electrons
-      ! n_e = Zave(1) * n_i(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
-      n_e = Zave(1) * n_i(1,1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic
+      ! n_e = sum_ni(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
+      n_e = sum_ni(1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic
       ! Diagonal matrix for Newton integration scheme
       ! db_dphi        = n_e / T_e ! Isotropic
-      db_dphi        =  Zave(1) * n_i(1,1) / (gamma_e * T_e) * &
+      db_dphi        =  sum_ni(1) / (gamma_e * T_e) * &
                         (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_dexp !Polytropic
 
       ! Check if the solution has converged