From 4511aede8c5fa9c653f617a1b885b244a49dfacd Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Wed, 13 Nov 2024 12:52:54 +0100
Subject: [PATCH] Nothing important

---
 plasmaExpansion.f90 | 20 +++++++++-----------
 1 file changed, 9 insertions(+), 11 deletions(-)

diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90
index 3cc155b..aea4b9e 100644
--- a/plasmaExpansion.f90
+++ b/plasmaExpansion.f90
@@ -99,7 +99,7 @@ program plasmaExpansion
 
   ! Set domain boundaries (non-dimensional units)
   r0 =  10.0e-6_dp / L_ref
-  rf =  50.0e-6_dp / L_ref
+  rf = 100.0e-6_dp / L_ref
   dr = 2.0e1_dp
   nr = nint((rf - r0) / dr) + 1
   dr =      (rf - r0) / float(nr-1)
@@ -109,7 +109,7 @@ program plasmaExpansion
   end do
 
   ! Set position to calculate cumulative sum of f (non-dimensional units)
-  rCum = 40.0e-6 / L_ref
+  rCum = 80.0e-6 / L_ref
 
   ! Index for cumulative sum
   rCum_index = minloc(abs(r - rCum), 1)
@@ -208,10 +208,9 @@ program plasmaExpansion
   Res = 0.0_dp
 
   ! Set boundary values
-  phi0 = 0.0_dp / phi_ref ! Dirichlet
+  phi0 = 0.0e0_dp / phi_ref ! Dirichlet
   phi(1) = phi0 ! Dirichlet
   ! phi0 = phi(1) ! Neumann
-  ! phiF = 0.0_dp ! Dirichlet
   allocate(f0(j0:nv))
   f0 = 0.0_dp
 
@@ -282,11 +281,6 @@ 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
@@ -294,8 +288,7 @@ program plasmaExpansion
       diag_high      =  1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * dr)
       diag(1)        =  1.0_dp ! Dirichlet
       diag_high(1)   =  0.0_dp ! Dirichlet
-      ! diag_high(1)   =  2.0_dp / dr**2 - db_dphi(1) ! Neumann
-      ! diag(nr)       =  2.0_dp ! Dirichlet
+      ! diag(nr)       =  1.0_dp ! Dirichlet
       ! diag_low(nr-1) =  0.0_dp ! Dirichlet
       diag_low(nr-1) =  2.0_dp / dr**2 - db_dphi(nr) ! Neumann
 
@@ -315,6 +308,11 @@ program plasmaExpansion
       phi = phi_old + 1.0e-1_dp*Res
       ! phi0=phi(1) ! Neumann
 
+      ! Calculate distribution of electrons
+      ! n_e = Zave * n_i ! Quasi-neutral
+      n_e = Zave(1) * n_i(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
+      ! 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
+
       ! Check if the solution has converged
       phiConv = maxval(abs(Res),1)
       if (phiConv < 1.0e-3_dp) then