I have a two-dimensional function that is periodic in both variables. We can then think of this function on the surface of a torus. Now imagine a line that endlessly spirals and wraps around the torus. I need to know the distribution of lengths that the line spends in regions where the function is above a set value. I will give some demo code and my current solution which works, but it does not scale well.

An example function T:

Clear[S, T];
S[f_, q_, i_, r_] := 
 Evaluate[RotationMatrix[i, {0, 0, 1}].RotationMatrix[
 f, {0, 1, 0}].RotationMatrix[q, {0, 0, 1}].{r, 0, 0}]
T[f_, q_, i_, r_, l_] := 
 ArcCos[(RotationMatrix[l, {0, 1, 0}].{1, 0, 0}).S[f, q, i, r]/Sqrt[
   S[f, q, i, r].S[f, q, i, r]]]

Simple color function to make a nice plot:

cf[z_] := Which[
z >= π/2, ColorData["Rainbow"][2 z/π - 1]
GrayLevel[2 z/π]

A line and some constants:

 Clear[ft, gt];
 ft[t_] := t;
 gt[t_] := -.2 Sqrt[2.] t;

 i = .2 π;
 r = 1.02;
 l = .4 π;

 len = 107.;
 n = 200;
 pts = {Mod[ft[#], 2 π], Mod[gt[#], 2 π]} & /@ 
  Subdivide[0, len, 10 n];
 arrows = Select[
 Partition[{Mod[ft[#], 2 π], Mod[gt[#], 2 π]} & /@ 
 Subdivide[0, len, n], 2], Norm[#[[2]] - #[[1]]] <  π &];

Make the plot:

 ContourPlot[T[f, q, i, r, l], {f, 0, 2 π}, {q, 0, 2 π},
  ColorFunctionScaling -> False,
  ColorFunction -> cf,
  Contours -> Subdivide[0, π, 20],
  PlotRange -> {0, π},
  Epilog -> {
    {Red, Thick, Arrow[#] & /@ arrows}}

Lines over the regions

If you were to follow one of the lines as it wraps around, sometimes it is in the colored regions and sometimes it is not. I need to know the distribution of the lengths of the line segments that are in the colored region.

My current method:

It is easy to plot the value of the function along the line:

tmax = 100;
Plot[T[ft[t], gt[t], i, r, l], {t, 0, tmax}]

function along line

To find all of the sections that are above $\frac{\pi}{2}$, I do the following:

 pic = Plot[T[ft[t], gt[t], i, r, l], {t, 0, tmax},
  PlotRange -> {π/2, All},
  Axes -> False

sec = Select[Cases[pic, Line[a_] :> a, ∞], 
  Length[#] > 3 &];

sec = {First[#], Last[#]} & /@ sec;

I make a plot of the function with the rage going from $\frac{\pi}{2}$ up to All then use Case[] to find the data points in the Line[] commands. And finally just take the endpoints of the Line[] paths.

Take care of the possible end cases if the function was above $\frac{\pi}{2}$ at t=0 or t=tmax

 ϵ = 10^-6;
 sec = Select[sec, 
   Chop[#[[1, 2]] - π/2, ϵ] == 0 && 
   Chop[#[[2, 2]] - π/2, ϵ] == 0 &];

Plot that shows it acuraly works:

 Plot[T[ft[t], gt[t], i, r, l], {t, 0, tmax},
  PlotStyle -> Thick,
  Epilog -> {
    Line[#] & /@ sec

Function with line segments

Finally: finding the distibutions of the segment lengths:

 ends = sec[[All, All, 1]];

  dT = #[[2]] - #[[1]] & /@ ends;



The issue with this is that as tmax grows without bound the plotting method starts to give bad results even when I increase PlotPoints and the amount of memory grows to the point of crashing the sub-kernels.

Are there better Mathematica ways of doing this? I have played with the Region* functions but I can never get them to work consistently.


The MeshFunctions capability of Plot[] is one route you might want to explore.

An aside first: your S and T functions can be implemented quite compactly:

S[f_, q_, i_, r_] := EulerMatrix[{i, f, q}].{r, 0, 0};
T[f_, q_, i_, r_, l_] :=
  ArcCos[Normalize[S[f, q, i, r]].RotationMatrix[l, {0, 1, 0}].{1, 0, 0}]

From there,

ft[t_] := t
gt[t_] := -Sqrt[2] t/5

tmax = 100;
plt = With[{i = .2 π, r = 1.02, l = .4 π},
           Plot[T[ft[t], gt[t], i, r, l], {t, 0, tmax},
                Mesh -> {{π/2}}, MeshFunctions -> {#2 &}]]

ends = Partition[Sort[Cases[Normal[plt], Point[{x_, y_}] :> x, ∞]], 2];
dt = -(Subtract @@@ ends);

and running Histogram[dT] ought to yield the same figure. All of this assumes the number of crossing points is always even, of course.

  • $\begingroup$ Thanks. I like the idea of using Mash and MeshFunctions which allows the explicit setting of the value in question. But one still needs to test if the endpoints are the beginning of a segment with values greater than the value in question. One could check the derivative at the endpoints. However, this method now requires Plot to generate all the data above and below the value in question and still does not scale well for very large tmax where the Plot function starts missing some of the oscillations. $\endgroup$ – c186282 Sep 24 '17 at 15:17
  • $\begingroup$ "...where the Plot function starts missing some of the oscillations..." - that's where you try increasing PlotPoints. $\endgroup$ – J. M.'s torpor Sep 24 '17 at 17:55
  • $\begingroup$ Yes, I have but around PlotPoints of 90000 things go sour. I need a method that "moves along the function" as opposed to "calculate everything up front and take what I need". If nothing obvious comes up I'll just write my own that uses the Pot method and subdivides the domain into reasonable sizes. $\endgroup$ – c186282 Sep 24 '17 at 18:10
  • $\begingroup$ "subdivides the domain into reasonable sizes." - I was going to suggest that next, but the bookkeeping involved in stitching together pieces of domains seems annoying. $\endgroup$ – J. M.'s torpor Sep 24 '17 at 18:19
  • $\begingroup$ Yes. which was the motivation for the question. I'm hoping for some Mathematica magic. $\endgroup$ – c186282 Sep 24 '17 at 20:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.