# How can I produce a branch cut in a surface plot based made with ListPlot3D?

I need to create some surface plots, based on a function that is computationally expensive (a few minutes per point). The function has some branch cuts.

The command Plot3D does not work because it just takes too long -- never finishes at any acceptable resolution. I can get enough data points (over a day or three), put them in a matrix, and use ListPlot3D to produce a plot. But I don't see how to make Mathematica recognize the branch cut and not connect the two sides at the discontinuity.

-
Have you considered using RegionFunction to exclude a narrow strip around the discontinuity? – Simon Woods Oct 30 '13 at 16:03
Hmm. Interesting idea. In my case, the branch cuts are not in any analytic form, so calculating them would be nontrivial. I did just think of another solution, making use of the fact that in my case, the branch cuts are closed loops. I could separate the data points "by hand" (values are greater or less than such and such) into two (or, in my case, three) separate matrices, and plot them with ListPlot3D on the same figure. Ad hoc, not elegant, and each figure would take some user effort, but I think doable. There should be a better way. – Richard Tasgal Oct 31 '13 at 7:41

One option is simply to create the ListPlot3D and post-process it to remove polygons which are too steep.

data = With[{g := GaussianFilter[RandomReal[{-1, 1}, {50, 50}], 25]}, ArcTan[g, g]];
plot = ListPlot3D[data, Mesh -> None] // Normal;
toosteep[Polygon[pts_, ___]] := Module[{z = pts[[All, 3]]}, Max[z] - Min[z] > 1]
DeleteCases[plot, _Polygon?toosteep, -1]


Another possibility is to use a RegionFunction to exclude regions where the function gradient is high. Here I use a bit of image processing to locate the branch cuts:

exclude = Image[data] ~GradientFilter~ 1 ~Binarize~ 1.5 ~Dilation~ 1


regionfunc = ListInterpolation[ImageData@exclude, InterpolationOrder -> 0];
ListPlot3D[data, Mesh -> None, RegionFunction -> Function[{x, y, z}, regionfunc[y, x] < 1]]


-
Thanks! I wonder if I can use the techniques in this library.wolfram.com/infocenter/MathSource/758 function to improve the accuracy of the location of the branch cuts (which in my case are smooth curves). (Apologies for not doing this myself and writing it here. I'm not that proficient in Mathematica. It will take me a few days to figure it out.) – Richard Tasgal Oct 31 '13 at 16:34
@RichardTasgal, you could, but this of course increases the number of evaluations of your expensive function as you iteratively refine the location of the cuts. This is why Plot3D takes so long - it is increasing the sampling density near the discontinuity. – Simon Woods Oct 31 '13 at 17:16
One other thing, if you are "not that proficient" in Mathematica it may be that your expensive function is coded inefficiently. What sort of computations are involved? – Simon Woods Oct 31 '13 at 17:18
The function is as follows: I have a 6th-order polynomial equation where the coefficients depend on a few variables, a0(k,a,b)+a1(k,a,b)y+a2(k,a,b)y^2+a3(k,a,b)y^3+a4(k,a,b)y^4+a5(k,a,b)y^5+a6(k,a,‌​b)y^6=0. The functions a_j(k,a,b) take a few tens of lines to define---involving functions of functions and a determinant---but they are algebraic expressions. Then I find the location k of the largest imaginary root of the 6 solutions y_j (Max,Imag,NSolve) over all k (NMaximize), holding a,b constant. I plot this over a,b. Compile doesn't help the speed. Maybe I'm missing a subtlety. – Richard Tasgal Nov 3 '13 at 8:16
@RichardTasgal, it sounds like its just a really big calculation. I'd be happy to have a look if you can post a link to the notebook, but it may be that several minutes per point is as fast as it gets. – Simon Woods Nov 3 '13 at 11:30