4
$\begingroup$
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]

Figure 1

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.

Figure2

  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}}]]

enter image description here

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

$\endgroup$
  • 1
    $\begingroup$ 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). $\endgroup$ – Michael E2 May 2 '16 at 23:18
4
$\begingroup$

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]]

enter image description here

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.

enter image description here

(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.

Addendum

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Qi Zhong May 2 '16 at 17:18
  • $\begingroup$ 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. $\endgroup$ – Qi Zhong May 2 '16 at 20:37
  • $\begingroup$ 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? $\endgroup$ – Qi Zhong May 2 '16 at 20:44
  • $\begingroup$ @QiZhong Nice solution. It can be made faster, and I shall provide an approach when I return late this evening. $\endgroup$ – bbgodfrey May 2 '16 at 21:53
  • $\begingroup$ @QiZhong As requested, I have addressed your two comments in my updated answer. Best wishes. $\endgroup$ – bbgodfrey May 3 '16 at 4:40

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.