Problems with Solve

Question

I am having a problem similar to the previous question here. Solve either takes forever (I assume infinity) to run or returns an error. I have tried various combinations of using Round, Chop, N and Real to neaten up output, but I can't quite pinpoint where the problem is. In my previous question it was a mixture of exact and inexact coefficients causing Solve to run indefinitely, but here I am not sure. Normally I would post an analogous problem to make it easier for those reading but as I'm not sure even roughly where the issue is this time, I'm afraid I will post the full problem.

(*Nearest neighbour vectors*)
d1 = {0, -1} - {rx, ry};
d2 = {Sqrt[3]/2, 1/2} - {rx, ry};
d3 = {-Sqrt[3]/2, 1/2} - {rx, ry};

(*The coefficients*)
a = Norm[d1]^-3; b = Norm[d2]^-3; c = Norm[d3]^-3;

(*modulus of fq*)
modfq[k_] := 
  Sqrt[a^2 + b^2 + c^2 + 2 a b Cos[k.(d1 - d2)] + 
    2 b c Cos[k.(d2 - d3)] + 2 a c Cos[k.(d1 - d3)]];

(*Differential of fq w.r.t x and y to find turning points and then find minimum*)
findMin[p_] := (rx = p[[1]]; ry = p[[2]];
  Min[modfq /@ (Re[{x, y} /. 
       Solve[a b (d2 - d1)[[1]] Sin[{x, y}.(d1 - d2)] + 
          b c (d3 - d2)[[1]] Sin[{x, y}.(d2 - d3)] + 
          a c (d3 - d1)[[1]] Sin[{x, y}.(d1 - d3)] == 
          a b (d2 - d1)[[2]] Sin[{x, y}.(d1 - d2)] + 
          b c (d3 - d2)[[2]] Sin[{x, y}.(d2 - d3)] + 
          a c (d3 - d1)[[2]] Sin[{x, y}.(d1 - d3)] == 0, {x, y}]])])

(*Array of positions*)
res = 2; (*<--Ultimately this will be a value of a few thousand*)
xArray = Range[-Sqrt[3]/2, Sqrt[3]/2, Sqrt[3]/res];
yArray = Range[3/4, -3/4, -3/(2 res)];
g[p_] := Transpose@{xArray, ConstantArray[p, Length[xArray]]};
rArray = g /@ yArray ;

(*Map findMin onto rArray*)
diagram = Map[findMin, rArray, {2}];
ArrayPlot[diagram, ColorFunction -> "Rainbow"]

In this form it runs indefinitely, and attempts at rounding, choping etc. just don't seem to be working for me. Any suggestions on what the problem is?

Many thanks for any help, and apologies for the big code dump.

Answer 1

Clear["Global`*"];
f[rx_, ry_] := Module[{d1, d2, d3, a, b, c, sol1, sol2, data},
    {d1, d2, d3} = # - {rx, ry} & /@ {{0, -1}, {Sqrt[3]/2, 1/2}, {-Sqrt[3]/2, 1/2}};
    {a, b, c} = Norm[#]^-3 & /@ {d1, d2, d3};
    sol1 = NSolve[a b (d2 - d1) Sin[{x, y}.(d1 - d2)] + 
       b c (d3 - d2) Sin[{x, y}.(d2 - d3)] + 
       a c (d3 - d1) Sin[{x, y}.(d1 - d3)] == {0, 0}, {x, y}][[All, All, 2]];
    sol2 = Re@Flatten[Table[sol1 /. {C[1] -> i, C[2] -> j}, {i, -1, 1}, {j, -1, 1}], 2];
    data = Sqrt[a^2 + b^2 + c^2 + 2 a b Cos[#.(d1 - d2)] + 
       2 b c Cos[#.(d2 - d3)] + 2 a c Cos[#.(d1 - d3)]] & /@ sol2;
    Min@Chop[data, 10^-7]
]
res = 20;
diagram = Monitor[Table[f[rx, ry], 
    {rx, -Sqrt[3]/2, Sqrt[3]/2, Sqrt[3]/res}, {ry, 3/4, -3/4, -3/(2 res)}], 
    N[(rx + Sqrt[3]/2)/Sqrt[3]]];
ArrayPlot[diagram, ColorFunction -> "Rainbow", ColorFunctionScaling -> False]