# ContourPlot of a discontinuous function

rainbow = {"Rainbow", {1.7, 3.4}};
color = ColorData[rainbow];
cf = Function[{z}, Directive[Opacity[0.8], ColorData[rainbow][z]]];
plot = {AspectRatio -> 1,
ImageSize -> {400, 400},
PlotRange -> {{-1, 1}, {-1, 1}},
FrameLabel -> {x, y},
LabelStyle -> Black,
FrameStyle -> Directive[Thickness[0.002], FontSize -> 14],
Mesh -> None, MeshStyle -> Gray,
InterpolationOrder -> 0,
ColorFunctionScaling -> False, ColorFunction -> cf};

g[x_, y_, z_] := Piecewise[{{(x - 1/2)^2 + y^2 + z^2 - 1/16, (x - 1/2)^2 + y^2 - 1/16 <= 0 && z >= 0}}, z];
(*t is a function of x1 and x2, and it is a discontinuous function as shown in the Figure 1*)
t[x1_, x2_] := Module[
{center, line, surface, intersection, interPoint, d, data},
center = {0, 0, 1};
line = Line[{{x1, x2, 0}, center}];
surface = ImplicitRegion[g[x, y, z] == 0., {x, y, z}];
intersection = MeshCoordinates@DiscretizeRegion@RegionIntersection[line, surface];
interPoint = Nearest[intersection, center][[1]];
d = EuclideanDistance[center, interPoint];
time = 2 d;
time
];
data = Table[t[x, y], {y, -1, 1.0, 0.05}, {x, -1, 1.0, 0.05}];
ListDensityPlot[data, DataRange -> {{-1, 1}, {-1, 1}}, Evaluate@plot]


You don't need to read Module carefully, t[x,y] is only a function of x and y and it is discontinuous at the region as shown in Figure.

Then I want to plot some contourplot lines. For example:

ContourPlot[t[x, y] == 2.3, {x, -1, 1}, {y, -1, 1}, PlotRange -> {{-1, 1}, {-1, 1}}, PlotPoints -> 20, MaxRecursion -> 1]


and it looks like Figure 2(a). But what I want is Figure2(b), because of the discontinuous region.

1. How to modify the code to plot Figure2(b)?
2. There is no problem when I run the ListDensityPlot , however, there is a warning when I run ContourPlot:

Nearest::near1: "MeshCoordinates[DiscretizeRegion[RegionIntersection[ImplicitRegion[Piecewise[{{<<2>>}},z]==0.,{x,y,z}],Line[{{x,y,0},{0,0,1}}]]]] is neither a list of real points nor a valid list of rules"

Why?

1. In the ListDensityPlot If the code is

data = Table[t[x, y], {x, -1, 1.0, 0.05}, {y, -1, 1.0, 0.05}];

The plot is rotated 90 degress with what it should be. I have to change the code to

data = Table[t[x, y], {y, -1, 1.0, 0.05}, {x, -1, 1.0, 0.05}];


Why?

Problem 2 and 3 are not important. They didn't have an effect on the results, I am just curious about it.

This is my solution:

Show[ContourPlot[t[x, y] == 1.8, {x, 0.3, 0.7}, {y, -0.2, 0.2}, ContourStyle -> color[1.8], PlotPoints -> 20, MaxRecursion -> 0,  PlotRange -> {{-1, 1}, {-1, 1}}],
ContourPlot[t[x, y] == 1.9, {x, 0.2, 0.8}, {y, -0.3, 0.3},  ContourStyle -> color[1.9], PlotPoints -> 20, MaxRecursion -> 0],
ContourPlot[t[x, y] == 2.0, {x, 0.2, 0.8}, {y, -0.3, 0.3}, ContourStyle -> color[2.0], PlotPoints -> 20, MaxRecursion -> 1, Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.1, {x, -0.33, 0.85}, {y, -0.33, 0.33},  ContourStyle -> color[2.1], PlotPoints -> 20, MaxRecursion -> 1,  Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.2, {x, -0.46, 0.9}, {y, -0.46, 0.46},  ContourStyle -> color[2.2], PlotPoints -> 20, MaxRecursion -> 0,  Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.3, {x, -0.57, 0.9}, {y, -0.57, 0.57}, ContourStyle -> color[2.3], PlotPoints -> 20, MaxRecursion -> 0,  Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.4, {x, -0.67, 0.9}, {y, -0.67, 0.67},  ContourStyle -> color[2.4], PlotPoints -> 20, MaxRecursion -> 0, Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.5, {x, -0.76, 0.9}, {y, -0.76, 0.76}, ContourStyle -> color[2.5], PlotPoints -> 20, MaxRecursion -> 0,  Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}],
ContourPlot[t[x, y] == 2.6, {x, -0.84, 0.84}, {y, -0.84, 0.84}, ContourStyle -> color[2.6], PlotPoints -> 20, MaxRecursion -> 0, Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}]]


I define the color=ColorData[{"Rainbow",{1.7,2.7}}]

• It's quite okay & normal to post an answer to your own question. (IMO, it's better than including an answer in the question itself). May 2, 2016 at 23:18

A simply way to obtain the desired plots is to excluding portions of the curves in the discontinuous region. Precisely where t is discontinuous does not seem well defined, so I have taken it to be approximately the almost circular region shown in the first plot in the question. Then,

nodiscontinuity = ContourPlot[{t[x, y] == 2.3, t[x, y] == 2.5, t[x, y] == 2.7},
{x, -1, 1}, {y, -1, 1}, PlotRange -> {{-1, 1}, {-1, 1}}, PlotPoints -> 20,
MaxRecursion -> 1,
RegionFunction -> Function[{x, y, z}, (x - .57)^2 + (y - .01)^2 > .08]]


Superimposing this plot on the result of

ListDensityPlot[data, DataRange -> {{-1, 1}, {-1, 1}}, Evaluate@plot];

Show[%, nodiscontinuity]


then shows that the curves indeed lie just outside the apparent discontinuous region.

(Note that plot (b) in the Question does not lie entirely outside the discontinuous region.)

In answer to the other two questions, the error messages occur because ContourPlot is calling t with symbolic input before calling it with numerical input. And, replacing

data = Table[t[x, y], {y, -1, 1.0, 0.05}, {x, -1, 1.0, 0.05}];


by

data = Table[t[x, y], {x, -1, 1.0, 0.05}, {y, -1, 1.0, 0.05}];


rotates the plot produced by ListDensityPlot, because the second line of code produces the Transpose of the array produced by the first line of code. This can be undone by

ListDensityPlot[data // Transpose, DataRange -> {{-1, 1}, {-1, 1}}, Evaluate@plot]


reproducing the original, not rotated plot.

In comments below, Qi Zhong asked whether his creative solution above could be computed more rapidly. This can be accomplished by calling ContourPlot fewer times. However, the 2.1 contour still must be computed separately, because it requires much higher resolution. On my PC, the final plot in the Question require 230 sec, but

Show[ContourPlot[t[x, y], {x, -1, 1}, {y, -1, 1}, PlotPoints -> 20,
Contours -> DeleteCases[Range[1.8, 2.6, 0.1], 2.1], ColorFunction -> (White &),
Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}},
ContourStyle -> (Directive[color[#], Thickness[.005], Opacity[1]] & /@
Range[1.8, 2.6, 0.1]), MaxRecursion -> 1],
ContourPlot[t[x, y] == 2.1, {x, -0.33, 0.85}, {y, -0.33, 0.33},
ContourStyle -> color[2.1], PlotPoints -> 20, MaxRecursion -> 1,
Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}}]]


requires only 106 sec.

• Yeah, I know this can work. The problem is I need to find the range of x manually. And in fact, I have several lines to plot, such as t[x, y] == 2.4,2.5... I will consider your method. But, is there any better one? May 2, 2016 at 17:18
• I found the exact function of the edge, Exclusions -> {{4 (x + 2)^2 == 19 (x^2 + y^2 + 1), x >= 3/8}} and it works. The third problem: normally, we should write {x,-1,1,0.1} before {y,-1,1,0.1} in Table, right? But when I do this, it seems the plot is not what I think it should be. Then I have to write y first. Is there any way to solve the second problem. May 2, 2016 at 20:37
• Another thing, it will take about >5 minutes to run the code I add above. Is it normal or is there any way to make it faster? May 2, 2016 at 20:44
• @QiZhong Nice solution. It can be made faster, and I shall provide an approach when I return late this evening. May 2, 2016 at 21:53