# Tracking solutions over a Brillouin zone

Thanks to @bbgodfrey I have a function that finds the two solutions to a set of equations bounded by a hexagonal region (corresponding to the Brillouin zone). It takes the angle {theta,phi} as the input and outputs the {{x1,y1},{x2,y2}} coordinate solutions. However there is a small issue as I demonstrate now:

If I keep theta=0.3 constant whilst varying phi over the entire range 0 to 2Pi at first glance it seems that everything is ok. The one solution that starts at approximately {2.4,0} travels in a clockwise closed continuous contour in the extended zone and returns to its starting point. In the Brillouin zone we see that this is equivalent to jumping from the middle corner A to B to C and back to A again. However note I have been plotting both solutions {x1,y1} and {x2,y2} fresh out of the function.

However if I choose either output [[1]] or [[2]] of the function to plot, things get messy. Whenever the function hits a Brillouin boundary the order in which it returns the solutions seems to reverse. i.e. Solutions {{x1,y1},{x2,y2}} are switched to {{x2,y2},{x1,y1}}. The figure below demonstrates this by only plotting the first entry of the output of the function.

So whenever the function hits a boundary wall it "switches" to the other solution (which is the negative of it). Now it's not the functions fault bless it, it doesn't have enough information to associate which solution is which. But then again I don't have quite enough intellect to give it this information which is why I'm asking for your help please. How can I force my function to return the solutions in the order that makes physical sense as demonstrated in the first figure? Or alternatively once I have a list of these solutions can I manipulate them to give the same effect?

Thanks very much for any and all help. Copied below is the function and the code to generate the plot:

DiracPoint[{theta_, phi_}] :=
Quiet[Module[{c1, c2, c3, x, y, cond1, cond2}, {c1, c2, c3} =
1 - 3 Sin[theta]^2 Cos[phi - 2 Pi (# - 1)/3]^2 & /@ {1, 2, 3};
Cases[If[0 <= ((c2 + c3)^2 - c1^2)/(4 c2 c3) <= 1,
cond1 = (c2 + c3) Cos[Sqrt[3] x/2] == -c1 Cos[3 y/2];
cond2 = (c2 - c3) Sin[Sqrt[3] x/2] == c1 Sin[3 y/2];
{x, y} /. Solve[cond1 && cond2, {x, y}], Null], {z1_, z2_} /;
Abs[z2] <= 2 Pi/3 && Abs[z1] <= (4 Pi/3 - Abs[z2])/Sqrt[3],
Infinity]]]

a = 2 Pi/(3 Sqrt[3]);
b = 2 Pi/3;
plot3 = ListPlot[{{2 a, 0}, {a, b}, {-a, b}, {-2 a,
0}, {-a, -b}, {a, -b}, {2 a, 0}}, PlotStyle -> Red,
Joined -> True];
angles = Table[phi, {phi, 0., 2 Pi, 0.01}];
array = DiracPoint[{0.3, #}] & /@ angles;
array = DeleteCases[array, {}];
array1 = #[[1]] & /@ array;
array2 = #[[2]] & /@ array;
plot1 = ListPlot[array1, Joined -> False];
plot2 = ListPlot[array2, Joined -> False];
Show[plot3, plot2, AspectRatio -> 1]