Jason Cole's blog is titled "Almost looks like work", and I have to agree - I could even go so far as to title it "You know you want to put off your real work and recreate this MATLAB project"
So your code was good, but I found myself wanting some sort of units system, otherwise you could have $v$ greater than $c$ (you did not, but I just wanted to make it explicit). The only parameter mentioned in the blog post is $\gamma$, which determines the velocity, so that's how I wrote the code. Also, in your code for the field, there are a couple of places where you take $\hat{\mathbf{n}}$ evaluted at $t_0$ instead of the retarded time. But you are correct that the main issue is the plotting scale.
We will set the speed of light and the radius of the circle to one, and measure the field in units of $q/(4 \pi \epsilon_0)$. This means the velocity is determined by the Lorentz factor, $\gamma$.
γ = 1.2;
w = Sqrt[γ^2 - 1]/γ;
tPeriod = 2 π/w;
xt[t_] := {Cos[w t], Sin[w t]};
β[t_] := w {-Sin[w t], Cos[w t]};
βprime[t_] := -w^2 xt[t];
{xmin, xmax} = {ymin, ymax} = {-8.05, 8.05};
I'm going to make the plots out of 2D lists, and I want to try and avoid evaluating the field at the same point in space where the electron is, so I offset the x
and y
grids by a small amount.
The retarded time must be calculated numerically. I wish I knew a faster or more reliable way to do this - the FindRoot
returns an error every once in a while, but not always at the same place. But it seems to give a good answer, and memoization will help later. The goal here is to make animations, so it is good to decide how many frames are needed. I want two full revolutions, and I tried using a timestep of 0.05 tPeriod
but found that to be too jerky so I went for 0.01 tPeriod
. You can speed things up by using a sparser time or spatial grid.
ClearAll@retime;
retime[t0_, x0_] :=
retime[t0, x0] =
tr /. FindRoot[t0 - Norm[x0 - xt[tr]] - tr, {tr, t0},
MaxIterations -> 1000];
So now we define the electric field and the radial component of the Poynting vector,
eField[t0_, x0_] :=
With[{
n = Chop@Normalize[x0 - xt[retime[t0, x0]]]~PadRight~3,
βvec = Chop@β[retime[t0, x0]]~PadRight~3,
βprimevec = Chop@βprime[retime[t0, x0]]~PadRight~3,
r = Norm[x0 - xt[retime[t0, x0]]]
},
Which[retime[t0, x0] < 0,
PadRight[(x0 - xt[0])/Norm[x0 - xt[0]]^3, 3]
, True,
(n - βvec)/(γ^2 r^2 (1 - βvec.n)^3) +
Cross[n, Cross[n - βvec, βprimevec]]/(
r (1 - βvec.n)^3)]];
poynting[t0_, x0_] := With[{
n = Chop@Normalize[x0 - xt[retime[t0, x0]]]~PadRight~3,
βvec = Chop@β[retime[t0, x0]]~PadRight~3,
βprimevec = Chop@βprime[retime[t0, x0]]~PadRight~3,
r = Norm[x0 - xt[retime[t0, x0]]]
},
Which[retime[t0, x0] < 0, 0,
True,
r^-2 Norm[Cross[n, Cross[n - βvec, βprimevec]]/(
r (1 - βvec.n)^3)]^2
]
];
It should be possible to save a bit of time above, by storing the field as two separate parts, the velocity and radiation fields, and then defining the Poynting vector as a function of the latter.
Now we generate lists for the retarded time, the field, and the Poynting vector,
Monitor[
Quiet[trlist1 =
Table[retime[t0 tPeriod, {x0, y0}], {t0, 0, 2, .04}, {y0, ymin,
ymax, .1}, {x0, xmin, xmax, .1}];], {x0, y0, t0}]
Monitor[elist1 =
Table[eField[t0 tPeriod, {x0, y0}], {t0, 0.00, 2, .01}, {y0, ymin,
ymax, .1}, {x0, xmin, xmax, .1}];, {x0, y0, t0}]
Monitor[plist1 =
Table[poynting[t0 tPeriod, {x0, y0}], {t0, 0.00, 2, .01}, {y0, ymin,
ymax, .1}, {x0, xmin, xmax, .1}];, {x0, y0, t0}]
The easiest to plot is the retarded time. No nonlinear scaling needed here, You just want to replace any negative values with 0, and then plot it using the MATLAB color map Parula, a visually appealing palette that is much better than the old Jet palette. I have this palette defined in a pastebin, which is what the first line below is,
<<"http://pastebin.com/raw.php?i=sqYFdrkY";
trplot[n_] :=
With[{rsdata =
Rescale[trlist1[[n]] /. {x_?Negative -> 0.0}, {0,
Max@trlist1}]},
Show[
ListDensityPlot[rsdata, PlotRange -> All,
ColorFunction -> ParulaCM,
DataRange -> {{xmin, xmax}, {ymin, ymax}}, Frame -> None],
Graphics@{Red, Point[xt[(n - 1) .01 tPeriod]]}
]];

Let's plot the x component of the electric field for one timestep (the 75th for no particular reason). Here it is using linear scaling with an automatic cutoff, the log scaling used in the blog, and an ArcSinh
scaling I use sometimes.
GraphicsRow[{ListDensityPlot[#, PlotRange -> Automatic,
ColorFunction -> ParulaCM, Frame -> None],
ListDensityPlot[Log@Abs@#, PlotRange -> All,
ColorFunction -> ParulaCM, Frame -> None],
ListDensityPlot[
Rescale[ArcSinh[10000 # /(Max@Abs@#)]/ArcSinh[10000], {-1,
1}] &@#, PlotRange -> All, ColorFunction -> ParulaCM, Frame -> None]} &@
elist1[[75,All,All,1]], ImageSize -> 700]

The log scale is the best visually, so we'll go with it. Now, when making an animation, it's very important to use the same scale with each frame. For these plots, since we sample the field on a grid, sometimes we have the electron very close to a grid point, and if you just let ListDensityPlot
choose the scale, your animation will have flashes of brightness, very offputting.
The rescaling values I have below I arrived at via trial and error - just trying to make the plots look best. If someone can think of a more programmatic way to do it, I'd be happy to hear it.
exloglist = Rescale[Log@Abs@elist1[[All,All,All,1]], {-9, 3}];
eyloglist = Rescale[Log@Abs@elist1[[All,All,All,2]], {-9, 3}];
ploglist = Rescale[Log@Abs@plist1, {-20, 6}];
trplot[n_] :=
With[{rsdata =
Rescale[trlist1[[n]] /. {x_?Negative -> 0.0}, {0, Max@trlist1}]},
Show[
ListDensityPlot[rsdata, PlotRange -> All,
ColorFunction -> ParulaCM,
DataRange -> {{xmin, xmax}, {ymin, ymax}}, Frame -> None],
Graphics@{Red, Point[xt[(n - 1) .01 tPeriod]]}
]]; explot[n_] :=
ListDensityPlot[exloglist[[n]], ColorFunction -> ParulaCM,
DataRange -> {{xmin, xmax}, {ymin, ymax}},
ColorFunctionScaling -> False, Frame -> False];
eyplot[n_] :=
ListDensityPlot[eyloglist[[n]], ColorFunction -> ParulaCM,
DataRange -> {{xmin, xmax}, {ymin, ymax}},
ColorFunctionScaling -> False, Frame -> False];
pplot[n_] :=
Show[ListDensityPlot[ConstantArray[0, {2, 2}],
DataRange -> {{xmin, xmax}, {ymin, ymax}},
ColorFunction -> ParulaCM, ColorFunctionScaling -> False,
Frame -> False],
ListDensityPlot[ploglist[[n]], ColorFunction -> ParulaCM,
DataRange -> {{xmin, xmax}, {ymin, ymax}},
ColorFunctionScaling -> False, Frame -> False]
];
gridplot[n_] := Grid[{{pplot[n], explot[n]}, {trplot[n], eyplot[n]}}];
gridplot[137]

Now to create the animation, you need to export the frames as image files and use another program to do the work. Mathematica is great, but for some things you should use other tools.
Quiet@CreateDirectory["syncimages"];
Quiet@CreateDirectory["syncimages/grd_1.2"];
Do[
img = gridplot[n];
Export[
"syncimages/grd_1.2/frame_" <> IntegerString[n, 10, 3] <> ".png",
img], {n, 76, Length@trlist1}];~Monitor~n
Then navigate to the folder "syncimages" in the command line and either create an mp4 video using ffmpeg,
ffmpeg -framerate 30 -i "grd_1.2/frame_%03d.png" -codec:v libx264 -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" -r 30 -pix_fmt yuv420p grd_1.2.mp4
or create an animated gif using ImageMagick
convert -delay 3 grd_1.2/* g12B.gif
This website won't let me link to a .gifv file, so here is the animation on imgur.
edit: Gonna try for the undulating example overnight :-)
ListDensityPlot
command? $\endgroup$