Simple Particle In Cell Code in Matlab

Posted on November 28th, 2011
Previous Article :: Next Article


This post includes the sample Particle-In-Cell (PIC) code that goes with our previous article on the electrostatic particle-in-cell method. It simulates the flow of uniform plasma past an obstruction – a charged plate in our case. This simple code is implemented in Matlab. If you don’t have Matlab, you can run the code using Octave, which is an open-source Matlab emulator. If you do end up getting Octave, as of the date of this article, the MinGW build available on Octave-Forge is out of date and contains a buggy contouring function. In order to use the latest version – and see the plots – you need to install Cygwin and add Octave as a Cygwin package…

File Organization

The sample PIC code is contained in two files: flow_around_plate.m and eval_2dpot_GS.m. The first file is the actual PIC code, while the second file holds the potential solver. The solver computes electric potential using the Finite Difference method. The code simulates the flow of solar wind (low density plasma) over a charged plate. The code is not particularly fast (it’s in Matlab after all!), but shouldn’t take more than few minutes to run on a standard PC. If you run the code through Octave, make sure to first execute more off to avoid pagination of the output.

Input Parameters

The simulation plasma environment is set on lines 24-28:

n0 = 1e12; %density in #/m^3
phi0 = 0; %reference potential
Te = 1; %electron temperature in eV
v_drift = 7000; %ion injection velocity, 7km/s
phi_p = -7; %plate potential

By changing the plate potential, you can alter the size of the wake that forms behind the plate. At sufficiently large values, ions are attracted back to the plate. This behavior is analogous to the plasma environment that forms around satellites utilizing electric propulsion. Slow moving ions are extracted out of the thruster plumes and can impact surfaces without a direct line of sight to the thruster.

The simulation domain size and the number of time steps is set starting with line 35:

nx = 16; %number of nodes in x direction
ny = 10; %number of nodes in y direction
ts = 200; %number of time steps

At the default settings, 200 steps is sufficient to reach steady state. The simulation time step is such that a particle traveling at the drift velocity will cross 10% of the cell length per time step. Since particles speed up as they approach the plate, this conservative time step keeps particles from traveling more than a cell per time step. This is a necessary constraint in all PIC codes, but in general you want to keep the particles traveling a shorter distance (I generally use 25% in my simulations). The number of particles generated at each time step is specified on line 39,

np_insert = (ny-1)*15; %insert 15 particles per cell

Plate dimension are set on lines 48 and 49:

box(1,:) = [floor(nx/3) floor(nx/3)+2]; %x range
box(2,:) = [1 floor(ny/2)]; %y range

For an explanation of the physics used in the code, please see the Particle In Cell page. Feel free to leave a comment if you have any questions.

plasma potential plasma density
Figure 1. Simulation results, plasma potential and ion density


Download the files here: flow_around_plate.m, eval_2dpot_GS.m

You may also be interested in a Finite Element Particle In Cell example.

94 comments to “Simple Particle In Cell Code in Matlab”

  1. November 28, 2011 at 8:20 am

    Please leave a comment if you have any questions. Thanks!

  2. Ozgur Tumuklu
    December 1, 2011 at 2:37 am

    Hi Lubos;
    in Matlab code your code the wall potential is given phi_p = -7V. Is there any special meaning under this or was it chosen randomly,
    and, what is the result of choosing positive wall voltage?
    Thank you..
    Best regards,

    • December 1, 2011 at 7:00 am

      Hi Ozgur,

      Yes, this is just some arbitrary value that produced nice looking plot. The potential is negative since the potential in the bulk plasma is zero and hence this attracts the ions towards the plate. If you chose a positive value, the ions will be repelled away.

  3. Ashley
    December 1, 2011 at 7:02 am

    Our professor assigned us an assignment on Particle in Cell simulation for his Plasma Engineering graduate class. We are supposed to use your code as a basis. After altering your code a bit, we basically moved and extended the wall/barrier to the right edge and changed a few parameters.

    Is there something in your code that generates a decrease in the density on the left hand side of your density plot? Since there’s a uniform injection of ions, shouldn’t there be a uniform density until it reaches the sheath region? We just can’t really understand why the left hand side is also going to zero since it’s technically not supposed to do that. In your density plot online, it seems that the density goes to zero, as well, so maybe we’re not wrong, I just don’t understand why.

    • December 1, 2011 at 7:03 am

      Hi Ashley, that has to do with the way particles are injected into the domain. The particles are born at a random “x” position in the first cell and also with a random velocity. So you basically end up with a non-uniform density in this first cell. You can avoid this by loading particles at x=0 and also doing the density update prior to the push. I left it this way since I didn’t think it made much difference on the solution at the plate location.

  4. Ozgur Tumuklu
    January 24, 2012 at 12:35 pm

    Hi Lubos,
    i try to simulate plasma flows by using Particle in Cell method. The situation is very similar to your matlab code.;(i.e. collisionless and quasi neutral plasma,and the model is hybrid). i have finished the code in non-rectangular domain. But i could not verified the results whether true or not. Can you help me about verification.
    Thank you very much for your attention.

    • January 24, 2012 at 2:09 pm

      Sure Ozgur, what did you have in mind? Some ways you can also try to verify your code is by testing the potential solver with some analytical function and also checking for energy conservation.

  5. Sergio Flores
    April 21, 2012 at 11:38 pm

    A PIC code in Matlab could be adapted to run in a computer cluster? and/or adapted to run in parallel?, Thanks

    • April 22, 2012 at 8:51 pm

      Hi Sergio, if I am not mistaken, Matlab now supports parallel processing. However, I am not familiar with it – we do most of our actual development either in C/C++ or Java. Matlab is a great prototyping language, but if you are after speed, you’ll be better off using a non-interpreted language.

  6. Amitavo
    April 25, 2012 at 9:25 am

    Hi…I am new in PIC..I want some help how to do assignment from particle position to grid position…how I can implement…

    Also how I may know how many particle are there in one particular cell of the whole geometry?

    I want to interpolate particle properties like velocities and charge to my grid (one dimensional only)

    As in each cell there are so many particles so each has it own contribution on particular node (let say k) depending on weight factor so I want to know how I can sum all those properties on a node (k and k+1).

    Why you have add 1.0 in the assignment part (fi = 1+part_x(p,1)/dh; )?
    Please elaborate the assignment part.

    Please help me…

    • April 27, 2012 at 5:28 am

      Amitavo, you basically want to compute the average velocity. This is no different from computing the regular average, with the difference that you are averaging fractional contributions. When you scatter data to the grid, you place a fraction of the velocity on the left node, and the rest on the right node, assuming 1D problem. So as an example, let’s say the velocity is 100 m/s. Perhaps the particle is really close to the right node, so the right node may get 90 m/s, while only 10 m/s is deposited to the left node. Next imagine that this cell contains 5 particles, and that they are all located at the same spatial position. You will end up with 450 m/s total velocity on the right node and 50 m/s on the left node. If you were to divide the velocity sum by 5, you will get that the velocity varies from 10 m/s on the left end to 90 m/s on the right end. This is clearly incorrect since all particles in the cell are at 100 m/s. To get the correct answer you need to divide the sum by the fractional particle count. This is why you want to also keep a “count” variable on the grid nodes to which you scatter the value of 1 (if all particles have equal weight) for each particle. If particles have a variable weight than you scatter the respective specific weight. You then divide each node value by the count, i.e. vel[i]=vel_sum[i]/count[i];. In this example, the counts would be 4.5 on right, and 0.5 on left, giving you 100 m/s on both the left and the right node.

      The second question about adding 1 is due to Matlab starting its indexing at 1 instead of 0. If your code is in C or Java, you don’t add this 1.

  7. Amitavo
    April 27, 2012 at 7:16 am


    • April 27, 2012 at 8:09 am

      It’s no different than scattering velocity, except that you are scattering particle weight. If all particles have the same weight, you simply scatter the value of 1.0 every time you scatter the velocity. If they have different weights, you scatter the particle weight, or the weight divided by a reference weight (to normalize, since weights are typically very large numbers).

      • Amitavo
        April 30, 2012 at 12:24 am

        Thanks a lot….

  8. Amitavo
    May 3, 2012 at 4:06 am

    I am facing problems in assingment process after pushing the particles.
    I am not able to implement the particle properties from the new position of the particles to the grid point which is shared by two cells.

    After the localization process for the new particle positions how I can implement the assignment operation?

    Please guide me as I am new in the field.

    As an example: Suppose we have 10 particles in (K,K+1) cell and 5 particles in (K-1,K) cell after pushing then I want to assign velocities and charge to K node due to all these particles. Please suggest me/ help me how can I perform this operation.


  9. John Mathai
    May 28, 2012 at 3:28 am

    I have found your writings very useful. I am learning plasma out of interest. I have a doubt: if I have a 2d plasma inside 4 walls of 0 Volt each and I want to plot the potential (calculated as in your tutorial, electrons as Boltzmann particle, ions : moving macroparticles), why do I get an inverted Gaussian type curve? I thought potential should be uniform across a plasma! Thanks.

    • May 28, 2012 at 7:40 am

      Hi John, that would be the case if the governing Poisson equation reduced to the zero forcing function Laplace form, $latex \nabla^2\phi=0$. However, that is not the case. Even though plasma is quasi-neutral, this condition does not hold in the vicinity of walls. Near the walls you will get what’s known as a sheath, where $latex n_e

  10. olatunji
    May 30, 2012 at 11:03 am

    Am working on simulation of DBD using particle in cell method.PLS Sir how will i go about it.

    • June 1, 2012 at 6:51 am

      That’s not a particularly easy problem. I must admit I do not have much background with modeling discharges (most of my work has been in plume modeling), however what I imagine you would need is a code that takes into account dielectric charging, surface deposition, as well as a model for breakdown. Have you seen the excellent book by Lieberman and Lichtenberg, Principles of Plasma Discharges and Materials Processing? It goes into details of various plasma discharges.

  11. John Mathai
    June 1, 2012 at 1:19 am

    Thanks for the explanation. If the electrons are also simulated as macroparticles like ions (with same specific weight but with opposite charge), what is to be done for the sheath formation in the simulation? Or, will it be inherent? How the time-step will be defined then?
    Secondly, in the Matlab code for flow around a plate, why the factor 0.5 is missing for x-component of velocity in the Particle Generation section of the code? It is there for y-component of velocity.

    • June 1, 2012 at 6:58 am

      The sheath will form self-consistently. It forms since the electrons are more mobile (lower mass -> higher velocity) and thus will be absorbed by the wall faster than the ions. This will in turn make the wall negative in respect to the bulk plasma (or the plasma potential will go up if the wall potential is held fixed). This potential gradient acts to reduce the electron flux. The wall potential will continue to go down until the electron current is balanced by the ion current. About your second question – good catch! I have no idea. It could very well be a typo. The way it’s written you will end up with anisotropic temperature. I wrote the code few years ago so I don’t quite remember what I was thinking back then.

      • John Mathai
        June 5, 2012 at 10:17 pm

        Thanks for your reply. Also, thank you for maintaining such an active website.

  12. John
    June 18, 2012 at 9:59 am

    If you had an unstructured triangular mesh instead of a mesh with quadrilateral elements, what would be a good algorithm for calculating the electric field at each of the mesh nodes. For the scatter and gather operations I was thinking you could use barycentric coordinates of a triangle.

    But once you have the potential at each of the nodes on an unstructured mesh, what would be a good way to calculate the E field at each node.

    • June 18, 2012 at 10:03 am

      This is really the foundation of the finite element method. Typically what you do is you define a basis element for each cell. These elements can be of different order. The derivative will be one order less. Quite typically, you may see linear elements used for potential which then give you zero-order (constant) electric field in the cell.

  13. John
    June 24, 2012 at 12:08 pm

    I tried implementing this PIC code using the finite element method solver from this blog post.

    Here are some results I got.

    However, in order to get a stable result I needed to set

    Te = 0.125; %electron temperature in eV

    instead of Te = 1.0; . When I do this I get a result that looks similar to what the original PIC code from this blog post generates with Te set to 0.125 .

    Do you know why the solver result might become unstable when Te is set higher than 0.125 ?

    I used the following in the MIT FEM code for f(x,y).

    % multiply Fe by f at centroid for load f(x,y): one-point quadrature!


    Does this look correct to you?

  14. mahnaz
    July 3, 2012 at 11:30 am

    Hi Dr. Brieda

    I have some question about ES-PIC code. Please help me.
    In your code, was electron and ion isothermal together? otherwise What is your scientific reason?
    How have been calculated time step (10^-12) for electrons?
    Why the drift velocity has been considered only for ions? note that 10^5 (order of thermal velocity of electrons) isn’t much larger than 10^3 (order of drift velocity)
    in calculation thermal velocity , vth=k*Te/m was used, in spite of temperature of electron was in ev not Kelvin!

    best regards

    • July 13, 2012 at 5:18 am

      Hi Mahzan, these assumptions are typically taken in the context of plume modeling. The reason is that including electrons as kinetic particles (or even a fluid) complicates the numerical model unnecessarily. If magnetic field is not present, the electron density follows the Boltzmann relationship.

      Thanks for catching that vth mistake, I have corrected it. I have also specified a different temperature for ions to indicate that generally LEO plasma will not be isothermal.

  15. john
    July 6, 2012 at 9:14 am

    In eval_2Dpot_GS.m you set the neumann boudary conditions around the border as follows.

    %set boundaries
    b(1:nx) = 0; %zero electric field on y=0;
    b(nn-nx+1:nn) = 0; %zero electric field on y=L;
    b(nx:nx:nn) = 0; %zero electric field on x=L;
    b(1:nx:nn) = phi0; %fixed potential on x=0;

    But in flow_around_plate.m you set the electric field around the border as follows.

    efx(1,:) = 2*(phi(1,:) – phi(2,:)); %forward difference on x=0
    efx(nx,:) = 2*(phi(nx-1,:) – phi(nx,:)); %backward difference on x=Lx
    efy(:,1) = 2*(phi(:,1) – phi(:,2)); %forward difference on y=0
    efy(:,ny) = 2*(phi(:,ny-1) – phi(:,ny)); %forward difference on y=Ly

    Should the normal electric field on the borders y=0, y=L, x=L just be set to zero to match the boundary conditions?

    efx(nx,:) = 0;
    efy(:,1) = 0;
    efy(:,ny) = 0;

  16. Maks
    October 18, 2012 at 7:10 pm

    Great example! Thank you.
    I downloaded all files and ran them at the Octave. I received strange results: at every time step I have 0 particles and blank graphic.
    Have you any suggestions about how I can solve my problem?

    • Maks
      October 18, 2012 at 7:23 pm

      I’m sorry for dumb question..
      I found where was a problem.

  17. sudhan
    November 17, 2012 at 6:47 pm

    I downloaded all files and ran them at the Matlab. At every time step I have 0 particles and blank graphic.

    • December 8, 2012 at 12:03 pm

      Hi Sudhan, thanks for catching this. The equation for thermal velocity must have gotten screwed up at some point. I fixed it. Try it again and let me know if it still doesn’t work.

  18. WEI Wei
    November 27, 2012 at 8:44 am

    Thanks for your sharing. I was looking for MPM on the website to solve some normal fluid problems. Then I got here. I saw some interesting results from MPM, for example: Interesting, haha.

    • December 8, 2012 at 11:04 am

      Thanks for the link!

  19. Sayan
    February 4, 2013 at 4:07 am

    Is it possible to incorporate magnetic field to this code??
    Another thing, How I can plot the velocity distribution of the particles from this code?

    • February 4, 2013 at 4:30 pm

      Yes, you can have a static magnetic field. See the article on the Boris method for details of implementing the particle pusher.

      You will generally need to write your own code to compute the VDF. Just create a number of bins and iterate over the particles and increment the bin counters based on their speed. Then you can normalize the bins if needed by the max value or the sum.

  20. Sayan
    February 7, 2013 at 3:55 am

    I have tried to plot the electric field using quiver command from the data efx and efy and got a result. Problem is that I am getting a swirl kind of structure in the electric field. which is not acceptable as i have also introduced two parallel plates instead of one. One I have kept at y=Ly and another at y=0. having width $ \delta $ x = 3

  21. February 11, 2013 at 11:40 pm

    Hi, first of all I would like to say this is a great blog. I had been reading and testing some of the codes you posted here.

    I am doing a thesis on an inertial electrostatic confinement reactor (basically have to model flow of plasma due to a static electric field). I have two concentric (mostly transparent) grids held at about 60KV potential. Gases are ionized at the outer grid and accelerate towards the central grid. In the center some undergo fusion. Ions will also undergo recirculations (if they don’t collide with the central grid).

    So I am wondering if it will be possible to alter the geometries of this code (perhaps 2D circular geometry), and have grids rather than plates? Would it be viable to model an IEC device using this method? I would greatly appreciate any help you can provide.

    Kind regards,
    Jimmy Zhan

    • April 10, 2013 at 4:12 am

      Hi Jimmy, you should check out our Starfish code (still in development). It is a general 2D solver that can be used for a number of different geometries.

  22. March 4, 2013 at 10:31 pm

    Great blog,I have been working on the codes posted here. The link you posted here was very useful.I am learning plasma out of interest. Thank you so much.

  23. Sonu Sen
    February 8, 2014 at 2:54 am

    How this particle in cell simulation may be used for “self-focusing of laser beam in plasmas”

  24. Prakash
    June 11, 2014 at 12:44 am

    This is really a very fascinating site on PIC. With the above example of 2D PIC, is it possible to calculate current recorded by the plate as a function of time ? Can we use the formula (charge hitting the plate in time dT) / dT to calculate current at a particular time instant ? Since it is a 2D simulation, what will the current mean physically ? Is this the current collected by a plate of given x and y dimensions but with UNIT DIMENSION along Z (i.e. 1 meter) ? Your thought on this will be very helpful. Thanks!

    • June 12, 2014 at 4:20 am

      Hi Prakash, yes, you would basically count the number of particles hitting the plate and average over some time. Since this is an ion-only simulation, you will also need to include the electron term, which for can be modeled analytically from the planar sheath equation. For a 2D XY simulation like this one, the current will be current per unit depth as you mentioned. If this were an axisymmetric simulation, then you could compute the actual surface area by considering the body of revolution.

      • Prakash
        June 12, 2014 at 5:05 am

        Thanks a lot. This is the only elaborately-explained resource on the web on Particle-In-Cell. Please keep it up.

  25. Helen
    July 18, 2014 at 5:42 am

    Hi,I would like to say this is a great blog. It helps me a lot in my final project.
    I have tried to understand your code. About Debye length calculation,I could not understand lD = sqrt(EPS0*Te/(n0*QE)); %Debye length. According to formula, I thought is lD = sqrt(EPS0*Te*K/(n0*QE*QE)). Can you explain it, please? Thank you!

    • July 18, 2014 at 7:35 pm

      Hi Helen, this is because in the above example, electron temperature Te is in electron volts. The conversion from Kelvin to electron volts is k*Te(K)/QE = Te(eV).

      • Helen
        July 21, 2014 at 2:16 pm

        Thanks,Lubos. Can i ask about u=(j-1)*nx+1? what is this for ?

  26. farid
    July 19, 2014 at 3:22 am

    can i use this code for higher densities (about 10^20 #/m^3) ?
    is this code suitable for tokamak plasma simulation?

  27. john
    August 11, 2014 at 12:09 am

    Hi Lubos,
    I’m interested to particle in cell LINAC (linear accelerator) application where it used protons charges. Is Your matlab code adapt for this application? However I have a series coupling resonant cavity where the proton beam through it. For example initial energy is 7 MeV and the cavity were designed to accelerate untill 10 MeV.

  28. hiral
    October 20, 2014 at 9:39 am

    Hi lubos,
    I faced one problem to understanding the 2 lines
    A(u,u-nx)=1/(dh*dh); %phi(i,j-1)
    A(u,u+nx)=1/(dh*dh); %phi(i,j+1)

    why you take ‘u’ constant for j ? and why you use ‘u-nx’ and ‘u+nx’?

    kindly waiting for your reply

    • October 23, 2014 at 6:44 am

      Hi Hiral, this is the “y” derivative part of the laplacian, $latex \nabla^2 \equiv (\partial^2/\partial x^2 + \partial^2/\partial y^2)$. The phi data is actually a 2d matrix. However, in order to perform a matrix-vector multiplication, this 2d data needs to be flattened into a 1D array. One way of doing this is to take all the values for j=0 (such as phi(1,0), phi(2,0),.. phi(nx,0)) followed by the values for j=1, etc… Using this scheme, a value corresponding to (i,j+1) will be located nx entries away from (i,j). If u corresponds to the value at (i,j), then (u+nx) will correspond to the value at (i,j+1). The actual computation is u=j*nx+i.

  29. Hiral
    November 13, 2014 at 4:43 am

    Hi Lubos,

    I am trying to make PIC for Lunar model using coulomb’s Law. But, I am confusing with output value of force. I am sending that matlab file on Email ID. Please, kindly inform me the solution.

  30. Hiral
    November 19, 2014 at 8:44 am

    hi Lubos

    Here, “phi” – is surface potential or electron potential? and also “den”- is electron density?

    Kindly waiting for your replay.

  31. zarana
    November 21, 2014 at 1:24 am

    hi Lubos,

    I want to simulate the same problem using circle instead of Rectangular plate..What changes I can do??

    Waiting for our reply..

    • December 4, 2014 at 12:26 pm

      Hi Zarana, you have basically two options. You can either write the code in 3D, or you can implement an axisymmetric (RZ) version of the code. The second is probably what you had in mind. Looking at the blog now, I just realized that I don’t think I have a post on an RZ PIC code. I will try to post it over the weekend.

  32. Daniel
    May 4, 2015 at 9:51 am

    Hi Lubos.
    please can u help me with the code to simulate particle (e.g Red Blood cell) trajectories in micro channel (e.g Deterministic lateral displacement).
    I will appreciate any assistance rendered.

    • May 7, 2015 at 9:29 am

      Hi Daniel, I sent you an email but replying here as well in case you didn’t get it. Can you email me more info on what you are trying to do?

  33. Ajay
    July 2, 2015 at 10:32 am

    Hi Lubos, How do those neumann and dirichlet boundary conditions in the code work.
    For example, in

    for i=1:nx
    u = i
    A(u,u) = -1/dh;
    A(u,u+nx) = 1/dh;

    How are the values for A(u,u) and A(u,u+nx) formulated?

    • July 2, 2015 at 11:09 am

      That just means -phi[i][j] + phi[i+1][j] = 0, or d(phi)/dx = 0, since on the boundary d(phi)/dx is estimated as (phi[i+1][j]-phi[i][j])/dx.

      • Ajay
        July 3, 2015 at 11:04 am

        But in the comments you mention that
        A(u,u) refers to phi(i,j), and A(u,u+nx) refers to phi(i,j+1). That would refer to the element right above y=0, not the one right next to phi(i,j)? Am I confusing something here?

        • Ajay
          July 3, 2015 at 11:08 am

          Actually, I think I figured out why it works that way šŸ™‚ Thanks you!

  34. Dan
    July 3, 2015 at 6:26 am

    I know nothing about PIC, but the instructor of a plasma related course told us to write a code of “PIC simulation of 1-D collision-less argon gas discharge”, including a model where the ions move and one that the ions don’t as our final project. A few different sets of data are required, like V-x, E-x, Ion density-x and so on. A lot of my classmates are as perplexed as me, as we have no prior experience with this type of work. Could you please be so kind and help me out?
    Also, I can’t seem to run the example codes on this page with MATLAB.

    • July 3, 2015 at 9:06 am

      Hi Dan, I responded to you through email. Check your spam folder if you don’t see it.

  35. Ajay
    July 6, 2015 at 10:22 am

    Hi Lubos,

    In the potential solver code, what exactly is the following doing (I get it that from your comments you are setting up the boundaries, but am having difficulty visualizing it, given the notation used).

    b(1:nx) = 0; //Rows 1 to 13 have value 0?
    b(nn-nx+1:nn) = 0; //Rows 157 to 169 have value 0?
    b(nx:nx:nn) = 0; // .. ?
    b(1:nx:nn) = phi0; // .. ?

    Thank you!

  36. Terry.Zhu
    July 9, 2015 at 1:06 am

    Hi Lubos, Iā€™m a rookie in PIC simulation
    The Leapfrog method is commonly used in PIC codes. The Leapfrog method is commonly used in PIC codes. The times at which velocity and positions caculated are offset from each other by half a time step. As such, the two quantities leap over each other.
    I want to know how to move V(0) backword half a time step to V(-dt/2). As such, the two quantities leap over each other.
    I want to know how to move V(0) backword half a time step to V(-dt/2).

  37. Terry.Zhu
    July 9, 2015 at 1:07 am

    Using backward-difference :V(0)-V(-dt/2)=E(0)*dt/2 ?

  38. Atefeh
    July 29, 2015 at 6:50 am

    Hi lubos
    which code of pic simulation is suitable for laser driven accelerator? can I buy it? how can I access that?
    Thank you

  39. Ajay
    August 5, 2015 at 2:44 pm

    Hi Lubos,

    Just a quick question. For the Right hand side of our defining equation you have -e/eps(n(i) – n(0)exp(phi/kTe)).

    In the code, in eval_2dpotGS you use
    b = b0 – n0*exp((x-phi0)/Te);

    where b0 = reshape(den) etc.

    Does this mean you are not considering the ion motion, and only the electrons?

    • August 23, 2015 at 10:16 pm

      Hi Ajay, the ion charge density is obtained from the particles, by scattering them to the grid. That’s what gives the “den” term. The n0*exp(..) term is the density of electrons from the Boltzmann relationship. Since charge density is number density times charge, you have rho = q*(den_i – den_e), assuming singly charged ions.

  40. Enrique Valderrama
    December 29, 2015 at 8:39 pm

    Thank you!!! these is a great work!
    many questions šŸ™‚

    what means 1D for the cell size? or what are the dimensions of the plate?
    why there is ion density inside of the plate? I guess is for the type of border condition, but still I dont understand, is not possible to put some kind of step function on that border, not just a point where is zero, but beyond as well?
    where the data get stored? there is a way to print in a file the density and potential matrices?
    if I want to have only ions, I put the temperature(or density) of the electrons very small?
    what is the difference between Ti and v_drift ?

    I am running octave-gui on windows 10, and so far so good!


    • December 30, 2015 at 10:08 am

      Enrique, let me try to answer your questions:

      1) lD is Debye length. The Poisson solver may have convergence problems if the cell size is bigger than this. This is because we don’t resolve non-neutrality if cells are greater than the Debye length. The box just spans 1/2 of the height.

      2) There isn’t. It’s a plotting issue. Density is stored on the nodes, and there is a node at the box outer surface. So you have one node with some finite value and another on the right with zero. The contour plotter fills the gap in between by linearly interpolating from the value down to zero.

      3) Sure, you can export the results to a file. This was just a simple Matlab example. If you for instance check out the more recent finite element PIC example, you will see that one exports to a file and Paraview is used for plotting.

      4) You could just set the n0 Boltzmann relationship term to zero. But this will likely result in very high potentials / the solver diverging / particles flying at crazy velocities.

      5) Vdrift is a constant velocity. Ti controls the ion temperature. Gas temperature is basically the spread in velocities. So if Ti was zero, all ions are born with v=v_drift. The non-zero Ti then adds some random component to the velocity.

      Hopefully these helped!

      • Enrique
        December 31, 2015 at 1:11 pm

        thank you that was very helpful

  41. Mike
    April 7, 2016 at 11:46 am

    Hi, thanks for your support in this interesting topic. I have a problem with the eval_2dpot_GS.m file,

    Error using eval_2dpot_GS (line 20)
    Not enough input arguments.
    The line mentioned is:
    x = reshape(phi, numel(phi),1);

    I tried writing phi0 instead of phi and pass it through, but then appears another mistake in the line
    x(i)=(b(i) – A(i,1:i-1)*x(1:i-1) – A(i,i+1:nn)*x(i+1:nn))/A(i,i);

    Problaby is for my version of matlab, (2014), but Im not sure. It seems that are kind of basic mistakes, but I dont have the expertise to solve it.

    thanks if you can help.

  42. Ajay
    April 27, 2016 at 1:45 pm

    Hi Lubos,

    I have been trying out different values for the parameters used in the code. Is there a lower limit to the values used for v_drift?

    For example, when I tried using a value of 990 m/s, while my graph for the electric potential was quite similar to that in the case of v = 7000 m/s, my density graph was quite messed up.

    I was wondering if this was an inherent limitation of the code?

    Thanks for your time.

    • April 27, 2016 at 10:38 pm

      Hi Ajay,

      I wonder if this is due to thermal velocity becoming larger than the drift velocity. In that case you may end up with newborn particles moving in the -x direction. Ideally the code would check for this and reflect them but that is not implemented as the code assumes that v_drift>v_th.

      • Ajay
        April 28, 2016 at 2:30 pm

        Hi Lubos,

        Why exactly are we requiring that v_drift > v_th. does this condition follow upon the nature of ion speeds particularly in the case of thrusters and such?

  43. Keiichi Tanaka
    October 25, 2016 at 11:28 pm


    I am just remaking flow_around_plate.m and eval_2dpot_GS.m. from ‘download’ site.
    But, following error has been happened.

    matrix unmatched

    error==> eval_2dpot_GS at 22
    b = b0 – n0*exp((x-phi0)/Te); %add boltzmann term for electrons

    error==>flow_around_plate at 136
    phi = eval_2dpot_GS(phi);

    How can I solve this problem.

    Keiichi Tanaka(26/Oct/2016)

    • October 31, 2016 at 10:58 am

      Hi Keiichi, what program are you using to run the example? I don’t have Matlab so I tested the code with Octave. I just ran it and it ran fine.

  44. souma
    November 28, 2016 at 1:24 am

    How can we follow the trajectory of a particle

    • November 30, 2016 at 9:23 am

      You can generate a trace file to which you would output position (and possibly velocities) of the same particle. If you are storing your particles in an array, you would do something like:

      p = random(num_particles)
      //main loop
      - update e field)
      - push particles
      - output pos[p], vel[p]

  45. MOSES
    January 3, 2017 at 9:49 am

    I am simulating on deuterium-tritium plasma heating due to nonlinear plasma interaction will pic be useful?. Also i had problem of assigning charge density on grid via weighting procedures(cic)how can this be done so that one solves the poison equation?

    • January 3, 2017 at 4:02 pm

      This really depends on the plasma density. Generally, PIC is more suitable to lower density plasmas. Using PIC for high density discharges (like in fusion) typically requires a prohibitively large number of particles and mesh cells.

      • MOSES
        January 4, 2017 at 5:19 am

        can you suggest to me a more realistic approach to this problem please

  46. Sukanya
    March 6, 2017 at 9:39 am

    I want to use pic code for the study of magnetic reconnection in astrophysical plasma. I have few questions about it. 1. Can I use your code in this context? 2. What are the boundary conditions and potentials applicable for it? 3. Can you suggest some simple examples of reconnection simulation in PIC?

  47. March 23, 2017 at 11:04 am

    Hello Sir,

    Can I know what are the interpolation techniques used for grid to particles and particle to grid interpolation?

  48. Madhu
    June 13, 2017 at 2:20 am

    In the Poisson solver, the boundary conditions are set in b – the RHS vector. Shouldn’t it be on x – the potential?

    I tried to replace the b by x in the boundaries in eval_2dpot_GS.m but get errors.

    Could you please clarify?

    • June 15, 2017 at 2:21 pm

      With a matrix solver, you solve for all rows in the “x” vector. In order to get fixed values, you set the corresponding coefficients in the A matrix and the RHS such that the appropriate equation reads: “1*x=b”. This then gives the fixed value “x=b”. Similarly for a Neumann boundary you would have something like “x[0] – x[1] = 0”, which implies that x[0]=x[1].

  49. Anirudh
    May 29, 2018 at 10:20 pm

    when i run the above code in octave it gives me an error
    “error: ‘eval_2dpot_GS’ undefined near line 44 column 11 ” ..

    And is there any way to compliment this way in c++ rather than matlab ?

  50. Aravind
    May 29, 2018 at 10:53 pm

    I am getting the following error using a online Octave compiler

    $octave -qf –no-window-system demo.m
    Solving potential for the first time. Please be patient, this could take a while.
    warning: function ./demo.m shadows a core library function
    error: ‘eval_2dpot_GS’ undefined near line 165 column 11
    error: execution exception in demo.m

  51. Valhalla
    July 13, 2018 at 3:25 pm

    Dear professor,
    I am now studying FEM-PIC code, it can running well with PLASMA_DEN=10e10.

    But, with PLASMA_DEN=10e14 or higher, there is error when running.

    With possible reasons, I do my best to read in details. But, still there is error.

    I wish getting your help.

    Thank you very much.

  52. MOSES
    July 16, 2018 at 3:18 am

    How would one deposit and interpolate charge on the grid in one dimension whereas your code is 2d.
    since I’m a new in using

  53. Jayanta
    October 1, 2019 at 6:58 pm

    Can you please suggest me how to reconfigure the code, to model an electron flow around an spherical object?

  54. Shams
    December 16, 2019 at 3:35 am

    I am a new user in MATLAB.

    Could you tell me, please, can I use MATLAB for LASER plasma interaction? Please, help me.

  55. amrit
    June 26, 2020 at 7:04 am

    plz help me. I am trying to solve 2d3v problem in matlb by using pic simulation which book should I follow. I want to learn from basic about pic simulation.

Leave a Reply

You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> <pre lang="" line="" escaped="" cssfile=""> In addition, you can use \( ...\) to include equations.