Merge pull request 'issue/edgesVolume' (#59 ) from issue/edgesVolume into development

Reviewed-on: https://plasmalab.aero.upm.es/git/git/JGonzalez/fpakc/pulls/59
Merge remote-tracking branch 'origin/issue/Maxwellian' into issue/edgesVolume
2026-02-26 19:21:52 +01:00 · 2026-02-04 19:59:59 +01:00 · 2026-02-04 13:52:24 +01:00 · 2026-02-03 10:15:12 +01:00 · 2026-02-02 14:34:47 +01:00 · 2026-02-02 14:31:58 +01:00
41 changed files with 1833 additions and 1528 deletions
--- a/CITATION.cff
+++ b/CITATION.cff
@ -0,0 +1,50 @@
 # This CITATION.cff file was generated with cffinit.
 # Visit https://bit.ly/cffinit to generate yours today!
 cff-version: 1.2.0
 title: Finite element Particle Kinetic Code
 message: >-
  If you use this software, please cite it using the
  metadata from this file.
 type: software
 authors:
  - given-names: Jorge
    family-names: Gonzalez
    email: jorge.gonzalez@upm.es
    affiliation: Universidad Politécnica de Madrid
    orcid: 'https://orcid.org/0000-0001-7905-5001'
 repository-code: 'https://gitlab.com/JorgeGonz/fpakc'
 abstract: >-
  Welcome to fpakc (Finite element PArticle Kinetic Code), a
  modern object oriented Fortran open-source code for
  particle simulations of plasma and gases. This code works
  by simulating charged and neutral particles, following
  their trajectories, collisions and boundary conditions
  imposed by the user.
  One of our aims is to make a code easy to maintain as well
  as easy to use by a variety of reserchers and students.
  This code is currenlty in very early steps of development.
  The code aims to be easy to maintain and easy to use,
  allowing its application from complex problems to easy
  examples that can be used, for example, as teaching
  exercises.
  Parallelization techniques such as OpenMP, MPI will be
  used to distribute the cpu load. We aim to make fpakc GPU
  compatible in the future.
  The codefpakc makes use of finite elements to generate
  meshes in complex geometries. Particle properties are
  deposited in the nodes and cells of the mesh. The
  electromagnetic field, with the boundary conditions
  imposed by the user, is solved also in this mesh.
 keywords:
  - particle-in-cell
  - plasma
  - finite elements
 license: GPL-3.0
 version: beta
 date-released: '2025-10-01'
--- a/data/collisions/IO_e-Kr.dat
+++ b/data/collisions/IO_e-Kr.dat
@ -0,0 +1,52 @@
 # D108525 "refs": {"B56": {"note": "CLM-R294 (1989)"}}
 # Relative energy (eV) cross section (m^2)
 1.40E+01 0
 1.62E+01 7.249E-21
 1.88E+01 1.199E-20
 2.18E+01 1.644E-20
 2.53E+01 2.1E-20
 2.94E+01 2.542E-20
 3.41E+01 2.937E-20
 3.95E+01 3.26E-20
 4.58E+01 3.499E-20
 5.32E+01 3.653E-20
 6.17E+01 3.726E-20
 7.15E+01 3.728E-20
 8.29E+01 3.671E-20
 9.62E+01 3.566E-20
 1.12E+02 3.426E-20
 1.29E+02 3.259E-20
 1.50E+02 3.075E-20
 1.74E+02 2.881E-20
 2.02E+02 2.682E-20
 2.34E+02 2.484E-20
 2.72E+02 2.289E-20
 3.15E+02 2.101E-20
 3.65E+02 1.922E-20
 4.24E+02 1.751E-20
 4.91E+02 1.592E-20
 5.70E+02 1.443E-20
 6.61E+02 1.305E-20
 7.67E+02 1.177E-20
 8.89E+02 1.06E-20
 1.03E+03 9.526E-21
 1.20E+03 8.547E-21
 1.39E+03 7.658E-21
 1.61E+03 6.851E-21
 1.87E+03 6.121E-21
 2.16E+03 5.462E-21
 2.51E+03 4.868E-21
 2.91E+03 4.334E-21
 3.38E+03 3.855E-21
 3.92E+03 3.426E-21
 4.54E+03 3.041E-21
 5.27E+03 2.698E-21
 6.11E+03 2.391E-21
 7.09E+03 2.118E-21
 8.22E+03 1.875E-21
 9.53E+03 1.658E-21
 1.11E+04 1.466E-21
 1.28E+04 1.295E-21
 1.49E+04 1.143E-21
 1.72E+04 1.009E-21
 2.00E+04 8.898E-22
--- a/data/collisions/IO_e-Xe.dat
+++ b/data/collisions/IO_e-Xe.dat
@ -0,0 +1,52 @@
 # EL cross sections extracted from PROGRAM MAGBOLTZ, VERSION 7.1 JUNE 2004 www.lxcat.net/Biagi-v7.1
 # Relative energy (eV) cross section (m^2)
 1.21E+01 0
 1.41E+01 3.923E-21
 1.64E+01 1.194E-20
 1.91E+01 2.1E-20
 2.22E+01 2.946E-20
 2.58E+01 3.65E-20
 3.00E+01 4.185E-20
 3.49E+01 4.552E-20
 4.06E+01 4.766E-20
 4.72E+01 4.85E-20
 5.49E+01 4.828E-20
 6.39E+01 5.031E-20
 7.43E+01 5.1E-20
 8.64E+01 5.1E-20
 1.01E+02 5.032E-20
 1.17E+02 4.906E-20
 1.36E+02 4.732E-20
 1.58E+02 4.521E-20
 1.84E+02 4.283E-20
 2.14E+02 4.029E-20
 2.49E+02 3.764E-20
 2.90E+02 3.497E-20
 3.37E+02 3.233E-20
 3.92E+02 2.975E-20
 4.56E+02 2.726E-20
 5.31E+02 2.489E-20
 6.17E+02 2.266E-20
 7.18E+02 2.056E-20
 8.35E+02 1.861E-20
 9.72E+02 1.68E-20
 1.13E+03 1.514E-20
 1.32E+03 1.361E-20
 1.53E+03 1.221E-20
 1.78E+03 1.094E-20
 2.07E+03 9.781E-21
 2.41E+03 8.735E-21
 2.80E+03 7.789E-21
 3.26E+03 6.938E-21
 3.79E+03 6.171E-21
 4.41E+03 5.484E-21
 5.13E+03 4.868E-21
 5.97E+03 4.316E-21
 6.94E+03 3.824E-21
 8.07E+03 3.385E-21
 9.39E+03 2.994E-21
 1.09E+04 2.646E-21
 1.27E+04 2.336E-21
 1.48E+04 2.062E-21
 1.72E+04 1.818E-21
 2.00E+04 1.602E-21
--- a/doc/logos/fpakc.pdf
+++ b/doc/logos/fpakc.pdf
--- a/doc/user-manual/.gitignore
+++ b/doc/user-manual/.gitignore
@ -6,6 +6,7 @@
 *.aux
 *.ps
 bibliography.bib.bak
 bibliography.bib.sav
 *.bbl
 *.blg
 *.out
--- a/doc/user-manual/figures/logos/PPartiC.eps
+++ b/doc/user-manual/figures/logos/PPartiC.eps
@ -1,177 +0,0 @@
 %!PS-Adobe-3.0 EPSF-3.0
 %%Creator: cairo 1.16.0 (https://cairographics.org)
 %%CreationDate: Wed Oct 07 12:01:13 2020
 %%Pages: 1
 %%DocumentData: Clean7Bit
 %%LanguageLevel: 2
 %%BoundingBox: 6 7 205 57
 %%EndComments
 %%BeginProlog
 50 dict begin
 /q { gsave } bind def
 /Q { grestore } bind def
 /cm { 6 array astore concat } bind def
 /w { setlinewidth } bind def
 /J { setlinecap } bind def
 /j { setlinejoin } bind def
 /M { setmiterlimit } bind def
 /d { setdash } bind def
 /m { moveto } bind def
 /l { lineto } bind def
 /c { curveto } bind def
 /h { closepath } bind def
 /re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
      0 exch rlineto 0 rlineto closepath } bind def
 /S { stroke } bind def
 /f { fill } bind def
 /f* { eofill } bind def
 /n { newpath } bind def
 /W { clip } bind def
 /W* { eoclip } bind def
 /BT { } bind def
 /ET { } bind def
 /BDC { mark 3 1 roll /BDC pdfmark } bind def
 /EMC { mark /EMC pdfmark } bind def
 /cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
 /Tj { show currentpoint cairo_store_point } bind def
 /TJ {
  {
    dup
    type /stringtype eq
    { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
  } forall
  currentpoint cairo_store_point
 } bind def
 /cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
    cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
 /Tf { pop /cairo_font exch def /cairo_font_matrix where
      { pop cairo_selectfont } if } bind def
 /Td { matrix translate cairo_font_matrix matrix concatmatrix dup
      /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
      /cairo_font where { pop cairo_selectfont } if } bind def
 /Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
      cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
 /g { setgray } bind def
 /rg { setrgbcolor } bind def
 /d1 { setcachedevice } bind def
 /cairo_data_source {
  CairoDataIndex CairoData length lt
    { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def }
    { () } ifelse
 } def
 /cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def
 /cairo_image { image cairo_flush_ascii85_file } def
 /cairo_imagemask { imagemask cairo_flush_ascii85_file } def
 %%EndProlog
 %%BeginSetup
 %%EndSetup
 %%Page: 1 1
 %%BeginPageSetup
 %%PageBoundingBox: 6 7 205 57
 %%EndPageSetup
 q 6 7 199 50 rectclip
 1 0 0 -1 0 71 cm q
 0 g
 49.766 27.305 m 49.766 29.898 49.078 31.813 47.703 33.039 c 46.328 34.27
 44.297 34.883 41.609 34.883 c 36.625 34.883 l 36.625 47.945 l 32.438 47.945
 l 32.438 14.352 l 41.516 14.352 l 44.234 14.352 46.285 14.953 47.672 16.148
 c 49.066 17.348 49.766 19.242 49.766 21.836 c h
 45.594 21.398 m 45.594 20.18 45.313 19.301 44.75 18.758 c 44.195 18.219
 43.313 17.945 42.094 17.945 c 36.625 17.945 l 36.625 31.289 l 42.094 31.289
 l 43.313 31.289 44.195 31.008 44.75 30.445 c 45.313 29.883 45.594 29 45.594
 27.789 c h
 45.594 21.398 m f
 75.641 27.305 m 75.641 29.898 74.953 31.813 73.578 33.039 c 72.203 34.27
 70.172 34.883 67.484 34.883 c 62.5 34.883 l 62.5 47.945 l 58.313 47.945
 l 58.313 14.352 l 67.391 14.352 l 70.109 14.352 72.16 14.953 73.547 16.148
 c 74.941 17.348 75.641 19.242 75.641 21.836 c h
 71.469 21.398 m 71.469 20.18 71.188 19.301 70.625 18.758 c 70.07 18.219
 69.188 17.945 67.969 17.945 c 62.5 17.945 l 62.5 31.289 l 67.969 31.289
 l 69.188 31.289 70.07 31.008 70.625 30.445 c 71.188 29.883 71.469 29 71.469
 27.789 c h
 71.469 21.398 m f
 99.891 47.945 m 99.473 47.945 99.016 47.82 98.516 47.57 c 98.023 47.313
 97.602 46.957 97.25 46.508 c 96.508 47.469 95.535 47.945 94.328 47.945 
 c 90.047 47.945 l 87.617 47.945 85.852 47.461 84.75 46.492 c 83.645 45.516
 83.094 43.836 83.094 41.461 c 83.094 40.508 l 83.094 35.551 85.426 33.07
 90.094 33.07 c 96.047 33.07 l 96.047 30.664 l 96.047 29.676 95.742 28.91
 95.141 28.367 c 94.535 27.816 93.641 27.539 92.453 27.539 c 85.016 27.539
 l 85.016 23.945 l 91.969 23.945 l 94.75 23.945 96.801 24.555 98.125 25.773
 c 99.457 26.984 100.125 28.938 100.125 31.633 c 100.125 43.148 l 100.125
 43.629 100.305 43.988 100.672 44.227 c 101.047 44.469 101.664 44.586 102.531
 44.586 c 102.531 47.945 l h
 93.844 44.352 m 94.801 44.352 95.406 44.129 95.656 43.68 c 95.914 43.234
 96.047 42.656 96.047 41.945 c 96.047 36.664 l 90.047 36.664 l 89.211 36.664
 88.523 36.941 87.984 37.492 c 87.441 38.035 87.172 38.719 87.172 39.539
 c 87.172 41.945 l 87.172 42.781 87.363 43.391 87.75 43.773 c 88.133 44.16
 88.742 44.352 89.578 44.352 c h
 93.844 44.352 m f
 114.047 44.352 m 114.047 27.539 l 110.547 27.539 l 110.547 23.945 l 117.891
 23.945 l 117.891 25.242 l 118.953 24.379 120.234 23.945 121.734 23.945 
 c 126.672 23.945 l 126.672 27.789 l 121.25 27.789 l 120.227 27.789 119.461
 28.117 118.953 28.773 c 118.441 29.43 118.172 30.094 118.141 30.758 c 118.141
 44.352 l 123.359 44.352 l 123.359 47.945 l 110.547 47.945 l 110.547 44.352
 l h
 114.047 44.352 m f
 146.359 47.945 m 144.723 47.945 143.473 47.422 142.609 46.367 c 141.754
 45.305 141.328 43.91 141.328 42.18 c 141.328 27.539 l 136.766 27.539 l 
 136.766 23.945 l 141.328 23.945 l 141.328 18.18 l 145.406 18.18 l 145.406
 23.945 l 152.406 23.945 l 152.406 27.539 l 145.406 27.539 l 145.406 42.43
 l 145.406 43.129 145.578 43.625 145.922 43.914 c 146.273 44.207 146.82 
 44.352 147.563 44.352 c 152.406 44.352 l 152.406 47.945 l h
 146.359 47.945 m f
 172.719 14.352 m 172.719 18.664 l 168.156 18.664 l 168.156 14.352 l h
 161.344 27.539 m 161.344 23.945 l 172.469 23.945 l 172.469 44.352 l 179.578
 44.352 l 179.578 47.945 l 168.391 47.945 l 168.391 27.539 l h
 161.344 27.539 m f
 197.25 47.945 m 194.469 47.945 192.332 47.27 190.844 45.914 c 189.352 44.551
 188.609 42.508 188.609 39.789 c 188.609 22.508 l 188.609 19.789 189.352
 17.754 190.844 16.398 c 192.332 15.035 194.469 14.352 197.25 14.352 c 204.063
 14.352 l 204.063 17.992 l 196.766 17.992 l 195.578 17.992 194.613 18.348
 193.875 19.055 c 193.145 19.754 192.781 20.598 192.781 21.586 c 192.781
 40.695 l 192.781 41.688 193.145 42.535 193.875 43.242 c 194.613 43.953 
 195.578 44.305 196.766 44.305 c 204.063 44.305 l 204.063 47.945 l h
 197.25 47.945 m f
 0.831373 0 0 rg
 17.277 53.707 m 17.277 56.629 14.91 58.996 11.988 58.996 c 9.066 58.996
 6.695 56.629 6.695 53.707 c 6.695 50.785 9.066 48.414 11.988 48.414 c 14.91
 48.414 17.277 50.785 17.277 53.707 c f
 0 0 0.501961 rg
 23.906 61.84 m 23.906 62.969 22.992 63.883 21.863 63.883 c 20.734 63.883
 19.82 62.969 19.82 61.84 c 19.82 60.715 20.734 59.801 21.863 59.801 c 22.992
 59.801 23.906 60.715 23.906 61.84 c f
 0.596078 g
 0.751181 w
 0 J
 0 j
 [] 0.0 d
 4 M q 1 0 0 1 0 0 cm
 16.941 50.758 m 29.465 47.813 l S Q
 27.848 47.258 m 30.453 47.578 l 28.262 49.023 l 28.523 48.41 28.355 47.699
 27.848 47.258 c h
 27.848 47.258 m f*
 0.137102 w
 1 j
 q -1 0.235294 -0.235294 -1 0 0 cm
 -15.851 -50.987 m -18.248 -51.872 l -15.849 -52.753 l -16.234 -52.23 -16.233
 -51.519 -15.851 -50.987 c h
 -15.851 -50.987 m S Q
 0.6 g
 0.751178 w
 0 j
 q 1 0 0 1 0 0 cm
 23.234 59.613 m 30.398 49.301 l S Q
 28.824 49.969 m 30.973 48.461 l 30.313 51 l 30.098 50.375 29.496 49.957
 28.824 49.969 c h
 28.824 49.969 m f*
 0.115667 w
 1 j
 q -0.694805 1 -1 -0.694805 0 0 cm
 20.193 -42.855 m 18.17 -43.597 l 20.191 -44.342 l 19.87 -43.904 19.87 -43.302
 20.193 -42.855 c h
 20.193 -42.855 m S Q
 Q Q
 showpage
 %%Trailer
 end
 %%EOF
--- a/doc/user-manual/figures/logos/PPartiC.svg
+++ b/doc/user-manual/figures/logos/PPartiC.svg
@ -1,159 +0,0 @@
 <?xml version="1.0" encoding="UTF-8" standalone="no"?>
 <svg
   xmlns:dc="http://purl.org/dc/elements/1.1/"
   xmlns:cc="http://creativecommons.org/ns#"
   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
   xmlns:svg="http://www.w3.org/2000/svg"
   xmlns="http://www.w3.org/2000/svg"
   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
   width="80mm"
   height="25mm"
   viewBox="0 0 80 25"
   version="1.1"
   id="svg8"
   inkscape:version="1.0 (4035a4fb49, 2020-05-01)"
   sodipodi:docname="PPartiC.svg">
  <defs
     id="defs2">
    <marker
       inkscape:stockid="Arrow2Mend"
       orient="auto"
       refY="0.0"
       refX="0.0"
       id="Arrow2Mend"
       style="overflow:visible;"
       inkscape:isstock="true">
      <path
         id="path881"
         style="fill-rule:evenodd;stroke-width:0.625;stroke-linejoin:round;stroke:#000000;stroke-opacity:1;fill:#000000;fill-opacity:1"
         d="M 8.7185878,4.0337352 L -2.2072895,0.016013256 L 8.7185884,-4.0017078 C 6.9730900,-1.6296469 6.9831476,1.6157441 8.7185878,4.0337352 z "
         transform="scale(0.6) rotate(180) translate(0,0)" />
    </marker>
    <marker
       inkscape:stockid="Arrow2Send"
       orient="auto"
       refY="0.0"
       refX="0.0"
       id="Arrow2Send"
       style="overflow:visible;"
       inkscape:isstock="true">
      <path
         id="path887"
         style="fill-rule:evenodd;stroke-width:0.625;stroke-linejoin:round;stroke:#989898;stroke-opacity:1;fill:#989898;fill-opacity:1"
         d="M 8.7185878,4.0337352 L -2.2072895,0.016013256 L 8.7185884,-4.0017078 C 6.9730900,-1.6296469 6.9831476,1.6157441 8.7185878,4.0337352 z "
         transform="scale(0.3) rotate(180) translate(-2.3,0)" />
    </marker>
    <marker
       inkscape:stockid="Arrow1Mend"
       orient="auto"
       refY="0.0"
       refX="0.0"
       id="Arrow1Mend"
       style="overflow:visible;"
       inkscape:isstock="true">
      <path
         id="path863"
         d="M 0.0,0.0 L 5.0,-5.0 L -12.5,0.0 L 5.0,5.0 L 0.0,0.0 z "
         style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt;stroke-opacity:1;fill:#000000;fill-opacity:1"
         transform="scale(0.4) rotate(180) translate(10,0)" />
    </marker>
    <marker
       inkscape:stockid="Arrow1Lend"
       orient="auto"
       refY="0.0"
       refX="0.0"
       id="Arrow1Lend"
       style="overflow:visible;"
       inkscape:isstock="true">
      <path
         id="path857"
         d="M 0.0,0.0 L 5.0,-5.0 L -12.5,0.0 L 5.0,5.0 L 0.0,0.0 z "
         style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt;stroke-opacity:1;fill:#000000;fill-opacity:1"
         transform="scale(0.8) rotate(180) translate(12.5,0)" />
    </marker>
    <marker
       inkscape:isstock="true"
       style="overflow:visible"
       id="Arrow2Send-8"
       refX="0"
       refY="0"
       orient="auto"
       inkscape:stockid="Arrow2Send">
      <path
         transform="matrix(-0.3,0,0,-0.3,0.69,0)"
         d="M 8.7185878,4.0337352 -2.2072895,0.01601326 8.7185884,-4.0017078 c -1.7454984,2.3720609 -1.7354408,5.6174519 -6e-7,8.035443 z"
         style="fill:#999999;fill-opacity:1;fill-rule:evenodd;stroke:#999999;stroke-width:0.625;stroke-linejoin:round;stroke-opacity:1"
         id="path887-3" />
    </marker>
  </defs>
  <sodipodi:namedview
     id="base"
     pagecolor="#ffffff"
     bordercolor="#666666"
     borderopacity="1.0"
     inkscape:pageopacity="0.0"
     inkscape:pageshadow="2"
     inkscape:zoom="2.8"
     inkscape:cx="-108.81501"
     inkscape:cy="51.644646"
     inkscape:document-units="mm"
     inkscape:current-layer="layer1"
     inkscape:document-rotation="0"
     showgrid="false"
     inkscape:window-width="2560"
     inkscape:window-height="1378"
     inkscape:window-x="-9"
     inkscape:window-y="26"
     inkscape:window-maximized="1" />
  <metadata
     id="metadata5">
    <rdf:RDF>
      <cc:Work
         rdf:about="">
        <dc:format>image/svg+xml</dc:format>
        <dc:type
           rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
        <dc:title></dc:title>
      </cc:Work>
    </rdf:RDF>
  </metadata>
  <g
     inkscape:label="Layer 1"
     inkscape:groupmode="layer"
     id="layer1">
    <text
       xml:space="preserve"
       style="font-style:normal;font-weight:normal;font-size:16.9333px;line-height:1.25;font-family:sans-serif;fill:#000000;fill-opacity:1;stroke:none;stroke-width:0.264583"
       x="9.921875"
       y="16.914433"
       id="text12"><tspan
         sodipodi:role="line"
         id="tspan10"
         x="9.921875"
         y="16.914433"
         style="font-style:italic;font-variant:normal;font-weight:normal;font-stretch:normal;font-size:16.9333px;font-family:'Share Tech Mono';-inkscape-font-specification:'Share Tech Mono Italic';stroke-width:0.264583">PPartiC</tspan></text>
    <circle
       style="fill:#d40000;stroke-width:0.264583"
       id="path14"
       cx="4.2286086"
       cy="18.946056"
       r="1.8662573" />
    <circle
       style="fill:#000080;stroke-width:0.576413"
       id="path16"
       cx="7.7130766"
       cy="21.816313"
       r="0.7205171" />
    <path
       style="fill:none;fill-opacity:1;stroke:#989898;stroke-width:0.265;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1;marker-end:url(#Arrow2Send)"
       d="M 5.9767484,17.906622 10.394345,16.867188"
       id="path18"
       sodipodi:nodetypes="cc" />
    <path
       sodipodi:nodetypes="cc"
       id="path18-6"
       d="m 8.195956,21.029626 2.527716,-3.638021"
       style="fill:none;stroke:#999999;stroke-width:0.264999;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1;marker-end:url(#Arrow2Send-8)" />
  </g>
 </svg>
--- a/doc/user-manual/figures/logos/fpakc.eps
+++ b/doc/user-manual/figures/logos/fpakc.eps
@ -1,294 +0,0 @@
 %!PS-Adobe-3.0 EPSF-3.0
 %%Creator: cairo 1.16.0 (https://cairographics.org)
 %%CreationDate: Sat Nov 21 21:23:19 2020
 %%Pages: 1
 %%DocumentData: Clean7Bit
 %%LanguageLevel: 2
 %%BoundingBox: 8 15 191 78
 %%EndComments
 %%BeginProlog
 50 dict begin
 /q { gsave } bind def
 /Q { grestore } bind def
 /cm { 6 array astore concat } bind def
 /w { setlinewidth } bind def
 /J { setlinecap } bind def
 /j { setlinejoin } bind def
 /M { setmiterlimit } bind def
 /d { setdash } bind def
 /m { moveto } bind def
 /l { lineto } bind def
 /c { curveto } bind def
 /h { closepath } bind def
 /re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
      0 exch rlineto 0 rlineto closepath } bind def
 /S { stroke } bind def
 /f { fill } bind def
 /f* { eofill } bind def
 /n { newpath } bind def
 /W { clip } bind def
 /W* { eoclip } bind def
 /BT { } bind def
 /ET { } bind def
 /BDC { mark 3 1 roll /BDC pdfmark } bind def
 /EMC { mark /EMC pdfmark } bind def
 /cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
 /Tj { show currentpoint cairo_store_point } bind def
 /TJ {
  {
    dup
    type /stringtype eq
    { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
  } forall
  currentpoint cairo_store_point
 } bind def
 /cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
    cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
 /Tf { pop /cairo_font exch def /cairo_font_matrix where
      { pop cairo_selectfont } if } bind def
 /Td { matrix translate cairo_font_matrix matrix concatmatrix dup
      /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
      /cairo_font where { pop cairo_selectfont } if } bind def
 /Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
      cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
 /g { setgray } bind def
 /rg { setrgbcolor } bind def
 /d1 { setcachedevice } bind def
 /cairo_data_source {
  CairoDataIndex CairoData length lt
    { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def }
    { () } ifelse
 } def
 /cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def
 /cairo_image { image cairo_flush_ascii85_file } def
 /cairo_imagemask { imagemask cairo_flush_ascii85_file } def
 %%EndProlog
 %%BeginSetup
 %%BeginResource: font f-0-0
 %!FontType1-1.1 f-0-0 1.0
 11 dict begin
 /FontName /f-0-0 def
 /PaintType 0 def
 /FontType 1 def
 /FontMatrix [0.001 0 0 0.001 0 0] readonly def
 /FontBBox {30 -240 640 735                                   } readonly def
 /Encoding 256 array
 0 1 255 {1 index exch /.notdef put} for
 dup 97 /a put
 dup 99 /c put
 dup 102 /f put
 dup 107 /k put
 dup 112 /p put
 readonly def
 currentdict end
 currentfile eexec
 f983ef0097ece636fb4a96c74d26ab84185f6dfa4a16a7a1c27bbe3f1156aea698df336d20b467
 b10e7f33846656653c5ac6962759d3056cbdb3190bac614b984bf5a132dc418192443014ba63de
 800d392b6fea026574bb2535fd7bb5338f35bf15a88ea328fdaa49670c7852e3d060f3c5d6b07f
 2ef6d0f22646c5d18e19a2ae3ee120390f6dd96f76dcf1e127de5e9299077a00c17c0d71e36e5b
 9d5ec58fceda57739a6a4214d4b79d6c48d2784b60c320323c7acddddf34db833cac0cf109f799
 69d114a330d372e5c978a66acc84e3fe5557f6240856a013ffaa0199444e5c5036f775eba4a5c5
 8cde66cf604b9aca2178431127b8a1ff7ed633a65c04600af5f573483112251caad907dcd8c61e
 23500065b1568be79a17d2379b63d656e8f7715cdfdf357f0e30d9ab91d113c88231d875d60897
 1dee27eb5da34a09d2f1b1a2e7ab8be1036393cd6ae53eff0d77f452e9bf45eccc9e4836ae77c5
 b74262fa43bfccfd9a147a18dcdae6e50cc518129a64f669a5fae69c8082dec571d7d01a4e9f05
 6d3e52de0ea004acdbd1b3bf9b19aa17f839a75365a77b7275442a967093ffdf1694a72f9978a2
 304ac331d401b48e3e62a3a92cd39516dd480f088980d1ad8993a1f6fefb07e0d86b6f0293bb41
 68ac465726267cacb7516a0e910fe67c2dbfef06d8b64a9811506650d32fa182a0adcab8e2e21e
 ca6d0dc81959c25ea2d3f7ccec13e0cb4a7ef88e97c36e74fa13010220d6835ebdcbabdb507d84
 239e5483e8a8b7a52d6e1ea4ea1f5e6bef4534710c4055265aaa86fb445f3b2fc62cfdd9e283d8
 8bd083d09f0971cde00f2031b58b304d5f647f02aabf7ba9062c33979cd391f692c72ee179b7a8
 16f9c9e668d20021bd2b6a0f0114898729c6228be2895a696aaef0ebbcc842e64d5e72cc1d9b75
 44314028987a238f8fc4c18a0db3546c9ea42194b6bbdc45587e36d605fe2b7608d9292ddc0c9b
 be3e420b36fd52f9aef97f13533e101f34d4f882848f4845a7a824da815a710abe11a1ec8363c1
 06daa18dffb5a451af7bf3b20d79a63c94e305050ea893b7a61cc8cfdb2e8a593a073c2021b298
 40863c70742ab1734e1d6811cf1927832da10f562f6895575b50044179588ddd9ec4c413f68c3e
 3063f9594dc94115af4d9d4d6259c2ebb5afa796131772de3d297a8cc04f7f10398acc9142b1aa
 2da9741ad314918ff1553dafe4751b4c0efdc9d6a549acbf1b3d209f6ebe8f6561d627f37bbce1
 7213b92bf332c27718ca9f868f1724cda0774ea4c3a5a2ba99509eb9128c456e5526f234dc3adc
 37ac61ead9dafba1b5d58a9443ceb92474535cd3515e9ce357420b230fe927e81f06b2363c70aa
 b6e00858a44972ad3f8759069235bba0b8ae2c65a59fe3ee5642f88a8550a765907eb4f9432ac4
 9e896114d0bc969bc2c14acd9a50c31e2095133b6b4fc11a1136dbba4b515eaabf0cba23ffe795
 3532a1fca89780a841f3a5fe2514d31fd6d41fbcb5b8caadd53c1fa7b06506963f37006269d0ae
 fc1d5d6bd7f6788544e01a77bdf35aaacdd41dd4fc16237759516c60ee9b57e7b56606e0fb5a25
 9e6cb3b2e22d3ac4db73c228fda1b327cace7cbe25594ea2a9445efeda7826604e3daf18ca977a
 9a2788cddcea95b5460b648c12bebc3c39302a07481fe18ceed4c9e12ac7d51f8104bd589cf9f6
 68371d3aa8510b5ff08c972b84fd4e48c356f9098b9d130085f59bc7e748ac59e9715060d29b38
 d828e479d8a85e0aaaa7b4ab34b547643b9c3417305b29be240ba9514cd4656079f26f378a6af9
 f7c694b1332a14baa203faa5c6745eea7c9cc0b3fc13f722ff941857d141db5665b32ada8b7831
 06f5c01361bd4b161450546ad8812d364b394658c4d3e1e2cf1daba58315787fbe299a27a387f2
 1545b5d4bc3022e6a6803e971698504ebc44c9f7c02256a5ae01c9fb3194baf05fb7c9cd862427
 95ca9cfd57ba76d6697fc842ee19ff35ce8b9d3299a185a4d28537058b5ea471a818f521183ffa
 4b171f078577707cfbb3b2c6de4659b6dbacb952a0213b0586bd5879096ec53540ad07463e5b25
 1a38d715de6a0bb6f92d3e20e1b25be654fe7aad2e057428b68e2d6830f2ee0ca384fb862813ba
 d48832ec89991fa5e6bf9d3135378f0000000000000000000000000000000000000000000000000000000000000000
 0000000000000000000000000000000000000000000000000000000000000000
 0000000000000000000000000000000000000000000000000000000000000000
 0000000000000000000000000000000000000000000000000000000000000000
 0000000000000000000000000000000000000000000000000000000000000000
 0000000000000000000000000000000000000000000000000000000000000000
 0000000000000000000000000000000000000000000000000000000000000000
 0000000000000000000000000000000000000000000000000000000000000000
 cleartomark
 %%EndResource
 %%EndSetup
 %%Page: 1 1
 %%BeginPageSetup
 %%PageBoundingBox: 8 15 191 78
 %%EndPageSetup
 q 8 15 183 63 rectclip
 1 0 0 -1 0 86 cm q
 0 g
 BT
 56.000126 0 0 -56.000126 52.419201 52.791732 Tm
 /f-0-0 1 Tf
 [(fpak)10(c)]TJ
 ET
 0.749999 w
 0 J
 0 j
 [] 0.0 d
 4 M q 1 0 0 1 0 0 cm
 8.996 27.238 m 20.621 9.359 l 34.656 9.363 l 34.621 14.746 l 25.941 18.195
 l 17.332 33.629 l 17.41 55.852 l 9.316 55.785 l h
 8.996 27.238 m S Q
 q 1 0 0 1 0 0 cm
 21.156 69.262 m 21.156 38.387 l 40.305 38.387 l 40.305 53.027 l 27.457 
 53.027 l 27.457 69.52 l h
 21.156 69.262 m S Q
 q 1 0 0 1 0 0 cm
 27.898 43.988 m 27.898 47.73 l 33.426 47.73 l 33.426 43.898 l h
 27.898 43.988 m S Q
 0.283465 w
 q 1 0 0 1 0 0 cm
 9.094 50.535 m 17.414 50.297 l S Q
 q 1 0 0 1 0 0 cm
 9.34 43.07 m 17.094 44.93 l S Q
 q 1 0 0 1 0 0 cm
 9.43 34.727 m 17.348 41.363 l S Q
 q 1 0 0 1 0 0 cm
 17.332 33.629 m 8.996 27.238 l S Q
 q 1 0 0 1 0 0 cm
 13.297 31.055 m 11.32 36.062 l 12.719 43.309 l 14.023 50.113 l 12.02 55.516
 l S Q
 q 1 0 0 1 0 0 cm
 17.332 33.629 m 12.203 21.801 l S Q
 q 1 0 0 1 0 0 cm
 21.246 26.223 m 12.203 21.801 l S Q
 q 1 0 0 1 0 0 cm
 15.379 16.828 m 21.246 26.223 l S Q
 q 1 0 0 1 0 0 cm
 23.562 22.152 m 15.379 16.828 l S Q
 q 1 0 0 1 0 0 cm
 25.941 18.195 m 15.379 16.828 l S Q
 q 1 0 0 1 0 0 cm
 13.297 31.055 m 14.234 26.23 l S Q
 q 1 0 0 1 0 0 cm
 14.234 26.23 m 21.246 26.223 l S Q
 q 1 0 0 1 0 0 cm
 25.941 18.195 m 17.918 13.906 l S Q
 q 1 0 0 1 0 0 cm
 20.621 9.359 m 25.941 18.195 l S Q
 q 1 0 0 1 0 0 cm
 34.656 9.363 m 25.941 18.195 l S Q
 q 1 0 0 1 0 0 cm
 25.941 18.195 m 26.527 9.582 l S Q
 1 0 0 rg
 23.523 67.027 m 23.523 67.027 23.523 67.031 23.523 67.031 c 23.52 67.031
 23.52 67.027 23.52 67.027 c 23.52 67.027 23.52 67.027 23.523 67.027 c h
 23.523 67.027 m f
 0 g
 0.749999 w
 q 1 0 0 1 0 0 cm
 23.523 67.027 m 23.523 67.027 23.523 67.031 23.523 67.031 c 23.52 67.031
 23.52 67.027 23.52 67.027 c 23.52 67.027 23.52 67.027 23.523 67.027 c h
 23.523 67.027 m S Q
 1 0 0 rg
 24.66 63.297 m 24.66 63.297 24.656 63.297 24.656 63.297 c 24.656 63.293
 24.66 63.293 24.66 63.293 c 24.66 63.293 24.66 63.293 24.66 63.297 c h
 24.66 63.297 m f
 0 g
 q 1 0 0 1 0 0 cm
 24.66 63.297 m 24.66 63.297 24.656 63.297 24.656 63.297 c 24.656 63.293
 24.66 63.293 24.66 63.293 c 24.66 63.293 24.66 63.293 24.66 63.297 c h
 24.66 63.297 m S Q
 0 0 1 rg
 23.797 64.039 m 23.797 64.273 23.609 64.465 23.371 64.465 c 23.137 64.465
 22.949 64.273 22.949 64.039 c 22.949 63.805 23.137 63.613 23.371 63.613
 c 23.609 63.613 23.797 63.805 23.797 64.039 c h
 23.797 64.039 m f
 1 0 0 rg
 37.844 42.566 m 37.844 42.945 37.535 43.254 37.156 43.254 c 36.773 43.254
 36.465 42.945 36.465 42.566 c 36.465 42.184 36.773 41.875 37.156 41.875
 c 37.535 41.875 37.844 42.184 37.844 42.566 c h
 37.844 42.566 m f
 0 g
 0.141732 w
 q 1 0 0 1 0 0 cm
 35.977 42.188 m 31.02 39.777 l S Q
 31.531 40.027 m 31.91 39.895 l 30.895 39.715 l 31.66 40.406 l h
 31.531 40.027 m f*
 0.0679747 w
 q 1 0.486446 -0.486446 1 0 0 cm
 41.243 19.965 m 41.497 19.708 l 40.605 19.963 l 41.496 20.22 l h
 41.243 19.965 m S Q
 0.141732 w
 q 1 0 0 1 0 0 cm
 29.828 39.73 m 23.582 45.012 l S Q
 24.016 44.645 m 24.047 44.246 l 23.473 45.102 l 24.414 44.68 l h
 24.016 44.645 m f*
 0.0577315 w
 q 1 -0.845211 0.845211 1 0 0 cm
 -8.002 37.881 m -7.787 37.664 l -8.544 37.88 l -7.787 38.098 l h
 -8.002 37.881 m S Q
 0.141732 w
 q 1 0 0 1 0 0 cm
 23.82 46.355 m 25.027 51.145 l S Q
 24.887 50.594 m 24.543 50.387 l 25.059 51.281 l 25.094 50.25 l h
 24.887 50.594 m f*
 0.0732958 w
 q -0.252177 -1 1 -0.252177 0 0 cm
 -53.469 11.403 m -53.193 11.129 l -54.156 11.402 l -53.195 11.679 l h
 -53.469 11.403 m S Q
 0.141732 w
 q 1 0 0 1 0 0 cm
 23.453 63.156 m 26.16 58.465 l S Q
 25.875 58.953 m 25.98 59.34 l 26.23 58.34 l 25.488 59.059 l h
 25.875 58.953 m f*
 0.0654792 w
 q -0.576787 1 -1 -0.576787 0 0 cm
 33.038 -44.931 m 33.282 -45.177 l 32.424 -44.932 l 33.284 -44.686 l h
 33.038 -44.931 m S Q
 0.141732 w
 q 1 0 0 1 0 0 cm
 26.215 57.066 m 26.281 52.215 l S Q
 26.273 52.785 m 26.555 53.07 l 26.285 52.074 l 25.984 53.062 l h
 26.273 52.785 m f*
 0.0755832 w
 q -0.0137977 1 -1 -0.0137977 0 0 cm
 52.413 -26.997 m 52.694 -27.282 l 51.702 -26.999 l 52.694 -26.711 l h
 52.413 -26.997 m S Q
 0.141732 w
 q 1 0 0 1 0 0 cm
 26.5 50.707 m 34.664 48.957 l S Q
 34.109 49.074 m 33.895 49.41 l 34.805 48.926 l 33.773 48.859 l h
 34.109 49.074 m f*
 0.0739077 w
 q -1 0.214603 -0.214603 -1 0 0 cm
 -22.54 -53.911 m -22.266 -54.188 l -23.235 -53.912 l -22.263 -53.637 l 
 h
 -22.54 -53.911 m S Q
 0.141732 w
 q 1 0 0 1 0 0 cm
 36.809 49.66 m 38.371 43.828 l S Q
 38.223 44.379 m 38.426 44.727 l 38.406 43.691 l 37.879 44.578 l h
 38.223 44.379 m f*
 0.0730145 w
 q -0.26796 1 -1 -0.26796 0 0 cm
 31.85 -46.757 m 32.123 -47.034 l 31.163 -46.757 l 32.122 -46.486 l h
 31.85 -46.757 m S Q
 Q Q
 showpage
 %%Trailer
 end
 %%EOF
--- a/doc/user-manual/figures/scatteringQuad.eps
+++ b/doc/user-manual/figures/scatteringQuad.eps
@ -1,303 +0,0 @@
 %!PS-Adobe-3.0 EPSF-3.0
 %%Creator: cairo 1.16.0 (https://cairographics.org)
 %%CreationDate: Thu Oct 08 15:41:58 2020
 %%Pages: 1
 %%DocumentData: Clean7Bit
 %%LanguageLevel: 2
 %%BoundingBox: 13 30 125 153
 %%EndComments
 %%BeginProlog
 50 dict begin
 /q { gsave } bind def
 /Q { grestore } bind def
 /cm { 6 array astore concat } bind def
 /w { setlinewidth } bind def
 /J { setlinecap } bind def
 /j { setlinejoin } bind def
 /M { setmiterlimit } bind def
 /d { setdash } bind def
 /m { moveto } bind def
 /l { lineto } bind def
 /c { curveto } bind def
 /h { closepath } bind def
 /re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
      0 exch rlineto 0 rlineto closepath } bind def
 /S { stroke } bind def
 /f { fill } bind def
 /f* { eofill } bind def
 /n { newpath } bind def
 /W { clip } bind def
 /W* { eoclip } bind def
 /BT { } bind def
 /ET { } bind def
 /BDC { mark 3 1 roll /BDC pdfmark } bind def
 /EMC { mark /EMC pdfmark } bind def
 /cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
 /Tj { show currentpoint cairo_store_point } bind def
 /TJ {
  {
    dup
    type /stringtype eq
    { show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
  } forall
  currentpoint cairo_store_point
 } bind def
 /cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
    cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
 /Tf { pop /cairo_font exch def /cairo_font_matrix where
      { pop cairo_selectfont } if } bind def
 /Td { matrix translate cairo_font_matrix matrix concatmatrix dup
      /cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
      /cairo_font where { pop cairo_selectfont } if } bind def
 /Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
      cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
 /g { setgray } bind def
 /rg { setrgbcolor } bind def
 /d1 { setcachedevice } bind def
 /cairo_data_source {
  CairoDataIndex CairoData length lt
    { CairoData CairoDataIndex get /CairoDataIndex CairoDataIndex 1 add def }
    { () } ifelse
 } def
 /cairo_flush_ascii85_file { cairo_ascii85_file status { cairo_ascii85_file flushfile } if } def
 /cairo_image { image cairo_flush_ascii85_file } def
 /cairo_imagemask { imagemask cairo_flush_ascii85_file } def
 %%EndProlog
 %%BeginSetup
 %%EndSetup
 %%Page: 1 1
 %%BeginPageSetup
 %%PageBoundingBox: 13 30 125 153
 %%EndPageSetup
 q 13 30 112 123 rectclip
 1 0 0 -1 0 171 cm q
 0 0.784314 0 rg
 20.836 84.285 52.273 43.371 re f
 0 g
 0.751178 w
 0 J
 0 j
 [] 0.0 d
 4 M q 1 0 0 1 0 0 cm
 20.836 84.285 52.273 43.371 re S Q
 1 0.627451 0 rg
 73.109 84.285 43.184 43.371 re f
 0 g
 q 1 0 0 1 0 0 cm
 73.109 84.285 43.184 43.371 re S Q
 0.0431373 0 1 rg
 73.109 32.199 43.184 52.086 re f
 0 g
 q 1 0 0 1 0 0 cm
 73.109 32.199 43.184 52.086 re S Q
 1 0 0 rg
 20.836 32.199 52.273 52.086 re f
 0 g
 q 1 0 0 1 0 0 cm
 20.836 32.199 52.273 52.086 re S Q
 q 1 0 0 1 0 0 cm
 20.836 32.199 95.457 95.457 re S Q
 0 0 1 rg
 22.902 127.77 m 22.902 128.914 21.973 129.844 20.824 129.844 c 19.68 129.844
 18.75 128.914 18.75 127.77 c 18.75 126.621 19.68 125.691 20.824 125.691
 c 21.973 125.691 22.902 126.621 22.902 127.77 c h
 22.902 127.77 m f
 q 1 0 0 1 0 0 cm
 22.902 127.77 m 22.902 128.914 21.973 129.844 20.824 129.844 c 19.68 129.844
 18.75 128.914 18.75 127.77 c 18.75 126.621 19.68 125.691 20.824 125.691
 c 21.973 125.691 22.902 126.621 22.902 127.77 c h
 22.902 127.77 m S Q
 1 0 0 rg
 118.059 127.23 m 118.059 128.379 117.129 129.309 115.98 129.309 c 114.836
 129.309 113.906 128.379 113.906 127.23 c 113.906 126.086 114.836 125.156
 115.98 125.156 c 117.129 125.156 118.059 126.086 118.059 127.23 c f
 q 1 0 0 1 0 0 cm
 118.059 127.23 m 118.059 128.379 117.129 129.309 115.98 129.309 c 114.836
 129.309 113.906 128.379 113.906 127.23 c 113.906 126.086 114.836 125.156
 115.98 125.156 c 117.129 125.156 118.059 126.086 118.059 127.23 c S Q
 0 0.784314 0 rg
 117.789 32.41 m 117.789 33.559 116.859 34.488 115.715 34.488 c 114.566 
 34.488 113.637 33.559 113.637 32.41 c 113.637 31.266 114.566 30.336 115.715
 30.336 c 116.859 30.336 117.789 31.266 117.789 32.41 c f
 q 1 0 0 1 0 0 cm
 117.789 32.41 m 117.789 33.559 116.859 34.488 115.715 34.488 c 114.566 
 34.488 113.637 33.559 113.637 32.41 c 113.637 31.266 114.566 30.336 115.715
 30.336 c 116.859 30.336 117.789 31.266 117.789 32.41 c S Q
 1 0.627451 0 rg
 23.238 32.41 m 23.238 33.559 22.309 34.488 21.16 34.488 c 20.016 34.488
 19.086 33.559 19.086 32.41 c 19.086 31.266 20.016 30.336 21.16 30.336 c
 22.309 30.336 23.238 31.266 23.238 32.41 c f
 q 1 0 0 1 0 0 cm
 23.238 32.41 m 23.238 33.559 22.309 34.488 21.16 34.488 c 20.016 34.488
 19.086 33.559 19.086 32.41 c 19.086 31.266 20.016 30.336 21.16 30.336 c
 22.309 30.336 23.238 31.266 23.238 32.41 c S Q
 0 g
 76.465 84.191 m 76.465 86.262 74.785 87.941 72.715 87.941 c 70.645 87.941
 68.965 86.262 68.965 84.191 c 68.965 82.121 70.645 80.441 72.715 80.441
 c 74.785 80.441 76.465 82.121 76.465 84.191 c h
 76.465 84.191 m f
 q 1 0 0 1 0 0 cm
 76.465 84.191 m 76.465 86.262 74.785 87.941 72.715 87.941 c 70.645 87.941
 68.965 86.262 68.965 84.191 c 68.965 82.121 70.645 80.441 72.715 80.441
 c 74.785 80.441 76.465 82.121 76.465 84.191 c h
 76.465 84.191 m S Q
 0.749999 w
 q 1 0 0 1 0 0 cm
 20.355 84.105 m 115.98 84.105 l S Q
 q 1 0 0 1 0 0 cm
 72.922 127.277 m 72.922 32.008 l S Q
 14.742 139.539 m 16.68 139.539 l 16.68 132.867 l 14.57 133.289 l 14.57 
 132.211 l 16.664 131.789 l 17.852 131.789 l 17.852 139.539 l 19.789 139.539
 l 19.789 140.539 l 14.742 140.539 l h
 14.742 139.539 m f
 120.676 136.695 m 124.816 136.695 l 124.816 137.695 l 119.254 137.695 l
 119.254 136.695 l 119.699 136.238 120.309 135.617 121.082 134.836 c 121.863
 134.047 122.352 133.535 122.551 133.305 c 122.934 132.891 123.199 132.535
 123.348 132.242 c 123.504 131.941 123.582 131.648 123.582 131.367 c 123.582
 130.898 123.414 130.52 123.082 130.227 c 122.758 129.938 122.336 129.789
 121.816 129.789 c 121.441 129.789 121.043 129.852 120.629 129.977 c 120.223
 130.102 119.785 130.301 119.316 130.57 c 119.316 129.367 l 119.793 129.18
 120.238 129.039 120.645 128.945 c 121.059 128.844 121.441 128.789 121.785
 128.789 c 122.691 128.789 123.414 129.02 123.957 129.477 c 124.496 129.926
 124.77 130.531 124.77 131.289 c 124.77 131.645 124.699 131.984 124.566 
 132.305 c 124.43 132.629 124.184 133.008 123.832 133.445 c 123.727 133.563
 123.414 133.891 122.895 134.43 c 122.371 134.973 121.633 135.727 120.676
 136.695 c h
 120.676 136.695 m f
 122.875 24.07 m 123.438 24.195 123.875 24.453 124.188 24.836 c 124.508 
 25.211 124.672 25.68 124.672 26.242 c 124.672 27.109 124.375 27.781 123.781
 28.258 c 123.188 28.727 122.344 28.961 121.25 28.961 c 120.883 28.961 120.504
 28.922 120.109 28.852 c 119.723 28.781 119.328 28.672 118.922 28.523 c 
 118.922 27.383 l 119.242 27.57 119.598 27.719 119.984 27.82 c 120.379 27.914
 120.789 27.961 121.219 27.961 c 121.957 27.961 122.52 27.816 122.906 27.523
 c 123.301 27.234 123.5 26.805 123.5 26.242 c 123.5 25.734 123.316 25.332
 122.953 25.039 c 122.586 24.75 122.086 24.602 121.453 24.602 c 120.422 
 24.602 l 120.422 23.633 l 121.5 23.633 l 122.07 23.633 122.516 23.52 122.828
 23.289 c 123.141 23.051 123.297 22.711 123.297 22.273 c 123.297 21.828 
 123.133 21.484 122.813 21.242 c 122.5 21.004 122.047 20.883 121.453 20.883
 c 121.117 20.883 120.766 20.922 120.391 20.992 c 120.023 21.055 119.617
 21.16 119.172 21.305 c 119.172 20.258 l 119.629 20.133 120.051 20.039 120.438
 19.977 c 120.832 19.914 121.203 19.883 121.547 19.883 c 122.453 19.883 
 123.164 20.086 123.688 20.492 c 124.207 20.898 124.469 21.453 124.469 22.148
 c 124.469 22.641 124.328 23.051 124.047 23.383 c 123.773 23.719 123.383
 23.945 122.875 24.07 c h
 122.875 24.07 m f
 17.223 19.746 m 14.238 24.418 l 17.223 24.418 l h
 16.91 18.715 m 18.41 18.715 l 18.41 24.418 l 19.66 24.418 l 19.66 25.402
 l 18.41 25.402 l 18.41 27.465 l 17.223 27.465 l 17.223 25.402 l 13.285 
 25.402 l 13.285 24.262 l h
 16.91 18.715 m f
 39.906 47.727 m 38.297 52.07 l 41.516 52.07 l h
 39.234 46.555 m 40.578 46.555 l 43.906 55.305 l 42.672 55.305 l 41.875 
 53.055 l 37.938 53.055 l 37.141 55.305 l 35.891 55.305 l h
 39.234 46.555 m f
 45.5 57.066 m 48.188 57.066 l 48.188 57.707 l 44.578 57.707 l 44.578 57.066
 l 44.867 56.766 45.266 56.359 45.766 55.848 c 46.266 55.34 46.582 55.012
 46.719 54.863 c 46.957 54.582 47.125 54.348 47.219 54.16 c 47.32 53.965
 47.375 53.777 47.375 53.598 c 47.375 53.297 47.266 53.051 47.047 52.863
 c 46.836 52.668 46.566 52.566 46.234 52.566 c 45.992 52.566 45.738 52.609
 45.469 52.691 c 45.195 52.777 44.91 52.902 44.609 53.066 c 44.609 52.301
 l 44.922 52.176 45.207 52.082 45.469 52.02 c 45.738 51.957 45.988 51.926
 46.219 51.926 c 46.801 51.926 47.27 52.074 47.625 52.363 c 47.977 52.656
 48.156 53.047 48.156 53.535 c 48.156 53.777 48.109 54 48.016 54.207 c 47.93
 54.418 47.773 54.66 47.547 54.941 c 47.484 55.016 47.281 55.23 46.938 55.582
 c 46.594 55.938 46.113 56.434 45.5 57.066 c h
 45.5 57.066 m f
 93.52 47.684 m 91.91 52.027 l 95.129 52.027 l h
 92.848 46.512 m 94.191 46.512 l 97.52 55.262 l 96.285 55.262 l 95.488 53.012
 l 91.551 53.012 l 90.754 55.262 l 89.504 55.262 l h
 92.848 46.512 m f
 98.582 57.02 m 99.832 57.02 l 99.832 52.676 l 98.473 52.957 l 98.473 52.254
 l 99.832 51.973 l 100.598 51.973 l 100.598 57.02 l 101.848 57.02 l 101.848
 57.66 l 98.582 57.66 l h
 98.582 57.02 m f
 45.98 98.063 m 44.371 102.406 l 47.59 102.406 l h
 45.309 96.891 m 46.652 96.891 l 49.98 105.641 l 48.746 105.641 l 47.949
 103.391 l 44.012 103.391 l 43.215 105.641 l 41.965 105.641 l h
 45.309 96.891 m f
 53.23 104.98 m 53.605 105.055 53.895 105.215 54.105 105.465 c 54.313 105.715
 54.418 106.023 54.418 106.387 c 54.418 106.949 54.219 107.387 53.824 107.699
 c 53.438 108.004 52.891 108.152 52.184 108.152 c 51.941 108.152 51.699 
 108.125 51.449 108.074 c 51.199 108.035 50.938 107.965 50.668 107.871 c 
 50.668 107.137 l 50.875 107.254 51.105 107.348 51.355 107.418 c 51.613 107.48
 51.887 107.512 52.168 107.512 c 52.645 107.512 53.012 107.418 53.262 107.23
 c 53.52 107.035 53.652 106.754 53.652 106.387 c 53.652 106.055 53.531 105.793
 53.293 105.605 c 53.051 105.418 52.723 105.324 52.309 105.324 c 51.652 
 105.324 l 51.652 104.684 l 52.34 104.684 l 52.723 104.684 53.016 104.613
 53.215 104.465 c 53.41 104.309 53.512 104.09 53.512 103.809 c 53.512 103.52
 53.406 103.293 53.199 103.137 c 52.988 102.98 52.691 102.902 52.309 102.902
 c 52.098 102.902 51.875 102.93 51.637 102.98 c 51.395 103.023 51.129 103.09
 50.84 103.184 c 50.84 102.496 l 51.129 102.414 51.402 102.355 51.652 102.324
 c 51.91 102.285 52.156 102.262 52.387 102.262 c 52.969 102.262 53.426 102.398
 53.762 102.668 c 54.105 102.93 54.277 103.285 54.277 103.73 c 54.277 104.043
 54.184 104.309 53.996 104.527 c 53.816 104.746 53.563 104.898 53.23 104.98
 c h
 53.23 104.98 m f
 93.707 97.875 m 92.098 102.219 l 95.316 102.219 l h
 93.035 96.703 m 94.379 96.703 l 97.707 105.453 l 96.473 105.453 l 95.676
 103.203 l 91.738 103.203 l 90.941 105.453 l 89.691 105.453 l h
 93.035 96.703 m f
 100.754 102.836 m 98.801 105.867 l 100.754 105.867 l h
 100.551 102.164 m 101.52 102.164 l 101.52 105.867 l 102.316 105.867 l 102.316
 106.508 l 101.52 106.508 l 101.52 107.852 l 100.754 107.852 l 100.754 106.508
 l 98.176 106.508 l 98.176 105.773 l h
 100.551 102.164 m f
 27.949 77.395 m 27.949 80.301 l 27.043 80.301 l 27.043 72.754 l 27.949 
 72.754 l 27.949 73.582 l 28.137 73.262 28.371 73.02 28.652 72.863 c 28.941
 72.707 29.293 72.629 29.699 72.629 c 30.363 72.629 30.902 72.895 31.309
 73.426 c 31.723 73.949 31.934 74.637 31.934 75.488 c 31.934 76.355 31.723
 77.051 31.309 77.582 c 30.902 78.105 30.363 78.363 29.699 78.363 c 29.293
 78.363 28.941 78.285 28.652 78.129 c 28.371 77.973 28.137 77.73 27.949 
 77.395 c h
 31.012 75.488 m 31.012 74.832 30.871 74.316 30.59 73.941 c 30.316 73.566
 29.949 73.379 29.48 73.379 c 29 73.379 28.625 73.566 28.355 73.941 c 28.082
 74.316 27.949 74.832 27.949 75.488 c 27.949 76.156 28.082 76.676 28.355
 77.051 c 28.625 77.426 29 77.613 29.48 77.613 c 29.949 77.613 30.316 77.426
 30.59 77.051 c 30.871 76.676 31.012 76.156 31.012 75.488 c h
 31.012 75.488 m f
 35.578 70.629 m 35.148 71.379 34.828 72.125 34.609 72.863 c 34.398 73.594
 34.297 74.332 34.297 75.082 c 34.297 75.832 34.398 76.578 34.609 77.316
 c 34.828 78.059 35.148 78.793 35.578 79.535 c 34.797 79.535 l 34.316 78.773
 33.953 78.027 33.703 77.285 c 33.461 76.547 33.344 75.813 33.344 75.082
 c 33.344 74.355 33.461 73.625 33.703 72.895 c 33.941 72.156 34.305 71.402
 34.797 70.629 c h
 35.578 70.629 m f
 41.875 72.754 m 39.906 75.41 l 41.984 78.223 l 40.922 78.223 l 39.328 76.066
 l 37.734 78.223 l 36.672 78.223 l 38.797 75.363 l 36.859 72.754 l 37.922
 72.754 l 39.375 74.707 l 40.813 72.754 l h
 41.875 72.754 m f
 43.121 79.676 m 44.168 79.676 l 44.168 76.066 l 43.027 76.301 l 43.027 
 75.707 l 44.152 75.488 l 44.793 75.488 l 44.793 79.676 l 45.84 79.676 l 
 45.84 80.223 l 43.121 80.223 l h
 43.121 79.676 m f
 47.621 76.988 m 48.652 76.988 l 48.652 77.816 l 47.855 79.379 l 47.215 
 79.379 l 47.621 77.816 l h
 47.621 76.988 m f
 55.117 72.754 m 53.148 75.41 l 55.227 78.223 l 54.164 78.223 l 52.57 76.066
 l 50.977 78.223 l 49.914 78.223 l 52.039 75.363 l 50.102 72.754 l 51.164
 72.754 l 52.617 74.707 l 54.055 72.754 l h
 55.117 72.754 m f
 56.801 79.676 m 59.035 79.676 l 59.035 80.223 l 56.02 80.223 l 56.02 79.676
 l 56.27 79.426 56.602 79.094 57.02 78.676 c 57.434 78.25 57.699 77.973 
 57.816 77.848 c 58.023 77.621 58.164 77.426 58.238 77.27 c 58.32 77.105 
 58.363 76.941 58.363 76.785 c 58.363 76.535 58.273 76.332 58.098 76.176 
 c 57.918 76.02 57.691 75.941 57.41 75.941 c 57.211 75.941 56.996 75.98 56.77
 76.051 c 56.551 76.113 56.316 76.219 56.066 76.363 c 56.066 75.707 l 56.316
 75.605 56.551 75.527 56.77 75.473 c 56.996 75.422 57.207 75.395 57.395 
 75.395 c 57.883 75.395 58.273 75.52 58.566 75.77 c 58.855 76.012 59.004 
 76.34 59.004 76.754 c 59.004 76.941 58.965 77.125 58.895 77.301 c 58.82 
 77.48 58.691 77.688 58.504 77.926 c 58.449 77.988 58.277 78.168 57.988 78.457
 c 57.707 78.75 57.309 79.156 56.801 79.676 c h
 56.801 79.676 m f
 60.484 70.629 m 61.266 70.629 l 61.754 71.402 62.117 72.156 62.359 72.895
 c 62.609 73.625 62.734 74.355 62.734 75.082 c 62.734 75.813 62.609 76.547
 62.359 77.285 c 62.117 78.027 61.754 78.773 61.266 79.535 c 60.484 79.535
 l 60.922 78.793 61.242 78.059 61.453 77.316 c 61.672 76.578 61.781 75.832
 61.781 75.082 c 61.781 74.332 61.672 73.594 61.453 72.863 c 61.242 72.125
 60.922 71.379 60.484 70.629 c h
 60.484 70.629 m f
 Q Q
 showpage
 %%Trailer
 end
 %%EOF
--- a/doc/user-manual/fpakc_UserManual.pdf
+++ b/doc/user-manual/fpakc_UserManual.pdf
--- a/doc/user-manual/fpakc_UserManual.tex
+++ b/doc/user-manual/fpakc_UserManual.tex
@ -1,5 +1,5 @@
 \documentclass[10pt,a4paper,twoside]{book}
-\usepackage[latin1]{inputenc}
+%\usepackage[latin1]{inputenc}
 \usepackage{amsmath}
 \usepackage{amsfonts}
 \usepackage{amssymb}
@ -72,7 +72,7 @@
      The \Gls{fpakc} is a simulation tool that models species in plasma (ions, electrons and neutrals) following the trajectories of macro-particles as they move and interact between them and the boundaries of the domain.
      Particles properties are scattered into a finite element mesh in 1, 2 or three dimensions, with the possibility to choose different geometries.
      The official repository can be found at: \url{https://gitlab.com/JorgeGonz/fpakc.git}.
-      The code is currently in very early steps of development and further improvements are expected very soon.
+      The code is currently in the very early steps of development and further refinements are expected very soon.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \section{Main Guidelines}
@ -86,11 +86,11 @@
        \item \acrshort{fpakc} is coded in a \textit{understandable} way.
              This means that the code is required to be written in a clear way that is easy to understand and maintain.
              Variables and procedure names need to be self-understanding.
-              This ease the process of fixing bugs and improving the codes by a large team of developers.
+              This eases the process of fixing bugs and improving the codes by a large team of developers.
              For more information, please refer to the \acrshort{fpakc} Coding Style document.
-        \item \acrshort{fpakc} requires to be ease to use.
+        \item \acrshort{fpakc} requires being ease to use.
              Input files are required to be in a \textit{human} format, meaning that the different options can be easily understood without constant reference to the user guide.
-              \acrshort{fpakc} is aimed to be used in a wide range of applications and by a variety of scientist: from very established ones to newcomers to the field and also students.
+              \acrshort{fpakc} is aimed to be used in a wide range of applications and by various scientists: from well-established ones to newcomers to the field and also students.
      \end{enumerate}
      These are foundation stones of \acrshort{fpakc} and its development and should always be followed, at least for the releases in the official repository.
@ -105,16 +105,16 @@
    \section{The Particle Method}
      \Gls{fpakc} uses macro-particles to simulate the dynamics of different plasma species (mainly ions, electrons and neutrals).
      These macro-particles could represent a large amount of real particles.
-      For now own, macro-particles will be referred as just particles by abusing of language.
+      For now own, macro-particles will be referred as just particles by abuse of language.
      During the initiation phase, the input and mesh file(s) are reading.
      If an initial distribution for a species is specified in the input file, particles to match that distribution are loaded into the cells.
      The general steps performed in each iteration are:
      \begin{enumerate}
        \item Firstly, new particles are introduced into the domain as specified in the input file.
-        \item Particles are then pushed accounting for possible acceleration by external forces.
+        \item Particles are then pushed, accounting for possible acceleration by external forces.
-              During this process, if a particle changes cell it is found using the connectivity between elements.
+              During this process, if a particle changes cell, it is found using the connectivity between elements.
-              If a particle encounters a boundary instead a new cell, the interaction between the boundary and the wall is computed.
+              If a particle encounters a boundary instead a new cell, the interaction between the boundary and the wall are computed.
              A particle may abandon the computational domain and is no longer accounted for.
        \item Next, collisions for the particles inside each cell are carried out.
              This may include different collision processes for each particle.
@ -124,10 +124,10 @@
        \item Finally, particle properties are scattered among the mesh nodes.
              These properties are density, momentum and the stress tensor.
        \item If requested, the electromagnetic field is computed.
-        \item If the number of iteration requires writing output files, it is done after all steps for the particles is completed.
+        \item If the number of iteration requires writing output files, it is done after all steps for the particles are completed.
      \end{enumerate}
-      \Gls{fpakc} has the capability to configure all the behavior of the simulation via the input file.
+      \Gls{fpakc} has the capability to configure all the behaviour of the simulation via the input file.
      Parameters as injection, the kind of pusher used for each species, boundary conditions or collisions are user-input parameters and will be described in Chap.~\ref{ch:input_file}.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -168,8 +168,8 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \section{Find new cell}      
      Once the position and velocity of the particle are updated, the new cell that contains the particle is searched.
-      This is done by a neighbor search, starting from the previous cell containing the particle.
+      This is done by a neighbour search, starting from the previous cell containing the particle.
-      In the process of finding the new cell, it is possible that a particle encounters a boundary.
+      In the process of finding the new cell, a particle might encounter a boundary.
      When the particle interacts with the boundary, the particle may continue its life in the simulation or might be eliminated from it.
      Once that the new cell is found or that the particle life has been terminated, the pushing is complete.
      If a secondary mesh is used for the Monte-Carlo Collision method, the new cell in that mesh in which the particle reside is also found by the same method, although no interaction with the boundaries is accounted for this step.
@ -178,7 +178,7 @@
    \section{Variable Weighting Scheme\label{sec:weightingScheme}}
      One of the issues in particle simulations, specially for axial-symmetrical cases, is that due to the disparate volume of cells, specially close to the axis, the statistics in some cells is usually poor.
      To try to fix that, the possibility to include a Variable Weighting Scheme in the simulations is available in \Gls{fpakc}.
-      These schemes detect when a particle change cells and split it if necessary to improve statistics.
+      These schemes detect when a particle changes cells and split it if necessary to improve statistics.
      The use of a Variable Weighting Scheme is defined by the user in the input file.
      Beware that this can increase the number of particles in the simulation and increase computational time.
@ -189,16 +189,16 @@
      \Gls{fpakc} distinguish between two types of interactions: \acrfull{mcc} and \acrfull{cs}.
      \acrshort{mcc} refers to the process in which two particles interact in short range.
      These processes include, but are not limited to: elastic collisions, ionization/recombination, charge-exchange, excitation/de-excitation\ldots
-      A secondary mesh, with cell sizes in the range of the mean-free path, can be used for this type of collisions.
+      A secondary mesh, with cell sizes in the range of the mean-free path, can be used for this type of collision.
      \acrshort{cs} refers to the large range interaction that a charged species suffer do to the charge of other particles.
      The interactions between the different species is defined by the user.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      \subsection{\acrlong{mcc}}
-        For each cell the maximum number of collisions between particle is computed.
+        For each cell, the maximum number of collisions between particle is computed.
        For each collision, a random pair of particles is chosen.
        A loop over all possible collisions for the pair of particles chosen is performed.
-        If a random number is above the probability of collision for that specific type, the collision take place.
+        If a random number is above the probability of collision for that specific type, the collision takes place.
        If not, the next type for the particle pair is checked.
        Below are described the type of collision process implemented in \acrshort{fpakc}:
@ -219,7 +219,7 @@
          \item Recombination.
                When an electron and an ion interact, there is a possibility for them to be recombined into a neutral particle.
-                The photon emitted by this process is not modelled yet.
+                The photons emitted by this process are not modelled yet.
        \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -238,7 +238,7 @@
      Once that the pushing is complete, the array of particles that remain inside the domain is copied to a new array.
      The new array containing only the particles inside the domain will be the one used in the next steps.
      In this section, particles are assigned to the list of particles inside each individual cell.
-      Unfortunately, this is done right now without parallelisation and is very CPU consuming.
+      Unfortunately, this is done right now without parallelization and is very CPU consuming.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \section{Probing}\label{sec:probing}
@ -250,21 +250,21 @@
      The user can decide the grid width and the number of points in each direction.
      The distribution function will be calculated and wrote with a time step decided by the user.
-      If a particle velocity resides outside of the velocity grid (in any direction), it wont be added to the tally of the distribution function.
+      If a particle velocity resides outside the velocity grid (in any direction), it will not be added to the tally of the distribution function.
-      Due to the limitation of only taking into account particles in the cell, and not neighbour particles, two probes for the same species at different positions but in the same cell will output the same results.
+      Due to the limitation of only considering particles in the cell, and not neighbour particles, two probes for the same species at different positions but in the same cell will output the same results.
-      A more advance method taking into account distance between the particles and the probe position as well as particles in neighbour cells could be implemented to improve the statistics of the distribution function.
+      A more advance method considering distance between the particles and the probe position as well as particles in neighbour cells could be implemented to improve the statistics of the distribution function.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \section{Scattering}
      The properties of each particle are deposited in the nodes from the containing cell.
-      This process depend on the cell type, but in general, each node receive a proportional part of the particle properties as a function of the particle position inside the cell.
+      This process depends on the cell type, but in general, each node receives a proportional part of the particle properties as a function of the particle position inside the cell.
-      Figure \ref{fig:scatteringQuad} shows how a particle at a generic position $p(x_1, x_2)$ inside the cell is scattered to the four nodes.
+      The figure \ref{fig:scatteringQuad} shows how a particle at a generic position $p(x_1, x_2)$ inside the cell is scattered to the four nodes.
      \begin{wrapfigure}{l}{0.4\textwidth}
        \centering
        \includegraphics{figures/scatteringQuad}
        \caption{\label{fig:scatteringQuad}Example of how a particle is weighted in a quadrilateral cell.}
      \end{wrapfigure}
-      Each node receives a proportional part of the area formed by dividing the cell in for rectangles using as an additional vertex the particle position.
+      Each node receives a proportional part of the area formed by dividing the cell in for rectangles, using as an additional vertex the particle position.
      These properties are dimensionless, but they are converted to the correct units once the output is printed.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -273,11 +273,11 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \section{Average scheme}
-      Particle-in-cell codes has an intrinsic statistical noise associated with them.
+      Particle-in-cell codes have an intrinsic statistical noise associated with them.
      Although this can be reduced by increasing the number of particles, this also increases the CPU requirements of the case.
      It is quite common that most cases reach a quasi-steady state after a number of iterations and time-average results can be obtained after to improve analysis, plotting and restarting the case using these time-average results as new species backgrounds.
-      Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the amount of iterations printed.
+      Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the number of iterations printed.
      \Gls{fpakc} implements a simple average scheme that, after a start time provided by the user, scores a mean and standard deviation of all the main species properties, and the electromagnetic field.
      This scheme is based on the Welford's online algorithm~\cite{welford1962note}.
@ -286,7 +286,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
  \chapter{Installation}
    \section{Required Packages}
-      In order to properly compile \gls{fpakc}, the following packages are required.
+      To properly compile \gls{fpakc}, the following packages are required.
      \subsection{Gfortran}
        The \Gls{opensource} free compiler \Gls{gfortran}\cite{gfortranURL} from GCC is the basic way to compile \acrshort{fpakc}.
        It is distributed with all GNU/Linux distributions.
@ -369,7 +369,7 @@ make
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \section{Case file}
      The required format for the case file is \Gls{json}.
-      \Gls{json} is a case-sensitive format, so input must be written with the correct capitalisation.
+      \Gls{json} is a case-sensitive format, so input must be written with the correct capitalization.
      The basic structure and options available for the case file are explained below.
      The order of the objects and variables is irrelevant, but the structure needs to be maintained.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -380,9 +380,9 @@ make
          \item \textbf{path}: Character.
                               Path for the output files. This path is also used to locate the mesh input file.
          \item \textbf{folder}: Character.
-                                 Base name of the folder in wich output files are placed.
+                                 Base name of the folder in which output files are placed.
                                 The date and time is appended to this name.
-                                 If none is provided, only the date and time is writted as the folder name.
+                                 If none is provided, only the date and time is written as the folder name.
          \item \textbf{triggerOutput}: Integer.
                                        Determines the number of iterations between writing output files for macroscopic quantities.
          \item \textbf{cpuTime}: Logical. 
@ -460,6 +460,10 @@ make
                                   \begin{itemize}
                                     \item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0. This has to be in ASCII format.
                                     \item \textbf{vtu}: \Gls{vtu} file format. This has to be in ASCII format.
                                     \item \textbf{text}: Plain text file format only intended for 1D cases.
                                                          This has to be in ASCII format and comma separated.
                                                          The first column represents the position and the second column the physical ID of the node.
                                                          Values have to be $1$ (left boundary), $2$ (right boundary), or $0$ (no boundary.)
                                   \end{itemize}
          \item \textbf{meshFile}: Character.
                                   Mesh filename.
@ -515,7 +519,7 @@ make
                                  \item \textbf{absorption}: Particle is eliminated from the domain.
                                                             The particle is first moved into the edge and its properties are scattered among the edge nodes.
                                  \item \textbf{transparent}: Particle abandon the numerical domain.
-                                  \item \textbf{wallTemperature}: Reflective wall with cosntant temperature that exchange heat with particles.
+                                  \item \textbf{wallTemperature}: Reflective wall with constant temperature that exchange heat with particles.
                                                                  Required parameters are:
                                                                  \begin{itemize}
                                                                    \item \textbf{temperature}: Real.
@ -526,8 +530,8 @@ make
                                                                                                 Specific heat capacity of the material.
                                                                  \end{itemize}
                                  \item \textbf{ionization}: Per each particle crossing the surface with this type of boundary, a number of ionization events are calculated.
-                                                             A pair of ion-electron is generated for each ionization event taking as a reference a neutral background.
+                                                             A pair of ion-electron is generated for each ionization event, taking as a reference a neutral background.
-                                                             Secondary electron is taken as same type as incident particle.
+                                                             The secondary electron is taken as the same type as the incident particle.
                                                             The available input is:
                                                             \begin{itemize}
                                                               \item \textbf{neutral}: Object.
@ -540,7 +544,7 @@ make
                                                                                         \item \textbf{mass}: Real.
                                                                                                              Units in $\unit{kg}$.
                                                                                                              Mass of neutral species.
-                                                                                                              If missing, the mass of the ion is ussed
+                                                                                                              If missing, the mass of the ion is used
                                                                                         \item \textbf{density}: Real.
                                                                                                                 Units in $\unit{m^{-3}}$.
                                                                                                                 Density of neutral background.
@ -558,18 +562,18 @@ make
                                                                                       \end{itemize}
                                                               \item \textbf{effectiveTime}: Real.
                                                                                             Units in $\unit{s}$.
-                                                                                             As the particle is no longer simulated once it crossed the boundary, this time represent the effective time in which the particle produces ionization processes in the neutral background.
+                                                                                             As the particle is no longer simulated once it crossed the boundary, this time represents the effective time in which the particle produces ionization processes in the neutral background.
                                                                                             Required parameter.
                                                               \item \textbf{energyThreashold}: Real.
                                                                                                Units in $\unit{eV}$.
                                                                                                Ionization energy threshold for the simulated process.
                                                                                                Required parameter.
                                                               \item \textbf{crossSection}: Character.
-                                                                                            Complete path to the cross section data for the ionization process.
+                                                                                            Complete path to the cross-section data for the ionization process.
                                                             \end{itemize}
                                  \item \textbf{axis}: Identifies the symmetry axis for 2D cylindrical simulations.
-                                                       If for some reason a particle interact with this axis, it is reflected.
+                                                       If , for some reason, a particle interacts with this axis, it is reflected.
                                \end{itemize}
        \end{itemize}
@ -585,18 +589,26 @@ make
                               Type of boundary.
                                Accepted values are:
                               \begin{itemize}
-                                 \item \textbf{dirichlet}: Elastic reflection of particles.
+                                 \item \textbf{dirichlet}: Constant value of electric potential on the surface.
                                 \item \textbf{dirichletTime}: Constant value of the electric potential with a time variable profile.
                                                               The value of \textbf{boundaryEM.potential} will be multiplied for the corresponding value in the file \textbf{boundaryEM.temporalProfile}.
                                \end{itemize}
-           \item \textbf{potential}: Real.
+          \item \textbf{potential}: Real.
-                                     Fixed potential for Dirichlet boundary condition.
+                                    Fixed potential for Dirichlet boundary condition.
          \item \textbf{physicalSurface}: Integer.
                                          Identification of the edge in the mesh file.
          \item \textbf{temporalProfile}: Character.
                                          Filename of the 2 column file containing the time variable profile.
                                          File must be located in \textbf{output.path}.
                                          The first column is the time in $\unit{s}$.
                                          The second column is the factor that will multiply the value of the boundary.
        \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      \subsection{inject}
        The array \textbf{inject} specifies the injection of particles from different surfaces.
-        The injection of particles need to be associated to a physicalSurface in the mesh file.
+        The injection of particles needs to be associated to a physicalSurface in the mesh file.
        Multiple injections can be associated to the same surface.
        \begin{itemize}
          \item \textbf{name}: Character.
@ -610,7 +622,9 @@ make
                                Available values are:
                                \begin{itemize}
                                  \item \textbf{A}: Ampere.
-                                  \item \textbf{sccm}: Standard cubic centimeter.
+                                  \item \textbf{Am2}: Ampere per square meter.
                                                      This value will be multiplied by the area of injection.
                                  \item \textbf{sccm}: Standard cubic centimetre.
                                  \item \textbf{part/s}: Particles (real) per second.
                                \end{itemize}
          \item \textbf{v}: Real.
@ -627,7 +641,7 @@ make
                                  \begin{itemize}
                                    \item \textbf{Maxwellian}: Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction.
                                    \item \textbf{Half-Maxwellian}: Half-Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction.
-                                                                    Only takes into account the positive part of the half-Maxwellian.
+                                                                    Only considers the positive part of the half-Maxwellian.
                                    \item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction.
                                  \end{itemize}
          \item \textbf{T}: Real.
@ -636,6 +650,11 @@ make
                            Temperature in each direction.
          \item \textbf{physicalSurface}: Integer.
                                          Identification of the edge in the mesh file.
          \item \textbf{particlesPerEdge}: Integer.
                                           Optional.
                                           Number of particles to be injected by each edge in the numerical domain.
                                           The weight of the particles for each edge will modified by the surface of the edge to ensure the right flux is injected.
                                           If no value is provided, the number of particles to inject per edge will be calculated with the species weight and the surface of the edge respect to the total one.
        \end{itemize}        
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      \subsection{reference}
@ -651,7 +670,7 @@ make
          \item \textbf{radius}: Real.
                                 Reference atomic radius in $\unit{m}$.
          \item \textbf{crossSection}: Real.
-                                       Reference cross section in $\unit{m^2}$.
+                                       Reference cross-section in $\unit{m^2}$.
                                       If this value is present, radius is ignored.
        \end{itemize}
@ -677,8 +696,8 @@ make
                                 Indicates the type of pusher used for each species:
                                 \begin{itemize}
                                   \item \textbf{Neutral}: Pushes a particle without any external force.
-                                   \item \textbf{Electrostatic}: Pushes a particle including the effect of the electrostatic field.
+                                   \item \textbf{Electrostatic}: Pushes a particle, including the effect of the electrostatic field.
-                                   \item \textbf{Electromagnetic}: Pushes particles accounting for the electromagnetic field.
+                                   \item \textbf{Electromagnetic}: Pushes a particle, accounting for the electromagnetic field.
                                 \end{itemize}
          \item \textbf{WeightingScheme}: Character.
                                          Indicates the variable weighting scheme to be used in the simulation.
@ -706,11 +725,15 @@ make
                                  \begin{itemize}
                                    \item \textbf{species}: Character.
                                                            Name of species as defined in the object \textbf{species}.
-                                    \item \textbf{file}:    Character.
+                                    \item \textbf{file}: Character.
-                                                            Output file from previous run used as an initial state for the species.
+                                                         Output file from previous run used as an initial state for the species.
-                                                            The file format must be the same as in \textbf{geometry.meshType}
+                                                         The file format must be the same as in \textbf{geometry.meshType}
-                                                            Initial particles are assumed to have a Maxwellian distribution.
+                                                         Initial particles are assumed to have a Maxwellian distribution.
-                                                            File must be located at \textbf{output.path}.
+                                                         File must be located in \textbf{output.path}.
                                    \item \textbf{particlesPerCell}: Integer.
                                                                     Optional.
                                                                     Initial number of particles per cell.
                                                                     If not, the number of particles per cell will be assigned based on the species weight and the cell volume.
                                  \end{itemize}
        \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -726,11 +749,11 @@ make
        \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      \subsection{interactions}\label{ssec:input_interactions}
-        This object determine the different interactions among species.
+        This object determines the different interactions among species.
        Acceptable values are:
        \begin{itemize}
          \item \textbf{folderCollisions}: Character.
-                                           Indicates the path to in which the cross section tables are allocated.
+                                           Indicates the path to in which the cross-section tables are allocated.
          \item \textbf{meshCollisions}: Character.
                                         Determines a specific mesh for \acrshort{mcc} processes.
                                         The file needs to be located in the folder \textbf{output.folder}.
@ -757,13 +780,18 @@ make
                                                                                     Accepted values are \textbf{elastic}, \textbf{chargeExchange}, \textbf{ionization} and \textbf{recombination}.
                                                                                     Please refer to Sec.~\ref{ssec:collisions} for a description of the different collision types.
                                                                \item \textbf{crossSection}: Character.
-                                                                                             File in \textbf{interactions.folderCollisions} that contains the cross section data as a 1D table of relative energy (in $\unit{eV}$) and cross section (in $\unit{m^-2}$).
+                                                                                             File in \textbf{interactions.folderCollisions} that contains the cross-section data as a 1D table of relative energy (in $\unit{eV}$) and cross-section (in $\unit{m^-2}$).
                                                                \item \textbf{energyThreshold}: Real.
                                                                                                Energy threshold of the collisional process in $\unit{eV}$.
                                                                                                Only valid for \textbf{ionization} and \textbf{recombination} processes.
                                                                \item \textbf{electron}: Character.
                                                                                         Name of species designed as electrons.
                                                                                         Only valid for \textbf{ionization} and \textbf{recombination} processes.
                                                                \item \textbf{electronSecondary}: Character.
                                                                                                  Optional.
                                                                                                  Name of species designed as secondary electrons.
                                                                                                  If none provided, \textbf{electron} is used.
                                                                                                  Only valid for \textbf{ionization}.
                                                              \end{itemize}
                                     \end{itemize}
          \item \textbf{Coulomb}: Array of objects.
@ -773,7 +801,7 @@ make
                                  \begin{itemize}
                                    \item \textbf{species\_i}, \textbf{species\_j}: Character.
                                                                                    Define the two species involved in the collision processes.
-                                                                                    Order is indiferent.
+                                                                                    Order is indifferent.
                                  \end{itemize}
        \end{itemize}
@ -799,9 +827,9 @@ make
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    \section{1D Emissive Cathode (1D\_Cathode)}
-      Emission from a 1D cathode in both, cartesian and radial coordinates.
+      Emission from a 1D cathode in both, Cartesian and radial coordinates.
-      Both cases insert the same amount of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field.
+      Both cases insert the same number of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field.
-      This case is useful to ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifying some parameters in the input file.
+      This case is useful to illustrate how \acrshort{fpakc} can deal with different geometries by just modifying some parameters in the input file.
      The same mesh file (\lstinline|mesh.msh|) is used for both cases.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- a/src/fpakc.f90
+++ b/src/fpakc.f90
@ -11,12 +11,12 @@ PROGRAM fpakc
  USE OMP_LIB
  IMPLICIT NONE
  ! t = time step
  INTEGER:: t
  ! arg1 = Input argument 1 (input file)
  CHARACTER(200):: arg1
  ! inputFile = path+name of input file
  CHARACTER(:), ALLOCATABLE:: inputFile
  ! generic integer for time step
  INTEGER:: t
  tStep = omp_get_wtime()
  !Gets the input file
@ -27,11 +27,18 @@ PROGRAM fpakc
  !Reads the json configuration file
  CALL readConfig(inputFile)
  !Create output folder and initial files
  CALL initOutput(inputFile)
  !Do '0' iteration
-  t = tInitial
+  timeStep = tInitial
  !$OMP PARALLEL DEFAULT(SHARED)
  !$OMP SINGLE
  ! Initial reset of probes
  CALL resetProbes()
  CALL verboseError("Initial scatter of particles...")
  !$OMP END SINGLE
  CALL doScatter()
@ -45,19 +52,21 @@ PROGRAM fpakc
  tStep = omp_get_wtime() - tStep
  !Output initial state
-  CALL doOutput(t)
+  CALL doOutput()
  CALL verboseError('Starting main loop...')
  !$OMP PARALLEL DEFAULT(SHARED) 
  DO t = tInitial + 1, tFinal
    !Insert new particles and push them
    !$OMP SINGLE
    tStep = omp_get_wtime()
    ! Update global time step index
    timeStep = t
    !Checks if a species needs to me moved in this iteration
-    CALL solver%updatePushSpecies(t)
+    CALL solver%updatePushSpecies()
    !Checks if probes need to be calculated this iteration
-    CALL resetProbes(t)
+    CALL resetProbes()
    tPush = omp_get_wtime()
    !$OMP END SINGLE
@ -75,7 +84,7 @@ PROGRAM fpakc
    !$OMP END SINGLE
    IF (doMCCollisions) THEN
-      CALL meshForMCC%doCollisions(t)
+      CALL meshForMCC%doCollisions()
    END IF
@ -120,12 +129,12 @@ PROGRAM fpakc
    !$OMP SINGLE
    tEMField = omp_get_wtime() - tEMField
-    CALL doAverage(t)
+    CALL doAverage()
    tStep = omp_get_wtime() - tStep
    !Output data
-    CALL doOutput(t)
+    CALL doOutput()
    !$OMP END SINGLE
  END DO
--- a/src/makefile
+++ b/src/makefile
@ -9,6 +9,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \
 			  	$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
 			  	$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
 				$(OBJDIR)/moduleMeshInputText.o $(OBJDIR)/moduleMeshOutputText.o \
 			  	$(OBJDIR)/moduleMesh3DCart.o \
 			  	$(OBJDIR)/moduleMesh2DCyl.o  \
 			  	$(OBJDIR)/moduleMesh2DCart.o \
--- a/src/modules/common/moduleCaseParam.f90
+++ b/src/modules/common/moduleCaseParam.f90
@ -2,8 +2,13 @@
 MODULE moduleCaseParam
  !Final and initial iterations
  INTEGER:: tFinal, tInitial = 0
  ! Global index of current iteration
  INTEGER:: timeStep
  ! Time step for all species
  REAL(8), ALLOCATABLE:: tau(:)
  ! Minimum time step
  REAL(8):: tauMin
  ! Time step for Monte-Carlo Collisions
  REAL(8):: tauColl
 END MODULE moduleCaseParam
--- a/src/modules/common/moduleRandom.f90
+++ b/src/modules/common/moduleRandom.f90
@ -40,31 +40,48 @@ MODULE moduleRandom
      INTEGER:: rnd
      REAL(8):: rnd01
-      rnd = 0.D0
+      rnd = 0
      CALL RANDOM_NUMBER(rnd01)
-      rnd = INT(REAL(b - a) * rnd01) + 1
+      rnd = a + FLOOR((b+1-a)*rnd01)
    END FUNCTION randomIntAB
    !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method
    function randomMaxwellian() result(rnd)
      USE moduleConstParam, only: pi
      implicit none
      real(8):: rnd
      real(8):: v1, v2, Rsquare
      v1 = 0.d0
      do while (v1 <= 0.d0)
        v1 = random()
      end do
      v2 = random()
      rnd = sqrt(-2.d0*log(v1))*cos(2*pi*v2)
    end function randomMaxwellian
    !Returns a random number in a Maxwellian distribution of mean 0 and width 1
-    FUNCTION randomMaxwellian() RESULT(rnd)
+    FUNCTION randomHalfMaxwellian() RESULT(rnd)
      USE moduleConstParam, ONLY: PI
      IMPLICIT NONE
      REAL(8):: rnd
-      REAL(8):: x, y
+      REAL(8):: x
      rnd = 0.D0
      x = 0.D0
      DO WHILE (x == 0.D0)
        CALL RANDOM_NUMBER(x)
      END DO
      CALL RANDOM_NUMBER(y)
-      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
+      rnd = DSQRT(-DLOG(x))
-    END FUNCTION randomMaxwellian
+    END FUNCTION randomHalfMaxwellian
    !Returns a random number weighted with the cumWeight array
    FUNCTION randomWeighted(cumWeight, sumWeight) RESULT(rnd)
@ -73,10 +90,21 @@ MODULE moduleRandom
      REAL(8), INTENT(in):: cumWeight(1:)
      REAL(8), INTENT(in):: sumWeight
      REAL(8):: rnd0b
-      INTEGER:: rnd
+      INTEGER:: rnd, i
-      rnd0b = random(0.D0, sumWeight)
+      rnd0b = random()
-      rnd = MINLOC(DABS(rnd0b - cumWeight), 1)
+      i = 1
      DO 
        IF (rnd0b <= cumWeight(i)/sumWeight) THEN
          rnd = i
          EXIT
        ELSE
          i = i +1
        END IF
      END DO
      ! rnd = MINLOC(DABS(rnd0b - cumWeight), 1)
    END FUNCTION randomWeighted
--- a/src/modules/common/moduleTable.f90
+++ b/src/modules/common/moduleTable.f90
@ -93,7 +93,7 @@ MODULE moduleTable
        f = self%fMax
      ELSE
-        i = MINLOC(x - self%x, 1)
+        i = MINLOC(ABS(x - self%x), 1)
        deltaX = x - self%x(i)
        IF (deltaX < 0 ) THEN 
          i = i - 1
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@ -84,20 +84,6 @@ MODULE moduleInput
      CALL readParallel(config)
      CALL checkStatus(config, "readParallel")
      !If everything is correct, creates the output folder
      CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
      !Copies input file to output folder
      CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
      !Copies particle mesh
      IF (mesh%dimen > 0) THEN
        CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
        IF (doubleMesh) THEN
          CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
        END IF
      END IF
    END SUBROUTINE readConfig
    !Checks the status of the JSON case file and, if failed, exits the execution.
@ -273,13 +259,17 @@ MODULE moduleInput
          !Read BC
          CALL readEMBoundary(config)
        CASE("ElectrostaticBoltzmann")
          !Read BC
          CALL readEMBoundary(config)
        CASE("ConstantB")
          !Read BC
          CALL readEMBoundary(config)
          !Read constant magnetic field
          DO i = 1, 3
-            WRITE(istring, '(i2)') i
+            WRITE(iString, '(i2)') i
-            CALL config%get(object // '.B(' // istring // ')', B(i), found)
+            CALL config%get(object // '.B(' // iString // ')', B(i), found)
            IF (.NOT. found) THEN
              CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
@ -322,7 +312,7 @@ MODULE moduleInput
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: object
      INTEGER:: nInitial
-      INTEGER:: i, j, p, e
+      INTEGER:: i, p, e
      CHARACTER(LEN=2):: iString
      CHARACTER(:), ALLOCATABLE:: spName
      INTEGER:: sp
@ -338,7 +328,8 @@ MODULE moduleInput
      REAL(8):: densityCen
      !Mean velocity and temperature at particle position
      REAL(8):: velocityXi(1:3), temperatureXi
-      INTEGER:: nNewPart = 0.D0
+      INTEGER:: nNewPart = 0
      REAL(8):: weight = 0.D0
      CLASS(meshCell), POINTER:: cell
      TYPE(particle), POINTER:: partNew
      REAL(8):: vTh
@ -357,6 +348,9 @@ MODULE moduleInput
          !Reads node values at the nodes
          filename = path // spFile
          CALL mesh%readInitial(filename, density, velocity, temperature)
          !Check if initial number of particles is given
          CALL config%get(object // '.particlesPerCell', nNewPart, found)
          !For each volume in the node, create corresponding particles
          DO e = 1, mesh%numCells
            !Scale variables
@ -369,7 +363,11 @@ MODULE moduleInput
            densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar)
            !Calculate number of particles
-            nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
+            IF (.NOT. found) THEN
              nNewPart = FLOOR(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
            END IF
            weight = densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / REAL(nNewPart)
            !Allocate new particles
            DO p = 1, nNewPart
@ -406,7 +404,7 @@ MODULE moduleInput
              partNew%n_in   = .TRUE.
-              partNew%weight = species(sp)%obj%weight
+              partNew%weight = weight
              !Assign particle to temporal list of particles
              CALL partInitial%add(partNew)
@ -634,7 +632,7 @@ MODULE moduleInput
      INTEGER:: i, k, ij
      INTEGER:: pt_i, pt_j
      REAL(8):: energyThreshold, energyBinding
-      CHARACTER(:), ALLOCATABLE:: electron
+      CHARACTER(:), ALLOCATABLE:: electron, electronSecondary
      INTEGER:: e
      CLASS(meshCell), POINTER:: cell
@ -711,8 +709,16 @@ MODULE moduleInput
                IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
                CALL config%get(object // '.electron', electron, found)
                IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
-                CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
+                CALL config%get(object // '.electronSecondary', electronSecondary, found)
-                                          crossSecFilePath, energyThreshold, electron)
+                IF (found) THEN
                  CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
                                            crossSecFilePath, energyThreshold, electron, electronSecondary)
                ELSE
                  CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
                                            crossSecFilePath, energyThreshold, electron)
                END IF
              CASE ('recombination')
                !Electorn impact ionization
@ -797,7 +803,7 @@ MODULE moduleInput
      TYPE(json_file), INTENT(inout):: config
      INTEGER:: i, s
-      CHARACTER(2):: istring, sString
+      CHARACTER(2):: iString, sString
      CHARACTER(:), ALLOCATABLE:: object, bType
      REAL(8):: Tw, cw !Wall temperature and specific heat
      !Neutral Properties
@ -805,16 +811,16 @@ MODULE moduleInput
      REAL(8), DIMENSION(:), ALLOCATABLE:: v0
      REAL(8):: effTime
      REAL(8):: eThreshold !Energy threshold
-      INTEGER:: speciesID
+      INTEGER:: speciesID, electronSecondaryID
-      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection
+      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, electronSecondary
      LOGICAL:: found
      INTEGER:: nTypes
      CALL config%info('boundary', found, n_children = nBoundary)
      ALLOCATE(boundary(1:nBoundary))
      DO i = 1, nBoundary
-        WRITE(istring, '(i2)') i
+        WRITE(iString, '(i2)') i
-        object = 'boundary(' // TRIM(istring) // ')'
+        object = 'boundary(' // TRIM(iString) // ')'
        boundary(i)%n = i
        CALL config%get(object // '.name',            boundary(i)%name,            found)
@ -823,62 +829,77 @@ MODULE moduleInput
        IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary')
        ALLOCATE(boundary(i)%bTypes(1:nSpecies))
        DO s = 1, nSpecies
-          WRITE(sString,'(i2)') s
+          associate(bound => boundary(i)%bTypes(s)%obj)
-          object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
+            WRITE(sString,'(i2)') s
-          CALL config%get(object // '.type', bType, found)
+            object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
-          SELECT CASE(bType)
+            CALL config%get(object // '.type', bType, found)
-          CASE('reflection')
+            SELECT CASE(bType)
-            ALLOCATE(boundaryReflection:: boundary(i)%bTypes(s)%obj)
+            CASE('reflection')
              ALLOCATE(boundaryReflection:: bound)
-          CASE('absorption')
+            CASE('absorption')
-            ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj)
+              ALLOCATE(boundaryAbsorption:: bound)
-          CASE('transparent')
+            CASE('transparent')
-            ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj)
+              ALLOCATE(boundaryTransparent:: bound)
-          CASE('ionization')
+            CASE('axis')
-            !Neutral parameters
+              ALLOCATE(boundaryAxis:: bound)
            CALL config%get(object // '.neutral.ion', speciesName, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
            speciesID = speciesName2Index(speciesName)
            CALL config%get(object // '.neutral.mass', m0, found)
            IF (.NOT. found) THEN
              m0 = species(s)%obj%m*m_ref
            END IF
            CALL config%get(object // '.neutral.density', n0, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
            CALL config%get(object // '.neutral.velocity', v0, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
            CALL config%get(object // '.neutral.temperature', T0, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
-            CALL config%get(object // '.effectiveTime', effTime, found)
+            CASE('wallTemperature')
-            IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')
+              CALL config%get(object // '.temperature', Tw, found)
              IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
              CALL config%get(object // '.specificHeat', cw, found)
              IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')
-            CALL config%get(object // '.energyThreshold', eThreshold, found)
+              CALL initWallTemperature(bound, Tw, cw)
            IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
-            CALL config%get(object // '.crossSection', crossSection, found)
+            CASE('ionization')
-            IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
+              !Neutral parameters
              CALL config%get(object // '.neutral.ion', speciesName, found)
              IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
              speciesID = speciesName2Index(speciesName)
              CALL config%get(object // '.neutral.mass', m0, found)
              IF (.NOT. found) THEN
                m0 = species(s)%obj%m*m_ref
              END IF
              CALL config%get(object // '.neutral.density', n0, found)
              IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
              CALL config%get(object // '.neutral.velocity', v0, found)
              IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
              CALL config%get(object // '.neutral.temperature', T0, found)
              IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')
-            CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
+              CALL config%get(object // '.effectiveTime', effTime, found)
-                                                           speciesID, effTime, crossSection, eThreshold)
+              IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')
-          CASE('wallTemperature')
+              CALL config%get(object // '.energyThreshold', eThreshold, found)
-            CALL config%get(object // '.temperature', Tw, found)
+              IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')
            IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
            CALL config%get(object // '.specificHeat', cw, found)
            IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')
-            CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw)
+              CALL config%get(object // '.crossSection', crossSection, found)
              IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
-          CASE('axis')
+              CALL config%get(object // '.electronSecondary', electronSecondary, found)
-            ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)
+              electronSecondaryID = speciesName2Index(electronSecondary)
              IF (found) THEN
                CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, &
                                                               speciesID, effTime, crossSection, eThreshold,electronSecondaryID)
-          CASE DEFAULT
+              ELSE
-            CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
+                CALL initIonization(bound, species(s)%obj%m, m0, n0, v0, T0, &
                                                               speciesID, effTime, crossSection, eThreshold)
-          END SELECT
+              END IF
            case('outflowAdaptive')
              allocate(boundaryOutflowAdaptive:: bound)
            CASE DEFAULT
              CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
            END SELECT
          end associate
        END DO
@ -895,6 +916,7 @@ MODULE moduleInput
      USE moduleMeshInputGmsh2, ONLY: initGmsh2
      USE moduleMeshInputVTU,   ONLY: initVTU
      USE moduleMeshInput0D,    ONLY: init0D
      USE moduleMeshInputText,  ONLY: initText
      USE moduleMesh3DCart
      USE moduleMesh2DCyl
      USE moduleMesh2DCart
@ -912,7 +934,6 @@ MODULE moduleInput
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
      REAL(8):: volume
      CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT
      object = 'geometry'
@ -954,9 +975,9 @@ MODULE moduleInput
        !Read the 0D mesh
        CALL mesh%readMesh(pathMeshParticle)
-        !Get the volumne
+        !Get the volume
        CALL config%get(object // '.volume', volume, found)
-        !Rescale the volumne
+        !Rescale the volume
        IF (found) THEN
          mesh%cells(1)%obj%volume = mesh%cells(1)%obj%volume*volume / Vol_ref
          mesh%nodes(1)%obj%v     = mesh%cells(1)%obj%volume
@ -1044,6 +1065,20 @@ MODULE moduleInput
          END IF
        case ("text")
          !Check if the geometry is right.
          if (mesh%dimen /= 1) then
            call criticalError("Text mesh is only allowed for 1D geometries", 'readGeometry')
          end if
          !Read the mesh
          call initText(mesh)
          if (doubleMesh) then
            call initText(meshColl)
          end if
        CASE DEFAULT
          CALL criticalError('Mesh format ' // meshFormat // ' not defined.', 'readGeometry')
@ -1090,13 +1125,13 @@ MODULE moduleInput
      TYPE(json_file), INTENT(inout):: config
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found
-      CHARACTER(2):: istring
+      CHARACTER(2):: iString
      INTEGER:: i
      CHARACTER(:), ALLOCATABLE:: speciesName
      REAL(8), ALLOCATABLE, DIMENSION(:):: r
      REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
      INTEGER, ALLOCATABLE, DIMENSION(:):: points
-      REAL(8):: timeStep
+      REAL(8):: everyTimeStep
      CALL config%info('output.probes', found, n_children = nProbes)
@ -1104,7 +1139,7 @@ MODULE moduleInput
      DO i = 1, nProbes
        WRITE(iString, '(I2)') i
-        object = 'output.probes(' // trim(istring) // ')'
+        object = 'output.probes(' // trim(iString) // ')'
        CALL config%get(object // '.species', speciesName, found)
        CALL config%get(object // '.position', r, found)
@ -1112,16 +1147,14 @@ MODULE moduleInput
        CALL config%get(object // '.velocity_2', v2, found)
        CALL config%get(object // '.velocity_3', v3, found)
        CALL config%get(object // '.points', points, found)
-        CALL config%get(object // '.timeStep', timeStep, found)
+        CALL config%get(object // '.timeStep', everyTimeStep, found)
        IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')
-        CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep)
+        CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep)
      END DO
      CALL resetProbes(tInitial)
    END SUBROUTINE readProbes
    SUBROUTINE readEMBoundary(config)
@ -1129,7 +1162,6 @@ MODULE moduleInput
      USE moduleOutput
      USE moduleErrors
      USE moduleEM
      USE moduleRefParam
      USE moduleSpecies
      USE json_module
      IMPLICIT NONE
@ -1137,34 +1169,72 @@ MODULE moduleInput
      TYPE(json_file), INTENT(inout):: config
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found
-      CHARACTER(2):: istring
+      CHARACTER(:), ALLOCATABLE:: typeEM
-      INTEGER:: i, e, s
+      REAL(8):: potential
      INTEGER:: physicalSurface
      CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath
      INTEGER:: b, s, n, ni
      CHARACTER(2):: bString
      INTEGER:: info
      EXTERNAL:: dgetrf
      CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
-      IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
+      IF (found) THEN
        ALLOCATE(boundaryEM(1:nBoundaryEM))
      END IF
-      DO i = 1, nBoundaryEM
+      DO b = 1, nBoundaryEM
-        WRITE(istring, '(I2)') i
+        WRITE(bString, '(I2)') b
-        object = 'boundaryEM(' // trim(istring) // ')'
+        object = 'boundaryEM(' // TRIM(bString) // ')'
-        CALL config%get(object // '.type', boundEM(i)%typeEM, found)
+        CALL config%get(object // '.type', typeEM, found)
-        SELECT CASE(boundEM(i)%typeEM)
+        SELECT CASE(typeEM)
        CASE ("dirichlet")
-          CALL config%get(object // '.potential', boundEM(i)%potential, found)
+          CALL config%get(object // '.potential', potential, found)
-          IF (.NOT. found) &
+          IF (.NOT. found) THEN
            CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
          boundEM(i)%potential = boundEM(i)%potential/Volt_ref
-          CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
+          END IF
-          IF (.NOT. found) &
+
-            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
+          CALL config%get(object // '.physicalSurface', physicalSurface, found)
          IF (.NOT. found) THEN
            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', &
                               'readEMBoundary')
          END IF
          CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential)
        CASE ("dirichletTime")
          CALL config%get(object // '.potential', potential, found)
          IF (.NOT. found) THEN
            CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', &
                               'readEMBoundary')
          END IF
          CALL config%get(object // '.temporalProfile', temporalProfile, found)
          IF (.NOT. found) THEN
            CALL criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', &
                               'readEMBoundary')
          END IF
          temporalProfilePath = path // temporalProfile
          CALL config%get(object // '.physicalSurface', physicalSurface, found)
          IF (.NOT. found) THEN
            CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', &
                               'readEMBoundary')
          END IF
          CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath)
        CASE DEFAULT
-          CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')
+          CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')
        END SELECT
@ -1183,18 +1253,28 @@ MODULE moduleInput
      END DO
-      IF (ALLOCATED(boundEM)) THEN
+      ! Modify K matrix due to boundary conditions
-        DO e = 1, mesh%numEdges
+      DO b = 1, nBoundaryEM
-          IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
+        SELECT TYPE(boundary => boundaryEM(b)%obj)
-            DO i = 1, nBoundaryEM
+        TYPE IS(boundaryEMDirichlet)
-              IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
+          DO n = 1, boundary%nNodes
-                CALL boundEM(i)%apply(mesh%edges(e)%obj)
+            ni = boundary%nodes(n)%obj%n
            mesh%K(ni, :)  = 0.D0
            mesh%K(ni, ni) = 1.D0
-              END IF
+          END DO
-            END DO
+
-          END IF
+        TYPE IS(boundaryEMDirichletTime)
-        END DO
+          DO n = 1, boundary%nNodes
-      END IF
+            ni = boundary%nodes(n)%obj%n
            mesh%K(ni, :)  = 0.D0
            mesh%K(ni, ni) = 1.D0
          END DO
        END SELECT
      END DO
      !Compute the PLU factorization of K once boundary conditions have been read
      CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
@ -1215,24 +1295,25 @@ MODULE moduleInput
      TYPE(json_file), INTENT(inout):: config
      INTEGER:: i
-      CHARACTER(2):: istring
+      CHARACTER(2):: iString
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: speciesName
      CHARACTER(:), ALLOCATABLE:: name
      REAL(8):: v
-      REAL(8), ALLOCATABLE:: T(:), normal(:)
+      REAL(8), ALLOCATABLE:: temperature(:), normal(:)
      REAL(8):: flow
      CHARACTER(:), ALLOCATABLE:: units
      INTEGER:: physicalSurface
      INTEGER:: particlesPerEdge
      INTEGER:: sp
      CALL config%info('inject', found, n_children = nInject)
      ALLOCATE(inject(1:nInject))
      nPartInj = 0
      DO i = 1, nInject
-        WRITE(istring, '(i2)') i
+        WRITE(iString, '(i2)') i
-        object = 'inject(' // trim(istring) // ')'
+        object = 'inject(' // trim(iString) // ')'
        !Find species
        CALL config%get(object // '.species', speciesName, found)
@ -1240,7 +1321,7 @@ MODULE moduleInput
        CALL config%get(object // '.name',            name,            found)
        CALL config%get(object // '.v',               v,               found)
-        CALL config%get(object // '.T',               T,               found)
+        CALL config%get(object // '.T',               temperature,     found)
        CALL config%get(object // '.n',               normal,          found)
        IF (.NOT. found) THEN
          ALLOCATE(normal(1:3))
@ -1249,8 +1330,10 @@ MODULE moduleInput
        CALL config%get(object // '.flow',            flow,            found)
        CALL config%get(object // '.units',           units,           found)
        CALL config%get(object // '.physicalSurface', physicalSurface, found)
        particlesPerEdge = 0
        CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
-        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
+        CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
        CALL readVelDistr(config, inject(i), object)
@ -1269,6 +1352,7 @@ MODULE moduleInput
      USE moduleCaseParam, ONLY: tauMin
      USE moduleMesh, ONLY: mesh
      USE moduleSpecies, ONLY: nSpecies
      USE moduleRefParam, ONLY: ti_ref
      IMPLICIT NONE
      TYPE(json_file), INTENT(inout):: config
@ -1282,8 +1366,11 @@ MODULE moduleInput
        CALL config%get('average.startTime', tStart, found)
        IF (found) THEN
-          tAverageStart = INT(tStart / tauMin)
+          tAverageStart = INT(tStart / ti_ref / tauMin)
        ELSE
          tAverageStart = 0
        END IF
        ALLOCATE(averageScheme(1:mesh%numNodes))
@ -1310,28 +1397,28 @@ MODULE moduleInput
      TYPE(injectGeneric), INTENT(inout):: inj
      CHARACTER(:), ALLOCATABLE, INTENT(in):: object
      INTEGER:: i
-      CHARACTER(2):: istring
+      CHARACTER(2):: iString
      CHARACTER(:), ALLOCATABLE:: fvType
      LOGICAL:: found
-      REAL(8):: v, T, m
+      REAL(8):: v, temperature, m
      !Reads species mass
      m = inj%species%m
      !Reads distribution functions for velocity
      DO i = 1, 3
-        WRITE(istring, '(i2)') i
+        WRITE(iString, '(i2)') i
-        CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found)
+        CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found)
-        IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // &
+        IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // &
                                            " found for " // object, 'readVelDistr')
        SELECT CASE(fvType)
        CASE ("Maxwellian")
-          T = inj%T(i)
+          temperature = inj%temperature(i)
-          CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
+          CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m)
        CASE ("Half-Maxwellian")
-          T = inj%T(i)
+          temperature = inj%temperature(i)
-          CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m)
+          CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m)
        CASE ("Delta")
          v = inj%vMod*inj%n(i)
@ -1372,5 +1459,37 @@ MODULE moduleInput
    END SUBROUTINE readParallel
    SUBROUTINE initOutput(inputFile)
      USE moduleRefParam
      USE moduleMesh,   ONLY: mesh, doubleMesh, pathMeshParticle, pathMeshColl
      USE moduleOutput, ONLY: path, folder
      IMPLICIT NONE
      CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
      INTEGER:: fileReference = 30
      !If everything is correct, creates the output folder
      CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
      !Copies input file to output folder
      CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
      !Copies particle mesh
      IF (mesh%dimen > 0) THEN
        CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
        IF (doubleMesh) THEN
          CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
        END IF
      END IF
      ! Write commit of fpakc
      CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpakc_commit.txt')
      ! Write file with reference values
      OPEN (fileReference, file=path // folder // '/' // 'reference.txt')
      WRITE(fileReference, "(7(1X,A20))")      'L_ref', 'v_ref', 'ti_ref', 'Vol_ref', 'EF_ref', 'Volt_ref', 'B_ref'
      WRITE(fileReference, "(7(1X,ES20.6E3))")  L_ref,   v_ref,   ti_ref,   Vol_ref,   EF_ref,   Volt_ref,   B_ref
      CLOSE(fileReference)
    END SUBROUTINE initOutput
 END MODULE moduleInput
--- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90
+++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90
@ -122,6 +122,8 @@ MODULE moduleMesh1DCart
      self%x = r1(1)
      self%surface = 1.D0
      self%normal = (/ 1.D0, 0.D0, 0.D0 /)
      !Boundary index
--- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90
+++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
@ -122,6 +122,8 @@ MODULE moduleMesh1DRad
      self%r = r1(1)
      self%surface = 1.D0
      self%normal = (/ 1.D0, 0.D0, 0.D0 /)
      !Boundary index
--- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90
+++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
@ -163,7 +163,7 @@ MODULE moduleMesh2DCart
      r2 = self%n2%getCoordinates()
      self%x  = (/r1(1), r2(1)/)
      self%y  = (/r1(2), r2(2)/)
-      self%weight = 1.D0
+      self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2)
      !Normal vector
      self%normal = (/ -(self%y(2)-self%y(1)), &
                         self%x(2)-self%x(1) , &
@ -318,6 +318,8 @@ MODULE moduleMesh2DCart
      INTEGER, INTENT(in):: nNodes
      REAL(8):: fPsi(1:nNodes)
      fPsi = 0.D0
      fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
                (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
                (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
@ -492,34 +494,36 @@ MODULE moduleMesh2DCart
    END FUNCTION insideQuad
-    !Transform physical coordinates to element coordinates
+    !Transform physical coordinates to element coordinates with a Taylor series
    PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 
      IMPLICIT NONE
      CLASS(meshCell2DCartQuad), INTENT(in):: self
      REAL(8), INTENT(in):: r(1:3)
      REAL(8):: Xi(1:3)
-      REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
+      REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2)
      REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
      REAL(8):: pDer(1:3, 1:3)
      REAL(8):: conv
      !Iterative newton method to transform coordinates
-      conv = 1.D0
+      conv  = 1.D0
-      XiO  = 0.D0
+      Xi0   = 0.D0
      Xi(3) = 0.D0
      DO WHILE(conv > 1.D-4)
-        dPsi = self%dPsi(XiO, 4)
+        fPsi           = self%fPsi(Xi0, 4)
-        pDer = self%partialDer(4, dPsi)
+        x0(1)          = dot_product(fPsi, self%x)
-        detJ = self%detJac(pDer)
+        x0(2)          = dot_product(fPsi, self%y)
-        invJ = self%invJac(pDer)
+        deltaR         = r(1:2) - x0
-        fPsi = self%fPsi(XiO, 4)
+        dPsi           = self%dPsi(Xi0, 4)
-        f = (/ DOT_PRODUCT(fPsi,self%x), &
+        pDer           = self%partialDer(4, dPsi)
-               DOT_PRODUCT(fPsi,self%y), &
+        detJ           = self%detJac(pDer)
-               0.D0 /) - r
+        pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-        Xi = XiO - MATMUL(invJ, f)/detJ
+        pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
-        conv = MAXVAL(DABS(Xi-XiO),1)
+        Xi(1:2)        = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ
-        XiO = Xi
+        conv           = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1)
        Xi0(1:2)       = Xi(1:2)
      END DO
@ -569,6 +573,7 @@ MODULE moduleMesh2DCart
      pDer = self%partialDer(4, dPsi)
      detJ = self%detJac(pDer)
      fPsi = self%fPsi(Xi, 4)
      !Compute total volume of the cell
      self%volume  = detJ*4.D0
      !Compute volume per node
@ -675,8 +680,8 @@ MODULE moduleMesh2DCart
      dPsi = 0.D0
-      dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /)
+      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
-      dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
+      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
    END FUNCTION dPsiTria
@ -762,6 +767,7 @@ MODULE moduleMesh2DCart
        pDer   = self%partialDer(3, dPsi)
        detJ   = self%detJac(pDer)
        invJ   = self%invJac(pDer)
        localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
      END DO
@ -820,19 +826,19 @@ MODULE moduleMesh2DCart
      CLASS(meshCell2DCartTria), INTENT(in):: self
      REAL(8), INTENT(in):: r(1:3)
      REAL(8):: Xi(1:3)
-      REAL(8):: deltaR(1:3)
+      REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2)
-      REAL(8):: dPsi(1:3, 1:3)
+      REAL(8):: dPsi(1:3,1:4)
      REAL(8):: pDer(1:3, 1:3)
      REAL(8):: invJ(1:3,  1:3), detJ
      !Direct method to convert coordinates
-      Xi     = 0.D0
+      Xi(3)          = 0.D0
-      deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /)
+      deltaR         = (/ r(1) - self%x(1), r(2) - self%y(1) /)
-      dPsi   = self%dPsi(Xi, 3)
+      dPsi           = self%dPsi(Xi, 3)
-      pDer   = self%partialDer(3, dPsi)
+      pDer           = self%partialDer(3, dPsi)
-      invJ   = self%invJac(pDer)
+      detJ           = self%detJac(pDer)
-      detJ   = self%detJac(pDer)
+      pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-      Xi     = MATMUL(invJ,deltaR)/detJ
+      pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
      Xi(1:2)        = MATMUL(pDerInv,deltaR)/detJ
    END FUNCTION phy2logTria
@ -907,8 +913,8 @@ MODULE moduleMesh2DCart
      invJ = 0.D0
-      invJ(1, 1:2) = (/  pDer(2,2), -pDer(1,2) /)
+      invJ(1, 1:2) = (/  pDer(2,2), -pDer(2,1) /)
-      invJ(2, 1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+      invJ(2, 1:2) = (/ -pDer(1,2),  pDer(1,1) /)
      invJ(3, 3)   = 1.D0
    END FUNCTION invJ2DCart
--- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
@ -144,6 +144,7 @@ MODULE moduleMesh2DCyl
      USE moduleSpecies
      USE moduleBoundary
      USE moduleErrors
      USE moduleConstParam, ONLY: PI
      IMPLICIT NONE
      CLASS(meshEdge2DCyl), INTENT(out):: self
@ -163,7 +164,15 @@ MODULE moduleMesh2DCyl
      r2 = self%n2%getCoordinates()
      self%z  = (/r1(1), r2(1)/)
      self%r  = (/r1(2), r2(2)/)
-      self%weight = r2(2)**2 - r1(2)**2
+      !Edge surface
      IF (self%z(2) /= self%z(1)) THEN
        self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1))
      ELSE
        self%surface = ABS(self%r(2)**2 - self%r(1)**2)
      END IF
      self%surface = self%surface * PI
      !Normal vector
      self%normal = (/ -(self%r(2)-self%r(1)), &
                         self%z(2)-self%z(1) , &
@ -223,21 +232,13 @@ MODULE moduleMesh2DCyl
      CLASS(meshEdge2DCyl), INTENT(in):: self
      REAL(8):: rnd
      REAL(8):: r(1:3)
-      REAL(8):: dr, dz
+      REAL(8):: p1(1:2), p2(1:2)
      rnd = random()
      dr  = self%r(2) - self%r(1)
      dz  = self%z(2) - self%z(1)
      IF (dr /= 0.D0) THEN
        r(2) = dr * DSQRT(rnd) + self%r(1)
        r(1) = dz * (r(2) - self%r(1))/dr + self%z(1)
      ELSE
        r(2) = self%r(1)
        r(1) = dz * rnd + self%z(1)
      END IF
      p1 = (/self%z(1), self%r(1) /)
      p2 = (/self%z(2), self%r(2) /)
      r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
      r(3)   =  0.D0
    END FUNCTION randPosEdge
@ -246,7 +247,6 @@ MODULE moduleMesh2DCyl
    !QUAD FUNCTIONS
    !Init element
    SUBROUTINE initCellQuad2DCyl(self, n, p, nodes)
      USE moduleRefParam
      IMPLICIT NONE
      CLASS(meshCell2DCylQuad), INTENT(out):: self
@ -326,6 +326,8 @@ MODULE moduleMesh2DCyl
      INTEGER, INTENT(in):: nNodes
      REAL(8):: fPsi(1:nNodes)
      fPsi = 0.D0
      fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
                (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
                (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
@ -496,7 +498,7 @@ MODULE moduleMesh2DCyl
    END FUNCTION elemFQuad
-    !Checks if Xi is inside the element
+    !Check if Xi is inside the element
    PURE FUNCTION insideQuad(Xi) RESULT(ins)
      IMPLICIT NONE
@ -508,34 +510,36 @@ MODULE moduleMesh2DCyl
    END FUNCTION insideQuad
-    !Transform physical coordinates to element coordinates
+    !Transform physical coordinates to element coordinates with a Taylor series
    PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 
      IMPLICIT NONE
      CLASS(meshCell2DCylQuad), INTENT(in):: self
      REAL(8), INTENT(in):: r(1:3)
      REAL(8):: Xi(1:3)
-      REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
+      REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2)
      REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
      REAL(8):: pDer(1:3, 1:3)
      REAL(8):: conv
      !Iterative newton method to transform coordinates
-      conv = 1.D0
+      conv  = 1.D0
-      XiO  = 0.D0
+      Xi0   = 0.D0
      Xi(3) = 0.D0
      DO WHILE(conv > 1.D-4)
-        dPsi = self%dPsi(XiO, 4)
+        fPsi           = self%fPsi(Xi0, 4)
-        pDer = self%partialDer(4, dPsi)
+        x0(1)          = dot_product(fPsi, self%z)
-        detJ = self%detJac(pDer)
+        x0(2)          = dot_product(fPsi, self%r)
-        invJ = self%invJac(pDer)
+        deltaR         = r(1:2) - x0
-        fPsi = self%fPsi(XiO, 4)
+        dPsi           = self%dPsi(Xi0, 4)
-        f = (/ DOT_PRODUCT(fPsi,self%z), &
+        pDer           = self%partialDer(4, dPsi)
-               DOT_PRODUCT(fPsi,self%r), &
+        detJ           = self%detJac(pDer)
-               0.D0 /) - r
+        pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-        Xi = XiO - MATMUL(invJ, f)/detJ
+        pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
-        conv = MAXVAL(DABS(Xi-XiO),1)
+        Xi(1:2)        = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ
-        XiO = Xi
+        conv           = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1)
        Xi0(1:2)       = Xi(1:2)
      END DO
@ -553,7 +557,7 @@ MODULE moduleMesh2DCyl
      XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
      nextInt = MAXLOC(XiArray,1)
-      !Selects the higher value of directions and searches in that direction
+      !Select the higher value of directions and searches in that direction
      NULLIFY(neighbourElement)
      SELECT CASE (nextInt)
      CASE (1)
@ -581,6 +585,7 @@ MODULE moduleMesh2DCyl
      REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3)
      self%volume  = 0.D0
      !2D 1 point Gauss Quad Integral
      Xi = 0.D0
      dPsi = self%dPsi(Xi, 4)
@ -589,18 +594,18 @@ MODULE moduleMesh2DCyl
      fPsi = self%fPsi(Xi, 4)
      r = DOT_PRODUCT(fPsi,self%r)
      !Computes total volume of the cell
-      self%volume  = r*detJ*PI8 !4*2*pi
+      self%volume  = r*detJ*PI8 !2*pi * 4 (weight of 1 point 2D-Gaussian integral)
-      !Computes volume per node
+      !Computes volume per node. Change the radius point to calculate the area to improve accuracy near the axis.
-      Xi = (/-5.D-1, -5.D-1, 0.D0/)
+      Xi = (/-5.D-1, -0.25D0, 0.D0/)
      r = self%gatherF(Xi, 4, self%r)
      self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8
-      Xi = (/ 5.D-1, -5.D-1, 0.D0/)
+      Xi = (/ 5.D-1, -0.25D0, 0.D0/)
      r = self%gatherF(Xi, 4, self%r)
      self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8
-      Xi = (/ 5.D-1,  5.D-1, 0.D0/)
+      Xi = (/ 5.D-1,  0.75D0, 0.D0/)
      r = self%gatherF(Xi, 4, self%r)
      self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8
-      Xi = (/-5.D-1,  5.D-1, 0.D0/)
+      Xi = (/-5.D-1,  0.75D0, 0.D0/)
      r = self%gatherF(Xi, 4, self%r)
      self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8
@ -702,8 +707,8 @@ MODULE moduleMesh2DCyl
      dPsi = 0.D0
-      dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /)
+      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
-      dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
+      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
    END FUNCTION dPsiTria
@ -855,19 +860,19 @@ MODULE moduleMesh2DCyl
      CLASS(meshCell2DCylTria), INTENT(in):: self
      REAL(8), INTENT(in):: r(1:3)
      REAL(8):: Xi(1:3)
-      REAL(8):: deltaR(1:3)
+      REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2)
-      REAL(8):: dPsi(1:3, 1:3)
+      REAL(8):: dPsi(1:3,1:4)
      REAL(8):: pDer(1:3, 1:3)
      REAL(8):: invJ(1:3,  1:3), detJ
      !Direct method to convert coordinates
-      Xi     = 0.D0
+      Xi(3)          = 0.D0
-      deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /)
+      deltaR         = (/ r(1) - self%z(1), r(2) - self%r(1) /)
-      dPsi   = self%dPsi(Xi, 3)
+      dPsi           = self%dPsi(Xi, 3)
-      pDer   = self%partialDer(3, dPsi)
+      pDer           = self%partialDer(3, dPsi)
-      invJ   = self%invJac(pDer)
+      detJ           = self%detJac(pDer)
-      detJ   = self%detJac(pDer)
+      pDerInv(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-      Xi     = MATMUL(invJ,deltaR)/detJ
+      pDerInv(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
      Xi(1:2)        = MATMUL(pDerInv,deltaR)/detJ
    END FUNCTION phy2logTria
@ -945,8 +950,8 @@ MODULE moduleMesh2DCyl
      invJ = 0.D0
-      invJ(1, 1:2) = (/  pDer(2,2), -pDer(1,2) /)
+      invJ(1, 1:2) = (/  pDer(2,2), -pDer(2,1) /)
-      invJ(2, 1:2) = (/ -pDer(2,1),  pDer(1,1) /)
+      invJ(2, 1:2) = (/ -pDer(1,2),  pDer(1,1) /)
      invJ(3, 3)   = 1.D0
    END FUNCTION invJ2DCyl
--- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90
+++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90
@ -109,6 +109,7 @@ MODULE moduleMesh3DCart
      USE moduleBoundary
      USE moduleErrors
      USE moduleMath
      USE moduleRefParam, ONLY: L_ref
      IMPLICIT NONE
      CLASS(meshEdge3DCartTria), INTENT(out):: self
@ -142,6 +143,8 @@ MODULE moduleMesh3DCart
      self%normal = crossProduct(vec1, vec2)
      self%normal = normalize(self%normal)
      self%surface = 1.D0/L_ref**2 !TODO: FIX THIS WHEN MOVING TO 3D
      !Boundary index
      self%boundary => boundary(bt)
      ALLOCATE(self%fBoundary(1:nSpecies))
--- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
+++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
@ -1,22 +1,22 @@
 MODULE moduleMeshOutput0D
  CONTAINS
-    SUBROUTINE printOutput0D(self, t)
+    SUBROUTINE printOutput0D(self)
      USE moduleMesh
      USE moduleRefParam
      USE moduleSpecies
      USE moduleOutput
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: i
      TYPE(outputFormat):: output
      CHARACTER(:), ALLOCATABLE:: fileName
      DO i = 1, nSpecies
        fileName='OUTPUT_' // species(i)%obj%name // '.dat'
-        IF (t == 0) THEN
+        IF (timeStep == 0) THEN
          OPEN(20, file = path // folder // '/' // fileName, action = 'write')
          WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", &
                                                                        "pressure (Pa)", "temperature (K)"
@ -27,14 +27,17 @@ MODULE moduleMeshOutput0D
        OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
        CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj)
-        WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature
+        WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density,    &
                                                                 output%velocity,   &
                                                                 output%pressure,   &
                                                                 output%temperature
        CLOSE(20)
      END DO
    END SUBROUTINE printOutput0D
-    SUBROUTINE printColl0D(self, t)
+    SUBROUTINE printColl0D(self)
      USE moduleMesh
      USE moduleRefParam
      USE moduleCaseParam
@ -43,12 +46,11 @@ MODULE moduleMeshOutput0D
      IMPLICIT NONE
      CLASS(meshGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: t
      CHARACTER(:), ALLOCATABLE:: fileName
      INTEGER:: k
      fileName='OUTPUT_Collisions.dat'
-      IF (t == tInitial) THEN
+      IF (timeStep == tInitial) THEN
        OPEN(20, file = path // folder // '/' // fileName, action = 'write')
        WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions"
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
@ -57,12 +59,12 @@ MODULE moduleMeshOutput0D
      END IF
      OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
-      WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
+      WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
      CLOSE(20)
    END SUBROUTINE printColl0D
-    SUBROUTINE printEM0D(self, t)
+    SUBROUTINE printEM0D(self)
      USE moduleMesh
      USE moduleRefParam
      USE moduleCaseParam
@ -70,7 +72,6 @@ MODULE moduleMeshOutput0D
      IMPLICIT NONE
      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
    END SUBROUTINE printEM0D
--- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
@ -108,6 +108,7 @@ MODULE moduleMeshInputGmsh2
      READ(10, *) totalNumElem
      !Count edges and volume elements
      numEdges = 0
      SELECT TYPE(self)
      TYPE IS(meshParticles)
        self%numEdges = 0
@ -328,7 +329,7 @@ MODULE moduleMeshInputGmsh2
      DO i = 1, numNodes
        !Reads the density
-        READ(10, *), e, density(i)
+        READ(10, *) e, density(i)
      END DO
@ -339,7 +340,7 @@ MODULE moduleMeshInputGmsh2
      DO i = 1, numNodes
        !Reads the velocity
-        READ(10, *), e, velocity(i, 1:3)
+        READ(10, *) e, velocity(i, 1:3)
      END DO
--- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2
    END SUBROUTINE writeGmsh2FooterElementData
    !Prints the scattered properties of particles into the nodes
-    SUBROUTINE printOutputGmsh2(self, t)
+    SUBROUTINE printOutputGmsh2(self)
      USE moduleMesh
      USE moduleRefParam
      USE moduleSpecies
      USE moduleOutput
      USE moduleMeshInoutCommon
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: n, i
      TYPE(outputFormat):: output(1:self%numNodes)
      REAL(8):: time
      CHARACTER(:), ALLOCATABLE:: fileName
-      time = DBLE(t)*tauMin*ti_ref
+      time = DBLE(timeStep)*tauMin*ti_ref
      DO i = 1, nSpecies
-        fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t)
+        fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (60, file = path // folder // '/' //  fileName)
        CALL writeGmsh2HeaderMesh(60)
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, time, 1, self%numNodes)
        DO n=1, self%numNodes
          CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
          WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
        END DO
        CALL writeGmsh2FooterNodeData(60)
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, time, 3, self%numNodes)
        DO n=1, self%numNodes
          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
        END DO
        CALL writeGmsh2FooterNodeData(60)
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, time, 1, self%numNodes)
        DO n=1, self%numNodes
          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
        END DO
        CALL writeGmsh2FooterNodeData(60)
-        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, time, 1, self%numNodes)
        DO n=1, self%numNodes
          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
        END DO
@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2
    END SUBROUTINE printOutputGmsh2
    !Prints the number of collisions into the volumes
-    SUBROUTINE printCollGmsh2(self, t)
+    SUBROUTINE printCollGmsh2(self)
      USE moduleMesh
      USE moduleRefParam
      USE moduleCaseParam
@ -145,7 +145,6 @@ MODULE moduleMeshOutputGmsh2
      IMPLICIT NONE
      CLASS(meshGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: numEdges
      INTEGER:: k, c
      INTEGER:: n
@ -167,9 +166,9 @@ MODULE moduleMeshOutputGmsh2
      END SELECT
      IF (collOutput) THEN
-        time = DBLE(t)*tauMin*ti_ref
+        time = DBLE(timeStep)*tauMin*ti_ref
-        fileName = formatFileName(prefix, 'Collisions', 'msh', t)
+        fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (60, file = path // folder // '/' //  fileName)
@ -179,7 +178,7 @@ MODULE moduleMeshOutputGmsh2
          DO c = 1, interactionMatrix(k)%amount
            WRITE(cString, "(I2)") c
            title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
-            CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells)
+            CALL writeGmsh2HeaderElementData(60, title, timeStep, time, 1, self%numCells)
            DO n=1, self%numCells
              WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c)
            END DO
@ -196,7 +195,7 @@ MODULE moduleMeshOutputGmsh2
    END SUBROUTINE printCollGmsh2
    !Prints the electrostatic EM properties into the nodes and volumes
-    SUBROUTINE printEMGmsh2(self, t)
+    SUBROUTINE printEMGmsh2(self)
      USE moduleMesh
      USE moduleRefParam
      USE moduleCaseParam
@ -205,7 +204,6 @@ MODULE moduleMeshOutputGmsh2
      IMPLICIT NONE
      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: n, e
      REAL(8):: time
      CHARACTER(:), ALLOCATABLE:: fileName
@ -214,27 +212,27 @@ MODULE moduleMeshOutputGmsh2
      Xi = (/ 0.D0, 0.D0, 0.D0 /)
      IF (emOutput) THEN
-        time = DBLE(t)*tauMin*ti_ref
+        time = DBLE(timeStep)*tauMin*ti_ref
-        fileName = formatFileName(prefix, 'EMField', 'msh', t)
+        fileName = formatFileName(prefix, 'EMField', 'msh', timeStep)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (20, file = path // folder // '/' //  fileName)
        CALL writeGmsh2HeaderMesh(20)
-        CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, time, 1, self%numNodes)
        DO n=1, self%numNodes
          WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
        END DO
        CALL writeGmsh2FooterNodeData(20)
-        CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells)
+        CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, time, 3, self%numCells)
        DO e=1, self%numCells
          WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref
        END DO
        CALL writeGmsh2FooterElementData(20)
-        CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, time, 3, self%numNodes)
        DO n=1, self%numNodes
          WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
        END DO
--- a/src/modules/mesh/inout/makefile
+++ b/src/modules/mesh/inout/makefile
@ -1,4 +1,4 @@
-all: vtu.o gmsh2.o 0D.o
+all: vtu.o gmsh2.o 0D.o text.o
 vtu.o: moduleMeshInoutCommon.o
 	$(MAKE) -C vtu all
@ -9,5 +9,8 @@ gmsh2.o:
 0D.o:
 	$(MAKE) -C 0D all
 text.o:
 	$(MAKE) -C text all
 %.o: %.f90
 	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
--- a/src/modules/mesh/inout/moduleMeshInoutCommon.f90
+++ b/src/modules/mesh/inout/moduleMeshInoutCommon.f90
@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon
  CHARACTER(LEN=4):: prefix = 'Step'
  CONTAINS
-    PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName)
+    PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName)
      USE moduleOutput
      IMPLICIT NONE
      CHARACTER(*), INTENT(in):: prefix, suffix, extension
-      INTEGER, INTENT(in), OPTIONAL:: t
+      INTEGER, INTENT(in), OPTIONAL:: timeStep
      CHARACTER (LEN=iterationDigits):: tString
      CHARACTER(:), ALLOCATABLE:: fileName
-      IF (PRESENT(t)) THEN
+      IF (PRESENT(timeStep)) THEN
-        WRITE(tString, iterationFormat) t
+        WRITE(tString, iterationFormat) timeStep
        fileName = prefix //  '_' // tString // '_' // suffix // '.' // extension
      ELSE
--- a/src/modules/mesh/inout/text/makefile
+++ b/src/modules/mesh/inout/text/makefile
@ -0,0 +1,7 @@
 all: moduleMeshInputText.o moduleMeshOutputText.o
 moduleMeshInputText.o: moduleMeshOutputText.o moduleMeshInputText.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 %.o: %.f90
 	$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@
--- a/src/modules/mesh/inout/text/moduleMeshInputText.f90
+++ b/src/modules/mesh/inout/text/moduleMeshInputText.f90
@ -0,0 +1,232 @@
 module moduleMeshInputText
  !The mesh is stored as a column-wise text file.
  !Aimed for simple geometries in 1D
  contains
    !Inits the text mesh
    subroutine initText(self)
      use moduleMesh
      use moduleMeshOutputText
      implicit none
      class(meshGeneric), intent(inout), target:: self 
      if (associated(meshForMCC,self)) then
        self%printColl => printCollText
      end if
      select type(self)
      type is (meshParticles)
        self%printOutput  => printOutputText
        self%printEM      => printEMText
        self%printAverage => printAverageText
        self%readInitial  => readInitialText
      end select
      self%readMesh => readText
    end subroutine initText
    !Reads the text mesh
    subroutine readText(self, filename)
      use moduleMesh
      use moduleMesh1DCart
      use moduleMesh1DRad
      use moduleErrors
      implicit none
      class(meshGeneric), intent(inout):: self
      character(:), allocatable, intent(in):: filename !Dummy file, not used
      integer:: fileID, reason
      character(len=256):: line
      integer:: nNodes
      real(8):: r(1:3) !dummy 3D coordinate
      integer:: physicalID
      integer:: n, c
      integer, allocatable:: p(:)
      integer:: bt
      fileID = 10
      open(fileID, file=trim(filename))
      !Skip header
      read(fileID, *)
      !Get number of nodes
      nNodes = 0
      do
        read(fileID, *, iostat=reason) line
        if (reason > 0) then
          call criticalError('Error reading mesh file', 'readText')
        else if (reason < 0) then
          exit
        else if (len(line) > 0) then
          nNodes = nNodes + 1
        end if
      end do
      if (nNodes == 0) then
        call criticalError('No nodes read in mesh file', 'readText')
      end if
      self%numNodes = nNodes
      allocate(self%nodes(1:self%numNodes))
      SELECT TYPE(self)
      TYPE IS(meshParticles)
        ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes))
        ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes))
        self%K = 0.D0
        self%IPIV = 0
      END SELECT
      self%numCells = nNodes - 1
      allocate(self%cells(1:self%numCells))
      select type(self)
      type is (meshParticles)
        self%numEdges = 2
        allocate(self%edges(1:self%numEdges))
      end select
      !Read the mesh now
      rewind(fileID)
      !Skip header
      read(fileID, *)
      !Allocate nodes and edges
      do n = 1, self%numNodes
        r = 0.D0
        read(fileID, *) r(1), physicalID
        select case(self%geometry)
        case("Cart")
          allocate(meshNode1DCart:: self%nodes(n)%obj)
        case("Rad")
          allocate(meshNode1DRad:: self%nodes(n)%obj)
        end select
        !Init nodes
        call self%nodes(n)%obj%init(n, r)
        !Allocate edges if required)
        select type(self)
        type is (meshParticles)
          if ((physicalID == 1) .or. (physicalID == 2)) then
            select case(self%geometry)
            case("Cart")
              allocate(meshEdge1DCart:: self%edges(physicalID)%obj)
            case("Rad")
              allocate(meshEdge1DRad:: self%edges(physicalID)%obj)
            end select
            allocate(p(1))
            p(1) = n
            bt = getBoundaryId(physicalID)
            call self%edges(physicalID)%obj%init(physicalID, p, physicalID, physicalID)
            deallocate(p)
          end if
        end select
      end do
      !Allocate cells
      n = 1
      allocate(p(1:2))
      do c = 1, self%numCells
        p(1) = n
        n = n + 1
        p(2) = n
        select case(self%geometry)
        case("Cart")
          allocate(meshCell1DCartSegm:: self%cells(c)%obj)
        case("Rad")
          allocate(meshCell1DRadSegm:: self%cells(c)%obj)
        end select
        call self%cells(c)%obj%init(c, p, self%nodes)
      end do
      deallocate(p)
      close(fileID)
      !Call mesh connectivity
      CALL self%connectMesh
    end subroutine readText
    subroutine readInitialText(filename, density, velocity, temperature)
      use moduleErrors
      implicit none
      character(:), allocatable, intent(in):: filename
      real(8), allocatable, intent(out), dimension(:):: density
      real(8), allocatable, intent(out), dimension(:,:):: velocity
      real(8), allocatable, intent(out), dimension(:):: temperature
      integer:: fileID, reason
      character(len=256):: line
      integer:: nNodes
      integer:: n
      fileID = 10
      open(fileID, file=trim(filename))
      do 
        read(fileID, *, iostat=reason) line
        if (reason > 0) then
          call criticalError('Error reading mesh file', 'readText')
        else if (reason < 0) then
          exit
        else if (len(line) > 0) then
          nNodes = nNodes + 1
        end if
      end do
      allocate(density(1:nNodes))
      allocate(velocity(1:nNodes, 1:3))
      allocate(temperature(1:nNodes))
      rewind(fileID)
      do n = 1, nNodes
        read(fileID, *) density(n), velocity(n, 1:3), temperature(n)
      end do
      close(fileID)
    end subroutine readInitialText
 end module moduleMeshInputText
--- a/src/modules/mesh/inout/text/moduleMeshOutputText.f90
+++ b/src/modules/mesh/inout/text/moduleMeshOutputText.f90
@ -0,0 +1,265 @@
 module moduleMeshOutputText
  contains
    subroutine writeSpeciesOutput(self, fileID, speciesIndex)
      use moduleMesh
      use moduleOutput
      use moduleRefParam, only: L_ref
      implicit none
      class(meshParticles), INTENT(in):: self
      integer, intent(in):: fileID
      integer, intent(in):: speciesIndex
      real(8):: r(1:3)
      type(outputFormat):: output
      integer:: n
      do n = 1, self%numNodes
        r = self%nodes(n)%obj%getCoordinates()
        call calculateOutput(self%nodes(n)%obj%output(speciesIndex), output, self%nodes(n)%obj%v, species(speciesIndex)%obj)
        write(fileID, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, output%density, output%velocity, output%temperature
      end do
    end subroutine writeSpeciesOutput
    subroutine writeCollOutput(self, fileID)
      use moduleMesh
      use moduleCollisions
      use moduleRefParam, only: L_ref
      implicit none
      class(meshGeneric), intent(in):: self
      integer, intent(in):: fileID
      integer:: n, k, c
      do n = 1, self%numCells
        write(fileID, '(I0)', advance='no') n
        do k = 1, nCollPairs
          do c = 1, interactionMatrix(k)%amount
            write(fileID, '(",",I0)', advance='no') self%cells(n)%obj%tallyColl(k)%tally(c)
          end do
        end do
        write(fileID, *)
      end do
    end subroutine writeCollOutput
    subroutine writeEMOutput(self, fileID)
      use moduleMesh
      use moduleRefParam, only: L_ref, Volt_ref, B_ref, EF_ref
      implicit none
      class(meshParticles), intent(in):: self
      integer, intent(in):: fileID
      integer:: n, c
      real(8):: r(1:3), Xi(1:3)
      do n = 1, self%numNodes
        r = self%nodes(n)%obj%getCoordinates()
        if (n == self%numNodes) then
          Xi = (/ 1.D0, 0.D0, 0.D0 /)
          c = self%numNodes - 1
        else
          Xi = (/ 0.D0, 0.D0, 0.D0 /)
          c = n
        end if
        associate(output => self%nodes(n)%obj%emData)
          write(fileID, '(7(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, &
                                                    output%phi*Volt_ref, &
                                                    self%cells(c)%obj%gatherElectricField(Xi)*EF_ref, &
                                                    output%B*B_ref
        end associate
      end do
    end subroutine writeEMOutput
    subroutine writeAverage(self, fileIDMean, &
                                  fileIDDeviation, &
                                  speciesIndex)
      use moduleMesh
      use moduleOutput
      use moduleAverage
      use moduleRefParam, only: L_ref
      implicit none
      class(meshParticles), intent(in):: self
      integer, intent(in):: fileIDMean, fileIDDeviation
      INTEGER, intent(in):: speciesIndex
      real(8):: r(1:3)
      type(outputFormat):: outputMean
      type(outputFormat):: outputDeviation
      integer:: n
      do n = 1, self%numNodes
        r = self%nodes(n)%obj%getCoordinates()
        call calculateOutput(averageScheme(n)%mean%output(speciesIndex),      outputMean,      &
                             self%nodes(n)%obj%v, species(speciesIndex)%obj)
        write(fileIDMean, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, outputMean%density, outputMean%velocity, outputMean%temperature
        call calculateOutput(averageScheme(n)%deviation%output(speciesIndex), outputDeviation, &
                             self%nodes(n)%obj%v, species(speciesIndex)%obj)
        write(fileIDDeviation, '(5(ES0.6E3,","),ES0.6E3)') r(1)*L_ref, outputDeviation%density, outputDeviation%velocity, outputDeviation%temperature
      end do
    end subroutine writeAverage
    subroutine printOutputText(self)
      use moduleMesh
      use moduleSpecies
      use moduleMeshInoutCommon
      use moduleCaseParam, ONLY: timeStep
      implicit none
      class(meshParticles), intent(in):: self
      INTEGER:: s, fileID
      character(:), allocatable:: fileName
      fileID = 60
      do s = 1, nSpecies
        fileName = formatFileName(prefix, species(s)%obj%name, 'csv', timeStep)
        write(*, "(6X,A15,A)") "Creating file: ", fileName
        open (fileID, file = path // folder // '/' //  fileName)
        write(fileID, '(5(A,","),A)') 'Position (m)', &
                                      'Density (m^-3)', &
                                      'Velocity (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', &
                                      'Temperature (K)'
        call writeSpeciesOutput(self, fileID, s)
        close(fileID)
      end do
    end subroutine printOutputText
    subroutine printCollText(self)
      use moduleMesh
      use moduleOutput
      use moduleMeshInoutCommon
      use moduleCaseParam, only: timeStep
      implicit none
      class(meshGeneric), intent(in):: self
      integer:: fileID
      character(:), allocatable:: fileName
      integer:: k, c
      character (len=2):: cString
      fileID = 62
      if (collOutput) then
        fileName = formatFileName(prefix, 'Collisions', 'csv', timeStep)
        write(*, "(6X,A15,A)") "Creating file: ", fileName
        open (fileID, file = path // folder // '/' //  fileName)
        write(fileID, '(A)', advance='no') "Cell"
        do k = 1, nCollPairs
          do c = 1, interactionMatrix(k)%amount
            write(cString, "(I2)") c
            write(fileID, '(",",A)', advance='no') 'Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
          end do
        end do
        write(fileID, *)
        call writeCollOutput(self, fileID)
        close(fileID)
      end if
    end subroutine printCollText
    subroutine printEMText(self)
      use moduleMesh
      use moduleMeshInoutCommon
      use moduleCaseParam, only: timeStep
      implicit none
      class(meshParticles), intent(in):: self
      integer:: fileID
      character(:), allocatable:: fileName
      fileID = 64
      if (emOutput) then
        fileName = formatFileName(prefix, 'EMField', 'csv', timeStep)
        write(*, "(6X,A15,A)") "Creating file: ", fileName
        open (fileID, file = path // folder // '/' //  fileName)
        write(fileID, '(8(A,","),A)') 'Position (m)', &
                                      'Potential (V)', &
                                      'Electric Field (V m^-1):0', 'Electric Field (V m^-1):1', 'Electric Field (V m^-1):2', &
                                      'Magnetic Field (T):0',      'Magnetic Field (T):1',      'Magnetic Field (T):2'
        call writeEMOutput(self, fileID)
        close(fileID)
      end if
    end subroutine printEMText
    subroutine printAverageText(self)
      use moduleMesh
      use moduleSpecies
      use moduleMeshInoutCommon
      implicit none
      class(meshParticles), intent(in):: self
      integer:: s
      integer:: fileIDMean, fileIDDeviation
      character(:), allocatable:: fileNameMean, fileNameDeviation
      fileIDMean = 66
      fileIDDeviation = 67
      do s = 1, nSpecies
        fileNameMean = formatFileName('Average_mean', species(s)%obj%name, 'csv', timeStep)
        write(*, "(6X,A15,A)") "Creating file: ", fileNameMean
        open (fileIDMean, file = path // folder // '/' //  fileNameMean)
        write(fileIDMean, '(5(A,","),A)') 'Position (m)', &
                                          'Density, mean (m^-3)', &
                                          'Velocity, mean (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', &
                                          'Temperature, mean (K)'
        fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv', timeStep)
        write(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation
        open (fileIDDeviation, file = path // folder // '/' //  fileNameDeviation)
        write(fileIDDeviation, '(5(A,","),A)') 'Position (m)', &
                                               'Density, deviation (m^-3)', &
                                               'Velocity, deviation (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', &
                                               'Temperature, deviation (K)'
        call writeAverage(self, fileIDMean, fileIDDeviation, s)
        close(fileIDMean)
        close(fileIDDeviation)
      end do
    end subroutine printAverageText
 end module moduleMeshOutputText
--- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU
      CLASS(meshGeneric), INTENT(inout):: self
      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
      REAL(8):: r(1:3) !3 generic coordinates
-      INTEGER:: fileID, error, found
+      INTEGER:: fileID
      CHARACTER(LEN=256):: line
      INTEGER:: numNodes, numElements, numEdges
      INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
@ -275,6 +275,7 @@ MODULE moduleMeshInputVTU
      END DO
      !Count the number of edges
      numEdges = 0
      SELECT CASE(self%dimen)
      CASE(3)
        !Edges are triangles, type 5 in VTK
@ -495,7 +496,7 @@ MODULE moduleMeshInputVTU
        END SELECT
      END DO
-        
+
      !Call mesh connectivity
      CALL self%connectMesh
@ -547,6 +548,8 @@ MODULE moduleMeshInputVTU
      CALL readDataBlock(fileID, numNodes, temperature)
      REWIND(fileID)
      close(fileID)
    END SUBROUTINE readInitialVTU
 END MODULE moduleMeshInputVTU
--- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
@ -11,7 +11,7 @@ MODULE moduleMeshOutputVTU
      WRITE(fileID,"(A)") '<?xml version="1.0"?>'
      WRITE(fileID,"(2X, A)") '<VTKFile type="UnstructuredGrid">'
-      WRITE(fileID,"(4X, A,ES20.6E3,A)") '<UnstructuredGrid>'
+      WRITE(fileID,"(4X, A)") '<UnstructuredGrid>'
      WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', nNodes, '" NumberOfCells="', nCells, '">'
    END SUBROUTINE writeHeader
@ -209,23 +209,22 @@ MODULE moduleMeshOutputVTU
      WRITE(fileID,"(8X,A)") '<CellData>'
      !Electric field
      WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="Electric Field (V m^-1)" NumberOfComponents="3">'
-      WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
+      WRITE(fileID,"(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
      WRITE(fileID,"(10X,A)") '</DataArray>'
      WRITE(fileID,"(8X,A)") '</CellData>'
    END SUBROUTINE writeEM
-    SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
+    SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection)
      USE moduleCaseParam
      USE moduleOutput
      USE moduleRefParam
      IMPLICIT NONE
      INTEGER:: fileID
      INTEGER, INTENT(in):: t
      CHARACTER(*):: fileNameStep, fileNameCollection
-      IF (t == tInitial) THEN
+      IF (timeStep == tInitial) THEN
        !Create collection file
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
        OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection)
@ -237,10 +236,11 @@ MODULE moduleMeshOutputVTU
      !Write iteration file in collection
      OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection, ACCESS='APPEND')
-      WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
+      WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") &
            '<DataSet timestep="', DBLE(timeStep)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
      !Close collection file
-      IF (t == tFinal) THEN
+      IF (timeStep == tFinal) THEN
        WRITE (fileID + 1, "(2X, A)") '</Collection>'
        WRITE (fileID + 1, "(A)") '</VTKFile>'
@ -307,22 +307,21 @@ MODULE moduleMeshOutputVTU
    END SUBROUTINE writeAverage
-    SUBROUTINE printOutputVTU(self,t)
+    SUBROUTINE printOutputVTU(self)
      USE moduleMesh
      USE moduleSpecies
      USE moduleMeshInoutCommon
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      CLASS(meshParticles), INTENT(in):: self
-      INTEGER, INTENT(in):: t
+      INTEGER:: i, fileID
      INTEGER:: n, i, fileID
      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
      TYPE(outputFormat):: output(1:self%numNodes)
      fileID = 60
      DO i = 1, nSpecies
-        fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t)
+        fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (fileID, file = path // folder // '/' //  fileName)
@ -338,29 +337,27 @@ MODULE moduleMeshOutputVTU
        !Write collection file for time plotting
        fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd')
-        CALL writeCollection(fileID, t, fileName, filenameCollection)
+        CALL writeCollection(fileID, fileName, filenameCollection)
      END DO
    END SUBROUTINE printOutputVTU
-    SUBROUTINE printCollVTU(self,t)
+    SUBROUTINE printCollVTU(self)
      USE moduleMesh
      USE moduleOutput
      USE moduleMeshInoutCommon
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      CLASS(meshGeneric), INTENT(in):: self
-      INTEGER, INTENT(in):: t
+      INTEGER:: fileID
      INTEGER:: n, i, fileID
      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
      CHARACTER (LEN=iterationDigits):: tstring
      TYPE(outputFormat):: output(1:self%numNodes)
      fileID = 62
      IF (collOutput) THEN
-        fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
+        fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (fileID, file = path // folder // '/' //  fileName)
@ -376,26 +373,26 @@ MODULE moduleMeshOutputVTU
        !Write collection file for time plotting
        fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd')
-        CALL writeCollection(fileID, t, fileName, filenameCollection)
+        CALL writeCollection(fileID, fileName, filenameCollection)
      END IF
    END SUBROUTINE printCollVTU
-    SUBROUTINE printEMVTU(self, t)
+    SUBROUTINE printEMVTU(self)
      USE moduleMesh
      USE moduleMeshInoutCommon
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: fileID
      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
      fileID = 64
      IF (emOutput) THEN
-        fileName = formatFileName(prefix, 'EMField', 'vtu', t)
+        fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (fileID, file = path // folder // '/' //  fileName)
@ -411,7 +408,7 @@ MODULE moduleMeshOutputVTU
        !Write collection file for time plotting
        fileNameCollection = formatFileName('Collection', 'EMField', 'pvd')
-        CALL writeCollection(fileID, t, fileName, filenameCollection)
+        CALL writeCollection(fileID, fileName, filenameCollection)
      END IF
@ -424,9 +421,8 @@ MODULE moduleMeshOutputVTU
      IMPLICIT NONE
      CLASS(meshParticles), INTENT(in):: self
-      INTEGER:: n, i, fileIDMean, fileIDDeviation
+      INTEGER:: i, fileIDMean, fileIDDeviation
      CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation
      TYPE(outputFormat):: output(1:self%numNodes)
      fileIDMean = 66
      fileIDDeviation = 67
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@ -59,6 +59,13 @@ MODULE moduleMesh
  END TYPE meshNodeCont
  ! Array of pointers to nodes.
  TYPE:: meshNodePointer
    CLASS(meshNode), POINTER:: obj
    CONTAINS
  END TYPE meshNodePointer
  !Type for array of boundary functions (one per species)
  TYPE, PUBLIC:: fBoundaryGeneric
    PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL()
@ -76,8 +83,8 @@ MODULE moduleMesh
    CLASS(meshCell), POINTER:: eColl => NULL()
    !Normal vector
    REAL(8):: normal(1:3)
-    !Weight for random injection of particles
+    ! Surface of edge
-    REAL(8):: weight = 1.D0
+    REAL(8):: surface = 0.D0
    !Pointer to boundary type
    TYPE(boundaryCont), POINTER:: boundary
    !Array of functions for boundary conditions
@ -337,10 +344,10 @@ MODULE moduleMesh
    !Array of cell elements
    TYPE(meshCellCont),  ALLOCATABLE:: cells(:)
    !PROCEDURES SPECIFIC OF FILE TYPE
-    PROCEDURE(readMesh_interface),     POINTER, PASS::   readMesh     => NULL()
+    PROCEDURE(readMesh_interface),    POINTER, PASS::   readMesh     => NULL()
-    PROCEDURE(readInitial_interface),  POINTER, NOPASS:: readInitial  => NULL()
+    PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial  => NULL()
-    PROCEDURE(connectMesh_interface),  POINTER, PASS::   connectMesh  => NULL()
+    PROCEDURE(connectMesh_interface), POINTER, PASS::   connectMesh  => NULL()
-    PROCEDURE(printColl_interface),    POINTER, PASS::   printColl    => NULL()
+    PROCEDURE(printColl_interface),   POINTER, PASS::   printColl    => NULL()
    CONTAINS
      !GENERIC PROCEDURES
      PROCEDURE, PASS:: doCollisions
@ -372,10 +379,9 @@ MODULE moduleMesh
    END SUBROUTINE connectMesh_interface
    !Prints number of collisions in each cell
-    SUBROUTINE printColl_interface(self, t)
+    SUBROUTINE printColl_interface(self)
      IMPORT meshGeneric
      CLASS(meshGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: t
    END SUBROUTINE printColl_interface
@ -403,18 +409,16 @@ MODULE moduleMesh
  ABSTRACT INTERFACE
    !Prints Species data
-    SUBROUTINE printOutput_interface(self, t)
+    SUBROUTINE printOutput_interface(self)
      IMPORT meshParticles
      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
    END SUBROUTINE printOutput_interface
    !Prints EM info
-    SUBROUTINE printEM_interface(self, t)
+    SUBROUTINE printEM_interface(self)
      IMPORT meshParticles
      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
    END SUBROUTINE printEM_interface
@ -495,28 +499,29 @@ MODULE moduleMesh
      IMPLICIT NONE
      CLASS(meshParticles), INTENT(inout):: self
-      INTEGER:: e
+      INTEGER:: c
      INTEGER:: nNodes
      INTEGER, ALLOCATABLE:: n(:)
      REAL(8), ALLOCATABLE:: localK(:,:)
      INTEGER:: i, j
-      DO e = 1, self%numCells
+      DO c = 1, self%numCells
-        nNodes = self%cells(e)%obj%nNodes
+        associate(nNodes => self%cells(c)%obj%nNodes)
-        ALLOCATE(n(1:nNodes))
+          ALLOCATE(n(1:nNodes))
-        ALLOCATE(localK(1:nNodes, 1:nNodes))
+          ALLOCATE(localK(1:nNodes, 1:nNodes))
-        n      = self%cells(e)%obj%getNodes(nNodes)
+          n      = self%cells(c)%obj%getNodes(nNodes)
-        localK = self%cells(e)%obj%elemK(nNodes)
+          localK = self%cells(c)%obj%elemK(nNodes)
-        DO i = 1, nNodes
+          DO i = 1, nNodes
-          DO j = 1, nNodes
+            DO j = 1, nNodes
-            self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
+              self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
            END DO
          END DO
-        END DO
+          DEALLOCATE(n, localK)
-        DEALLOCATE(n, localK)
+        end associate
      END DO
@ -613,6 +618,7 @@ MODULE moduleMesh
      INTEGER:: sp
      INTEGER:: i
      CLASS(meshNode), POINTER:: node
      REAL(8):: pFraction !Particle fraction
      cellNodes = self%getNodes(nNodes)
      fPsi = self%fPsi(part%Xi, nNodes)
@ -623,10 +629,11 @@ MODULE moduleMesh
      DO i = 1, nNodes
        node => mesh%nodes(cellNodes(i))%obj
        pFraction = fPsi(i)*part%weight
        CALL OMP_SET_LOCK(node%lock)
-        node%output(sp)%den          = node%output(sp)%den          + part%weight*fPsi(i)
+        node%output(sp)%den          = node%output(sp)%den          + pFraction
-        node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*fPsi(i)*part%v(:)
+        node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + pFraction*part%v(:)
-        node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*tensorS
+        node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS
        CALL OMP_UNSET_LOCK(node%lock)
      END DO
@ -787,7 +794,7 @@ MODULE moduleMesh
    END FUNCTION findCellBrute
    !Computes collisions in element
-    SUBROUTINE doCollisions(self, t)
+    SUBROUTINE doCollisions(self)
      USE moduleCollisions
      USE moduleSpecies
      USE moduleList
@ -795,10 +802,10 @@ MODULE moduleMesh
      USE moduleRandom
      USE moduleOutput
      USE moduleMath
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      CLASS(meshGeneric), INTENT(inout), TARGET:: self
      INTEGER, INTENT(in):: t
      INTEGER:: e
      CLASS(meshCell), POINTER:: cell
      INTEGER:: k, i, j
@ -814,7 +821,7 @@ MODULE moduleMesh
      REAL(8):: rnd_real !Random number for collision
      INTEGER:: rnd_int  !Random number for collision
-      IF (MOD(t, everyColl) == 0) THEN
+      IF (MOD(timeStep, everyColl) == 0) THEN
        !Collisions need to be performed in this iteration
        !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
        DO e=1, self%numCells
@ -911,7 +918,9 @@ MODULE moduleMesh
                    !Loop over collisions
                    DO c = 1, interactionMatrix(k)%amount
                      IF (rnd_real <= probabilityColl(c)) THEN
                        !$OMP CRITICAL
                        CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
                        !$OMP END CRITICAL
                        !If collisions are gonna be output, count the collision
                        IF (collOutput) THEN
@ -1021,6 +1030,9 @@ MODULE moduleMesh
          ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3))
          ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3))
          ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount))
          deltaV_ij = 0.D0
          p_ij      = 0.D0
          mass_ij   = 0.D0
          !Loop over particles of species_i
          partTemp => cell%listPart_in(i)%head
          p = 1
@ -1105,6 +1117,9 @@ MODULE moduleMesh
            ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3))
            ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3))
            ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount))
            deltaV_ji = 0.D0
            p_ji      = 0.D0
            mass_ji   = 0.D0
            !Loop over particles of species_j
            partTemp => cell%listPart_in(j)%head
            p = 1
--- a/src/modules/mesh/moduleMeshBoundary.f90
+++ b/src/modules/mesh/moduleMeshBoundary.f90
@ -77,6 +77,20 @@ MODULE moduleMeshBoundary
    END SUBROUTINE transparent
    !Symmetry axis. Reflects particles.
    !Although this function should never be called, it is set as a reflective boundary
    !to properly deal with possible particles reaching a corner and selecting this boundary.
    SUBROUTINE symmetryAxis(edge, part)
      USE moduleSpecies
      IMPLICIT NONE
      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part
      CALL reflection(edge, part)
    END SUBROUTINE symmetryAxis
    !Wall with temperature
    SUBROUTINE wallTemperature(edge, part)
      USE moduleSpecies
@ -147,7 +161,13 @@ MODULE moduleMeshBoundary
            ALLOCATE(newElectron)
            ALLOCATE(newIon)
-            newElectron%species => part%species
+            IF (ASSOCIATED(bound%electronSecondary)) THEN
              newElectron%species => bound%electronSecondary
            ELSE
              newElectron%species => part%species
            END IF
            newIon%species      => bound%species
            newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0))
@ -198,19 +218,27 @@ MODULE moduleMeshBoundary
    END SUBROUTINE ionization
-    !Symmetry axis. Reflects particles.
+    subroutine outflowAdaptive(edge, part)
-    !Although this function should never be called, it is set as a reflective boundary
+      use moduleRandom
-    !to properly deal with possible particles reaching a corner and selecting this boundary.
+      implicit none
    SUBROUTINE symmetryAxis(edge, part)
      USE moduleSpecies
      IMPLICIT NONE
-      CLASS(meshEdge), INTENT(inout):: edge
+      class(meshEdge), intent(inout):: edge
-      CLASS(particle), INTENT(inout):: part
+      class(particle), intent(inout):: part
-      CALL reflection(edge, part)
+      select type(bound => edge%boundary%bTypes(part%species%n)%obj)
      type is(boundaryOutflowAdaptive)
-    END SUBROUTINE symmetryAxis
+        if (random() < 0.844d0) then
          call reflection(edge, part)
        else
          call transparent(edge, part)
        end if
      end select
    end subroutine outflowAdaptive
    !Points the boundary function to specific type
    SUBROUTINE pointBoundaryFunction(edge, s)
@ -230,14 +258,17 @@ MODULE moduleMeshBoundary
      TYPE IS(boundaryTransparent)
        edge%fBoundary(s)%apply => transparent
      TYPE IS(boundaryAxis)
        edge%fBoundary(s)%apply => symmetryAxis
      TYPE IS(boundaryWallTemperature)
        edge%fBoundary(s)%apply => wallTemperature
      TYPE IS(boundaryIonization)
        edge%fBoundary(s)%apply => ionization
-      TYPE IS(boundaryAxis)
+      type is(boundaryOutflowAdaptive)
-        edge%fBoundary(s)%apply => symmetryAxis
+        edge%fBoundary(s)%apply => outflowAdaptive
      CLASS DEFAULT
        CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')
--- a/src/modules/moduleBoundary.f90
+++ b/src/modules/moduleBoundary.f90
@ -26,6 +26,12 @@ MODULE moduleBoundary
  END TYPE boundaryTransparent
  !Symmetry axis
  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
    CONTAINS
  END TYPE boundaryAxis
  !Wall Temperature boundary
  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature
    !Thermal velocity of the wall: square root(Wall temperature X specific heat)
@ -38,6 +44,7 @@ MODULE moduleBoundary
  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
    REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
    CLASS(speciesGeneric), POINTER:: species !Ion species
    CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron
    TYPE(table1D):: crossSection
    REAL(8):: effectiveTime
    REAL(8):: eThreshold
@ -46,11 +53,13 @@ MODULE moduleBoundary
  END TYPE boundaryIonization
-  !Symmetry axis
+  !Boundary for quasi-neutral outflow adjusting reflection coefficient
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
+  type, public, extends(boundaryGeneric):: boundaryOutflowAdaptive
-    CONTAINS
+    real(8):: outflowCurrent
    real(8):: reflectionFraction
    contains
-  END TYPE boundaryAxis
+  end type boundaryOutflowAdaptive
  !Wrapper for boundary types (one per species)
  TYPE:: bTypesCont
@ -103,17 +112,19 @@ MODULE moduleBoundary
    END SUBROUTINE initWallTemperature
-    SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, speciesID, effTime, crossSection, eThreshold)
+    SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
      USE moduleRefParam
      USE moduleSpecies
      USE moduleCaseParam
      USE moduleConstParam
      USE moduleErrors
      IMPLICIT NONE
      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
      REAL(8), INTENT(in):: me !Electron mass
      REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties
-      INTEGER:: speciesID
+      INTEGER, INTENT(in):: ion
      INTEGER, OPTIONAL, INTENT(in):: electronSecondary
      REAL(8):: effTime
      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
      REAL(8), INTENT(in):: eThreshold
@ -126,7 +137,22 @@ MODULE moduleBoundary
        boundary%n0 = n0 * Vol_ref
        boundary%v0 = v0 / v_ref
        boundary%vTh = DSQRT(kb*T0/m0)/v_ref
-        boundary%species => species(speciesID)%obj
+        boundary%species => species(ion)%obj
        IF (PRESENT(electronSecondary)) THEN
          SELECT TYPE(sp => species(electronSecondary)%obj)
          TYPE IS(speciesCharged)
            boundary%electronSecondary => sp
          CLASS DEFAULT
            CALL criticalError("Species " // sp%name // " chosen for " // &
                               "secondary electron is not a charged species", 'initIonization')
          END SELECT
        ELSE
          boundary%electronSecondary => NULL()
        END IF
        boundary%effectiveTime = effTime / ti_ref
        CALL boundary%crossSection%init(crossSection)
        CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
--- a/src/modules/moduleCollisions.f90
+++ b/src/modules/moduleCollisions.f90
@ -43,7 +43,8 @@ MODULE moduleCollisions
  TYPE, EXTENDS(collisionBinary):: collisionBinaryIonization
    REAL(8):: eThreshold !Minimum energy (non-dimensional units) required for ionization
    REAL(8):: deltaV !Change in velocity due to exchange of eThreshold
-    CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons
+    CLASS(speciesCharged), POINTER:: electron          !Pointer to species considerer as electrons
    CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron
    CONTAINS
      PROCEDURE, PASS:: collide => collideBinaryIonization
@ -241,7 +242,7 @@ MODULE moduleCollisions
    !ELECTRON IMPACT IONIZATION
    !Inits electron impact ionization
-    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron)
+    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron, electronSecondary)
      USE moduleTable
      USE moduleRefParam
      USE moduleConstParam
@ -253,7 +254,8 @@ MODULE moduleCollisions
      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
      REAL(8), INTENT(in):: energyThreshold
      CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
-      INTEGER:: electronIndex
+      CHARACTER(:), ALLOCATABLE, OPTIONAL, INTENT(in):: electronSecondary
      INTEGER:: electronIndex, electronSecondaryIndex
      ALLOCATE(collisionBinaryIonization:: collision)
@ -278,10 +280,27 @@ MODULE moduleCollisions
        CLASS DEFAULT
          CALL criticalError("Species " // sp%name // " chosen for " // &
-                             "secondary electron is not a charged species", 'initBinaryIonization')
+                             "impacting electron is not a charged species", 'initBinaryIonization')
        END SELECT
        IF (PRESENT(electronSecondary)) THEN
          electronSecondaryIndex = speciesName2Index(electronSecondary)
          SELECT TYPE(sp => species(electronSecondaryIndex)%obj)
          TYPE IS(speciesCharged)
            collision%electronSecondary => sp
          CLASS DEFAULT
            CALL criticalError("Species " // sp%name // " chosen for " // &
                               "secondary electron is not a charged species", 'initBinaryIonization')
          END SELECT
        ELSE
          collision%electronSecondary => NULL()
        END IF
        !momentum change per ionization process
        collision%deltaV = sqrt(collision%eThreshold / collision%electron%m)
@ -336,6 +355,12 @@ MODULE moduleCollisions
        !Copy basic information from primary electron
        newElectron = electron
        !If secondary electron species indicates, convert
        IF (ASSOCIATED(self%electronSecondary)) THEN
          newElectron%species => self%electronSecondary
        END IF
        !Secondary electorn gains energy from ionization
        newElectron%v       =  vChange
@ -362,7 +387,7 @@ MODULE moduleCollisions
          CALL sp%ionize(neutral)
        CLASS DEFAULT
-          ! CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
+          CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
          RETURN
        END SELECT
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@ -54,15 +54,16 @@ MODULE moduleInject
    INTEGER:: id
    CHARACTER(:), ALLOCATABLE:: name
    REAL(8):: vMod !Velocity (module)
-    REAL(8):: T(1:3) !Temperature
+    REAL(8):: temperature(1:3) !Temperature
    REAL(8):: n(1:3) !Direction of injection
    LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
    INTEGER:: nParticles !Number of particles to introduce each time step
    CLASS(speciesGeneric), POINTER:: species !Species of injection
    INTEGER:: nEdges
    INTEGER, ALLOCATABLE:: edges(:) !Array with edges
-    REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability
+    INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge
-    REAL(8):: sumWeight
+    REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge
    REAL(8):: surface ! Total surface of injection
    TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
    CONTAINS
      PROCEDURE, PASS:: init         => initInject
@ -75,7 +76,7 @@ MODULE moduleInject
  CONTAINS
    !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
+    SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
      USE moduleMesh
      USE moduleRefParam
      USE moduleConstParam
@ -86,49 +87,29 @@ MODULE moduleInject
      CLASS(injectGeneric), INTENT(inout):: self
      INTEGER, INTENT(in):: i
-      REAL(8), INTENT(in):: v, n(1:3), T(1:3)
+      REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
-      INTEGER, INTENT(in):: sp, physicalSurface
+      INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge
      REAL(8):: tauInject
      REAL(8), INTENT(in):: flow
      CHARACTER(:), ALLOCATABLE, INTENT(in):: units
      INTEGER:: e, et
      INTEGER:: phSurface(1:mesh%numEdges)
      INTEGER:: nVolColl
      REAL(8):: fluxPerStep = 0.D0
-      self%id         =  i
+      self%id          =  i
-      self%vMod       =  v / v_ref
+      self%vMod        =  v / v_ref
-      self%n          =  n / NORM2(n)
+      self%n           =  n / NORM2(n)
-      self%T          =  T / T_ref
+      self%temperature =  temperature / T_ref
      self%species    => species(sp)%obj
      tauInject = tau(self%species%n)
      SELECT CASE(units)
      CASE ("sccm")
        !Standard cubic centimeter per minute
        self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight)
      CASE ("A")
        !Input current in Ampers
        self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight))
      CASE ("part/s")
        !Input current in Ampers
        self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight)
      CASE DEFAULT
        CALL criticalError("No support for units: " // units, 'initInject')
      END SELECT
      !Scale particles for different species steps
      IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
      !Gets the edge elements from which particles are injected
      DO e = 1, mesh%numEdges
        phSurface(e) = mesh%edges(e)%obj%physicalSurface
      END DO
      self%nEdges = COUNT(phSurface == physicalSurface)
-      ALLOCATE(inject(i)%edges(1:self%nEdges))
+      ALLOCATE(self%edges(1:self%nEdges))
      ALLOCATE(self%particlesPerEdge(1:self%nEdges))
      ALLOCATE(self%weightPerEdge(1:self%nEdges))
      et = 0
      DO e=1, mesh%numEdges
        IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
@ -160,15 +141,82 @@ MODULE moduleInject
      END DO
-      !Calculates cumulative probability
+      !Calculates total area
-      ALLOCATE(self%cumWeight(1:self%nEdges))
+      self%surface = 0.D0
-      et = 1
+      DO et = 1, self%nEdges
-      self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight
+        self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface
      DO et = 2, self%nEdges
        self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1)
      END DO
-      self%sumWeight = self%cumWeight(self%nEdges)
+
      ! Information about species and flux
      self%species => species(sp)%obj
      tauInject    = tau(self%species%n)
      ! Convert units
      SELECT CASE(units)
      CASE ("sccm")
        !Standard cubic centimeter per minute
        fluxPerStep = flow*sccm2atomPerS
      CASE ("A")
        !Current in Ampers
        SELECT TYPE(sp => self%species)
        CLASS IS(speciesCharged)
          fluxPerStep = flow/(qe*abs(sp%q))
        CLASS DEFAULT
          call criticalError('Attempted to assign a flux in "A" to a species without charge.', 'initInject')
        END SELECT
      CASE ("Am2")
        !Input current in Ampers per square meter
        SELECT TYPE(sp => self%species)
        CLASS IS(speciesCharged)
          fluxPerStep = flow*self%surface*L_ref**2/(qe*abs(sp%q))
        CLASS DEFAULT
          call criticalError('Attempted to assign a flux in "Am2" to a species without charge.', 'initInject')
        END SELECT
      CASE ("part/s")
        !Input current in Ampers
        fluxPerStep = flow
      CASE DEFAULT
        CALL criticalError("No support for units: " // units, 'initInject')
      END SELECT
      fluxPerStep = fluxPerStep * tauInject * ti_ref / self%surface
      !Assign particles per edge
      IF (particlesPerEdge > 0) THEN
        ! Particles per edge defined by the user
        self%particlesPerEdge = particlesPerEdge
        DO et = 1, self%nEdges
          self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge)
        END DO
        self%nParticles = SUM(self%particlesPerEdge)
      ELSE
        ! No particles assigned per edge, use the species weight
        self%weightPerEdge = self%species%weight
        DO et = 1, self%nEdges
          self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight))
        END DO
        self%nParticles = SUM(self%particlesPerEdge)
        !Rescale weight to match flux
        self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles))
      END IF
      !Scale particles for different species steps
      IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
    END SUBROUTINE initInject
@ -203,23 +251,23 @@ MODULE moduleInject
    END SUBROUTINE doInjects
-    SUBROUTINE initVelDistMaxwellian(velDist, T, m)
+    SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
      IMPLICIT NONE
      CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: T, m
+      REAL(8), INTENT(in):: temperature, m
-      velDist = velDistMaxwellian(vTh = DSQRT(T/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
    END SUBROUTINE initVelDistMaxwellian
-    SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m)
+    SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m)
      IMPLICIT NONE
      CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: T, m
+      REAL(8), INTENT(in):: temperature, m
-      velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m))
+      velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
    END SUBROUTINE initVelDistHalfMaxwellian
@ -241,7 +289,7 @@ MODULE moduleInject
      REAL(8):: v
      v = 0.D0
-      v = self%vTh*randomMaxwellian()
+      v = self%vTh*randomMaxwellian()/sqrt(2.d0)
    END FUNCTION randomVelMaxwellian
@ -254,10 +302,7 @@ MODULE moduleInject
      REAL(8):: v
      v = 0.D0
-      DO WHILE (v <= 0.D0)
+      v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0)
        v = self%vTh*randomMaxwellian()
      END DO
    END FUNCTION randomVelHalfMaxwellian
@ -282,9 +327,8 @@ MODULE moduleInject
      IMPLICIT NONE
      CLASS(injectGeneric), INTENT(in):: self
-      INTEGER:: randomX
+      INTEGER, SAVE:: nMin
-      INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
+      INTEGER:: i, e
      INTEGER:: i
      INTEGER:: n, sp
      CLASS(meshEdge), POINTER:: randomEdge
      REAL(8):: direction(1:3)
@ -299,59 +343,65 @@ MODULE moduleInject
        END IF
      END DO
      nMin = nMin + 1
      nMax = nMin + self%nParticles - 1
      !Assign weight to particle.
      partInj(nMin:nMax)%weight = self%species%weight
      !Particle is considered to be outside the domain
      partInj(nMin:nMax)%n_in = .FALSE.
      !$OMP END SINGLE
      !$OMP DO
-      DO n = nMin, nMax
+      DO e = 1, self%nEdges
-        randomX = randomWeighted(self%cumWeight, self%sumWeight)
+        ! Select edge for injection
        randomEdge => mesh%edges(self%edges(e))%obj
        ! Inject particles in edge
        DO i = 1, self%particlesPerEdge(e)
          ! Index in the global partInj array
          n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i
          !Particle is considered to be outside the domain
          partInj(n)%n_in = .FALSE.
          !Random position in edge
          partInj(n)%r = randomEdge%randPos()
          !Assign weight to particle.
          partInj(n)%weight = self%weightPerEdge(e)
          !Volume associated to the edge:
          IF (ASSOCIATED(randomEdge%e1)) THEN
            partInj(n)%cell = randomEdge%e1%n
-        randomEdge => mesh%edges(self%edges(randomX))%obj
+          ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
-        !Random position in edge
+            partInj(n)%cell = randomEdge%e2%n
        partInj(n)%r = randomEdge%randPos()
        !Volume associated to the edge:
        IF (ASSOCIATED(randomEdge%e1)) THEN
          partInj(n)%cell = randomEdge%e1%n
-        ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
+          ELSE
-          partInj(n)%cell = randomEdge%e2%n
+            CALL criticalError("No Volume associated to edge", 'addParticles')
-        ELSE
+          END IF
-          CALL criticalError("No Volume associated to edge", 'addParticles') 
+          partInj(n)%cellColl = randomEdge%eColl%n
          sp = self%species%n
-        END IF
+          !Assign particle type
-        partInj(n)%cellColl = randomEdge%eColl%n
+          partInj(n)%species => self%species
        sp = self%species%n
-        !Assign particle type
+          if (all(self%n == 0.D0)) then
-        partInj(n)%species => self%species
+            direction =  randomEdge%normal
-        direction = self%n
+          else
            direction = self%n
-        partInj(n)%v = 0.D0
+          end if
-        partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
+          partInj(n)%v = 0.D0
                                                 self%v(2)%obj%randomVel(), &
                                                 self%v(3)%obj%randomVel() /)
-        !If velocity is not in the right direction, invert it
+          do while(dot_product(partInj(n)%v, direction) <= 0.d0)
-        IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN
+            partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
-          partInj(n)%v = - partInj(n)%v
+                                                     self%v(2)%obj%randomVel(), &
                                                     self%v(3)%obj%randomVel() /)
          end do
        END IF
-        !Obtain natural coordinates of particle in cell
+          !Obtain natural coordinates of particle in cell
-        partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
+          partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
-        !Push new particle with the minimum time step
+          !Push new particle with the minimum time step
-        CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
+          CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
-        !Assign cell to new particle
+          !Assign cell to new particle
-        CALL solver%updateParticleCell(partInj(n))
+          CALL solver%updateParticleCell(partInj(n))
        END DO
      END DO
      !$OMP END DO
--- a/src/modules/moduleProbe.f90
+++ b/src/modules/moduleProbe.f90
@ -27,7 +27,7 @@ MODULE moduleProbe
  CONTAINS
    !Functions for probeDistFunc type
-    SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep)
+    SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep)
      USE moduleCaseParam
      USE moduleRefParam
      USE moduleSpecies
@ -41,7 +41,7 @@ MODULE moduleProbe
      REAL(8), INTENT(in):: r(1:3)
      REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2)
      INTEGER, INTENT(in):: points(1:3)
-      REAL(8), INTENT(in):: timeStep
+      REAL(8), INTENT(in):: everyTimeStep
      INTEGER:: sp, i
      REAL(8):: dv(1:3)
@ -91,17 +91,17 @@ MODULE moduleProbe
                      1:self%nv(3)))
      !Number of iterations between output
-      IF (timeStep == 0.D0) THEN
+      IF (everyTimeStep == 0.D0) THEN
        self%every = 1
      ELSE
-        self%every = NINT(timeStep/ tauMin / ti_ref)
+        self%every = NINT(everyTimeStep/ tauMin / ti_ref)
      END IF
      !Maximum radius
      !TODO: Make this an input parameter
-      self%maxR = 1.D0
+      self%maxR = 1.D-2/L_ref
      !Init the probe lock
      CALL OMP_INIT_LOCK(self%lock)
@ -148,7 +148,7 @@ MODULE moduleProbe
        deltaR = NORM2(self%r - part%r)
        !Only include particle if it is inside the maximum radius
-        IF (deltaR < self%maxR) THEN
+        ! IF (deltaR < self%maxR) THEN
          !find lower index for all dimensions
          CALL self%findLowerIndex(part%v, i, j, k, inside)
@ -162,40 +162,40 @@ MODULE moduleProbe
            fk  = self%vk(k+1) - part%v(3)
            fk1 = part%v(3)    - self%vk(k)
-            ! weight = part%weight * DEXP(deltaR/self%maxR)
+            weight = part%weight * DEXP(-deltaR/self%maxR)
-            weight = part%weight
+            ! weight = part%weight
            !Lock the probe
            CALL OMP_SET_LOCK(self%lock)
            !Assign particle weight to distribution function
-            self%f(i  , j  , k  ) = fi  * fj  * fk  * weight
+            self%f(i  , j  , k  ) = self%f(i  , j  , k  ) + fi  * fj  * fk  * weight
-            self%f(i+1, j  , k  ) = fi1 * fj  * fk  * weight
+            self%f(i+1, j  , k  ) = self%f(i+1, j  , k  ) + fi1 * fj  * fk  * weight
-            self%f(i  , j+1, k  ) = fi  * fj1 * fk  * weight
+            self%f(i  , j+1, k  ) = self%f(i  , j+1, k  ) + fi  * fj1 * fk  * weight
-            self%f(i+1, j+1, k  ) = fi1 * fj1 * fk  * weight
+            self%f(i+1, j+1, k  ) = self%f(i+1, j+1, k  ) + fi1 * fj1 * fk  * weight
-            self%f(i  , j  , k+1) = fi  * fj  * fk1 * weight
+            self%f(i  , j  , k+1) = self%f(i  , j  , k+1) + fi  * fj  * fk1 * weight
-            self%f(i+1, j  , k+1) = fi1 * fj  * fk1 * weight
+            self%f(i+1, j  , k+1) = self%f(i+1, j  , k+1) + fi1 * fj  * fk1 * weight
-            self%f(i  , j+1, k+1) = fi  * fj1 * fk1 * weight
+            self%f(i  , j+1, k+1) = self%f(i  , j+1, k+1) + fi  * fj1 * fk1 * weight
-            self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * weight
+            self%f(i+1, j+1, k+1) = self%f(i+1, j+1, k+1) + fi1 * fj1 * fk1 * weight
            !Unlock the probe
            CALL OMP_UNSET_LOCK(self%lock)
          END IF
-        END IF
+        ! END IF
      END IF
    END SUBROUTINE calculate
-    SUBROUTINE output(self, t)
+    SUBROUTINE output(self)
      USE moduleOutput
      USE moduleRefParam
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      CLASS(probeDistFunc), INTENT(inout):: self
      INTEGER, INTENT(in):: t
      CHARACTER (LEN=iterationDigits):: tstring
      CHARACTER (LEN=3):: pstring
      CHARACTER(:), ALLOCATABLE:: filename
@ -204,14 +204,14 @@ MODULE moduleProbe
      !Divide by the velocity cube volume
      self%f = self%f * self%dvInv
-      WRITE(tstring, iterationFormat) t
+      WRITE(tstring, iterationFormat) timeStep
      WRITE(pstring, "(I3.3)") self%id
      fileName='Probe_' // tstring// '_f_' // pstring // '.dat'
      WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
      OPEN (10, file = path // folder // '/' //  fileName)
      WRITE(10, "(A1, 1X, A)")               "# ", self%species%name
      WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m"
-      WRITE(10, "(A6, ES15.6E3, A2)")    "# t = ", REAL(t)*tauMin*ti_ref, " s"
+      WRITE(10, "(A6, ES15.6E3, A2)")    "# t = ", REAL(timeStep)*tauMin*ti_ref, " s"
      WRITE(10, "(A1, A19, 3(A20))")     "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f"
      DO i = 1, self%nv(1)
        DO j = 1, self%nv(2)
@ -252,15 +252,15 @@ MODULE moduleProbe
    END SUBROUTINE doProbes
-    SUBROUTINE outputProbes(t)
+    SUBROUTINE outputProbes()
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      INTEGER, INTENT(in):: t
      INTEGER:: i
      DO i = 1, nProbes
        IF (probe(i)%update) THEN
-          CALL probe(i)%output(t)
+          CALL probe(i)%output()
        END IF
@ -268,15 +268,15 @@ MODULE moduleProbe
    END SUBROUTINE outputProbes
-    SUBROUTINE resetProbes(t)
+    SUBROUTINE resetProbes()
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      INTEGER, INTENT(in):: t
      INTEGER:: i
      DO i = 1, nProbes
        probe(i)%f = 0.D0
-        probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0
+        probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0
      END DO
--- a/src/modules/output/moduleOutput.f90
+++ b/src/modules/output/moduleOutput.f90
@ -22,7 +22,6 @@ MODULE moduleOutput
  !Type for EM data in node
  TYPE emNode
    CHARACTER(:), ALLOCATABLE:: type
    REAL(8):: phi
    REAL(8):: B(1:3)
@ -160,12 +159,12 @@ MODULE moduleOutput
    END SUBROUTINE calculateOutput
-    SUBROUTINE printTime(t, first)
+    SUBROUTINE printTime(first)
      USE moduleSpecies
      USE moduleCompTime
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
      INTEGER, INTENT(in):: t
      LOGICAL, INTENT(in), OPTIONAL:: first
      CHARACTER(:), ALLOCATABLE:: fileName
@ -187,7 +186,7 @@ MODULE moduleOutput
        OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
-        WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
+        WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
        CLOSE(20)
--- a/src/modules/solver/electromagnetic/moduleEM.f90
+++ b/src/modules/solver/electromagnetic/moduleEM.f90
@ -1,55 +1,202 @@
 !Module to solve the electromagnetic field
 MODULE moduleEM
  USE moduleMesh
  USE moduleTable
  IMPLICIT NONE
-  TYPE:: boundaryEM
+  ! Generic type for electromagnetic boundary conditions
-    CHARACTER(:), ALLOCATABLE:: typeEM
+  TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric
-    INTEGER:: physicalSurface
+    INTEGER:: nNodes
    TYPE(meshNodePointer), ALLOCATABLE:: nodes(:)
    CONTAINS
      PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply
  END TYPE boundaryEMGeneric
  ABSTRACT INTERFACE
    ! Apply boundary condition to the load vector for the Poission equation
    SUBROUTINE applyEM_interface(self, vectorF)
      IMPORT boundaryEMGeneric
      CLASS(boundaryEMGeneric), INTENT(in):: self
      REAL(8), INTENT(inout):: vectorF(:)
    END SUBROUTINE applyEM_interface
  END INTERFACE
  TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet
    REAL(8):: potential
    CONTAINS
-      PROCEDURE, PASS:: apply
+      ! boundaryEMGeneric DEFERRED PROCEDURES
      PROCEDURE, PASS:: apply => applyDirichlet
-  END TYPE boundaryEM
+  END TYPE boundaryEMDirichlet
  TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
    REAL(8):: potential
    TYPE(table1D):: temporalProfile
    CONTAINS
      ! boundaryEMGeneric DEFERRED PROCEDURES
      PROCEDURE, PASS:: apply => applyDirichletTime
  END TYPE boundaryEMDirichletTime
  ! Container for boundary conditions
  TYPE:: boundaryEMCont
    CLASS(boundaryEMGeneric), ALLOCATABLE:: obj
  END TYPE boundaryEMCont
  INTEGER:: nBoundaryEM
-  TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)
+  TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:)
  !Information of charge and reference parameters for rho vector
  REAL(8), ALLOCATABLE:: qSpecies(:)
  CONTAINS
-    !Apply boundary conditions to the K matrix for Poisson's equation
+    SUBROUTINE findNodes(self, physicalSurface)
    SUBROUTINE apply(self, edge)
      USE moduleMesh
      IMPLICIT NONE
-      CLASS(boundaryEM), INTENT(in):: self
+      CLASS(boundaryEMGeneric), INTENT(inout):: self
-      CLASS(meshEdge):: edge
+      INTEGER, INTENT(in):: physicalSurface
-      INTEGER:: nNodes
+      CLASS(meshEdge), POINTER:: edge
-      INTEGER, ALLOCATABLE:: nodes(:)
+      INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:)
-      INTEGER:: n
+      INTEGER:: nNodes, nodesNew
      INTEGER:: e, n
-      nNodes = edge%nNodes
+      !Temporal array to hold nodes
-      nodes = edge%getNodes(nNodes)
+      ALLOCATE(nodes(0))
-      DO n = 1, nNodes
+      ! Loop thorugh the edges and identify those that are part of the boundary
-        SELECT CASE(self%typeEM)
+      DO e = 1, mesh%numEdges
-        CASE ("dirichlet")
+        edge => mesh%edges(e)%obj
-          mesh%K(nodes(n), :)        = 0.D0
+        IF (edge%physicalSurface == physicalSurface) THEN
-          mesh%K(nodes(n), nodes(n)) = 1.D0
+          ! Edge is of the right boundary index
-        
+          ! Get nodes in the edge
-          mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
+          nNodes = edge%nNodes
-          mesh%nodes(nodes(n))%obj%emData%phi  = self%potential
+          nodesEdge = edge%getNodes(nNodes)
          ! Collect all nodes that are not already in the temporal array
          DO n = 1, nNodes
            IF (ANY(nodes == nodesEdge(n))) THEN
              ! Node already in array, skip
              CYCLE
-        END SELECT
+            ELSE
              ! If not, add element to array of nodes
              nodes = [nodes, nodesEdge(n)]
            END IF
          END DO
        END IF
      END DO
-    END SUBROUTINE
+      ! Point boundary to nodes
      nNodes = SIZE(nodes)
      ALLOCATE(self%nodes(nNodes))
      self%nNodes = nNodes
      DO n = 1, nNodes
        self%nodes(n)%obj => mesh%nodes(nodes(n))%obj
      END DO
    END SUBROUTINE findNodes
    ! Initialize Dirichlet boundary condition
    SUBROUTINE initDirichlet(self, physicalSurface, potential)
      USE moduleRefParam, ONLY: Volt_ref
      IMPLICIT NONE
      CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
      INTEGER, INTENT(in):: physicalSurface
      REAL(8), INTENT(in):: potential
      ! Allocate boundary edge
      ALLOCATE(boundaryEMDirichlet:: self)
      SELECT TYPE(self)
      TYPE IS(boundaryEMDirichlet)
        self%potential = potential / Volt_ref
        CALL findNodes(self, physicalSurface)
      END SELECT
    END SUBROUTINE initDirichlet
    ! Initialize Dirichlet boundary condition
    SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile)
      USE moduleRefParam, ONLY: Volt_ref, ti_ref
      IMPLICIT NONE
      CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
      INTEGER, INTENT(in):: physicalSurface
      REAL(8), INTENT(in):: potential
      CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile
      ! Allocate boundary edge
      ALLOCATE(boundaryEMDirichletTime:: self)
      SELECT TYPE(self)
      TYPE IS(boundaryEMDirichletTime)
        self%potential = potential / Volt_ref
        CALL findNodes(self, physicalSurface)
        CALL self%temporalProfile%init(temporalProfile)
        CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0)
      END SELECT
    END SUBROUTINE initDirichletTime
    !Apply Dirichlet boundary condition to the poisson equation
    SUBROUTINE applyDirichlet(self, vectorF)
      USE moduleMesh
      IMPLICIT NONE
      CLASS(boundaryEMDirichlet), INTENT(in):: self
      REAL(8), INTENT(inout):: vectorF(:)
      INTEGER:: n, ni
      DO n = 1, self%nNodes
        self%nodes(n)%obj%emData%phi  = self%potential
        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
      END DO
    END SUBROUTINE applyDirichlet
    !Apply Dirichlet boundary condition with time temporal profile
    SUBROUTINE applyDirichletTime(self, vectorF)
      USE moduleMesh
      USE moduleCaseParam, ONLY: timeStep, tauMin
      IMPLICIT NONE
      CLASS(boundaryEMDirichletTime), INTENT(in):: self
      REAL(8), INTENT(inout):: vectorF(:)
      REAL(8):: timeFactor
      INTEGER:: n, ni
      timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)
      DO n = 1, self%nNodes
        self%nodes(n)%obj%emData%phi  = self%potential * timeFactor
        vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
      END DO
    END SUBROUTINE applyDirichletTime
    !Assemble the source vector based on the charge density to solve Poisson's equation
-    SUBROUTINE assembleSourceVector(vectorF)
+    SUBROUTINE assembleSourceVector(vectorF, n_e)
      USE moduleMesh
      USE moduleRefParam
      IMPLICIT NONE
@ -58,8 +205,9 @@ MODULE moduleEM
      REAL(8), ALLOCATABLE:: localF(:)
      INTEGER, ALLOCATABLE:: nodes(:)
      REAL(8), ALLOCATABLE:: rho(:)
      REAL(8), INTENT(in), OPTIONAL:: n_e(1:mesh%numNodes)
      INTEGER:: nNodes
-      INTEGER:: e, i, ni
+      INTEGER:: e, i, ni, b
      CLASS(meshNode), POINTER:: node
      !$OMP SINGLE
@ -77,6 +225,10 @@ MODULE moduleEM
          ni = nodes(i)
          node => mesh%nodes(ni)%obj
          rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref))
          IF (PRESENT(n_e)) THEN
            rho(i) = rho(i) - n_e(i)
          END IF
        END DO
@ -97,18 +249,12 @@ MODULE moduleEM
      !$OMP END DO
      !Apply boundary conditions
-      !$OMP DO
+      !$OMP SINGLE
-      DO i = 1, mesh%numNodes
+      do b = 1, nBoundaryEM
-        node => mesh%nodes(i)%obj
+        call boundaryEM(b)%obj%apply(vectorF)
-        SELECT CASE(node%emData%type)
+      end do
-        CASE ("dirichlet")
+      !$OMP END SINGLE
          vectorF(i) = node%emData%phi
        END SELECT
      END DO
      !$OMP END DO
    END SUBROUTINE assembleSourceVector
@ -156,4 +302,86 @@ MODULE moduleEM
    END SUBROUTINE solveElecField
    FUNCTION BoltzmannElectron(phi, n) RESULT(n_e)
      USE moduleRefParam
      USE moduleConstParam
      IMPLICIT NONE
      INTEGER, INTENT(in):: n
      REAL(8), INTENT(in):: phi(1:n)
      REAL(8):: n_e(1:n)
      REAL(8):: n_e0 = 1.0D16, phi_0 = -500.0D0, T_e = 11604.0
      INTEGER:: i
      n_e = n_e0 / n_ref * exp(qe * (phi*Volt_ref - phi_0) / (kb * T_e))
      RETURN
    END FUNCTION BoltzmannElectron
    SUBROUTINE solveElecFieldBoltzmann()
      USE moduleMesh
      USE moduleErrors
      IMPLICIT NONE
      INTEGER, SAVE:: INFO
      INTEGER:: n
      REAL(8), ALLOCATABLE, SAVE:: tempF(:)
      REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:), phi(:)
      INTEGER:: k
      EXTERNAL:: dgetrs
      !$OMP SINGLE
      ALLOCATE(tempF(1:mesh%numNodes))
      ALLOCATE(n_e(1:mesh%numNodes))
      ALLOCATE(phi_old(1:mesh%numNodes))
      ALLOCATE(phi(1:mesh%numNodes))
      !$OMP END SINGLE
      !$OMP DO
      DO n = 1, mesh%numNodes
        phi_old(n) = mesh%nodes(n)%obj%emData%phi
      END DO
      !$OMP END DO
      !$OMP SINGLE
      DO k = 1, 100
        n_e = BoltzmannElectron(phi_old, mesh%numNodes)
        CALL assembleSourceVector(tempF, n_e)
        CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
                                mesh%IPIV, tempF,  mesh%numNodes, info)
        phi = tempF
        PRINT *, MAXVAL(n_e), MINVAL(n_e)
        PRINT *, MAXVAL(phi), MINVAL(phi)
        PRINT*, k, "diff = ", MAXVAL(ABS(phi - phi_old))
        phi_old = phi
      END DO
      !$OMP END SINGLE
      IF (info == 0) THEN
        !Suscessful resolution of Poission equation
        !$OMP DO
        DO n = 1, mesh%numNodes
          mesh%nodes(n)%obj%emData%phi = phi_old(n)
        END DO
        !$OMP END DO
      ELSE
        !$OMP SINGLE
        CALL criticalError('Poisson equation failed', 'solveElecFieldBoltzmann')
        !$OMP END SINGLE
      END IF
      !$OMP SINGLE
      DEALLOCATE(tempF, n_e, phi_old, phi)
      !$OMP END SINGLE
    END SUBROUTINE solveElecFieldBoltzmann
 END MODULE moduleEM
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@ -138,6 +138,9 @@ MODULE moduleSolver
      CASE('Electrostatic','ConstantB')
        self%solveEM => solveElecField
      CASE('ElectrostaticBoltzmann')
        self%solveEM => solveElecFieldBoltzmann
      END SELECT
    END SUBROUTINE initEM
@ -491,52 +494,46 @@ MODULE moduleSolver
    END SUBROUTINE updateParticleCell
    !Update the information about if a species needs to be moved this iteration
-    SUBROUTINE updatePushSpecies(self, t)
+    SUBROUTINE updatePushSpecies(self)
      USE moduleSpecies
      USE moduleCaseparam, ONLY: timeStep
      IMPLICIT NONE
      CLASS(solverGeneric), INTENT(inout):: self
      INTEGER, INTENT(in):: t
      INTEGER:: s
      DO s=1, nSpecies
-        self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0
+        self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0
      END DO
    END SUBROUTINE updatePushSpecies
    !Output the different data and information
-    SUBROUTINE doOutput(t)
+    SUBROUTINE doOutput()
      USE moduleMesh
      USE moduleOutput
      USE moduleSpecies
      USE moduleCompTime
      USE moduleProbe
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE
-      INTEGER, INTENT(in):: t
+      CALL outputProbes()
      IF (t == tInitial) THEN
        CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpack_commit.txt')
      END IF
      CALL outputProbes(t)
      counterOutput = counterOutput + 1
      IF (counterOutput >= triggerOutput .OR. &
-          t == tFinal .OR. t == tInitial) THEN
+          timeStep == tFinal .OR. timeStep == tInitial) THEN
        !Resets output counter
        counterOutput=0
-        CALL mesh%printOutput(t)
+        CALL mesh%printOutput()
-        IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
+        IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl()
-        CALL mesh%printEM(t)
+        CALL mesh%printEM()
-        WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal
+        WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal
        WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
-        IF (t == 0) THEN
+        IF (timeStep == 0) THEN
          WRITE(*, "(5X,A21,F8.1,A2)")   "     init time: ",    1.D3*tStep, "ms"
        ELSE
@ -554,34 +551,32 @@ MODULE moduleSolver
      counterCPUTime = counterCPUTime + 1
      IF (counterCPUTime >= triggerCPUTime .OR. &
-          t == tFinal .OR. t == tInitial) THEN
+          timeStep == tFinal .OR. timeStep == tInitial) THEN
        !Reset CPU Time counter
        counterCPUTime = 0
-        CALL printTime(t, t == 0)
+        CALL printTime(timeStep == 0)
      END IF
      !Output average values
-      IF (useAverage .AND. t == tFinal) THEN
+      IF (useAverage .AND. timeStep == tFinal) THEN
        CALL mesh%printAverage()
      END IF
    END SUBROUTINE doOutput
-    SUBROUTINE doAverage(t)
+    SUBROUTINE doAverage()
      USE moduleAverage
      USE moduleMesh
      IMPLICIT NONE
      INTEGER, INTENT(in):: t
      INTEGER:: tAverage, n
      IF (useAverage) THEN
-        tAverage = t - tAverageStart
+        tAverage = timeStep - tAverageStart
        IF (tAverage == 1) THEN
          !First iteration in which average scheme is used
Author	SHA1	Message	Date
Jorge Gonzalez	796176b4b4	Merge pull request 'issue/edgesVolume' (#59 ) from issue/edgesVolume into development Reviewed-on: https://plasmalab.aero.upm.es/git/git/JGonzalez/fpakc/pulls/59	2026-02-26 19:21:52 +01:00
Jorge Gonzalez - Galactica	7208b2e0da	Merge remote-tracking branch 'origin/issue/Maxwellian' into issue/edgesVolume	2026-02-04 19:59:59 +01:00
JGonzalez	23bba31005	Fixed an issue with the random Maxwellian subroutines	2026-02-04 13:52:24 +01:00
JGonzalez	159f2e7620	Trying partial reflection	2026-02-03 10:15:12 +01:00
Jorge Gonzalez	34b5bdbb52	Merge pull request 'Wrong surface in edges' (#58 ) from issue/edgesVolume into development Reviewed-on: https://plasmalab.aero.upm.es/git/git/JGonzalez/fpakc/pulls/58	2026-02-02 14:34:47 +01:00
Jorge Gonzalez - Galactica	8d5cb6a516	Wrong surface in edges Some edges were calculating the area surface incorrectly, which was leading to wrong densities.	2026-02-02 14:31:58 +01:00
Jorge Gonzalez	4747673d0a	Bones of boundary	2026-01-28 14:56:38 +01:00
Jorge Gonzalez	6a901a0cac	Merge pull request 'feature/meshText' (#57 ) from feature/meshText into development Reviewed-on: https://plasmalab.aero.upm.es/git/git/JGonzalez/fpakc/pulls/57	2026-01-23 15:27:27 +01:00
JGonzalez	56a17300f3	Change in manual	2026-01-23 15:26:23 +01:00
JGonzalez	0c0017b84a	txs auto checkin	2026-01-23 15:21:26 +01:00
JGonzalez	1b1a574edc	Read initial Now we can read an initial condition as we do with other formats. Only documentation left.	2026-01-23 15:04:39 +01:00
JGonzalez	d34e72d2e8	Collisions done Thanks to the advance statement the output of collisions is done.	2026-01-23 14:51:49 +01:00
JGonzalez	46b9f263ea	Average done Now we can output averages in this format.	2026-01-23 14:26:39 +01:00
JGonzalez	8c8c6409f6	Output for EM The output for the EM field properties is done.	2026-01-23 12:16:52 +01:00
JGonzalez	27158c7c1d	Output for species First output done. Now it is mostly Copy Paste. There are so many things I need to reform in other parts... Learning new things is awful.	2026-01-22 18:44:25 +01:00
JGonzalez	3083e20ff7	Subroutine for reading mesh It is now complete and working. I also added some minor improvements to moduleMesh.f90 using 'associate' to practice. I noticed there are a lot of things (allocating K, for example) that are common for all meshes and should be moved to a general module.	2026-01-22 14:06:20 +01:00
JGonzalez	7d4f4b98c3	Implementing input subroutines	2026-01-19 15:37:31 +01:00
JGonzalez	9f9bacca7c	Skeleton implementation of text mesh for simple 1D cases.	2026-01-19 14:48:06 +01:00
Jorge Gonzalez	d5bcf8ce50	Merge branch 'issue/triangles' into 'development' Fixed issue with volume in triangles See merge request JorgeGonz/fpakc!56	2025-10-11 12:01:28 +00:00
JGonzalez	55e062a9ef	Fixed issue with volume in triangles The right value in 2D Cartesian is used for calculating the volume.	2025-10-11 14:00:16 +02:00
Jorge Gonzalez	c9b551458d	Merge branch 'citationFile' into 'development' Add citation file See merge request JorgeGonz/fpakc!55	2025-09-24 17:34:24 +00:00
JGonzalez	9d789ef7a6	Add citation file I think it follows the standard.	2025-09-24 19:33:25 +02:00
Jorge Gonzalez	b1f21a2c02	Merge branch 'IEPC2025' into 'development' Merge IEPC2025 See merge request JorgeGonz/fpakc!54	2025-09-23 16:42:04 +00:00
JGonzalez	dff9a87f0d	Implemenint injecting particles without direction I was almost sure this was implemented in the past, but it was not working. Now, if n = 0 or if n is not provided, particles are injected with the normal to the surface.	2025-08-08 19:27:27 +02:00
JGonzalez	5166a650d2	Data for Kr, just to test	2025-08-06 10:59:03 +02:00
JGonzalez	3d7b1ce476	Fixed! So it seems that rectangles and triangles are now working properly. I have also checked the phy2log routine, now it seems a bit more complicated, but it is much clearer. Maybe in the future is worth rethinking to improve speed (specially for quad elements)	2025-08-03 22:14:19 +02:00
JGonzalez	102fd013f3	Issue with triangles fixed Now they give the right electric field. I have to change 2DCyl. However, there was some insonsistency between the change of coordinates in phy2log and the Jacobian for the K matrix. I fixed it putting a transpose() in phy2log, but I don't like that solution. I need to review the basic procedure of phy2log.	2025-08-03 20:46:12 +02:00
JGonzalez	7e193b9fa8	Minor changes, no improvement made yet	2025-08-03 15:32:55 +02:00
Jorge Gonzalez - Galactica	d86c8480f3	Fixed issue with some velocities. Still, at some point I need to rething all the injection thing.	2025-08-02 16:51:00 +02:00
JGonzalez	4b040c35c3	Fixes to random variables After reading some works and reviewing what I had, I've done some corrections to how the randomb velicities in Maxwellian distributions are calculated. These should be correct now.	2025-08-02 13:25:48 +02:00
Jorge Gonzalez - Galactica	78cb9a2453	No reflection of particles at injection, that should be a boundary condition	2025-08-02 10:31:06 +02:00
JGonzalez	d1e73297eb	Adjust flux when no particlesPerEdge is used This does not affects the cases of the IEPC, but I am also doing other stuff.	2025-07-27 18:17:01 +02:00
JGonzalez	8e531ede08	Vol_ref was the right answer	2025-07-27 17:16:57 +02:00
JGonzalez	76c5af89b2	Merge branch 'IEPC2025' of gitlab.com:JorgeGonz/fpakc into IEPC2025	2025-07-27 17:15:36 +02:00
JGonzalez	7f73b69dc2	Fix injection Half-Maxwellian distribution should inject particles correctly	2025-07-27 17:14:38 +02:00
Jorge Gonzalez - Galactica	69215ef66d	Change in calculating ionization I don't know why I normalizing density n_0 by Vol_ref and not n_ref	2025-07-22 19:52:39 +02:00
JGonzalez	a2f9957f32	I am dumb The Poisson equation was not working because I didn't finish implementing the new type of BCs. Dirichlet is probably untested. I should stop doing shitty developments and no testing.	2025-07-18 16:31:52 +02:00
JGonzalez	d28dd16c2e	Average fix and data for Xe	2025-07-17 18:34:11 +02:00
Jorge Gonzalez	221de46734	Merge branch 'development' into feature/BoltzmannElectrons	2024-10-13 14:54:34 +02:00
Jorge Gonzalez	2268a97e23	Merge branch 'mergeDevleopment' into 'development' Issue with injecting current See merge request JorgeGonz/fpakc!53	2024-10-13 12:41:03 +00:00
Jorge Gonzalez	3d8d912722	Merge branch 'development' of gitlab.com:JorgeGonz/fpakc into development	2024-10-13 13:35:34 +02:00
Jorge Gonzalez	2af10acd70	Issue with injecting current Values were not right in 1D geometry. Fixed.	2024-10-13 13:32:57 +02:00
Jorge Gonzalez	98ee3e9c9c	Still not working, but will be fixed I have the solution in the plasma expansion code, but I need to do other stuff.	2024-09-30 17:06:25 +02:00
Jorge Gonzalez	0679fa58bc	Merge branch 'feature/temporalDirichlet' into 'development' Time variable Dirichlet condition See merge request JorgeGonz/fpakc!52	2024-07-13 11:01:10 +00:00
JGonzalez	6f185c4188	Organizing things Move the array of pointers to the nodes to moduleMesh.	2024-07-13 12:35:42 +02:00
JGonzalez	e4dfba45f8	Manual updated New dirichletTime condition is documented.	2024-07-13 12:13:39 +02:00
JGonzalez	2d4b405fb1	Functionality added Now we have a new boundary condition that can change the value of the potential in a surface based on a file.	2024-07-13 12:06:41 +02:00
JGonzalez	10dee05922	NOT WORKING: Compilation okay, but not Dirichlet BC The code compiles but the right BC is not being applied to the vectorF. I'll check this tomorrow.	2024-07-12 23:30:35 +02:00
JGonzalez	ac27725940	Big one... I should've commited before, but I wanted to make things compile. The big change is that I've added a global time step so the parameter does not need to be passed in each function. This is useful as we are moving towards using time profiles for boundary conditions and injection of particles (not in this branch, but in the future and the procedure will be quite similar)	2024-07-12 23:08:19 +02:00
JGonzalez	49025a6965	Starting changes Planning the new way to do BC in the EM field solver. Probably I have to change how things are read, but I don't think this is going to affect the input file.	2024-07-12 19:21:00 +02:00
JGonzalez	f0a27c0529	More comments So if the source vector is being updated every time step, it might be "easy" to implement these things.	2024-07-12 13:17:02 +02:00
JGonzalez	abedb79b16	Some comments Just some comments on how I am going to make the desired changes (have a Dirichlet boundary condition for the electric potential that changes with time). This might be a good opportunity to rework the boundary conditions in the electrostatic field and include other things like a Newmann boundary condition. We will see.	2024-07-12 11:02:26 +02:00
Jorge Gonzalez	065bb1d13e	Merge branch 'issue/injection2DCyl' into 'development' Rework of injection of particles with a special focus in 2DCyl to ensure an homogeneous distribution. See merge request JorgeGonz/fpakc!51	2024-07-11 16:51:42 +00:00
Jorge Gonzalez	fb9d1ad837	Update moduleMesh2DCyl.f90 This is not needed.	2024-07-11 16:44:29 +00:00
JGonzalez	152ae6b097	No longer needed The correction in the node volume is no longer needed as now things are being calculated right with the last change. Still, at some point I should review the calculation of the node volume in 2DCyl.	2024-07-11 18:19:26 +02:00
JGonzalez	72dd3c9073	I think I have it So the radius for the volume integral in the nodes has to be set with a 1/4 3/4 ratio to match the change in volume as the radius changes along the volume. This has a bigger impact close to the edge. Still unsure if this is the "right" thing to do, but at least it works okay (it seems).	2024-07-11 17:40:57 +02:00
Jorge Gonzalez	bc8f205287	Clean-up! Trying to reduce warnings and unused variables in the code. This should not be in this branch.	2024-07-11 15:55:02 +02:00
Jorge Gonzalez	d710cc5e66	Probes fixed! This should've been commited in another branch, but I am taking advantage of this to clean up the code a bit.	2024-07-11 15:52:13 +02:00
Jorge Gonzalez	03f3b81811	Manual updated The manual has been modified to account for the new particlePerEdge option and to indicate that now when the flux is giving by a flux like in the Am2 units, the real surface of the inject is used to scale it.	2024-07-11 15:31:01 +02:00
Jorge Gonzalez	6d62807fb4	Remove testing 'print' statements 2DCyl still seems to have some problems but waiting for simulation.	2024-07-11 14:50:14 +02:00
Jorge Gonzalez	96c563c146	Finally, some progress I rewrote how particles are injected. Now the particles per edge and its weight are calculated in the initialization. There is the possibility for the user to select the particles per edge. TODO: Write documentation for new feature. TODO: Test in 2DCyl	2024-07-11 14:39:56 +02:00
Jorge Gonzalez	e23fc2fc2c	Small progress I made some small changes to how things are calculated. I have also discovered that the issue with different density when changing injection is not related with the node volume but with the way injection is carried out. When loading particles from a file, all provide the same density regardless the cell (node) volume. I am doing testing in 2DCart as it is easier to set up.	2024-07-11 11:21:38 +02:00
JGonzalez	9d961bb85f	Removing weight from edges This parameter is no longer needed. Also removing the cumulative sum of weights from the injection.	2024-07-10 21:57:22 +02:00
JGonzalez	cb92462f36	New injection based on surface to all geometries. WARNING: 3DCart still not working (too tired to calculate things and I'm not ussing it...)	2024-07-10 21:55:45 +02:00
JGonzalez	d608639e77	Forgot to save Sorry.	2024-07-09 22:05:16 +02:00
JGonzalez	667a2ecd93	So, no idea... Basically things do not work. I've added a correction to the node volume in the axis which gives okays results but still this is not perfect. I need to find a better way to do things. Also, I've noticed that the density changes with the size of the cells, which should not happen! I'vw to check this issue.	2024-07-09 21:57:32 +02:00
JGonzalez	11831a973d	Still not working Trying to have a very simple volume per node assuming a rectangle and the density at the axis it higher than it should (kinda like when using the more accurate volume calculation). This is still weird. I also suspect that the size of the first cell in the axis will also affect this...	2024-07-09 21:25:30 +02:00
JGonzalez	b36f9c2615	Shifting towards constant number of particles per edge So now each edge has the same number of particles and the weight of each particle is calculated based on the surface of each edge compared to the total one. Only in 2DCyl, still to extend to other geometries. Not perfect constant density, but the issue might be the node volume.	2024-07-09 17:49:42 +02:00
JGonzalez	fa23f9481a	Issue with bibliography Backup files are properly ignored now.	2024-07-07 14:40:26 +02:00
JGonzalez	1f2ec8d82f	New option for initial distribution of species The number of particles per cell can be defined when giving an initial distribution fora species. If not, the typical method of using the species weight is used. This is particularly useful for cylindrical coordinates in which very little particles might end up in the axis if a constant weight is used.	2024-07-07 14:37:34 +02:00
JGonzalez	5bc064d018	Some parallel issues Forgot to check Gmsh2 format with new changes. Cartesian coordinates were not calculating things properly.	2024-07-07 11:53:56 +02:00
JGonzalez	626e970d82	Some progress Fixed an issue with random integer numbers. Cylindrical coordinates are not perfect yet: - Box (cylinder) with initial constant density loses particles at r = 0 - Injection density still low in r = 0	2024-07-06 19:14:44 +02:00
JGonzalez	6b5ac16e4b	Still working on it No uniform density yet...	2024-07-06 10:12:03 +02:00
JGonzalez	b972120ed5	Switching to variable particle weight I have to change the injection of particles. Each edge will receive a similar number of particles and their weight will change to have a constant density based on the geometry. Still testing.	2024-07-04 10:56:13 +02:00
JGonzalez	6b78ae3738	Added an adhoc factor to correct injection I think that the volume of the nodes is not being well calculated, maybe we need a better volume calculation for this, using multiple points (as it is done for K)	2024-07-03 22:04:28 +02:00
Jorge Gonzalez	0471517526	Merge branch 'issue/injection' into 'development' Change in injection to achieve uniform density of particles. See merge request JorgeGonz/fpakc!50	2024-06-30 08:50:25 +00:00
Jorge Gonzalez	e277fe6ddb	PI is no needed here	2024-06-30 08:48:51 +00:00
JGonzalez	59a322a4c7	Clean up Fixing calculation of node volumes.	2024-06-30 10:46:05 +02:00
JGonzalez	cd7bf66bd8	A workaround The random position for edges in the axis is corrected so that there is a more uniform charge density in the axis. Still, things are not perfect and this is something to really look into in the future.	2024-06-30 10:36:36 +02:00
JGonzalez	4cadfe5367	Seems are a bit better There is still less density in the axis. I don't find a reason why. There must be a modification to the weight...	2024-06-29 22:22:10 +02:00
JGonzalez	5386114d15	Cylindrical injection working better Seems things are a bit better. Still, more cases are needed and still not perfectly uniform...	2024-06-29 14:58:48 +02:00
JGonzalez	6389c8ba2d	Quick because food Cartesian fixed now	2024-06-27 12:08:08 +02:00
Jorge Gonzalez	0ce921a48c	Merge branch 'issue/quadLocation' into 'development' Issue with calculating coordinates in quads See merge request JorgeGonz/fpakc!49	2024-06-26 13:14:22 +00:00
JGonzalez	c6470819e8	Issue with calculating coordinates in quads The third coordinate (unused) was causing some errors when it was becomming too large.	2024-06-26 15:11:01 +02:00
JGonzalez	e4f7987f90	Trying to solve Still I don't understand this basic thing...	2024-05-19 16:45:03 +02:00
JGonzalez	a3bdf8230a	Implementation of Boltzmann electrons Still not working, just saving code.	2024-05-19 10:55:20 +02:00
Jorge Gonzalez	62b8873786	Merge branch 'feature/Amps_per_m2' into 'development' Implementing injection with current density See merge request JorgeGonz/fpakc!48	2024-03-29 12:09:30 +00:00
Jorge Gonzalez	d86b3a3417	Implementing injection with current density WARNING: This current denstiy will be multiplied by the reference length, no the surface area that is being used for injection! New units in the injection of particles 'Am2' to inject a density current. Manual has been modified accordingly. Reference parameters are now also printed in the case folder.	2024-03-28 09:45:46 +01:00
Jorge Gonzalez	a3d7b38e3b	Merge branch 'improve/secondaryElectronIonization' into 'development' Improve ionization See merge request JorgeGonz/fpakc!47	2023-11-24 09:34:30 +00:00
Jorge Gonzalez	e41b448ef8	Different species for secondary electrons The option to have a different species than the impacting electron for secondary electrons from ionization is introduced.	2023-11-24 10:30:50 +01:00
Jorge Gonzalez	4585390b50	Fixing issue Fixing an issue with reading tables led me to other issues with collisions that I think are fixed right now. I am testing with the 1D ionization model for ALPHIE and things seems to be working properly.	2023-11-21 09:53:36 +01:00