How to draw symmetrical parts of a circle

Question

Let's consider the following simple case:

V = 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
h = 0.3;
rad = 6*h;

S1 = ContourPlot[V, {x, -2, 2}, {y, -2, 2}, Contours -> {h}, 
ContourShading -> False, ContourStyle -> {Darker[Green], Thick}, 
PlotPoints -> 100, PerformanceGoal :> "Speed", AspectRatio -> 1];
S2 = ContourPlot[x^2 + y^2, {x, -2, 2}, {y, -2, 2}, Contours -> {rad},
ContourShading -> False, ContourStyle -> {Blue, Thick, Dashed}, 
PlotPoints -> 100, PerformanceGoal :> "Speed", AspectRatio -> 1];
P0 = Show[{S1, S2}, FrameLabel -> {"x", "y"}, RotateLabel -> False, 
FrameStyle -> Directive[FontSize -> 20, FontFamily -> "Helvetica"], 
PlotRange -> 1.5, ImageSize -> 550]

which gives

We see, that the blue dashed circle intersects the green contour at the three openings (exits). I would like to draw the symmetrical parts of the circle inside these exits. By the term symmetrical, I mean with reverse curvature. I added by hand and with red color in the plot the symmetrical parts.

Any ideas how to obtain this?

EDIT

A copy of the above code (in version 9) can be found here:code.

Answer 1

First we need to find the intersection points. It's more convenient to find them as angles around the circle, rather than $(x,y)$ coordinates, otherwise we might not be sure what order they're in.

r = Sqrt[rad];
angles = θ /. Solve[{(V /. {x -> r Cos[θ], y -> r Sin[θ]}) == h, 0 <= θ < 2 π}, θ]

{1.32755, 1.81404, 3.42195, 3.90843, 5.51634, 6.00283}

By inspection we see that the arcs we want are from the first to the second, then the third to the fourth, and the fifth to the sixth point. With a little geometry, we can define a function to draw the reflected arc between any two points on a circle.

reflectedArc[{x_, y_}, r_, {θ1_, θ2_}] :=
 Circle[{x, y} + r {Cos[θ1] + Cos[θ2], Sin[θ1] + Sin[θ2]}, r, {π + θ1, π + θ2}]

We apply this onto our pairs of points and produce the final figure.

arcs = Graphics[{Red, reflectedArc[{0, 0}, r, #] & /@ Partition[angles, 2]}];
Show[{S1, S2, arcs}, FrameLabel -> {"x", "y"}, RotateLabel -> False, 
 FrameStyle -> Directive[FontSize -> 20, FontFamily -> "Helvetica"], 
 PlotRange -> 1.5, ImageSize -> 550]

(By the way, I must point out that you could have defined S2 much more easily as Graphics[{Blue, Thick, Dashed, Circle[{0, 0}, Sqrt[rad]]}].)

Answer 2

This solution is different from the other solutions in that it does not use the fact that the blue curve is a circle. Instead, we search the interior of the blue curve for points which will lie on the boundary when reflected.

I recast your problem statement as

f1[x_, y_] := 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
f2[x_, y_] := x^2 + y^2;
h = 0.3; rad = 6*h;

S1 = ContourPlot[f1 @@ {x, y}, {x, -2, 2}, {y, -2, 2}, 
   Contours -> {h}, ContourShading -> False, 
   ContourStyle -> {Darker[Green], Thick}, PlotPoints -> 100, 
   PerformanceGoal :> "Speed", AspectRatio -> 1];
S2 = ContourPlot[f2 @@ {x, y}, {x, -2, 2}, {y, -2, 2}, 
   Contours -> {rad}, ContourShading -> False, 
   ContourStyle -> {Blue, Thick, Dashed}, PlotPoints -> 100, 
   PerformanceGoal :> "Speed", AspectRatio -> 1];

Then I found the points of intersection, and lines to be reflected across. This is a weak point in the solution, as it depends on manually identifying the lines to reflect across.

pts = {x, y} /. Solve[{f2[x, y] == rad, f1[x, y] == h}, {x, y}, Reals]; // Quiet
lines = Partition[pts, 2];

Define reflections geometrically (Projection is built-in, that was new to me)

Rejection[u_, v_] := u - Projection[u, v];
Reflection[{a_, b_}, c_] := a + Projection[c - a, b - a] - Rejection[c - a, b - a]

Then use my Reflection and the RegionFunction option

S3 = ContourPlot[f2 @@ Reflection[#, {x, y}], {x, -2, 2}, {y, -2, 2}, 
   RegionFunction -> Function[{x, y}, f2[x, y] < rad], 
   Contours -> {rad}, ContourShading -> False, 
   ContourStyle -> {Red, Thick, Dashed}, PlotPoints -> 100, 
   PerformanceGoal :> "Speed", AspectRatio -> 1] & /@ lines;

And finally combine plots

P0 = Show[{S1, S2, S3}, FrameLabel -> {"x", "y"}, 
  RotateLabel -> False, 
  FrameStyle -> Directive[FontSize -> 20, FontFamily -> "Helvetica"], 
  PlotRange -> 1.5, ImageSize -> 550]

Here is a picture with your function, as well as one with f2[x_,y_]:=x^4+y^6

Answer 3

Let me join this party.

I propose a purely graphical approach (it looks like quite compact):

f1[x_, y_] := 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
f2[x_, y_] := x^2 + y^2;
h = 0.3; rad = 6 h;

Show[Normal@ContourPlot[f2[x, y] == rad, {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
    ContourStyle -> {Red,Thick}, RegionFunction -> (f1[#1, #2] < h &), MaxRecursion -> 4] /.
     Line[{p1_, pp___, p2_}] :> GeometricTransformation[Line[{p1, pp, p2}], 
       ReflectionTransform[Cross[p2 - p1], p2]], 
 ContourPlot[f1[x, y] == h, {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
  ContourStyle -> {Darker[Green], Thick}]]

Here Normal converts GraphicsComplex to Lines. Then I reflect this lines around ending points.

Answer 4

We can try as the following.

V = 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
h = 0.3; rad = 6*h; center = {0, 0};
pts = {x, y} /.NSolve[{x^2 + y^2 == 6*h, V == h}, {x, y},Reals];(*find intersections*)
midpts = (Mean /@ Partition[pts, 2]);(*Find mid points the black ones*)
dist =Norm /@ midpts;(*Distance from center*)
Off[Divide::infy];(*equations of straight line joining center and mid points*)
eqs = (y - #[[2]] == (x - #[[1]]) First@Ratios@(center - #) & /@ midpts) /. 
   a_ == ComplexInfinity -> (x == 0) /.ComplexInfinity == b_ -> (y == 0);
(*pattern check for straight lines x=0 and y=0*)
Off[NSolve::ratnz];(*Find equation of the reverse circles*)res = 
 MapThread[Total[({x, y}-First@(Select[{x, y} /.Chop@NSolve[{#1, Norm[#2 - {x, y}] == #3},
   {x, y},Reals], # != center &]))^2] == rad &, {eqs, midpts,dist}];

We find collinear point on the line joining the center and the mid point (in black) lying at the same distance from the mid points as the distance from the midpoint to the center. In res we have the equations of the reverse circles centered at those points. You can ContourPlot these circle equations now on a right x-domain to get the right segment. Red and black points are drawn on you last Show command using Epilog.

reverse = 
  MapThread[
   ContourPlot[#1, Evaluate@Flatten@{x, #2[[1]]}, {y, -2, 2}, 
     Contours -> {rad}, ContourShading -> False, 
     ContourStyle -> {Red, Dashed}, PlotPoints -> 100, 
     PerformanceGoal :> "Speed", AspectRatio -> 1] &, {res,Transpose /@ Partition[pts, 2]}];

Answer 5

This is not particularly clever but works:

V = 1/2*(x^2 + y^2) + x^2*y - 1/3*y^3;
cp = ContourPlot[{V == 0.3, x^2 + y^2 == 1.8}, {x, -2, 2}, {y, -2, 2}];
sol = NSolve[{V == 0.3, x^2 + y^2 == 1.8}, {x, y}];
roots = {x, y} /. sol
rem[u_] := Module[{f}, f[x_] := (f[Reverse@x] = Sequence[]; x); f /@ u];
ans = rem[Nearest[roots, #, 2] & /@ roots];
cent[{{p_, q_}, {r_, s_}}] := 
 First@Cases[{x, y} /. 
    Chop[NSolve[{(x - p)^2 + (y - q)^2 == 
        1.8, (x - r)^2 + (y - s)^2 == 1.8}, {x, y}]], Except[{0, 0}]];
vec = {#, cent[#]} & /@ ans;

arcs = {#[[2]], 
    Arg[Complex @@@ (Map[Function[x, x - #[[2]]], #[[1]]])]} & /@ vec;
Show[cp, Graphics[Circle[#[[1]], Sqrt[1.8], #[[2]]] & /@ arcs]]

Essentially, just find the intersections, second circle centre and arcs then compose. Styling could be altered.

Answer 6

This is like Rahul's solution only imagine that we didn't realize that we could use trigonometry to find the origin of the reflected arc analytically, and we also didn't realize polar coordinates were more suited for the task hehe.

(* Find the points of intersection. *)
sol = Solve[{
    V == h,
    x^2 + y^2 == rad
    }, {x, y}, Reals];

(* Find three sets with two points in each. Group by how close they \
are to each other. *)
pts = FindClusters[{x, y} /. sol, 3];

(* Find the origin of the other circle that intersects each set of \
coordinates *)
origos = {x, y} /. First@FindInstance[{
                Norm[{x, y} - #[[1]]] == Norm[{x, y} - #[[2]]],
                Norm[{x, y} - #[[1]]] == Sqrt[rad],
                Norm[{x, y}] > Sqrt[rad]
                }, {x, y}, Reals] & /@ pts;

(* Create data list: {origin,intersection 1,intersection 2} *)

data = MapThread[{#, #2[[1]], #2[[2]]} &, {origos, pts}];

(* Generate the circle segments *)
circles = Circle[
          #[[1]], Norm[#[[1]] - #[[2]]], {
       ArcTan @@ (#[[2]] - #[[1]]),
       ArcTan @@ (#[[3]] - #[[1]])
      }
          ] & /@ data;

(* Draw. *)
Show[P0, Graphics[{Red, circles}] ]