From a0b73194d268bff1f7212ec98088e9d480b77719 Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Tue, 12 Nov 2024 11:37:40 +0100
Subject: [PATCH] Very small progress, E at r0 still high

---
 plasmaExpansion.f90 | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90
index a4f492b..3cc155b 100644
--- a/plasmaExpansion.f90
+++ b/plasmaExpansion.f90
@@ -165,8 +165,8 @@ program plasmaExpansion
   u_i     = 0.0_dp
   E_i     = 0.0_dp
   T_i     = 0.0_dp
-  n_e     = 1.0_dp
-  T_e     = 1.0_dp
+  n_e     = 0.0_dp
+  T_e     = 0.0_dp
   Zave    = 0.0_dp
   phi     = 0.0_dp
   phi_old = 0.0_dp
@@ -209,6 +209,7 @@ program plasmaExpansion
 
   ! Set boundary values
   phi0 = 0.0_dp / phi_ref ! Dirichlet
+  phi(1) = phi0 ! Dirichlet
   ! phi0 = phi(1) ! Neumann
   ! phiF = 0.0_dp ! Dirichlet
   allocate(f0(j0:nv))
@@ -232,6 +233,7 @@ program plasmaExpansion
     ! 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)
     f0 = f0 * n_bc / (sum(f0)*dv)
+    T_e = Temp_bc
 
     ! Boundary conditions
     ! r = r0, v>0
@@ -280,6 +282,11 @@ program plasmaExpansion
       ! Store previous value
       phi_old = phi
 
+      ! Calculate distribution of electrons
+      ! n_e = Zave * n_i ! Quasi-neutral
+      n_e = Zave(1) * n_i(1) * exp((phi_old - phi0) / T_e) ! Isothermal (Boltzmann)
+      ! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi_old-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working
+
       ! Diagonal matrix for Newton integration scheme
       db_dphi        = n_e / T_e ! Isotropic
       diag           = -2.0_dp / dr**2 - db_dphi
@@ -306,12 +313,7 @@ program plasmaExpansion
       ! Iterate system
       call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info)
       phi = phi_old + 1.0e-1_dp*Res
-      phi0=phi(1) ! Neumann
-
-      ! Calculate distribution of electrons
-      T_e = Temp_bc
-      n_e = Zave(1) * n_i(1) * exp((phi - phi0) / T_e)
-      ! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working
+      ! phi0=phi(1) ! Neumann
 
       ! Check if the solution has converged
       phiConv = maxval(abs(Res),1)