# 3D Parametric Plotting for N Charged Particles in a Nonuniform Electric Field

I have a code (very bottom of post) which plots a static electric field as it passed through some metal sheet with an aperture in it.

I seek to observe the motion of some N charged particles (point charges) randomly dispersed near the aperture-field (seen above). I have created a set of equations for the N particles in x[t], y[t] and z[t] as they sit in the electric field, and can solve the set of differential equations which exist for each particle.

When it comes time to plot the motion of the particles and the electric field/aperture (seen above), I am running into some problems in the last few lines of the code

frames = Table[
Show[v, ParametricPlot3D[
Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j,
numbodies}], {t, 0, tf}, PlotRange -> Automatic,
Axes -> False],
Graphics3D[
Table[{Hue[.35],
Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.0005]}, {j,
numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}];
ListAnimate[frames]


which produces these errors:

Here is the full code:

Clear["Global*"];
Needs["NDSolveFEM"]

q = 1.60217733*10^-19*10;(*Net ion Charge*)

R = Import["https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=1"];
RegionDimension[R];
M = BoundaryMeshRegion[MeshCoordinates[R], MeshCells[R, 2]];
RegionDimension[M];
Volume[M];

r = RegionDifference[
RegionDifference[
RegionDifference[Cuboid[{0, 0, -0.5}, {2, 2, 0.5}], M],
Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]],
Cuboid[{0, 0, -0.5}, {2, 2, -0.4}]];
ToElementMesh[r]["Wireframe"];
pol = -1;

V0 = 4000;
sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] == -pol*
V0/2, (0.4 <= z <= 0.5) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
pol*V0/2, (0.0071 <= z <= 0.0072) && (0 <= y <= 2) && (0 <= x <=
2)], DirichletCondition[
V[x, y, z] ==
0, (0 <= z <= 0.0070) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
0, (-0.5 <= z <= -0.4) && (0 <= y <= 2) && (0 <= x <= 2)]},
V, {x, y, z} \[Element] r]

electricField[x_, y_, z_] = -Grad[sol[x, y, z], {x, y, z}];
v = Show[VectorPlot3D[
electricField[x, y, z], {x, 0.5, 1}, {y, 0.5,
1}, {z, -0.25, -0.00001}, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PerformanceGoal -> "Quality",
VectorScale -> 0.05, VectorPoints -> 7], M];

eforce1 = q*electricField1[x, y, z];

vecForce =
Show[VectorPlot3D[
q*electricField[x, y, z], {x, 0.5, 1}, {y, 0.5,
1}, {z, -0.25, -0.00001}, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PerformanceGoal -> "Quality",
VectorScale -> 0.05, VectorPoints -> 7], M]

mass = 6.52*10^-8;(*particle mass in kg/m^3*)
numbodies = 3;
vel0 = Table[
Partition[{RandomReal[{-0.0001, 0.0001}],
RandomReal[{-0.0001, 0.0001}], RandomReal[{-0.0001, 0.0001}]}, 1],
numbodies]
pos0 = Table[
Partition[{RandomReal[{0, 2}], RandomReal[{0, 2}],
RandomReal[{0.00001, 1}]}, 1], numbodies]

eqs = Table[{x1[j]''[t] == 1/mass*eforce1[[1]] /. {x -> x1[j][t]},
y1[j]''[t] == 1/mass*eforce1[[1]] /. {y -> y1[j][t]},
z1[j]''[t] == 1/mass*eforce1[[1]] /. {z -> z1[j][t]},
x1[j][0] == pos0[[j, 1]],
y1[j][0] == pos0[[j, 2]],
z1[j][0] == pos0[[j, 3]],
x1[j]'[0] == vel0[[j, 1]],
y1[j]'[0] == vel0[[j, 2]],
z1[j]'[0] == vel0[[j, 3]]},
{j, numbodies}];
vars = Flatten[Table[{x1[j], y1[j], z1[j]}, {j, numbodies}]]

event = Table[{WhenEvent[
x1[j][t] == 0, {x1[j]'[t] -> 0, y1[j]'[t] -> 0, z1[j]'[t] -> 0}],
WhenEvent[
y1[j][t] == 0, {x1[j]'[t] -> 0, y1[j]'[t] -> 0, z1[j]'[t] -> 0}],
WhenEvent[
z1[j][t] == 0, {x1[j]'[t] -> 0, y1[j]'[t] -> 0,
z1[j]'[t] -> 0}]} /. j -> i, {i, numbodies}];

tfin = 15
sol1 = NDSolve[{eqs, event}, vars, {t, 0, tfin}][[1]]

plotXZ = ContourPlot[sol[x, 0.75, z], {x, 0, 2}, {z, -0.4, 0.1},
ContourShading -> Automatic, ColorFunction -> "Rainbow",
Contours -> 100];

frames = Table[
Show[v, ParametricPlot3D[
Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j,
numbodies}], {t, 0, tf}, PlotRange -> Automatic,
Axes -> False],
Graphics3D[
Table[{Hue[.35],
Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.0005]}, {j,
numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}];
ListAnimate[frames]
$$$$

• Your initial conditions are lists with 1 element, but numbers are expected. E.g. instead of x1[j][0] == pos0[[j, 1, ]] write x1[j][0] == pos0[[j, 1, 1]] Jun 28, 2022 at 18:17
• @DanielHuber That was it! And of course it was something simple that I missed lol. Thank you.
– Zach
Jun 28, 2022 at 18:38
• @DanielHuber How do I mark a comment as the answer? I am new and apparently can't figure it out.
– Zach
Jun 28, 2022 at 18:50
• You can not. I will write an answer. Jun 28, 2022 at 18:52