"Plasma Physics via Computer Simulation" by C.K. Birdsall, A.B Langdon covers almost everything you need to start PIC simulations.

If you already have implemented PIC, you do something like this:

   double s = x/dx;

   int         i = floor(s);

   s -= i;

  rho[i] += q*(1 - s)/dx;

  rho[i+1] += q*s/dx;

To do zeroth order (NGP), you simply replace the inner piece with:

  int i = floor(x/dx);

  rho[i] += q/dx;

Thank you, sir for your help. I have one more doubt.

When we implement the NGP scheme whether the plot of instability change?? Initially when I used CIC the plot was perfectly coming, now the plot looks different

You may check few things to test your code:

1 Check if your code converges to the same answer when you decrease time and spatial steps.

2 Total energy conservation.

3. Growth rate compared to the theoretical one.

4. Growth rate of the most unstable mode your system supports (if you use periodic boundary conditions, that is modes with λ=L/N, N=1, 2, .., L - length of your system).

5. Check if the spectrum of your instability is well resolved; i.e, you have many unstable modes, not just the λ=L and λ=L/2. For example, you should take L = 30*λ, where λ is the wavelength of the most unstable mode.

6. For the case of L >> λ, you should excite the initial seed by hand to check 1, 2, 3. Otherwise you will have many modes developing, and the growth rate will be a mixed average. For PIC and λ/h∼10 you should have very good agreement with the theory (of order of 10^-4 or much better; in relativistic code I did catch the perfect discrepancy due to relativistic mass change for c = 0.05, which is about (1 - γ) ∼ 10^-3 => your accuracy should be better).

The book (Russian translation) has chapter 5 (section 5.6) dedicated exactly to your problem. If numeration is the same - just grok the book :-).

Thank You Sir... thanks a lot