Synchrotron Radiation and ListDensityPlot

Question

I'm trying to reproduce the following plots:

As described by this blog entry: Synchrotron radiation

As you can see in the linked article, the task is essentially to make a DensityPlot of a vector field. This vector field depends on the position of the observer $(x,y)$ and on the time $(t_r)$. This time has the peculiarity that it must be computed numerically from the real time using a simple equation:

$$t_r=t_0-\frac{ x_0-x(t_r)}{c}$$

Using this we can evaluate the electric vector field as:

$$\mathbf{E}(\mathbf{r},t)=\frac{q}{4\pi\varepsilon_0}\left[\frac{\hat{\mathbf{n}}-\vec{\beta}}{\gamma^2R^2(1-\vec{\beta}\mathbf{\cdot}\hat{\mathbf{n}})^3}+\frac{\hat{\mathbf{n}}\times[(\hat{\mathbf{n}}-\vec{\beta})\times\dot{\vec{\beta}}\,]}{c\,R\,(1-\vec{\beta}\mathbf{\cdot}\hat{\mathbf{n}})^3}\right]_{\mathrm{retarded}} \qquad \qquad (2)$$

In this expression, $\mathbf{\beta}$ is the velocity of the particle divided by the speed of light, $\gamma$ is the relativistic factor of the particle, R is the distance from the particle to the observer and $\mathbf{n}$ is the unit vector from the particle to the observer. The retarded label means that all the quantities are evaluated as a function of the retarded time. For example:

$$\beta=\beta(t_r(t))$$

where the function $t_r(t)$ comes from resolving the first equation in this text.

---

Here is the code I'm trying to use to produce the plots.

First, we define a function to calculate the retarded time from the normal time:

tr[t0_, x0_, y0_] := 
FindRoot[
    t0 - Norm[{x0, y0} - {x[tr], y[tr]}]/10 - tr, {tr, t0}][[1, 
   2]]

Next, we define the trajectory of the particle

R=5
x[t_] := R Cos[t]
y[t_] := R Sin[t]

We define auxiliar functions for $\beta$ and $n$:

n[x0_, y0_] := n[x0_, y0_, t_] := ({x[t], y[t], 0} - {x0, y0, 0})/
 Norm[{x[t], y[t], 0} - {x0, y0, 0}]
β[t0_] := {x'[t0], y'[t0], 0}

And finally, the electric field:

Efield[γ_, x0_, y0_,t0_] := ((n[x0, y0, t0] - β[
   tr[t0, x0, 
    y0]])/(γ^2 Norm[{x0, y0} - {x[tr[t0, x0, y0]], 
      y[tr[t0, x0, y0]]}]^2 (1 - β[tr[t0, x0, y0]].n[x0, y0,
        t0])^3) + Cross[n[x0, y0, t0], (  
  Cross[     (n[x0, y0, t0] - β[tr[t0, x0, y0]]), β'[
    tr[t0, x0, y0]]          ]     )]/(10 Norm[{x0, 
     y0} - {x[tr[t0, x0, y0]], 
     y[tr[t0, x0, y0]]}] (1 - β[tr[t0, x0, y0]].n[x0, y0, 
       t0])^3))

When I try to use Efield1, Efield[0] to make a ListDensityPlot the (wrong) following result:

So the question is the following:

Is there some way to generate the original plots using ListDensityPlot?

Maybe the problem is the scaling of the z-axis in the DensityPlot, but what I have tried so far does not work. For scaling I pass the following options to ListDensityPlot:

 ColorFunction -> Function[z, ColorData["DeepSeaColors"][z/10]],
 ColorFunctionScaling -> False

---

Edit 1:

1.Corrected a typo in the code for n[x0_, y0_] (changed Abs to Norm).

2.The code of the DensityPlot is:

   data = Flatten[Quiet@Table[{x0, y0,Efield[1, x0, y0, 34 \[Pi]/3 (*For example*)][[1]]}, {x0, -40, 40,1}, {y0, -40, 40, 1}],1]
  ListDensityPlot[data, InterpolationOrder -> 2, PlotRange -> All, 
 ColorFunction -> Function[z, ColorData["DeepSeaColors"][z/10]], 
 ColorFunctionScaling -> False]

---

Edit 2

It seems that is a problem with DensityPlot. The blue lines in the plots I wanted have much much lower values that the rest of the plot. So the problem is how to rescale these lines with the appropriate color function. Using ArcTanh as a scaling function shows the "dipole" lines, but no the pretty blue lines.

Answer 1

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 :-)