Basically I am generating a grid in space (so a vector where v[[1]]=position 1 in x) and a random number of particles in space (so a vector where p[[1]]= position of particle 1 on space).

Afterwards, I want to see in which spaces of the grid the particles are and do somethings.

The Do is the specific part of the code I want to optimize.

rho = Compile[{{xp, _Real, 1}, {xg, _Real, 
1}, {carga, _Real}, {np, _Real}, {ng, _Real}, {dx, _Real}},


n = ng + 1;
  ρ = Table[carga*np/(dx*ng), n];
  ρ[[ng + 1]] = 0;
  
  Do[If[Abs[xp[[i]] - xg[[j]]] <= 
     dx, ρ[[j]] = ρ[[j]] - 
      carga*(1 - (Abs[xp[[i]] - xg[[j]]])/dx)/dx; r += 1, 0] , {i, 1, 
    np}, {j, 1, ng + 1}];
  ρ[[1]] = ρ[[1]] + ρ[[ng + 1]];]

Basically I am worried that I am thinking too much in terms of other languages and maybe mathematica has a better way to do this instead of using

If[Abs[xp[[i]] - xg[[j]]] <= 
     dx

That manually makes a lot of operations

Basically I am generating a grid in space (so a vector where v[[1]]=position 1 in x) and a random number of particles in space (so a vector where p[[1]]= position of particle 1 on space).

Afterwards, I want to see in which spaces of the grid the particles are and do somethings.

The Do is the specific part of the code I want to optimize.

rho = Compile[{{xp, _Real, 1}, {xg, _Real, 
1}, {carga, _Real}, {np, _Real}, {ng, _Real}, {dx, _Real}},


n = ng + 1;
  ρ = Table[carga*np/(dx*ng), n];
  ρ[[ng + 1]] = 0;
  
  Do[If[Abs[xp[[i]] - xg[[j]]] <= 
     dx, ρ[[j]] = ρ[[j]] - 
      carga*(1 - (Abs[xp[[i]] - xg[[j]]])/dx)/dx; r += 1, 0] , {i, 1, 
    np}, {j, 1, ng + 1}];
  ρ[[1]] = ρ[[1]] + ρ[[ng + 1]];]

Basically I am worried that I am thinking too much in terms of other languages and maybe mathematica has a better way to do this instead of using

If[Abs[xp[[i]] - xg[[j]]] <= 
     dx

That manually makes a lot of operations

Here is a part of the code I am using

eps = 8.85*10^(-12);
q = 80;
m = 50;
L = 2;
nparticles = 2000;
ngrid = 200;
vth = 2;
u = 1;

dx = L/ngrid;

x = RandomVariate[UniformDistribution[{0, L}], nparticles];
r = 0;

xgrid = Range[0, L, dx] // N;
rho[x, xgrid, q, nparticles, ngrid, dx]; // AbsoluteTiming
r