From 49edcd0df9ee73586a99a386deabf8d13c6d717c Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Fri, 18 Oct 2024 09:51:15 +0200
Subject: [PATCH] Cleaning

Trying to do polytropic electrons.
---
 plasmaExpansion.f90 | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90
index ee875a5..ab18b34 100644
--- a/plasmaExpansion.f90
+++ b/plasmaExpansion.f90
@@ -70,7 +70,7 @@ program plasmaExpansion
   real(dp), allocatable, dimension(:):: b
   integer:: info
 
-  real(dp), allocatable, dimension(:):: phi, phi_old, E
+  real(dp), allocatable, dimension(:):: phi, phi_old, E, db_dphi
   real(dp):: phiConv
   real(dp):: phi0
   real(dp):: T_e
@@ -131,7 +131,7 @@ program plasmaExpansion
   end if
 
   t0 = 0.0_dp
-  tf = 1.0e-6_dp / t_ref
+  tf = 2.0e-6_dp / t_ref
   ! tf = 1.0e1_dp * (rf - r0) / c_s
   dt = 1.0e-2_dp*dr/c_s
   nt = nint((tf - t0) / dt) 
@@ -176,10 +176,12 @@ program plasmaExpansion
   ! Allocate matrix for Poisson equation
   allocate(diag(1:nr), diag_low(1:nr-1), diag_high(1:nr-1))
   allocate(b(1:nr))
+  allocate(db_dphi(1:nr))
   diag      = 0.0_dp
   diag_low  = 0.0_dp
   diag_high = 0.0_dp
   b         = 0.0_dp
+  db_dphi   = 0.0_dp
   diag           = -2.0_dp / dr**2
   diag_low       =  1.0_dp / dr**2 - 1.0_dp / (r(2:nr)   * 2.0_dp * dr)
   diag_high      =  1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * 2.0_dp * dr)
@@ -278,7 +280,8 @@ program plasmaExpansion
       phi_old = phi
 
       ! Diagonal matrix for Newton integration scheme
-      diag           = -2.0_dp / dr**2 - n_e / T_e
+      db_dphi        = n_e / T_e ! Isotropic
+      diag           = -2.0_dp / dr**2 - db_dphi
       diag_low       =  1.0_dp / dr**2 - 1.0_dp / (r(2:nr)   * 2.0_dp * dr)
       diag_high      =  1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * 2.0_dp * dr)
       diag(1)        =  1.0_dp ! Dirichlet
@@ -286,7 +289,7 @@ program plasmaExpansion
       ! diag_high(1)   =  2.0_dp / dr**2 - n_e(1) ! Neumann
       ! diag(nr)       =  1.0_dp ! Dirichlet
       ! diag_low(nr-1) =  0.0_dp ! Dirichlet
-      diag_low(nr-1) =  2.0_dp / dr**2 - n_e(nr) / T_e ! Neumann
+      diag_low(nr-1) =  2.0_dp / dr**2 - db_dphi(nr) ! Neumann
 
       ! Calculate charge density
       b = -(Zave*n_i - n_e)
@@ -307,7 +310,7 @@ program plasmaExpansion
       ! Calculate distribution of electrons
       T_e = 1.0_dp
       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_i(1))**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working
+      ! 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)