# How can I reliably plot ellipse-field streamlines?

I am trying to replicate some standard structures in the study of ellipse fields in Mathematica, and I'm having much more trouble than is really reasonable, so I'd like to put this on the floor here in the hope that someone here can point out some simpler method that I'm currently not seeing.

More specifically, I would like to have simple and reliable ways to produce streamline plots like the following,

from the paper

Geometry of phase and polarization singularities, illustrated by edge diffraction and the tides. M.V. Berry. M.S. Soskin, M.V. Vasnetsov (Eds.), Second International Conference on Singular Optics, Proc. SPIE, 4403, p. 1 (2001). Author eprint.

Basically, the deal is the following: you have some complex-valued vector field, say, $$\mathbf F(x,y) = (1+iby, i(1-x))$$ (where $b$ is a parameter, the main interesting points of which are $b\in\{-1,1,3\}$), and you see it as encoding an ellipse at each point much in the same way as for an electromagnetic wave (i.e. the ellipse is traced out by $\mathrm{Re}(\mathbf F(x,y)e^{-it})$ for $t\in[0,2\pi]$); you then look for the streamlines of the major and minor axis of these ellipses.

The major and minor axis of each ellipse can be obtained as in this question as the real and imaginary parts of $$\mathbf A+i\mathbf B = \frac{\sqrt{\mathbf E^*\cdot \mathbf E^*}}{\left|\sqrt{\mathbf E^*\cdot \mathbf E^*}\right|}\mathbf E,$$ and as a naive start, one can do fairly well by putting these in explicitly and asking for a StreamPlot:

Block[{F, A, B, R = 0.6},
F[x_, y_] = {1 + I b y, I (1 - x)};
A[x_, y_] =
1/Abs[Sqrt[F[x, y]\[Conjugate].F[x, y]\[Conjugate]]] Re[
Sqrt[F[x, y]\[Conjugate].F[x, y]\[Conjugate]] F[x, y]];
B[x_, y_] =
1/Abs[Sqrt[F[x, y]\[Conjugate].F[x, y]\[Conjugate]]] Im[
Sqrt[F[x, y]\[Conjugate].F[x, y]\[Conjugate]] F[x, y]];
Table[
StreamPlot[
A[x, y]
, {x, -R, R}, {y, -R, R}
, ImageSize -> 400
]
, {b, {-1, 1, 5}}]
]


However, I have a couple of problems here:

• I cannot control where the streamlines start and stop, and the defaults do some pretty ugly spacing and they start and eliminate streamlines in mid-flight, which I would like to avoid. For some reason, if I try to give it an explicit StreamPoints specification, Mathematica just hangs forever.

• The stream plots have a discontinuity on the negative $x$ axis, where the major axis $\mathbf A(x,y)$ discontinuously changes sign. This discontinuity is both:

• spurious, because ellipses are symmetric and it doesn't matter which way you choose the axis, but also
• unavoidable, because of the way these ellipse points have been set up: if you track the ellipse around the origin, it will come back to itself while acquiring a 180° rotation; this won't affect the ellipse but it means that $\mathbf A(x,y)$ must change direction at some point.

I would like to have a solid construction that allows me to get past both of these limitations, but after several attacks I can't really get there. I'm OK with NDSolveing my way through this, but it keeps choking on the discontinuity and with the accumulation point at the origin.

(Also, I am wary of solutions that place explicit stock in the symmetry of the problem, by e.g. splitting up the plane into different regions. It's not a deal killer, but I would like to be able to handle cases where the branch cut falls along some non-symmetrically-placed line.)

Does anyone see how this can be done?

• For a different way of visualizing general polarization states, see also Visualizing vector spherical harmonics - that's not an answer to your question, unfortunately. – Jens Aug 22 '18 at 17:08

# Use NDSolve to solve for y(x)

Perhaps this serves as a starting point. Shown below is the case when b = -1.

Each streamline is obtained from NDSolve. Initial conditions in the form of $y(x_n) = y_n$ are chosen with a simple lattice, points (shown as red dots).

Not all systems can be solved, at least when $y(x_n < 0) = 0$, so the outputs from NDSolve that are not results of successful solution are filtered out by Cases. But then not all InterpolatingFunction are valid either; those that aren't have equal domain endpoints. They are filtered out by Select.

Clear[r, b, f, ff, x, y, z, vectorfield, points, plots, sols];

r = 6/10;
b = -1;
f[x_, y_] := {1 + I b y, I (1 - x)};
ff[x_, y_] := Sqrt[Conjugate[f[x, y]] . Conjugate[f[x, y]]];
z[x_, y_] := f[x, y] ff[x, y] / Abs[ff[x, y]];

vectorfield = Re @ z[x, y[x]];
points = Tuples[Subdivide[-r, r, 5], 2];

sols = NDSolveValue[
{y'[x] == Divide @@ Reverse @ vectorfield, y[First @ #] == Last @ #},
y,
{x, -r, r}
] & /@ points;

plots = Composition[
Select[#, Unequal @@ First[#["Domain"]] &] &,
Cases[#, _InterpolatingFunction] &
] @ sols;

Show[
{
ListPlot[points, PlotStyle -> Red],
Plot[#[x], Evaluate @ {x, Sequence @@ First[#["Domain"]]}] & /@ plots
},
Axes -> False,
Frame -> True,
AspectRatio -> 1,
PlotRange -> {{-r, r}, {-r, r}},
]


The problem really is how to construct points such that you have beautifully spaced streamlines.

# Specify Better Sample Points

Since this is an aesthetic issue, human intervention is going to be necessary, on a case-by-case basis.

For example, with

points = Composition[
Map[First],
Cases[#, _Point] &,
RegionIntersection /@ # &,
Subsets[#, {2}] &
] @ {
InfiniteLine[{{-r, 0}, {r, r}}],
InfiniteLine[{{-r, 0}, {r, -r}}],
Sequence @@ Table[InfiniteLine[{{n, -r}, {n, r}}], {n, -r, r, 1/10}]
}


we can obtain this:

Now I can produce other plots too. When b = 3 (monstar):

When b = 1 (lemon):

# Downsides

Note though that the approach finds $y(t)$ through $y(x(t))$, which is a workaround for not knowing the bounds of $t$ within the domain and range of interest. One downside is that if a particular streamline isn't a function $y(x)$ in a domain of interest, NDSolve only goes as far as it is a function and stops. For example, with

vectorfield = {-1 - x^2 + y, 1 + x - y^2} /. y -> y[x]
points = {{-2, 1/2}}


we get this:

In such cases, we have to pick up from the terminal point somehow, in order to continue the streamline towards its accumulation point or the frame edge.

For example, after obtaining plots[[1]] above, we could slightly shift the terminal point (in an appropriate direction), feed it back to points:

terminus = Last @ Last @ plots[[1]]["Domain"]
points = Append[{terminus, plots[[1]][terminus] + 1/1000}] @ points


and solve again to obtain

It's easy to specify points that build continuous streamlines for the OP's problem. This is not easy in general.

Another related downside of the approach is that it doesn't work if some parts of the streamlines within the frame of interest is parallel to the y-axis.

# Solve for x(t) and y(t) instead

To address both downsides, we could switch to the conventional way, i.e. solving for $y(t)$ and $x(t)$. We may designate $t = 0$ to be where our sample points are. But because we don't know the appropriate bounds of $t$ to specify for NDSolve, we may try some small range at first, e.g. $t \in [-1/10, 1/10]$.

vectorfield = Re @ z[x[t], y[t]];
points = Composition[
Map[First],
Cases[#, _Point] &,
RegionIntersection /@ # &,
Subsets[#, {2}] &
] @ {
InfiniteLine[{{-r, 0}, {r, r}}],
InfiniteLine[{{-r, 0}, {r, -r}}],
Sequence @@ Table[InfiniteLine[{{n, -r}, {n, r}}], {n, -r, r, 1/10}]
};
sols = NDSolveValue[
{
x'[t] == First @ vectorfield,
y'[t] == Last @ vectorfield,
x[0] == First @ #,
y[0] == Last @ #
},
{x, y},
{t, -1/10, 1/10}
] & /@ points;
plots = Cases[sols, {_InterpolatingFunction, _InterpolatingFunction}];
Show[
{
ListPlot[points, PlotStyle -> Red],
ParametricPlot[
{#[[1]][t], #[[2]][t]},
Evaluate @ {t, Sequence @@ First[#[[1]]["Domain"]]}
] & /@ plots
},
Axes -> False,
Frame -> True,
AspectRatio -> 1,
PlotRange -> {{-r, r}, {-r, r}},
]


from which we obtain

which tells us that perhaps $t \in [-1, 1]$ (i.e. putting {t, -1, 1} in NDSolve) should suffice to fill up the frame. And indeed, it leads to the previous plot above.

To illustrate that the downsides have been eliminated, consider a wholly vertical vector field

vectorfield = {0, 1}
points = Table[{n, 0}, {n, -r, r}]


from which we obtain

The approach should be robust enough for most well-behaved vector fields now. For example,

r = 3;
vf = {y, -y + x - x^3};
vectorfield = vf /. {x -> x[t], y -> y[t]};
points = Composition[
Fold[If[EuclideanDistance[Last[#1], #2] < 1, #1, Append[#1, #2]] &],
Level[#, {-2}] &,
MeshPrimitives[#, 0] &,
DiscretizeRegion,
ImplicitRegion[Norm[#] == 3, {x, y}] &
] @ vf;
sols = NDSolveValue[
{
x'[t] == First @ vectorfield,
y'[t] == Last @ vectorfield,
x[0] == First @ #,
y[0] == Last @ #
},
{x, y},
{t, -10, 10}
] & /@ points;
plots = Cases[sols, {_InterpolatingFunction, _InterpolatingFunction}];
Show[
{
StreamPlot[vf, {x, -r, r}, {y, -r, r}, StreamStyle -> LightGray],
ListPlot[points, PlotStyle -> Red],
ParametricPlot[
{#[[1]][t], #[[2]][t]},
Evaluate @ {t, Sequence @@ First[#[[1]]["Domain"]]}
] & /@ plots
},
Axes -> False,
Frame -> True,
AspectRatio -> 1,
PlotRange -> {{-r, r}, {-r, r}},