Major changes in the ccpla script

pietromandracci · Aug 15, 2022 · 101b18c · 101b18c
1 parent 7c1deb1
commit 101b18c
Show file tree

Hide file tree

Showing 35 changed files with 2,596 additions and 1,504 deletions.
diff --git a/changes.rst b/changes.rst
@@ -1,3 +1,21 @@
+=============
+Version 0.2.1
+=============
+
+* Major changes in the script ccpla
+  - ion-neutral scattering is now calculated in the center of mass frame of reference
+  - more data are saved to file during the simulation, including potential distribution and electric current
+
+* Changes in the script ccpla_analysis
+  - it is now possible to plot other data, such as the  potential distribution
+
+=============
+Version 0.1.0
+=============
+
+* Major changes in the script ccpla
+  - electron impact excitation processes have been added to the simulation
+  - some changes in the GUI have made to correct some erros that did happen when reloading configuration files from the GUI  
 
 =============
 Version 0.0.1

diff --git a/doc/plasmapro/ccpla_manual.html b/doc/plasmapro/ccpla_manual.html
@@ -3,7 +3,7 @@
 <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
 <head>
 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
-<meta name="generator" content="Docutils 0.16: http://docutils.sourceforge.net/" />
+<meta name="generator" content="Docutils 0.17.1: http://docutils.sourceforge.net/" />
 <title>CCPLA: Capacitively Coupled PLAsma simulation</title>
 <style type="text/css">
 

diff --git a/pysica/analysis/univariate.py b/pysica/analysis/univariate.py
@@ -983,6 +983,7 @@ def plot_lag(self, lag=1, interface=PLOT_INTERFACE, symbol=PLOT_SYMBOL,
             return (status, message)
 
         return (status, message)
+
 
     def plot_autocorrelation(self, start=0, step=1, n_points=None,
                              alternate_method=False,

diff --git a/pysica/fortran/fmathematics.cpython-310-x86_64-linux-gnu.so b/pysica/fortran/fmathematics.cpython-310-x86_64-linux-gnu.so
diff --git a/pysica/fortran/fmathematics.cpython-39-x86_64-linux-gnu.so b/pysica/fortran/fmathematics.cpython-39-x86_64-linux-gnu.so
diff --git a/pysica/functions/pdf.py b/pysica/functions/pdf.py
@@ -173,17 +173,17 @@ def pdf_margenau_energy(e, B, M, E, w, l):
 
 
 def pdf_maxwell_speed(v, m, kt):
-    """ Maxwell probability distribution function for the particle velocity
+    """ Maxwell probability distribution function for the particle speed
 
-        f(e) = sqrt(2/pi * (m/kt)**3) * v**2 * exp(-m*v**2 / (2*kt))
+        f(v) = sqrt(2/pi * (m/kt)**3) * v**2 * exp(-m*v**2 / (2*kt))
 
         Parameters
         ----------
  
         v:   speed
         m:   mass of the particle
         kt:  temperature in energy units 
-             (must be compatible to the ones used for v)
+             (must be compatible to the ones used for v and m)
 
         Returns
         -------
@@ -192,3 +192,45 @@ def pdf_maxwell_speed(v, m, kt):
     """
 
     return sqrt( 2/pi * (m/kt)**3 ) * v*v * exp(- m*v*v / (2*kt) )
+
+
+def pdf_maxwell_reduced(x):
+    """ Maxwell probability distribution function for the adimensional particle speed
+
+        f(x) = 4/sqrt(pi) * x**2 * exp(-x**2)
+
+        where x = 2/sqrt(pi) v / v_mean
+
+        Parameters
+        ----------
+ 
+        x:     adimensional speed x = 2/sqrt(pi) v / v_mean
+        Returns
+        -------
+
+        pdf_maxwell_reduced: value of the pdf
+    """
+
+    return 4 / sqrt(pi) * x*x * exp(- x*x )
+
+
+def pdf_maxwell_velocity_component(vx, m, kt):
+    """ Maxwell probability distribution function for a component of the particle velocity vector
+
+        f(vx) = (m / (2*pi*kt))**1/2 * exp(-m*v**2 / (2*kt))
+
+        Parameters
+        ----------
+ 
+        vx:  component of the velocity vector
+        m:   mass of the particle
+        kt:  temperature in energy units 
+             (must be compatible to the ones used for v and m)
+
+        Returns
+        -------
+
+        pdf_maxwell_speed: value of the pdf
+    """
+
+    return sqrt(m / (2*pi*kt)) * exp(- m*vx*vx / (2*kt) )
diff --git a/pysica/functions/physics.py b/pysica/functions/physics.py
@@ -35,6 +35,24 @@
 # | Physical functions |
 # +--------------------+
 
+def mean_speed_maxwell(temperature, mass):
+    """ Calculates the mean speed of the molecules in a gas at thermal equilibrium 
+
+        Parameters
+        ----------
+
+        temperature:    gas temperature / K
+        mass:           molecular mass / atomic units
+
+        Returns
+        -------
+
+        mean speed / m s**-1
+
+    """
+
+    return numpy.sqrt(8 * K_BOLTZMANN * temperature / (numpy.pi * mass * ATOMIC_UNIT_MASS))
+
 
 def number_density(pressure, temperature):
     """ Calculates the number density of a gas at given pressure and temperature.
@@ -51,7 +69,7 @@ def number_density(pressure, temperature):
         number density: numer density of gas molecules / m**-3
     """
 
-    return pressure / (K_BOLTZMANN*temperature)
+    return pressure / (K_BOLTZMANN * temperature)
 
 
 def pressure_conversion(pressure_in, unit_in, unit_out):

diff --git a/pysica/functions/random_pdf.py b/pysica/functions/random_pdf.py
@@ -0,0 +1,82 @@
+# COPYRIGHT 2020 by Pietro Mandracci
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+""" PYthon tools for SImulation and CAlculus: random generation utilities.
+
+    This module contains some functions to generate random number with specific distributions.
+
+    Documentation is also available in the docstrings.
+"""
+
+#import math
+import numpy
+from random import random
+
+from ..constants import *
+from .pdf import *
+
+def random_maxwell_velocity(v_mean, modulus=False):
+    """ Generate the x,y,z components of a velocity vector following Maxwell distribution
+
+        Parameters
+        ----------
+
+        v_mean:   mean value of the velocity module
+        modulus:  if True, the modulus of the velocity vector is returned also
+
+        Returns
+        -------
+
+        if modulus==False: (vx, vy, vz):   components of the random velocity vector
+        if modulus==False: (vx, vy, vz,v): components and modulus of the random velocity vector
+    """
+
+    v = random_maxwell_speed(v_mean)
+
+    costheta = 2 * random() - 1
+    sintheta = numpy.sqrt(1 - costheta * costheta)
+    phi      = 2 * numpy.pi * random()
+    vx       = v * sintheta * numpy.cos(phi)
+    vy       = v * sintheta * numpy.sin(phi)
+    vz       = v * costheta
+
+    if modulus: return (vx, vy, vz, v)
+    else:       return (vx, vy, vz)      
+
+
+def random_maxwell_speed(v_mean):
+    """ Generate a random speed according to the Maxwell distribution.
+
+        Parameters
+        ----------
+        
+        v_mean: mean value of the velocity module       
+
+        Returns
+        -------
+        
+        random speed value
+    """
+
+    xmax = 4.0
+    ymax = 4.0 / (numpy.sqrt(numpy.pi) * numpy.e)
+    f = 0
+    y = 1
+    while(f < y):
+        x = random() * xmax
+        y = random() * ymax
+        f = pdf_maxwell_reduced(x)
+
+    return sqrt(pi) / 2 * v_mean * x
diff --git a/pysica/parameters.py b/pysica/parameters.py
@@ -24,6 +24,8 @@
 PLOT_INTERFACE  = 'pylab'
 PLOT_SYMBOL     = '+'
 PLOT_COLOR      = 'red'
+PLOT_SYMBOLS    = [ '.',     'o',    '+',    'x',     '^',    'v',       '>',   '<',      'D'      ]
+N_PLOT_SYMBOLS  = len(PLOT_SYMBOLS)
 PLOT_COLORS     = [ 'black', 'grey', 'blue', 'green', 'cyan', 'magenta', 'red', 'orange', 'yellow' ]
 N_PLOT_COLORS   = len(PLOT_COLORS)
 FILL_COLOR      = 'grey'

diff --git a/pysica/plasmapro/ccpla.py b/pysica/plasmapro/ccpla.py
@@ -94,7 +94,7 @@ def exit_ccpla(message="Exiting program", gui=False, error=False):
 
 def print_message(message, gui=False, error=False):
     if gui:
-        create_message_window(message=message, error=error)
+        create_message_window(message, error)
     else:
         if error: message = 'ERROR: ' + message
         print('\n' + message + '\n')
@@ -168,7 +168,7 @@ def print_message(message, gui=False, error=False):
                   default=debug_level_python, help=help_string)
 
 # Debug level for Fotran code
-help_string = "Fortran debug level [0..2] (default="+str(debug_level_fortran)+")"
+help_string = "Fortran debug level [0..3] (default="+str(debug_level_fortran)+")"
 parser.add_option("-D", "--debug-level-fortran", action="store", type="int", dest="debug_lev_for",
                   default=debug_level_fortran, help=help_string)        
 
@@ -197,8 +197,8 @@ def print_message(message, gui=False, error=False):
     exit_ccpla('ERROR: Python debug level must be in range 0..2' + EOL,
                gui=cl_options.gui_mode, error=True)
 
-if (cl_options.debug_lev_for not in list(range(3))):
-    exit_ccpla('ERROR: Fortran debug level must be in range 0..2' + EOL,
+if (cl_options.debug_lev_for not in list(range(4))):
+    exit_ccpla('ERROR: Fortran debug level must be in range 0..3' + EOL,
                gui=cl_options.gui_mode, error=True)
 
 if cl_options.gui_mode:
@@ -268,7 +268,6 @@ def print_message(message, gui=False, error=False):
 if (not cl_options.gui_mode): print(GPL_MESSAGE)
 
 if (not cl_options.batch_mode): wait_input()
-
 
 
 # +-------------------------+
@@ -304,7 +303,7 @@ def print_message(message, gui=False, error=False):
 # | Load gas properties |
 # +---------------------+
 
-if (cl_options.verbosity > 0): print('\nDefining neutrals ensamble from file \"' + FILENAME_NEUTRALS + '\" ...')
+if (cl_options.verbosity > 0): print('\nCreating neutrals ensamble from file \"' + FILENAME_NEUTRALS + '\" ...')
 neutrals = target_particles.TargetParticles( parameters.N_sigma,
                                              parameters.N_sigma_ions,
                                              parameters.T_neutrals,
@@ -345,7 +344,7 @@ def print_message(message, gui=False, error=False):
 
 # Define the ensamble of charged particles moving in the plasma: electrons and ions
 if (cl_options.verbosity > 0): print('\nDefining charged particles ensambles ...')
-charges = moving_particles.MovingParticles(neutrals.types+1, parameters.Nmax_particles, START_WEIGHT, parameters.rescale_factor)
+charges = moving_particles.MovingParticles(neutrals.types+1, parameters.Nmax_particles, parameters.start_weight, parameters.rescale_factor)
 
 (status, message) = initialize_ensambles(charges, neutrals, parameters, cl_options)
 if (status!=0): exit_ccpla(message + EOL, gui=cl_options.gui_mode, error=True)
@@ -391,9 +390,12 @@ def print_message(message, gui=False, error=False):
 # +--------------------------+
 
 if (parameters.save_delay > 0):
-    charges.initialize_savefiles(parameters.filename_stat_ele, parameters.filename_distrib_ele, parameters.filename_distrib_ion,
+    charges.initialize_savefiles(parameters.filename_stat_ele,
+                                 parameters.filename_distrib_ele,
+                                 parameters.filename_distrib_ion,
                                  append=False, sep='\t', ext=EXT)
     neutrals.initialize_savefile(parameters.filename_stat_neu, append=False, sep='\t', ext=EXT)
+    ccp.initialize_savefiles(parameters.filename_I, parameters.filename_V, append=False, sep='\t', ext=EXT)    
 
 
 # +----------------+
@@ -444,6 +446,7 @@ def print_message(message, gui=False, error=False):
         if ( (i_save_data >= parameters.save_delay) and (charges.n_active(0) > 0) ):
             charges.save_data_to_files()
             neutrals.save_data_to_files(charges.time)
+            ccp.save_data_to_files(charges.time)
             i_save_data = 0
         i_save_data += 1