# Handling Self-Intersecting Regions

I am trying to test various parametric regions of a function for membership of the origin (point {0,0}). I generate the boundary by creating four parametric boundary curves, and then use RegionMember to tell whether the origin is inside the region (indicating a root of the function). However, the function often overlaps itself, as is shown in the result from my code below. The MeshBoundaryRegion function below gives an error of the form

BoundaryMeshRegion::binsect: The boundary curves self-intersect or cross each other in...


Consequently, RegionMember also throws an error. I also tried creating the region as a polygon, but I found that RegionMember also did not work with the resultant region. (It simply did not evaluate.) The code is below.

l = 1;
j = 4;
i = 2;
wint = Re[zGrid[[j, i + 0]]] Abs[Ws[snorm]];
wfin = Re[zGrid[[j, i + 1]]] Abs[Ws[snorm]];
yint = Im[zGrid[[j + 0, i]]] Abs[Ws[snorm]];
yfin = Im[zGrid[[j + 1, i]]] Abs[Ws[snorm]];
wres = 36;
yres = 36;
winc = (wfin - wint)/(wres - 1);
yinc = (yfin - yint)/(yres - 1);
w = Table[wint + (i - 1) winc, {i, wres}];
y = Table[yint + (i - 1) yinc, {i, yres}];

paraTab1 = Table[{Re[f[k[[l]], w[[i]], yint]], Im[f[k[[l]], w[[i]], yint]]}, {i, wres}];
paraTab2 = Table[{Re[f[k[[l]], wfin, y[[i]]]], Im[f[k[[l]], wfin, y[[i]]]]}, {i, yres}];
paraTab3 = Table[{Re[f[k[[l]], w[[wres - i + 1]], yfin]], Im[f[k[[l]], w[[wres - i + 1]], yfin]]}, {i, wres}];
paraTab4 = Table[{Re[f[k[[l]], wint, y[[wres - i + 1]]]], Im[f[k[[l]], wint, y[[yres - i + 1]]]]}, {i, yres}];

paraTabFullEx = Join[paraTab1, paraTab2, paraTab3, paraTab4]/norm;

SetOptions[Graphics,
ImageSize -> Large,
AxesOrigin -> {0, 0},
GridLines -> {{0}, {0}},
LabelStyle -> Directive[Black, FontFamily -> "Cambria", FontSize -> 20]
];

paraTabEx = Delete[paraTabFullEx, {{wres}, {wres + yres}, {2 wres + yres}, {2 wres + 2 yres}}];

Rg = BoundaryMeshRegion[paraTabEx, Line[Join[Table[n, {n, Length[paraTabEx]}], {1}]]];

paraTabLinEx = Line[Join[paraTabEx, {paraTab1[[1]]/norm}]];
Graphics[paraTabLinEx]

poly = Polygon[paraTabEx];
Graphics[poly]


Another, perhaps more visually comprehensible, example is below. (Please disregard the badly formatted tick labels.)

Notice that in this second example, the function does not pass through the origin; the small inner loop is not part of the region.

Is there a way to break the space traced out of the function into overlapping and non-overlapping components? Alternatively, is there a better approach that I should consider? Increasing the resolution of the function space sample has proven prohibitively computationally expensive.

Thank you

To actually replicate these results, the full code below is needed

(* Fundamental Constants *)
\[Mu]0 = 4 \[Pi] 10^-7; (* H m^-1 *)
c = 299792458; (* m s^-1 *)
ep0 = 1/(\[Mu]0 c^2); (* F m^-1 *)
eEl = (1.602176664*10^-19); (* C *)
evtj = 1.602177164 10^-19; (* J *)
me = 9.1093864 10^-31; (* kg *)
mH = 1.67262264 10^-27; (* kg *)
ZH = 1; (* C *)

(* Control Parameters *)
s = {"e", "i"};

ionSpec = "H";

mr = 1836;

ni = 1*10^18;(* m^-3 *)

TeeV = 185;(* eV *)
Tre = 1;
TieV = 50;(* eV *)
Tri = 10;

B = 0.01464; (* T *)

(* Derived Constants *)
Ns = Length[s];

mi = ToExpression["m" <> ionSpec] (mr/(mH/me));
Zi = ToExpression["Z" <> ionSpec];

qe = -eEl;
qi = Zi eEl;

ne = Zi ni;(* m^-3 *)

Te = TeeV evtj;(* J *)
Ti = TieV evtj;(* J *)

Ts[s_] := ToExpression["T" <> s];
Trs[s_] := ToExpression["Tr" <> s];
ms[s_] := ToExpression["m" <> s];
qs[s_] := ToExpression["q" <> s];
ns[s_] := ToExpression["n" <> s];
trm[s_] := ToExpression[s <> "trm"];

Ws[s_] := (qs[s] B)/ms[s]; (* s^-1 *)

ws[s_] := Sqrt[(ns[s] qs[s]^2)/(ep0 ms[s])]; (* s^-1 *)

vtsl[s_] := Sqrt[Ts[s]/ms[s]];(* m s^-1 *)

ys[s_] := Trs[s] - 1;

(* Functions *)
z[k_, w_, m_, s_] := (w - m Ws[s])/(k vtsl[s]);

W[z_] := 1 - Sqrt[\[Pi]/2] (Erfi[z/Sqrt[2]] - I) z E^(-z^2/2);

Z[z_] := -Sqrt[(\[Pi]/2)] (Erfi[z/Sqrt[2]] - I) E^(-z^2/2);

(* Dielectric Response Function Components *)

IT := ( {
{1, 0, 0},
{0, 1, 0},
{0, 0, 0}
} );

Ab[k_, w_, s_] := ( {
{ys[s], 0, 0},
{0, ys[s], 0},
{0, 0, (z[k, w, 0, s])^2 W[z[k, w, 0, s]]}
} );

Bb = 1/2 ( {
{1, I, 0},
{-I, 1, 0},
{0, 0, 0}
} );

eps[k_, w_, mMax_, s_] := -(ws[s]/w)^2 (Ab[k, w, s] + Bb (ys[s] z[k, w, +1, s] + z[k, w, 0, s]) Z[z[k, w, +1, s]] + Bb\[Conjugate] (ys[s] z[k, w, -1, s] + z[k, w, 0, s]) Z[z[k, w, -1, s]]);

ep[k_, w_] := IdentityMatrix[3] - \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$Ns$$]$$eps[k, w, trm[s[\([i]$$]], s[$$[i]$$]]\)\);
delta[k_, w_] := ep[k, w] - ((k c)/w)^2 IT;
f[k_, w_, y_] := Det[delta[+k, w + I y]] Det[delta[-k, w + I y]];

snorm = "i";
knorm = c/\[Omega]s["i"];

kint = (015/100)/knorm;
kfin = (300/100)/knorm;
kpow = 4;
kres = 2^kpow;
kinc = (kfin - kint)/(kres - 1);
k = Table[kint + (i - 1) kinc, {i, kres}];

lowfr = 1/10;

xxmin = lowfr;
xxmax = (2/2);
xxres = 16;
xxinc = (xxmax - xxmin)/(xxres - 1);
xx = Table[xxmin + (i - 1) xxinc, {i, xxres}];

yymin = -(1/10);
yymax = (30/100);
yyres = 16;
yyinc = (yymax - yymin)/(yyres - 1);
yy = Table[ymin + (i - 1) yyinc, {i, yyres}];

zGrid = Table[xx[[i]] + I yy[[j]], {j, yyres}, {i, xxres}];


Use l = 1 in the first block of code. The ordered pair j = 4 and i = 2 produce the other plot above. Please let me know if the code above does not run properly.

• You've included the code, and made your question clear, that is fantastic. But I'm missing a few definitions, without which I can't run it. zGrid, Ws, snorm, f, k, and norm` are all still blue after I run the code. – Jason B. Mar 24 '16 at 8:59
• Welcome! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Mar 24 '16 at 9:02
• JasonB, thank you so much for trying to get it running. Unfortunately, many of those quantities I can't give simple definitions here for you. It would be far more helpful if I could simply provide my whole .nb file. – S2167 Mar 24 '16 at 11:24
• The second block of code should provide the needed definitions. – S2167 Mar 24 '16 at 12:14
• maybe forget region functions and just directly implement the winding rule? – george2079 Mar 24 '16 at 12:37