Suppose that I have a list of coordinates {x,y} following some geometry in space, for example
TRIANGLE = Triangle[{{0, 0}, {1, 1}, {2, 0}}];
pts = RandomPoint[TRIANGLE, 3000];
And, I want to construct a matrix in which follow some rule
First we create the matrix with zeros (I would like to know if its really necessary to create it)
MATRIX = Table[0, {i, 1, Length[pts]}, {j, 1, Length[pts]}];
Or using Array
MATRIX = Array[0 &, {Length[pts], Length[pts]}];
And then, fill the matrix with the following condition (Is there a better way to construct this part?)
Table[
If[
EuclideanDistance[pts[[i]], pts[[j]]] ==
EuclideanDistance[pts[[1]], pts[[2]]](*or equal to other condition*),
MATRIX[[i, j]] = r
,
Nothing
]
, {i, 1, Length[pts]}, {j, 1, Length[pts]}
];
Edit1: (Timing check for 1000x1000 matrices)
The method above gives: 4.26563
And
MAT = SparseArray[{{i_, j_} /;
EuclideanDistance[pts[[i]], pts[[j]]] ==
EuclideanDistance[pts[[1]], pts[[2]]] -> r}, {Length[pts],
Length[pts]}];
Gives: 4.40625
I know this example is not so good because pts follows a random distribution, but I'm using it because I want to create big matrices like that.
I would like to know if there is another way to do it, in view of, for larger matrices it takes a lot of computational time.
Edit2: For a realistic geometrical system (The real problem!)
I want to construct a bilayered system in which consist of two layers interacting with each other. In the following code, I proceed with only one layer because the other will be analogous with its respectives lists.
So, we are going to check in both layers separated when the points are near by a distance of 1.
m = 1.;
n = m + 1;
t=1;
(*Basis vectors*)
a1 = {Sqrt[3]/2, -1/2}*Sqrt[3];
a2 = {Sqrt[3]/2, 1/2}*Sqrt[3];
(*Unit cell vectors*)
t1 = m*a1 + n*a2;
t2 = (n + m)*a1 - m*a2;
(*k vector (Eigenvalues will deppend upon it)*)
k = {kx, ky};
(*This function is responsible to distribute the hexagon points in x,y plane*)
pts[x_, y_] :=
Block[{j, k},
Flatten[Table[{{Sqrt[3] j, 1 k} Sqrt[3]}, {j, 0, x}, {k, 0, y}],
2]];
hexagon = {{0, 1/2} Sqrt[3], {Sqrt[3]/6, 1} Sqrt[
3], {Sqrt[3]/2, 1} Sqrt[3], {(2 Sqrt[3])/3, 1/2} Sqrt[
3], {Sqrt[3]/2, 0} Sqrt[3], {Sqrt[3]/6, 0} Sqrt[3]};
(*Angle of rotation between the layers*)
\[Theta] = ArcCos[1/2*(n^2 + 4 n*m + m^2)/(n^2 + n*m + m^2)];
rot = RotationTransform[N[\[Theta]], {0, 0}];
(*----Constructing the layers----*)
LAYER1 = TranslationTransform[# - hexagon[[4]]][hexagon] & /@
pts[(t1 + t2)[[1]]/2, -(t1 + t2).RotationMatrix[-60 Degree][[2]]/2];
LAYER2 = rot /@ TranslationTransform[# - hexagon[[3]]][hexagon] & /@
pts[(t1 + t2)[[1]]/2, -(t1 + t2).RotationMatrix[-60 Degree][[2]]/2];
(*--NOW I'LL PROCEED ONLY WITH THE ROTATED PART--*)
L2 = Flatten[LAYER2, 1];
(*---ROTATED---*)
ROT = DeleteDuplicates[Table[
If[
RegionMember[
ConvexHullMesh[({t1, t2, t1 + t2, {0, 0}} +
0.1).RotationMatrix[-60 Degree]], L2[[i]]]
, L2[[i]], Nothing]
, {i, 1, Length[L2]}]];
(*--OUTSIDE THE UNIT CELL TWISTED--*)
OUTREGIONTWISTED =
DeleteDuplicates[
Table[If[
RegionMember[
ConvexHullMesh[({{-2, 2}, t1 + {2, 2}, (t1 + t2) + {2, -2},
t2 + {-2, -2}, {-2, 2}} + 0.1).RotationMatrix[-60 Degree]],
L2[[u]]], L2[[u]], Nothing], {u, 1, Length[L2]}]];
(*--FUNCTION TO REMOVE THE DUPLICATED POINTS INSIDE OUTREGIONTWISTED--*)
SIZOR[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
o = Ordering[c];
x = 1 - 2 UnitStep[-1 - Length[b] + o];
x = FoldList[Max[#, 0] + #2 &, x];
x[[o]] = x;
Pick[c, x, -1]]
(*Points to construct the periodic system*)
PERIODICTWISTED = SIZOR[OUTREGIONTWISTED, ROT];
(*TWISTED*)
TWISTEDMATRIX= Table[0, {i, 1, Length[ROT]}, {j, 1, Length[ROT]}];
(*-----------ELEMENTS RELATED TO THE UNITCELL AND OUTSIDE-------------*)
Table[
If[
(*----CHECK THE FIRST DISTANCE BETWEEN THE POINTS INSIDE THE UNIT \
CELL AND OUTSITE WITH THAT WE GOT THE "j" POSITION IN MATRIX---------*)
EuclideanDistance[PERIODICTWISTED[[h]], ROT[[j]]] == 1 &&
(*----------ANOTHER PERIODIC CONTITION TO OBTAIN THE "i" POSITION \
IN MATRIX---------*)
PERIODICTWISTED[[h]] - (nn t1 + mm t2) == ROT[[i]],
(*--------THEN THE ELEMENT POSITION----------*)
TWISTEDMATRIX[[i, j]] = t*Exp[I*k.(ROT[[j]] - PERIODICTWISTED[[h]])];
,
Nothing]
, {nn, -1, 1}, {mm, -1, 1}, {i, 1, Length[ROT]}, {j, 1,
Length[ROT]}, {h, 1, Length[PERIODICTWISTED]}];
(*------------ELEMENTS RELATED WITH INSIDE THE UNIT CELL------------*)
(*HERE I CAN USE WHAT kglr SUGGESTED*)
(*MAT1=t*(1-Unitize[DistanceMatrix[ROT]-1])*)
(*AND THEN DO MAT1+TWISTEDMATRIX*)
Flatten[Table[If[EuclideanDistance[ROT[[i]], ROT[[j]]] == 1,
TWISTEDMATRIX[[i, j]] = t
, Nothing], {i, 1, Length[ROT]}, {j, 1, Length[ROT]}], 1];
Then, the idea is to combine all the matrices doing ArrayFlatten[{{TWISTEDMATRIX,SOME_INTERACTION_MATRIX},{SOME_INTERACTION_MATRIX,NOTTWISTEDMATRIX}}], and solve the Eigenvalues for kx and ky, for example from 0 to 2Pi to make a surfaceplot (ListPlot3D) or the other option is doing as follows
MAT[kx_,ky_]=Table[TWISTEDMATRIX[[i,j]],{i,1,Length[TWISTEDMATRIX]},{j,1,Length[TWISTEDMATRIX]}];
ClearAll[\[CapitalGamma], K, M, b1, b2, T1, T2, B1, B2];
T2 = Flatten[ArrayReshape[t2, {1, 3}]];
T1 = Flatten[ArrayReshape[t1, {1, 3}]];
B2 = 2 Pi (Cross[{0, 0, 1}, T1]/T2.Cross[{0, 0, 1}, T1]);
B1 = 2 Pi (Cross[T2, {0, 0, 1}]/T1.Cross[T2, {0, 0, 1}]);
(*--Reciprocal Vectors--*)
b2 = {B2[[1]], B2[[2]]};
b1 = {B1[[1]], B1[[2]]};
\[CapitalGamma] = {0, 0};
M = b1/2;
K = (b1 - b2)/3;
\[CapitalGamma]x = \[CapitalGamma][[1]];
\[CapitalGamma]y = \[CapitalGamma][[2]];
Kx = K[[1]];
Ky = K[[2]];
Mx = M[[1]];
My = M[[2]];
N1 = 500;
datK\[CapitalGamma] = Module[{k, kx, ky},
Transpose[Table[
k = -nn*
Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2]/N1;
kx = k*(Kx - \[CapitalGamma]x)/
Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2];
ky = k*(Ky - \[CapitalGamma]y)/
Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2]; \
({k + 0.92, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn,
1, N1}]]];
dat\[CapitalGamma]M = Module[{k, kx, ky}, Transpose[Table[
k = nn*
Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2]/N1;
kx = \[CapitalGamma]x +
k*(Mx - \[CapitalGamma]x)/
Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2];
ky = \[CapitalGamma]y +
k*(My - \[CapitalGamma]y)/
Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2]; \
({Sqrt[(\[CapitalGamma]x - Kx)^2 + (\[CapitalGamma]y - Ky)^2] +
k, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn,
1, N1}]]];
datMK = Module[{k, kx, ky},
Transpose[Table[
kx = Mx + k*(Kx - Mx)/Sqrt[(Kx - Mx)^2 + (Ky - My)^2];
ky = My + k*(Ky - My)/Sqrt[(Kx - Mx)^2 + (Ky - My)^2];
k = nn*Sqrt[(Kx - Mx)^2 + (Ky - My)^2]/
N1; ({Sqrt[(\[CapitalGamma]x - Kx)^2 + (\[CapitalGamma]y -
Ky)^2] +
Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2] \
+ k, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 1,
N1}]]];
ListLinePlot[Union[datK\[CapitalGamma], dat\[CapitalGamma]M, datMK]]
Edit3: Implementing the Henrik suggestion
The following code is what I got above together with what Henrik have developed, for both cases, for example for m=5 I get imaginary eigenvalues, as can be seen in Listplot corresponding to blank spaces. First of all, our matrix is supposed to be Hermitian, so, only real eigenvalues are expected. My guess is that some elements in upper triangular part are numerically different from the botom triangular part of MAT[kx,ky].
m = 5.; (*for example*)
n = m + 1;
a1 = {Sqrt[3]/2, -1/2}*Sqrt[3];
a2 = {Sqrt[3]/2, 1/2}*Sqrt[3];
r = n*a1 + m*a2;
t1 = m*a1 + n*a2;
t2 = (n + m)*a1 - m*a2;
pts[x_, y_] :=
Block[{j, k},
Flatten[Table[{{Sqrt[3] j, 1 k} Sqrt[3]}, {j, 0, x}, {k, 0, y}],
2]];
\[Theta] = ArcCos[1/2*(n^2 + 4 n*m + m^2)/(n^2 + n*m + m^2)];
hexagon = {{0, 1/2} Sqrt[3], {Sqrt[3]/6, 1} Sqrt[
3], {Sqrt[3]/2, 1} Sqrt[3], {(2 Sqrt[3])/3, 1/2} Sqrt[
3], {Sqrt[3]/2, 0} Sqrt[3], {Sqrt[3]/6, 0} Sqrt[3]};
rot = RotationTransform[N[\[Theta]], {0, 0}];
LAYER1 = TranslationTransform[# - hexagon[[4]]][hexagon] & /@
pts[(t1 + t2)[[1]]/2, -(t1 + t2).RotationMatrix[-60 Degree][[2]]/2];
LAYER2 = rot /@ TranslationTransform[# - hexagon[[3]]][hexagon] & /@
pts[(t1 + t2)[[1]]/2, -(t1 + t2).RotationMatrix[-60 Degree][[2]]/2];
L2 = Flatten[LAYER2, 1];
ROT = DeleteDuplicates[Table[
If[
RegionMember[
ConvexHullMesh[({t1, t2, t1 + t2, {0, 0}} +
0.1).RotationMatrix[-60 Degree]], L2[[i]]]
, L2[[i]], Nothing]
, {i, 1, Length[L2]}]];
OUTREGIONTWISTED =
DeleteDuplicates[
Table[If[
RegionMember[
ConvexHullMesh[({{-2, 2}, t1 + {2, 2}, (t1 + t2) + {2, -2},
t2 + {-2, -2}, {-2, 2}} + 0.1).RotationMatrix[-60 Degree]],
L2[[u]]], L2[[u]], Nothing], {u, 1, Length[L2]}]];
SIZOR[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
o = Ordering[c];
x = 1 - 2 UnitStep[-1 - Length[b] + o];
x = FoldList[Max[#, 0] + #2 &, x];
x[[o]] = x;
Pick[c, x, -1]]
PERIODICTWISTED = SIZOR[OUTREGIONTWISTED, ROT];
R = Developer`ToPackedArray[ROT];
P = Developer`ToPackedArray[PERIODICTWISTED];
\[Epsilon] = 0.00001;
A = ConstantArray[0., {Length[R], Length[R]}];
Do[If[Abs[EuclideanDistance[R[[i]], R[[j]]] - 1] <= \[Epsilon],
A[[i, j]] += t Exp[I k.(R[[j]] - R[[i]])];], {i, 1, Length[R]}, {j,
1, Length[R]}];
Do[If[Abs[EuclideanDistance[P[[h]], R[[j]]] - 1] <= \[Epsilon],
Do[If[Norm[P[[h]] - (R[[i]] + (nn t1 + mm t2))] <= \[Epsilon],
A[[i, j]] += t Exp[-I k.(P[[h]] - R[[j]])];];, {nn, -1,
1}, {mm, -1, 1}, {i, 1, Length[R]}]];, {j, 1, Length[R]}, {h, 1,
Length[P]}];
Or in alternative way
(*Make sure that `SparseArray` assemble is additive.*)
SetSystemOptions[
"SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}];
(*Create NearestFunction for R.*)
Rnf = Nearest[R -> "Index"];
(*Find all pairs {i,j} satisfying \
Abs[EuclideanDistance[R[[i]],R[[j]]]-1]\[LessEqual]\[Epsilon].*)
{ilist1, jlist1} =
Transpose[
Join @@ Map[Thread[{First[#], Rest[#]}] &,
Rnf[R, {\[Infinity], 1 + \[Epsilon]}]]];
(*Use the {i,j} pairs to assemble the matrix.*)
A1 = SparseArray[
Transpose[{ilist1, jlist1}] ->
t Exp[I (R[[jlist1]] - R[[ilist1]]).k], {Length[R], Length[R]}];
stencil = Flatten[Table[(nn t1 + mm t2), {nn, -1, 1}, {mm, -1, 1}], 1];
(*Displacing the points in R by each vector in stencil.*)
T = Flatten[Outer[Plus, R, stencil, 1], 1];
(*Create NearestFunction for T.*)
Tnf = Nearest[T -> "Index"];
(*Find all pairs {j,h} satisfying \
Abs[EuclideanDistance[P[[h]],R[[j]]]-1]\[LessEqual]\[Epsilon].*)
{jcandidates, hcandidates} =
Transpose[
Join @@ MapIndexed[Thread[{#1, First[#2]}] &,
Rnf[P, {\[Infinity], 1 + \[Epsilon]}]]];
(*For each pair {j,h} above,search for the corresponding i satisfying \
Norm[P[[h]]-(R[[i]]+(nn t1+mm t2))]\[LessEqual]\[Epsilon] by using \
the NearestFunction of T.*)
icandidates =
Quotient[Tnf[P[[hcandidates]], {\[Infinity], \[Epsilon]}] - 1,
Length[stencil]] + 1;
(*Create all {i,j,h} triplets.*)
{ilist2, jlist2, hlist2} =
Transpose[
Join @@ MapThread[
Thread[{#1, #2, #3}] &, {icandidates, jcandidates, hcandidates}]];
(*Use these triplets to assemble the matrix.*)
A2 = SparseArray[
Transpose[{ilist2, jlist2}] ->
t Exp[-I (P[[hlist2]] - R[[jlist2]]).k], {Length[R], Length[R]}];
A = A1 + A2;
Now, we proceed to the calculation of the Eigenvalues as showed above
MAT[kx_, ky_] = Table[A[[i, j]], {i, 1, Length[A]}, {j, 1, Length[A]}];
ClearAll[\[CapitalGamma], K, M, b1, b2, T1, T2, B1, B2];
T2 = Flatten[ArrayReshape[t2, {1, 3}]];
T1 = Flatten[ArrayReshape[t1, {1, 3}]];
B2 = 2 Pi (Cross[{0, 0, 1}, T1]/T2.Cross[{0, 0, 1}, T1]);
B1 = 2 Pi (Cross[T2, {0, 0, 1}]/T1.Cross[T2, {0, 0, 1}]);
(*--Reciprocal Vectors--*)
b2 = {B2[[1]], B2[[2]]};
b1 = {B1[[1]], B1[[2]]};
\[CapitalGamma] = {0, 0};
M = b1/2;
K = (b1 - b2)/3;
\[CapitalGamma]x = \[CapitalGamma][[1]];
\[CapitalGamma]y = \[CapitalGamma][[2]];
Kx = K[[1]];
Ky = K[[2]];
Mx = M[[1]];
My = M[[2]];
N1 = 500;
datK\[CapitalGamma] =
Module[{k, kx, ky},
Transpose[
ParallelTable[
k = -nn*Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \
\[CapitalGamma]y)^2]/N1;
kx =
k*(Kx - \[CapitalGamma]x)/
Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2];
ky =
k*(Ky - \[CapitalGamma]y)/
Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2]; \
({k + 0.92, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 1, N1}]]];
dat\[CapitalGamma]M =
Module[{k, kx, ky},
Transpose[
ParallelTable[
k = nn*Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \
\[CapitalGamma]y)^2]/N1;
kx = \[CapitalGamma]x +
k*(Mx - \[CapitalGamma]x)/
Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2];
ky = \[CapitalGamma]y +
k*(My - \[CapitalGamma]y)/
Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2]; \
({Sqrt[(\[CapitalGamma]x - Kx)^2 + (\[CapitalGamma]y - Ky)^2] +
k, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 1, N1}]]];
datMK = Module[{k, kx, ky},
Transpose[
ParallelTable[
kx = Mx + k*(Kx - Mx)/Sqrt[(Kx - Mx)^2 + (Ky - My)^2];
ky = My + k*(Ky - My)/Sqrt[(Kx - Mx)^2 + (Ky - My)^2];
k = nn*
Sqrt[(Kx - Mx)^2 + (Ky - My)^2]/
N1; ({Sqrt[(\[CapitalGamma]x - Kx)^2 + (\[CapitalGamma]y -
Ky)^2] +
Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2] \
+ k, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 1, N1}]]];
ListLinePlot[Union[datK\[CapitalGamma], dat\[CapitalGamma]M, datMK]]
SparseArray
$\endgroup$MAT1 =r (1 - Unitize[DistanceMatrix[pts] - EuclideanDistance @@ pts[[;; 2]]])
? $\endgroup$EuclideanDistance @@ pts[[;; 2]]]
with0.5
$\endgroup$