Phonon Modes for 2D Lattice Vibrations

Clear[glConstMinSquaredDistanceFromOrigin,glConstMinAngleBetweenLatticeVectors,glConstLocDefault,glConstDefaultMass,glConstKDefaults,glConstKMin,glConstKMax,glConstMMax,glConstWindowHalfWidth,glConstSpringColors,glConstPrimaryDisplaySize,glConstTextDisplayOffsetHack,glConstNeighboringCellIndexes,glConstDynPlotRange,glConstSpringColorsByN];

glConstMinSquaredDistanceFromOrigin = 0.1 ;
glConstMinAngleBetweenLatticeVectors=Pi/6;
glConstLocDefault={{0.1,1.3},{1.9,0.8}};
glConstDefaultMass=20;
glConstKDefaults={0.25,0.25,0.25,0.25,0.25};
glConstKMin=0.05;
glConstKMax=1;
glConstMMax=30;
glConstMaxNumberOfMasses = 4 ;
glConstWindowHalfWidth=3;
glConstSpringColors=Darker@{Orange,Green,Purple,Cyan,Yellow};
glConstPrimaryDisplaySize=.9{380,380};
glConstTextDisplayOffsetHack={0.25,0.25};
glConstNeighboringCellIndexes={{1,0},{0,1},{1,1},{1,-1}};
glConstDynPlotRange = {{-glConstWindowHalfWidth/2, glConstWindowHalfWidth}, {-glConstWindowHalfWidth/2, glConstWindowHalfWidth}} ;

(*
{{-1,0}->RGBColor[...],{0,-1}->RGBColor[...], ...}
*)
glConstSpringColorsByN = DeleteDuplicates@Flatten[Table[s Append[glConstNeighboringCellIndexes,{0,0}][[n]] -> glConstSpringColors[[n]], {s, -1, 1, 2}, {n, 5} ],1] ;

(*
Future possible maintainance TODOs:
({{implement a min separation for the equilibrium mass positions too (as done for the lattice vectors). 
Trickier, since that should be based on the separations in the origin cell, 
and not the obsolute separations of the locators.
Even if that is done, that doesn't mean that the scaling won't allow for overlap in the dynamic tab.}, {would be cool to make qLoc a pair of Locators that:
 trace out constrained positions along the 2 Pi scaled reciprocal lattice vectors.
(returning the distance from the origin of each of those locators)}, {a message box on user error:
- when an angle or distance from origin change has driven a reset of the lattice vectors.
- when too many of the locators have been deleted and a forced reset to the defaults occurs.}, {Can the frequency qPlaneLoc LocatorPane look more like Slider2D in the dynamics tab for consistency?}, {Am using a calculated value for the max freq for the freq plot intersection plane. 
Hardcoded this to magic # that's out of the graphic (23), but that could be wrong for some params}, {Scaling of the masses ended up different for the dynPlot vs. parameters tab. 
Hacked the dynPlot scaling for now, to be a bit closer.}})
*)

ClearAll[checkAgainstRange] ;
checkAgainstRange::usage = "This is to deal with InputForm Manipulators, that allow non-numeric values to be passed, or values that exceed the range specified in the Manipulator." ;
checkAgainstRange[v_,default_,lowerLimit_, upperLimit_, typeFunc_ : NumberQ] := 
Module[{result},
result = If [ typeFunc[v],v, default] ;

result = If [ result <= upperLimit, result,default ] ;
result = If [ result >= lowerLimit, result,default ] ;

result
] ;

ClearAll[springPoints] ;
(* Based on my ListLinePlot answer posted in: http://mathematica.stackexchange.com/a/37228/10 *)
springPoints::usage = "springPoints[ {point1, point2}, numberOfTurns, height, fractionToDrawAsLinesAtEnds ].  Example:

ListLinePlot[springPoints[{{1,2},{3,5}}], AspectRatio->Automatic, PlotStyle -> Darker[ Purple ] ]" ;

springPoints[ a12_List, n_Integer:8,h_:.05, f_: 0.1 ] :=Module[{a1, a2, n1, springDiff, nd, r, r1 },
{a1, a2} = a12 ;
n1 = Norm[a1] ;
springDiff = a2 - a1 ;
nd = Norm[springDiff] ;
r = RotationMatrix[ArcTan @@  springDiff ] ;
r1 = r . {n1, 0} ;

{
Table[ a1 -r1 + r . { n1 + nd f + t (1 - 2f) nd, h Sin[ 2 Pi n t]}, {t, 0, 1, 0.01 } ],
Table[ a1 -r1 + r . { n1 + nd f + (1 - 2f) nd + t f nd, 0}, {t, 0, 1, 0.01 } ],
Table[ a1 -r1 + r . { n1 +t f nd, 0}, {t, 0, 1, 0.01 } ]
}
] ;


ClearAll[indexLabel] ;
indexLabel::usage = "(Subscript[k, i])(or other character) in italic. indexLabel['k', 1]" ;
indexLabel = Subscript[Style[#1,Italic], #2] & ;


ClearAll[kLable] ;
kLable::usage = "SubscriptBox[k,i\](or other character) in italic and colored by spring index. kLable[k]" ;
kLable = Style[ indexLabel["k", #], FontColor->glConstSpringColors[[#]] ] & ;


ClearAll[ massColors ] ;
massColors := ( Darker[{ Blue, Green, Purple, Red, Orange }][[Mod[#, 5] + 1] ] & );


ClearAll[ massLabel ] ;
massLabel := Style[indexLabel["m", #], massColors[#]] & ;

ClearAll[calcLineTable];
calcLineTable[ n_List, b_List, i_List ] :=Module[{first, second},
{first, second} = i ;
Table[
Line[
{-n[[first]]b[[first]] + j b[[second]],n[[first]]b[[first]] + j b[[second]]} 
] 
, {j, -n[[second]], n[[second]]}
]
] ;


ClearAll[calcReciprocalBasis];
calcReciprocalBasis::usage = "Return a reciprocal frame basis for a set of vectors.  This doesn't include the 2 π scaling that is common in solid state physics.  Example, displaying the desired Kronicker delta behaviour:

b = {{2,1},{-1/4, 2}} ;
g = calcReciprocalBasis[ b ] ;

{g[[1]].loc[[1]],
g[[1]].loc[[2]],
g[[2]].loc[[1]],
g[[2]].loc[[2]]}
" ;
calcReciprocalBasis[loc_List] := Inverse[ Transpose[ loc ] ] ;


ClearAll[nMassesFromLocArray];
nMassesFromLocArray[uLoc_] := ((Dimensions[uLoc] // First) -2) ;


ClearAll[adjustedmoSel];
adjustedmoSel[rNumM_,mOrigInd_,m1Ind_]:=If[rNumM>2,
If[ mOrigInd==m1Ind,
First[DeleteCases[Range[rNumM],m1Ind]],
mOrigInd 
],
First@Delete[Range[2],m1Ind]
];


ClearAll[constructKArrayElements];
constructKArrayElements[ i_Integer, j_Integer  ] := Module[{a},
a = Flatten[Table[{i, j, s glConstNeighboringCellIndexes[[n]]} -> glConstKDefaults[[n]], {s, -1, 1, 2}, {n, 4} ], 1] ;

If [ i < j, Append[a,{i,j,{0,0}} -> glConstKDefaults[[5]]], a ] 
] ;


ClearAll[constructKArray];
constructKArray[ r_Integer ] := Flatten[Table[constructKArrayElements[i,j], {i, r}, {j,r}], 2] ;


ClearAll[returnAlteredKArrayElements];
returnAlteredKArrayElements[kA_List,m1Ind_Integer,m2Ind_Integer,ni_Integer,v_]:=(kA/.{Rule[{m1Ind,m2Ind,glConstNeighboringCellIndexes[[ni]]},_]:>Rule[{m1Ind,m2Ind,glConstNeighboringCellIndexes[[ni]]},v],Rule[{m1Ind,m2Ind,-glConstNeighboringCellIndexes[[ni]]},_]:>Rule[{m1Ind,m2Ind,-glConstNeighboringCellIndexes[[ni]]},v]});


ClearAll[returnAlteredKarrayOriginElement];
returnAlteredKarrayOriginElement[kA_List,m1Ind_Integer,moInd_Integer,v_]:=Module[{m1oSet},
m1oSet=Append[Sort[{m1Ind,moInd}],{0,0}];

kA/.(Rule[m1oSet,_]:>Rule[m1oSet,v])
];


ClearAll[selectVFromKarray];
(*To refresh k1,k2,k3,k4 after changing m1Sel or m2Sel.*)
selectVFromKarray[kA_List,m1Ind_Integer,m2Ind_Integer,ni_Integer]:=({m1Ind,m2Ind,glConstNeighboringCellIndexes[[ni]]}/.kA);


ClearAll[selectOriginVFromKarray];
(*To refresh k5 after changing m1Sel and moSel*)
selectOriginVFromKarray[kA_List,m1Ind_Integer,moInd_Integer]:=(Append[Sort[{m1Ind,moInd}],{0,0}]/.kA);


ClearAll[returnKarrayForAdditionalMassIndex];
returnKarrayForAdditionalMassIndex[ kA_List,nmNew_Integer ] := Module[{k2},
k2 = Flatten[
(constructKArrayElements[#[[1]],#[[2]]] &)/@ (Select[Flatten[Table[ {i,j}, {i, nmNew}, {j, nmNew}], 1], Max[#] == nmNew &]), 2] ;

(*AppendTo[kArray, k2]*)
Flatten[{kA, k2}, 1] 
] ;


ClearAll[resetLatticeVectorsIfAngleTooSmall];
resetLatticeVectorsIfAngleTooSmall[ uLoc_List,minAngle_] := Module[{t},
t = Abs[ArcCos[Normalize[uLoc[[1]]] . Normalize[uLoc[[2]]]]] ;
t = If [ t > Pi/2, Pi-t, t] ;

If[ t < minAngle,
glConstLocDefault, {uLoc[[1]], uLoc[[2]]}]
] ;


ClearAll[projOpU];
projOpU[v_List] := {{v[[1]]^2, v[[1]]v[[2]]},{v[[1]]v[[2]], v[[2]]^2}} ;


projOp::usage = "given an input vector (OverVector[v]) = {(Subscript[v, x]), (Subscript[v, y])}, compute the projection matrix operator along the unit vector in that direction.

   projOp[{1, 0}] // MatrixForm = \!\(\*TagBox[
RowBox[{"(", "", GridBox[{
{"1", "0"},
{"0", "0"}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, "RowsIndexed" -> {}},
GridBoxSpacings->{"Columns" -> {
Offset[0.27999999999999997`], {
Offset[0.7]}, 
Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {
Offset[0.2], {
Offset[0.4]}, 
Offset[0.2]}, "RowsIndexed" -> {}}], "", ")"}],
Function[BoxForm`e$, 
MatrixForm[BoxForm`e$]]]\)
projOp[{0, 1}] // MatrixForm = \!\(\*TagBox[
RowBox[{"(", "", GridBox[{
{"0", "0"},
{"0", "1"}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, "RowsIndexed" -> {}},
GridBoxSpacings->{"Columns" -> {
Offset[0.27999999999999997`], {
Offset[0.7]}, 
Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {
Offset[0.2], {
Offset[0.4]}, 
Offset[0.2]}, "RowsIndexed" -> {}}], "", ")"}],
Function[BoxForm`e$, 
MatrixForm[BoxForm`e$]]]\)
projOp[{a,b}] // MatrixForm = \!\(\*Cell[TextData[Cell[BoxData[
RowBox[{"(", "", GridBox[{
{
FractionBox[
SuperscriptBox["a", "2"], 
RowBox[{
SuperscriptBox["a", "2"], "+", 
SuperscriptBox["b", "2"]}]], 
FractionBox[
RowBox[{"a", " ", "b"}], 
RowBox[{
SuperscriptBox["a", "2"], "+", 
SuperscriptBox["b", "2"]}]]},
{
FractionBox[
RowBox[{"a", " ", "b"}], 
RowBox[{
SuperscriptBox["a", "2"], "+", 
SuperscriptBox["b", "2"]}]], 
FractionBox[
SuperscriptBox["b", "2"], 
RowBox[{
SuperscriptBox["a", "2"], "+", 
SuperscriptBox["b", "2"]}]]}
},
GridBoxAlignment->{"Columns" -> {{Center}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, "RowsIndexed" -> {}},
GridBoxSpacings->{"Columns" -> {
Offset[0.27999999999999997`], {
Offset[0.7]}, 
Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {
Offset[0.2], {
Offset[0.4]}, 
Offset[0.2]}, "RowsIndexed" -> {}}], "", ")"}]],
CellChangeTimes->{{3.5990911441851797`*^9, 3.599091144641206*^9}, 3.5990913830438423`*^9, 3.5990915191646276`*^9, 3.599091555028679*^9, 3.599091735722014*^9, {3.599091916353346*^9, 3.599091928988068*^9}}]]]\)
" ;
ClearAll[projOp];
projOp[v_List] := projOpU[v]/(v. v) ;


ClearAll[ relativeDifferences] ;
relativeDifferences::usage = "Example: relativeDifferences[ {aLoc, bLoc}, {mLoc1, mLoc2, ...}, kArray[[1]] /. Rule -> List ].

Here, the list of the form: (kArray[[1]] /. Rule -> List) // First = {1,1,{-1,0}} = {m1index, m2index, {nIndex1, nIndex2}},

 is a tuple for the spring constant between m1 in the origin cell to m2 in the lattice cell {n1,n2} ;
" ;
relativeDifferences[ r_List,mp_List, {mIndex1_Integer, mIndex2_Integer, n_List}] := Module[{d, pOrigin, pOther, dn},
pOrigin = mp[[ mIndex1 ]] ;
pOther = mp[[ mIndex2 ]]+ n . r  ;
d = pOther - pOrigin ;
dn = d .d ;

{
dn, pOrigin, pOther,(projOpU[d]/dn  (*// MatrixForm*))
}
] ;


ClearAll[calculateCouplings];
calculateCouplings::usage = "Returns a pair of sets for each origin m_i coupling:

With:
(Subscript[(OverVector[p]), k]) for the position of the (k^th) mass in the origin cell, 
j for the index of the neighbouring cell mass, 
(OverVector[m]) for the cell index of the neighbouring cell,

this function returns elements of the form:
{k, j, (OverVector[m]), (Subscript[K, (OverVector[n]), k; (OverVector[m]), j]), ((Subscript[(OverVector[a]), (OverVector[m]), j; (OverVector[0]), k])^2), (Subscript[(OverVector[p]), k]), (Subscript[(OverVector[r]), (OverVector[m])]) + (Subscript[(OverVector[p]), j]), (Subscript[Proj, (Subscript[(OverVector[a]), (OverVector[m]), j; (OverVector[0]), k])])}

There are two sets of lists returned:
First has all the origin cell couplings for that mass.
Second has all the nn couplings for each mass sorted by distance from mass.  Can use that to select up to num-interactions." ;
calculateCouplings[ld_List, kA_List] := Module[{t, t1, t2, r, mp, numberOfMasses},

{r, mp, numberOfMasses} = {"latticeBasis", "mPosFirstCell", "numberOfMasses"}  /. ld ;

t = Append[#, relativeDifferences[r, mp, First[#]] ]&/@ (kA/. Rule -> List) ;
t = Flatten[{#[[1]], {#[[2]]}, #[[3]]}, 1] &/@ t ;

t1 = Table[
Sort[(Select[ t, (#[[1]] == i) && (#[[3]] ≠ {0, 0}) &]), #1[[5]] < #2[[5]] &]
, {i, numberOfMasses}
] ;

t2 = Select[ t,#[[3]] == {0, 0} &] ;
(* the rest of the permutations: *)
t2 = Flatten[{t2, Flatten[{{#[[2]], #[[1]]}, Drop[#, 2]}, 1] &/@ t2}, 1]  ;
t2 = Table[
Sort[(Select[ t2, (#[[1]] == i) &]), #1[[5]] < #2[[5]] &]
, {i, numberOfMasses}
] ;

{t2, t1}
] ;

ClearAll[spScaled];
spScaled[{s1_List,s2_List}, scaleKValue_, referenceSepSquared_] := Module[{ds},
ds = s2 - s1 ;
springPoints[ {s1,s2},Ceiling[12 scaleKValue Sqrt[(ds.ds)/referenceSepSquared]] ]
] ;


ClearAll[oneTermOfdynamicsMatrix];
oneTermOfdynamicsMatrix::usage = "\!\(TraditionalForm\`\(grab\\ data\\ \(from : \\ rawArrayOfMasses\), \\ parameterLocatorDependencies, 
\(\(couplingData[\([_]\)]\)[\([_]\)]\)[\([_]\)]\\  = \\ {k, \\ j, \\ \*OverscriptBox[\(m\), \(⇀\)], \\ \*SubscriptBox[\(K\), \(\*OverscriptBox[\(n\), \(⇀\)], k; \*OverscriptBox[\(m\), \(⇀\)], j\)], \\ \*SuperscriptBox[SubscriptBox[OverscriptBox[\(a\), \(⇀\)], \(\*OverscriptBox[\(m\), \(⇀\)], j; \*OverscriptBox[\(0\), \(⇀\)], k\)], \(2\)], \\ \*SubscriptBox[OverscriptBox[\(p\), \(⇀\)], \(k\)], \\ \*SubscriptBox[OverscriptBox[\(r\), \(⇀\)], OverscriptBox[\(m\), \(⇀\)]]\\  + \\ \*SubscriptBox[OverscriptBox[\(p\), \(⇀\)], \(j\)], \\ \*SubscriptBox[\(Proj\), SubscriptBox[OverscriptBox[\(a\), \(⇀\)], \(\*OverscriptBox[\(m\), \(⇀\)], j; \*OverscriptBox[\(0\), \(⇀\)], k\)]]}

ex : \\ {1, 2, {\(-1\), 0}, 0.5`, 1.`, {0.9`, 0.7`}, {0.30000000000000004`, \(-0.10000000000000009`\)}, {{0.36`, 0.48`}, {0.48`, 0.6400000000000001`}}}


to\\ form\\ equations\\ of\\ the\\ \(form : 
\*SuperscriptBox[\(ω\), \(2\)]\\ \(\*SubscriptBox[OverscriptBox[\(ϵ\), \(⇀\)], \(k\)](\*OverscriptBox[\(q\), \(⇀\)])\)\)\\  = \\ \*UnderscriptBox[\(∑\), \(\*OverscriptBox[\(0\), \(⇀\)], k\\  ≠ \\ \*OverscriptBox[\(m\), \(⇀\)], j\)]\*SubscriptBox[\(K\), \(\*OverscriptBox[\(0\), \(⇀\)], k; \*OverscriptBox[\(m\), \(⇀\)], j\)] \*SubscriptBox[\(Proj\), SubscriptBox[OverscriptBox[\(a\), \(^\)], \(\*OverscriptBox[\(n\), \(⇀\)], k; \*OverscriptBox[\(m\), \(⇀\)], j\)]]\\ \((\*FractionBox[\(\*SubscriptBox[OverscriptBox[\(ϵ\), \(⇀\)], \(k\)] \((\*OverscriptBox[\(q\), \(⇀\)])\)\), SubscriptBox[\(m\), \(k\)]] - \\ \*FractionBox[\(\*SubscriptBox[OverscriptBox[\(ϵ\), \(⇀\)], \(j\)] \((\*OverscriptBox[\(q\), \(⇀\)])\)\), SqrtBox[\(\*SubscriptBox[\(m\), \(k\)] \*SubscriptBox[\(m\), \(j\)]\)]] \*SuperscriptBox[\(e\), \(I\\ \((\*SubscriptBox[OverscriptBox[\(r\), \(⇀\)], OverscriptBox[\(m\), \(⇀\)]]\\  . \\ \*OverscriptBox[\(q\), \(⇀\)]\\ )\)\)])\)

Note\\ that\\ qVec\\ here\\ is\\ intended\\ to\\ be\\ an\\ undefined\\ global . \\ \\ Is\\ there\\ a\\ way\\ to\\ do\\ that\\ \(\(programatically\)\(?\)\)\)\)" ;
oneTermOfdynamicsMatrix[{k_Integer, j_Integer, mVec_List, kS_, _,_,_,proj_List}, r_List, mk_, epsilon_List, mA_List] := (kS proj  .( epsilon[[k]]/mk - epsilon[[j]]E^(I((mVec . r). qVec))/Sqrt[ mk mA[[j]] ])  )  ;


ClearAll[dynamicsMatrix];
dynamicsMatrix[ld_List, cd_List] := Module[{cdFlatAndSorted, r, res, epsilon, m,mScalarArray, numberOfMasses},

cdFlatAndSorted = SplitBy[ Sort@Flatten[cd, 2] , First ] ;
{r,mScalarArray, numberOfMasses} = {"latticeBasis", "mScalarArray", "numberOfMasses"}  /. ld ; 

(*Here (Subscript[(OverVector[ϵ]), k]) (=epsilon[[k]]) is represented as a 2x(2N) matrix with the vector factored out.  For example with N = 3, 
(Subscript[(OverVector[ϵ]), 1])= (({{1, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}}))

This will be pre-multiplied by a 2x2 projection matrix.
The whole thing should be logically treated as if it is post-multiplied by

(({{(Subscript[ϵ, 1,1])}, {(Subscript[ϵ, 1,2])}, {(Subscript[ϵ, 2,1])}, {(Subscript[ϵ, 2,2])}, {(Subscript[ϵ, 3,1])}, {(Subscript[ϵ, 3,2])}}))
*)
epsilon = Partition[IdentityMatrix[ 2 numberOfMasses ],2] ;

res = Partition[Table[ Total@(oneTermOfdynamicsMatrix[#, r, mScalarArray[[i]], epsilon, mScalarArray ] &/@(cdFlatAndSorted[[i]])), {i, numberOfMasses}] // Flatten, 2 numberOfMasses] ;
res /. qVec -> # & 


(*;cdFlatAndSorted[[1]]*)
(*;mScalarArray[[1]]*)
(*;oneTermOfdynamicsMatrix[#, r, mScalarArray[[1]], epsilon, mScalarArray ] &/@(cdFlatAndSorted[[1]])*)
(*;res*)
(*;epsilon*)
(*;oneTermOfdynamicsMatrix[ parameterLocatorDependencies, couplingData ][#, r,mScalarArray[[1]], epsilon,mScalarArray ] &/@(cd[[1]])*)
] ;



ClearAll[calcDynamics];
calcDynamics::usage = "a list of 2 N elements {(Subscript[ω, k]),\!(OverVector[ϵ])}, one for each of the system eigenvalues, when evaluated at qLoc." ;
(*Example:

Module[{ld, cd, basisAndMassPoints, kArray, matrix, dynamics, qPoint, omegaIndex},
basisAndMassPoints = { {0.6,1}, {1,-0.4}, {0.3,3} } ;
ld = locDependent[ basisAndMassPoints, 1,{10} ]  ;
kArray = constructKArray[ 1 ] ;
cd = calculateCouplings[ ld, kArray ] ;
qPoint = ("qMaxVec" /. ld)/2 ;
matrix = dynamicsMatrix[ ld, cd ] ;
dynamics = calcDynamics[ matrix[ qLoc ] ] // First // TableForm
]*)
calcDynamics[mAtQ_List] :=({Sqrt[#[[1]]], #[[2]]}&/@(Eigensystem[ mAtQ ] // Transpose))  ;


ClearAll[frequencyPlotCalc];
frequencyPlotCalc::usage = "frequencyPlotCalc[ parameterLocatorDependencies, matrix, freqMeshSize]" ;
frequencyPlotCalc[ld_List,m_,meshSz_Integer] := Module[{eigTable2D, eigTable1D, qMax, omegaPointList2D, omegaPointList1D, range, numberOfMasses},

{qMax, numberOfMasses} = {"qMaxVec","numberOfMasses"} /. ld ;

eigTable2D = Flatten[Re[Table[  {{qx, qy} ,Eigenvalues[ m[ {qx, qy}  ] // N ]}, {qx, -qMax[[1]]/2, qMax[[1]]/2, qMax[[1]]/ meshSz}, {qy, -qMax[[2]]/2, qMax[[2]]/2, qMax[[2]]/ meshSz} ] ],1] ;

omegaPointList2D[nn_] :=  Flatten[{#[[1]],Sqrt[#[[2]]][[nn]]}]&/@ eigTable2D  ;

range = ((2 numberOfMasses) // Range) ;

ListPlot3D[ 
omegaPointList2D[#] &/@ range, 
PlotRange -> Full , 
ImageSize->glConstPrimaryDisplaySize,
AxesLabel->{"(Subscript[q, x])", "(Subscript[q, y])"},
PerformanceGoal->"Quality"
]
] ;

 

ClearAll[showFreqPlot];
showFreqPlot::usage = "showFreqPlot[ qPlaneLoc, freqPlot ]" ;
showFreqPlot[qpl_List, freqPlotData_] := Module[{aa, bb},

{aa, bb} =  qpl ;

Show[{freqPlotData,
Graphics3D[ Polygon[{{aa[[1]], aa[[2]],0}, {bb[[1]], bb[[2]],0},  {bb[[1]], bb[[2]],23},{aa[[1]], aa[[2]],23}}]]
}]
] ;



ClearAll[pointsTable];
pointsTable::usage = "Adds a position within a cell to a set of vectors to the origins of all the cells of interest.  Example in findLatticeCells::usage." ;
pointsTable[ massPosition_List, latticeOrigins_List ] := ((# + massPosition) &/@ latticeOrigins) ;



ClearAll[ findLatticeCells ] ;
findLatticeCells::usage = "All the cell index pairs that lie within a rectangular region.

Example: 
Module[{a, b, nAndPoints, cellOrigins, sol},
a = {1/2,2} ;
b = {2,-1/3} ;

nAndPoints = findLatticeCells[{a,b}, {{-3, 5},{-1,6}}]  ;
cellOrigins = (nAndPoints // Transpose // Last) ;

Show[{
ListPlot[cellOrigins, PlotStyle->{Darker@Red, PointSize[Large]},
Epilog-> {Darker@Blue,Arrow[{{0,0},a}], Darker@Green, Arrow[{{0,0},b}]} ]
,
ListPlot[pointsTable[{0.3, 0.5}, cellOrigins], PlotStyle->{Darker@Purple, PointSize[Medium]} ]
} ]
]
" ;
findLatticeCells[latticeVectors_List,limits_List] := Module[{a, b,sol, oneSolution},
{a, b} = latticeVectors ;

(* Rationalize to fix: to: http://mathematica.stackexchange.com/questions/6055/how-to-get-rid-of-warnings-when-using-solve-on-an-equation-with-inexact-coeffici 
*)
oneSolution[o_] := {n,m} /.Solve[
Rationalize[
(limits[[1]][[1]] <  (n a + m b + o) . {1,0}) &&
( (n a + m b + o) . {1,0} < limits[[1]][[2]]) &&

(limits[[1]][[1]] <  (n a + m b +o) . {0,1}) &&
( (n a + m b +o) . {0,1} < limits[[1]][[2]]) ]
, {n,m}, Integers] ;

sol = Flatten[oneSolution[#] &/@ {{0,0},a, b, a+b} ,1] // Sort // DeleteDuplicates ;

{#,a #[[1]] + b#[[2]]} &/@ sol 
] ;


ClearAll[locDependent];
locDependent::usage = "Locator dependent calculations (i.e. based on the mass positions and the unit cell basis vectors)

Example:

ld = locDependent[ { {1/2,1}, {1,1/2}, {1/2,3} }, 1, {10} ] ;
" ;
locDependent[ uLoc_List,numMasses_Integer, rma_List ]  := Module[{m,latticeBasis, numberLatticeLinesToDisplay,reciprocalBasis,mObliqueComponents, mPosFirstCell, reciprocalNorms, allLatticeCellOriginsInDynPlotRange, cellOrigins},
latticeBasis =Take[uLoc, 2] ;
m = Drop[uLoc, 2] ;
numberLatticeLinesToDisplay = (Ceiling[  Abs[glConstWindowHalfWidth/ latticeBasis[[#]][[#]]]] & /@ Range[2]) ;
reciprocalBasis = calcReciprocalBasis[ latticeBasis ] ;
reciprocalNorms = (Norm[reciprocalBasis[[#]]] & /@ Range[2]) ;

mObliqueComponents = Table[m[[ i ]] . reciprocalBasis[[ j ]], {i, numMasses}, {j, 2}] ;
mPosFirstCell = (m[[#]] - Floor[mObliqueComponents[[#]]] . latticeBasis ) & /@ Range[numMasses] ;
allLatticeCellOriginsInDynPlotRange = findLatticeCells[latticeBasis,glConstDynPlotRange] // Transpose ;
cellOrigins = (allLatticeCellOriginsInDynPlotRange // Last) ;

{
"mScalarArray" -> rma,
"numberOfMasses" -> numMasses,
"latticeBasis" -> latticeBasis,
"latticeNorms"->(Norm[latticeBasis[[#]]] & /@ Range[2]),
"latticeUnitVectors"->(Normalize[latticeBasis[[#]]] & /@ Range[2]),
"numberLatticeLinesToDisplay"-> numberLatticeLinesToDisplay,
"latticeCellOriginIndexes" -> (allLatticeCellOriginsInDynPlotRange // First),
"reciprocalBasis"-> reciprocalBasis,
"reciprocalNorms"-> reciprocalNorms,
"qMaxVec" -> (2 Pi reciprocalNorms),
"mObliqueComponents"-> mObliqueComponents,
"mPosFirstCell"-> mPosFirstCell,
"pointsDataTable"-> ((pointsTable[mPosFirstCell[[#]], cellOrigins]) &/@ Range[numMasses]),
"lineTable" -> (calcLineTable[ numberLatticeLinesToDisplay, latticeBasis, # ] & /@ Permutations[{1,2}])
}
] ;


ClearAll[plotSprings];
plotSprings::usage = "Example:

Module[{ld, cd, basisAndMassPoints, kArray},
basisAndMassPoints = { {0.6,1}, {1,-0.4}, {0.3,3} } ;
ld = locDependent[ basisAndMassPoints, 1,{10} ]  ;
kArray = constructKArray[ 1 ] ;
cd = calculateCouplings[ ld, kArray ] ;
plotSprings[basisAndMassPoints, ld, cd, 1, 1, 1 ] 
]
" ;
plotSprings[uLoc_List, ld_List,cd_List,m1Ind_Integer, m2Ind_Integer, moInd_Integer] := Module[{aLoc, bLoc,mLoc,pointsList,latticeBasis,reciprocalBasis,pointsDataTable, numberLatticeLinesToDisplay, lines, cOrigin, cN, g, mScalarArray, numberOfMasses, scaleLoc},

latticeBasis = Take[ uLoc, 2 ] ;
aLoc = uLoc[[1]] ;
bLoc = uLoc[[2]] ;
scaleLoc = aLoc.aLoc  + bLoc .bLoc ;

{mLoc,reciprocalBasis,pointsDataTable,numberLatticeLinesToDisplay, lines, mScalarArray, numberOfMasses}={"mPosFirstCell","reciprocalBasis","pointsDataTable","numberLatticeLinesToDisplay", "lineTable", "mScalarArray", "numberOfMasses"}  /. ld ;

pointsList[n_Integer] :={massColors[n],
,PointSize[Sqrt[mScalarArray[[n]]/glConstMMax/500]]
,Point[ # ] & /@ pointsDataTable[[n]]
,Text[massLabel[ n],mLoc[[n]] + glConstTextDisplayOffsetHack]
} ;

cOrigin = Select[ cd[[1]][[m1Ind]], Take[#, 2] == {m1Ind, moInd} & ] ;
cN = Select[ cd[[2]][[m1Ind]], Take[#, 2] =={m1Ind, m2Ind}& ] ;

g = ListLinePlot[
spScaled[Take[#,{6,7}], #[[4]]/glConstKMax, scaleLoc],
AspectRatio->Automatic, 
PlotStyle -> (#[[3]] /. glConstSpringColorsByN) ] &/@ cN ;

If [ numberOfMasses > 1,AppendTo[g,
ListLinePlot[
spScaled[Take[#,{6,7}], #[[4]]/glConstKMax, scaleLoc],
AspectRatio->Automatic, 
PlotStyle -> (#[[3]] /. glConstSpringColorsByN) ] &/@ cOrigin], ] ;

Show[{
Graphics[
Flatten[{
{
 lines
,Blue,Arrow[{{0,0}, reciprocalBasis[[#]]}] & /@ Range[2]
,Thick,Arrowheads[0.05]
,Red,Arrow[{{0,0}, latticeBasis[[#]]}] & /@ Range[2]
,Text[OverVector[Style["a",Italic]],aLoc/2 + glConstTextDisplayOffsetHack]
,Text[OverVector[Style["b",Italic]],bLoc/2 + glConstTextDisplayOffsetHack]
(* This polygon (with doubling of PlotRange -> glConstDynPlotRange below) was to visually test the Solve logic used to find all the lattice cells contained or just outside of the region. *)
(*,Lighter@Lighter@Gray
,Opacity[0.5]
,Polygon[
({{glConstDynPlotRange[[1]][[1]],glConstDynPlotRange[[2]][[1]]},{glConstDynPlotRange[[1]][[2]],glConstDynPlotRange[[2]][[1]]},{glConstDynPlotRange[[1]][[2]],glConstDynPlotRange[[2]][[2]]},{glConstDynPlotRange[[1]][[1]],glConstDynPlotRange[[2]][[2]]}} )
 ] *)
},
pointsList[#] &/@ Range[numberOfMasses]
}]

,PlotRange -> (*2*) glConstDynPlotRange
,ImageSize->glConstPrimaryDisplaySize
],
g
} ] 
] ;

ClearAll[dynamicsGraphicsTable];
dynamicsGraphicsTable::usage = "Calculate a time indexed table of positions for each of the masses" ;
(*
Example:

Module[{ld, cd, basisAndMassPoints, kArray, matrix, dynamics, sc, qPoint, omegaIndex,  lines,mScalarArray, numberOfMasses, pointsData,dynAtOmegaIndex, dynMesh,dynTable},
basisAndMassPoints = { {0.6,1}, {1,-0.4}, {0.3,3} } ;
ld = locDependent[ basisAndMassPoints, 1,{10} ]  ;
kArray = constructKArray[ 1 ] ;
cd = calculateCouplings[ ld, kArray ] ;
sc = 0.25 ;
qPoint = ("qMaxVec" /. ld)/2 ;
matrix = dynamicsMatrix[ ld, cd ] ;
dynamics = calcDynamics[ matrix[ qPoint ] ] ;
omegaIndex = 1 ;
dynAtOmegaIndex = dynamics[[omegaIndex]] ;
dynMesh = 30 ;
dynTable = dynamicsGraphicsTable[ld, dynAtOmegaIndex, sc, qPoint, dynMesh ] ;

(*With[ {n=1},
showDynPlot[ dynTable, n ] ]*)
Animate[ 
showDynPlot[ dynTable, n ], {n, 1, dynMesh, 1}]
]*)
dynamicsGraphicsTable[ld_List, dynAtOmegaIndex_List, sc_, qPoint_List, dynMesh_Integer ] :=Module[{pointsDataTable,e, omega, lines, nu, a,b, mScalarArray, numberOfMasses, latticeCellOriginIndexes, mPosFirstCell,pointsData},

{{a, b},pointsDataTable,lines, mScalarArray, numberOfMasses, latticeCellOriginIndexes, mPosFirstCell}={"latticeBasis","pointsDataTable", "lineTable", "mScalarArray", "numberOfMasses", "latticeCellOriginIndexes", "mPosFirstCell"}  /. ld ;

{omega, e} = dynAtOmegaIndex ;
e = Partition[e,2] ;

nu = 2 Pi If[ omega == 0, 1, 1/omega] ;

pointsData := Table[
a #[[1]] + b #[[2]] 
+ mPosFirstCell[[k]] 
+ sc Re[ e[[ k ]] E^(I( qPoint . (a #[[1]] + b #[[2]] )  - nu tau )) ] 
&/@ latticeCellOriginIndexes,
{k, numberOfMasses},{tau,0,1,1/dynMesh}
] ;

{ pointsData, mScalarArray,numberOfMasses,lines,dynMesh }
] ;

ClearAll[showDynPlot];
showDynPlot::usage = "showDynPlot[ dynamicsGraphicsTable[ld, dynAtOmegaIndex, sc, qPoint, dynMesh ] ]" ;
showDynPlot[{ pointsData_List, mScalarArray_List,numberOfMasses_Integer,lines_List, dynMesh_Integer }, n_ ] := Show[Flatten[
{
Table[
ListPlot[ pointsData[[k]][[n]]
,PlotRange -> glConstDynPlotRange
,ImageSize->glConstPrimaryDisplaySize
, PlotStyle->Directive[PointSize[Sqrt[mScalarArray[[k]]/glConstMMax/350]],massColors[k]]
], {k,numberOfMasses}]
,
{
Graphics[{
lines
}]
}
}, 1]
] ;

Manipulate[
tick;
(* This If[] is effectively my Initialization :> {} block.  Done with a variable so that copies of the output cell won't re-execute their Initialization cells.  This allows "interesting" settings to be maintained per-copy of the output cell, but still provide an Initialization block for the primary output cell. *)
 If [ bInitDefaults,(
u = Flatten[ {glConstLocDefault, mLocDefault}, 1] ;

parameterLocatorDependencies = locDependent[ u, rawNumberOfMasses,rawArrayOfMasses ] ;
kArray = constructKArray[ 1 ] ;
Module[{n},n = 2 ;While[n≤rawNumberOfMasses,kArray = returnKarrayForAdditionalMassIndex[kArray,n] ; n++] ] ;

qLoc = ("qMaxVec" /. parameterLocatorDependencies)/2 ;
qPlaneLoc = {-qLoc,qLoc} ;
couplingData = calculateCouplings[ parameterLocatorDependencies, kArray ] ;

(* bRefreshMatrix: dependent on: parameterLocatorDependencies, couplingData *)
matrix = dynamicsMatrix[ parameterLocatorDependencies, couplingData ] ;

(* bRefreshDynamics: dependent on matrix, qLoc, *)
dynamics = calcDynamics[ matrix[ qLoc ] ] ;

(* bRefreshDynPlot: dependent on dynamics, omegaIndex, scale *)
dynPlot = dynamicsGraphicsTable[ parameterLocatorDependencies, dynamics[[omegaIndex]], scale, qLoc, dynMeshSize ] ;

(* bRefreshFreqPlot: dependent on matrix, freqMeshSize *)
freqPlot = frequencyPlotCalc[ parameterLocatorDependencies, matrix, freqMeshSize] ;

bInitDefaults = False ;
),] ;



(* This is the main panel display *)
Dynamic@Switch[ tabNumber,
dynTab,
(
If [ bRefreshMatrix,
(matrix = dynamicsMatrix[ parameterLocatorDependencies, couplingData ] ;bRefreshMatrix = False ;bRefreshDynamics = True ; bRefreshFreqPlot = True ),];

If[ bRefreshDynamics,
(dynamics = calcDynamics[ matrix[ qLoc ] ] ;bRefreshDynamics = False ;bRefreshDynPlot = True),];

If[ bRefreshDynPlot,
(dynPlot = dynamicsGraphicsTable[ parameterLocatorDependencies, dynamics[[omegaIndex]], scale, qLoc,dynMeshSize ] ;bRefreshDynPlot = False),] ;

showDynPlot[ dynPlot, tauSliderValue ]
),
freqTab,
(
If [ bRefreshMatrix,
(matrix = dynamicsMatrix[ parameterLocatorDependencies, couplingData ] ;bRefreshMatrix = False ;bRefreshDynamics = True ; bRefreshFreqPlot = True ),];

If[bRefreshFreqPlot,
(freqPlot = frequencyPlotCalc[ parameterLocatorDependencies, matrix, freqMeshSize] ;
bRefreshFreqPlot = False ), ] ;

showFreqPlot[ qPlaneLoc, freqPlot ]
),
_,
LocatorPane[ Dynamic[u,(
u = If [ (Dimensions[#] // First)<3, 
Flatten[ {glConstLocDefault, mLocDefault}, 1], 

(*cap the number of Locators at glConstMaxNumberOfMasses (4-5), so the selection Pane doesn't start truncating the 'to neighbouring' SetterBar *)
If[ (Dimensions[#] // First) > glConstMaxNumberOfMasses + 2, Take[#,glConstMaxNumberOfMasses], #] ] ;

Module[{rm,km},
rm = First@Dimensions[rawArrayOfMasses] ;

If [rm > nMassesFromLocArray[u],
(*Shrink rawArrayOfMasses *)
rawArrayOfMasses = Take[ rawArrayOfMasses, nMassesFromLocArray[u] ],

(*Grow rawArrayOfMasses if required*)
While[rm< nMassesFromLocArray[u],AppendTo[rawArrayOfMasses, glConstDefaultMass] ; rm++ ]
] ;

km = kArray[[All,1]][[All,1]] // Max ;

If [km > nMassesFromLocArray[u],
(*Shrink kArray *)
kArray = Select[kArray, (Max[{#[[1]], #[[2]]}] ≤  nMassesFromLocArray[u]&)],
(*Grow kArray if required*)
While[km< nMassesFromLocArray[u],km++ ;kArray = returnKarrayForAdditionalMassIndex[kArray,km] ] 
]
] ;

(* Don't allow the lattice vector end points to be too close to the origin *)
u[[1]] = If[ u[[1]]. u[[1]] < glConstMinSquaredDistanceFromOrigin, glConstLocDefault[[1]], u[[1]] ] ;
u[[2]] = If[ u[[2]]. u[[2]] < glConstMinSquaredDistanceFromOrigin, glConstLocDefault[[2]], u[[2]] ] ;

(* Don't allow the angle between lattice vectors get too small *)
{u[[1]], u[[2]]} = resetLatticeVectorsIfAngleTooSmall[ u,  glConstMinAngleBetweenLatticeVectors ] ;



rawNumberOfMasses = nMassesFromLocArray[u] ;

(*These are in case the number of locators were changed and we have a mass selected that is now invalid.*)
m1Sel = If [ m1Sel >rawNumberOfMasses, 1, m1Sel] ;
m2Sel = If [ m2Sel >rawNumberOfMasses, 1, m2Sel] ;
moSel = If [ moSel >rawNumberOfMasses, 1, moSel ] ;
moSel = adjustedmoSel[ rawNumberOfMasses, moSel, m1Sel ] ;

omegaIndex = If [ omegaIndex >(2rawNumberOfMasses), 1, omegaIndex] ;

(*Adjust these for Slider display in case m1Sel, m2Sel, or moSel were changed just above.*)
k1 = selectVFromKarray[kArray, m1Sel, m2Sel, 1] ;
k2 = selectVFromKarray[kArray, m1Sel, m2Sel,2] ;
k3 = selectVFromKarray[kArray, m1Sel, m2Sel,3] ;
k4 = selectVFromKarray[kArray, m1Sel, m2Sel,4] ;
k5 = selectOriginVFromKarray[ kArray, m1Sel, moSel ] ;

parameterLocatorDependencies = locDependent[ u, rawNumberOfMasses,rawArrayOfMasses ] ;

(* Reset in case a,or b changed, and the last selected value is now out of range. *)
qLoc = ("qMaxVec"/.parameterLocatorDependencies )/2  ;
qPlaneLoc = {-qLoc,qLoc} ;

couplingData = calculateCouplings[ parameterLocatorDependencies, kArray ] ;

bRefreshMatrix = True ;

tick=Not[tick])&],
(* Why doesn't Alt-click to remove existing Locator refresh this plot?  Workaround: move one of the other locators to refresh *)
plotSprings[u, parameterLocatorDependencies,couplingData,m1Sel, m2Sel, moSel] ,
LocatorAutoCreate->True, 
ContinuousAction->False
] (*LocatorPane*)
] (*Switch*)

(* Here is the declaration of the dynamics, the controls for which are in this first "variable" declaration. *)
,
TabView[
{
"dynamics" ->  Column[ { 
Dynamic@Row[{
OverVector[Style["q",Italic, 10]],
" = ",
Module[{qm},

qm = ("qMaxVec"  /.parameterLocatorDependencies) ;
Slider2D[
Dynamic[qLoc,(
qLoc=#;
bRefreshDynamics = True ;

tick=Not[tick])&], (*{{0.01,0.01}, qm}*) ({-#,#}/2 & @ qm)
,ContinuousAction->False
, ImageSize -> (150 #/Max[#] &) @ qm
]
], 
Style[" 2 π ", 10],  Style[Dynamic[(NumberForm[(qLoc/(2 Pi) // MatrixForm), {3,2}] )], 10]
}],

Style[
Row[{ Text@"number of time intervals ",
Manipulator[Dynamic[dynMeshSize,(
dynMeshSize = checkAgainstRange[#, 14, 4, 30, IntegerQ];
bRefreshDynPlot = True ;
tick=Not[tick])&],{4,30,1},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}]," ", Dynamic[dynMeshSize] 
}],
10],

Style[
Row[{Text@"time, normalized to one period, ",Style["t",Italic],"/",Style["T",Italic]," = ", Dynamic[N[(tauSliderValue-1)/dynMeshSize]]}]
,10],

Row[{
Dynamic@Manipulator[Dynamic[tauSliderValue,(
tauSliderValue = # ;
tick=Not[tick])&],{1,dynMeshSize+1,1},ImageSize->Tiny,ContinuousAction->True,Appearance->"Open",
AppearanceElements -> {"StepLeftButton","PlayPauseButton","StepRightButton","FasterSlowerButtons","DirectionButton"
}
] 
}
(*, ImageSize->{200,60}*)
]
,
Style[
Row[{Text@"oscillation scale factor:"}]
,10],

Style[
Row[{Manipulator[Dynamic[scale,(
scale=checkAgainstRange[#, 0.05, 0.05, 1];
bRefreshDynPlot = True ;
tick=Not[tick])&],{0.05,1},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}],
" ",  Dynamic[(NumberForm[scale, {3,2}] )]
}],
10]
,


Style[Row[{
Text@"angular frequency ω(",Style["q",Italic],Text@"), selection: ",
Dynamic@SetterBar[Dynamic[omegaIndex,(
omegaIndex=#;
bRefreshDynPlot = True ;

tick=Not[tick])&], Range[2rawNumberOfMasses] ]
} ],10]


}] (* Column *)
,"ω((OverVector[q]))" ->  Column[ {


Style[
Row[{ Text@"number of mesh intervals ",
Manipulator[Dynamic[freqMeshSize,(
freqMeshSize=checkAgainstRange[#, 12, 2, 30, EvenQ];
bRefreshFreqPlot = True ;
tick=Not[tick])&],{2,30,2},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}]," ", Dynamic[freqMeshSize] 
}]
,10]

,

LocatorPane[Dynamic[qPlaneLoc(*, (qPlaneLoc = #; tick=Not[tick] &)*)],
Dynamic[
Graphics[{
Lighter@Gray,
Polygon[
({{-#[[1]],-#[[2]]},{#[[1]],-#[[2]]},{#[[1]],#[[2]]},{-#[[1]],#[[2]]}} )
 ],
Black,Line[qPlaneLoc]
}
, PlotRange -> {{-#[[1]], #[[1]]},{-#[[2]], #[[2]]}} /2
,ImageSize ->  200 #/Max[#]
] & @ ("qMaxVec" /.parameterLocatorDependencies)
]
,
ContinuousAction->False
]
} ]
, "parameters"->Grid[ {
{
Style[
Row[{
Text["mass: "], 
Dynamic@If [ rawNumberOfMasses > 1,
SetterBar[
Dynamic[m1Sel,(
m1Sel=#;
massValue = ("mScalarArray" /. parameterLocatorDependencies)[[m1Sel]];

k1 = selectVFromKarray[kArray, m1Sel, m2Sel, 1] ;
k2 = selectVFromKarray[kArray, m1Sel, m2Sel,2] ;
k3 = selectVFromKarray[kArray, m1Sel, m2Sel,3] ;
k4 = selectVFromKarray[kArray, m1Sel, m2Sel,4] ;

moSel = adjustedmoSel[ rawNumberOfMasses, moSel, m1Sel ] ;
k5 = selectOriginVFromKarray[ kArray, m1Sel, moSel ] ;

tick=Not[tick])&], 
(# -> massLabel[ # ]&/@ Range[rawNumberOfMasses])
],
massLabel[ 1]
]
}]
,10],
Style[
Row[{
Manipulator[
Dynamic[massValue,(
massValue = checkAgainstRange[#, 0.25, 0.25, glConstMMax] ;
rawArrayOfMasses[[m1Sel]]=massValue ;
parameterLocatorDependencies = (parameterLocatorDependencies /. (Rule["mScalarArray", _] :> Rule["mScalarArray", rawArrayOfMasses]) ) ;
bRefreshMatrix = True ;
tick=Not[tick])&]
,{0.25,Dynamic@glConstMMax},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}
]
," ",Dynamic[NumberForm[massValue,{2,1}]]
}],
10]
}
,{Text@"\n", ""}
, {
Text[Style["coupling to neighboring: " ,10]], 
Style[
Dynamic@If [ rawNumberOfMasses > 1,
SetterBar[
Dynamic[m2Sel,(
m2Sel=#;

k1 = selectVFromKarray[kArray, m1Sel, m2Sel, 1] ;
k2 = selectVFromKarray[kArray, m1Sel, m2Sel,2] ;
k3 = selectVFromKarray[kArray, m1Sel, m2Sel,3] ;
k4 = selectVFromKarray[kArray, m1Sel, m2Sel,4] ;
tick=Not[tick])&], 
(# -> massLabel[ #] &/@ Range[rawNumberOfMasses])
],massLabel[ 1]
],
10]
}
, {


Style[
Row[{ Text@"horizontal: ",kLable[1], " || ", OverVector[Style["a",Italic]]," "}]
,10]
,
Style[
Row[{
Manipulator[
Dynamic[k1,(
k1=checkAgainstRange[#, glConstKMin, glConstKMin, glConstKMax] ;
kArray = returnAlteredKArrayElements[ kArray, m1Sel, m2Sel, 1, k1 ] ;
couplingData = calculateCouplings[ parameterLocatorDependencies, kArray ] ;
bRefreshMatrix = True ;
tick=Not[tick])&]
,{Dynamic@glConstKMin,Dynamic@glConstKMax},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}
]
,
" ",Dynamic[NumberForm[k1,{3,2}]]
}],
10]
}
, {


Style[
Row[{ Text@"vertical: ",kLable[2], " || ", OverVector[Style["b",Italic]]," "}]
,10]

,



Style[
Row[{Manipulator[
Dynamic[k2,(
k2=checkAgainstRange[#, glConstKMin, glConstKMin, glConstKMax];
kArray = returnAlteredKArrayElements[ kArray, m1Sel, m2Sel,  2, k2 ] ;
couplingData = calculateCouplings[ parameterLocatorDependencies, kArray ] ;
bRefreshMatrix = True ;
tick=Not[tick])&]
,{Dynamic@glConstKMin,Dynamic@glConstKMax},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}
]
," ",Dynamic[NumberForm[k2,{3,2}]]
}] ,
10]
}
, {


Style[
Row[{ Text@"diagonal: ",kLable[3], " || (",  OverVector[Style["b",Italic]], " + ", OverVector[Style["a",Italic]], ") "}]
,10]
,
Style[
Row[{Manipulator[
Dynamic[k3,(
k3=checkAgainstRange[#, glConstKMin, glConstKMin, glConstKMax];
kArray = returnAlteredKArrayElements[ kArray, m1Sel, m2Sel,  3, k3 ] ;
couplingData = calculateCouplings[ parameterLocatorDependencies, kArray ] ;
bRefreshMatrix = True ;
tick=Not[tick])&]
,{Dynamic@glConstKMin,Dynamic@glConstKMax},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}
]
," ",Dynamic[NumberForm[k3,{3,2}]]
}]
,10]
}
, {

Style[
Row[{ Text@"diagonal: ",kLable[4], " || (",  OverVector[Style["b",Italic]], " - ", OverVector[Style["a",Italic]], ") "}]
,10]
,
Style[
Row[{Manipulator[
Dynamic[k4,(
k4=checkAgainstRange[#, glConstKMin, glConstKMin, glConstKMax];
kArray = returnAlteredKArrayElements[ kArray, m1Sel, m2Sel,  4, k4 ] ;
couplingData = calculateCouplings[ parameterLocatorDependencies, kArray ] ;
bRefreshMatrix = True ;
tick=Not[tick])&]
,{Dynamic@glConstKMin,Dynamic@glConstKMax},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}
]
," ",Dynamic[NumberForm[k4,{3,2}]]
}],
10]
}

,{Text@"\n", ""}
,
{
Dynamic@If[ rawNumberOfMasses == 1
,""
, Style[
Row[{ Text@"coupling from ", Dynamic@massLabel[ m1Sel],Text@" to: "}]
, 10]
],
Dynamic@Switch[rawNumberOfMasses,
1,"",
2, Style[Dynamic@massLabel[ First@Delete[ Range[2], m1Sel] ], 10],
_, Dynamic@SetterBar[
Dynamic[moSel,(
moSel=#;
k5 = selectOriginVFromKarray[ kArray, m1Sel, moSel ] ;
tick=Not[tick])&], 
(# -> massLabel[ #] &/@ (DeleteCases[Range[rawNumberOfMasses],m1Sel]))
]
]
},
{
Dynamic@If[ rawNumberOfMasses == 1,"",Style[
Row[{kLable[5], ": "}], 10]
],
Dynamic@If[ rawNumberOfMasses == 1,"",
Style[
Row[{Dynamic@Manipulator[
Dynamic[k5,(
k5=checkAgainstRange[#, glConstKMin, glConstKMin, glConstKMax];
kArray = returnAlteredKarrayOriginElement[ kArray, m1Sel, moSel, k5 ] ;
couplingData = calculateCouplings[ parameterLocatorDependencies, kArray ] ;
bRefreshMatrix = True ;
tick=Not[tick])&]
,{Dynamic@glConstKMin,Dynamic@glConstKMax},ImageSize->Tiny,ContinuousAction->False,AppearanceElements->{"InputField"}
]
," ",Dynamic@NumberForm[ k5,{3,2}]
}],
10]
]
}
} ] (*Grid*)
}
, Dynamic[tabNumber, (tabNumber = # ;tick=Not[tick])&]
, ImageSize->{250,350} 
] 

(* Localized variable declarations using ControlType -> None Manipulators. *)
,{{tick,False},None}
,{{qPlaneLoc,{}}, None}
,{{qLoc,{}},None}
,{{tauSliderValue,1},None}
,{{omegaIndex, 1}, None}

,{{scale, 0.2}, None}
,{{freqMeshSize,12},None}
,{{dynMeshSize,14},None}

,{{matrix, {}}, None}
,{{dynamics, {}}, None}
,{{dynPlot, {}}, None}
,{{freqPlot, {}}, None}
,{{bInitDefaults, True}, None}
,{{bRefreshMatrix, False}, None}
,{{bRefreshDynamics, False}, None}
,{{bRefreshDynPlot, False}, None}
,{{bRefreshFreqPlot, False}, None}

,{{tabNumber,3},None}
,{{dynTab, 1}, None}
,{{freqTab, 2}, None}
,{{parametersTab, 3}, None}

(* The number of default masses are coupled by these next four variables. *)
(*,{{rawNumberOfMasses, 3}, None}
,{{mLocDefault, {{0.9,0.7},{0.4, 1.2}, {1.3, 1.4}}}, None}
,{{rawArrayOfMasses, {10, 15, 20}}, None}*)
,{{rawNumberOfMasses, 1}, None}
,{{mLocDefault, {{0.9,0.7}(*,{0.4, 1.2}, {1.3, 1.4}*)}}, None}
,{{rawArrayOfMasses, {10(*, 15, 20*)}}, None}
,{{moSel, 2}, None} (* pick anything that isn't m1Sel *)

,{{m1Sel, 1}, None}
,{{m2Sel, 1}, None}

,{{kArray, {}}, None} 
,{{u,{}}, None}
,{{parameterLocatorDependencies, {}}, None}
,{{couplingData, {}}, None}
,{{glConstSpringColorsByN, {}}, None}

,{{k1, 0.25},None}
,{{k2, 0.25},None}
,{{k3, 0.25},None}
,{{k4, 0.25},None}
,{{k5, 0.25},None}

,TrackedSymbols:>{tick}
,ControlPlacement->Left
,SaveDefinitions->True
]

A lattice of atoms can be modeled as harmonic oscillators, with forces proportional to the displacements of the atoms from equilibrium positions. The simplest such model introduces coupling only between nearest-neighbor atoms. In this Demonstration, a lattice cell containing one to five atoms is modeled, with nearest-neighbor harmonic coupling to the masses in each nearby cell. Normal mode solutions to these equations of motion are plotted. Controls are provided to alter the coupling "spring constants" and other free parameters, as well as controls to select from the reciprocal space vectors and angular frequencies associated with the normal mode solutions. A time control is also provided to display changes of the lattice through one period of the lattice vibration. A plot of the dispersion relation, showing the angular velocities associated with each reciprocal vector, is also provided.

The positions of masses within a two-dimensional periodic array of cells can be described by summing the lattice vector (Subscript[(OverVector[r]), (OverVector[n]) = ((Subscript[n, 1]), (Subscript[n, 2]))]) =(Subscript[n, 1])(OverVector[a]) +(Subscript[n, 2 ])(OverVector[b]), representing the origin of each of the lattice cells, and a relative vector to the position of each of the masses. With (Subscript[(OverVector[p]), k]) representing the equilibrium position of the (k^th) mass in cell (Subscript[(OverVector[r]), (OverVector[n])]), the position of that mass is (Subscript[(OverVector[r]), (OverVector[n])])+(Subscript[(OverVector[p]), k]).

Let (Subscript[(OverVector[a]), (OverVector[n]),k;(OverVector[m]),j])=(Subscript[(OverVector[r]), (OverVector[n])])+(Subscript[(OverVector[p]), k])-(Subscript[(OverVector[r]), (OverVector[m])])-(Subscript[(OverVector[p]), j]), with direction (Subscript[(OverHat[a]), (OverVector[n]),k;(OverVector[m]),j]), represent the equilibrium separation of the (k^th) mass in cell (Subscript[(OverVector[r]), (OverVector[n])]) from the (j^th) mass in cell (Subscript[(OverVector[r]), (OverVector[m])]). If the harmonic coupling between these masses has magnitude (Subscript[K, (OverVector[n]),k;(OverVector[m]),j]), then the system of equations describing the vector displacement (Subscript[(OverVector[u]), n,k]) for the (k^th) mass in unit cell (OverVector[n]) from the equilibrium position is given by

(Subscript[m, k]) (Subscript[(Overscript[(OverVector[u]), ..]), (OverVector[n]),k])=-(∑_((OverVector[n]),k ≠(OverVector[m]),j))(Subscript[K, (OverVector[n]),k;(OverVector[m]),j])(Subscript[Proj, (Subscript[(OverHat[a]), (OverVector[n]),k;(OverVector[m]),j])])((Subscript[(OverVector[u]), (OverVector[n]),k])-(Subscript[(OverVector[u]), (OverVector[m]),j])).

In general, we have one such equation for each (OverVector[n]), k pair. A trial solution of the form: (Subscript[(OverVector[u]), (OverVector[n]), k ])(t)=(((Subscript[(OverVector[ϵ]), k])((OverVector[q])))/(Sqrt[(Subscript[m, k])])) (e^(I ((Subscript[(OverVector[r]), (OverVector[n])]). (OverVector[q]) -ω t))) can be used to decouple this system, resulting in a single equation for each (k^th) mass of the form

(ω^2) (Subscript[(OverVector[ϵ]), k])((OverVector[q])) = (∑_((OverVector[n]),k ≠ (OverVector[m]),j))(Subscript[K, (OverVector[n]),k;(OverVector[m]),j])(Subscript[Proj, (Subscript[(OverHat[a]), (OverVector[n]),k;(OverVector[m]),j])]) ((((Subscript[(OverVector[ϵ]), k])((OverVector[q])))/(Subscript[m, k]))- (((Subscript[(OverVector[ϵ]), j])((OverVector[q])))/(Sqrt[(Subscript[m, k])(Subscript[m, j])]))(e^(I ((Subscript[(OverVector[r]), (OverVector[m])]) - (Subscript[(OverVector[r]), (OverVector[n])])). (OverVector[q])))).

This describes all the steady-state lattice vibrations, the normal modes, that can be modeled by this trial solution. Here (OverVector[q]) is a vector in reciprocal space, effectively parameterizing the angular velocity ω=ω((OverVector[q])). The vector (Subscript[(({{(Subscript[(OverVector[ϵ]), k])((OverVector[q]))}})), k]) is an eigenvector of the equations of motion of the system for this assumed solution, where (ω^2) are the eigenvalues of this system. For an N-atom basis, there are 2 N such (ω^2) eigenvalues per (OverVector[q]) point, each resulting in a different characteristic motion.

The rank of the resulting eigenvalue problem depends on the number of masses per unit cell, but the complexity of the matrix expression depends on the number of neighboring interactions that are considered. For example, given lattice vectors (OverVector[a]), (OverVector[b]), diagonals (OverVector[r])=(OverVector[a]) + (OverVector[b]), (OverVector[s])=(OverVector[a]) - (OverVector[b]) , and a one-atom basis, where each unit cell contains a single mass coupled with harmonic oscillator forces between only nearest neighbors, the normal mode solutions follow from the solution of the eigenvalue problem

(({{(ω^2), 0}, {0, (ω^2)}})) (OverVector[ϵ]) =(4/m)( (Subscript[k, 1 ])(OverHat[a]) ((OverHat[a])^T) (sin^2)( (OverVector[a]) . (OverVector[q])/2)+(Subscript[k, 2 ])(OverHat[b]) (OverHat[b])ᵀ(sin^2)( (OverVector[b]) . (OverVector[q])/2)+(Subscript[k, 3])(OverHat[r]) (OverHat[r])ᵀ(sin^2)(((OverVector[b]) +(OverVector[a]) ). (OverVector[q])/2)+(Subscript[k, 4])(OverHat[s]) (OverHat[s])ᵀ(sin^2)(( (OverVector[b])-(OverVector[a])). (OverVector[q])/2 ))(OverVector[ϵ]).

Controls are provided to display the dynamics associated with each of the characteristic angular frequencies ω, for given reciprocal vector values (OverVector[q]).

Three tabs are provided in this Demonstration. The primary tab displays the dynamics of the solution for particular pair of ((OverVector[q]),ω((OverVector[q]))) values. In that tab, selecting "Play" for the time control animates the lattice vibrations. A scaling control is provided to alter the initial magnitude of the eigenvectors, tantamount to picking the initial time boundary value constraints. Note that it is possible to select physically unrealistic scaling factors that allow for collisions that are not modeled by this system.

A second tab provides the dispersion relation, the dependence of angular velocity ω((OverVector[q])) on all (OverVector[q]) points.

Finally, a parameters tab provides controls for the spring constants (Subscript[K, (OverVector[n]),k;(OverVector[m]),j]), the primitive unit cell lattice vectors (OverVector[a]),(OverVector[b]), and the positions of the masses (Subscript[(OverVector[p]), k]) within each unit cell of the lattice. Additional mass position locators, up to four total, may be added or deleted by Alt-clicking on the lattice cell in the desired location. For mass locator removal, it is necessary to adjust one of the other locators so that the removal takes effect. Note that it is possible to select equilibrium mass positions that are too close, leading to physically unrealistic dynamics, such as masses passing through each other.

The total number of interactions, even when those interactions are restricted to just the neighboring lattice cells, increases quickly as additional masses per cell are added. For example, with three masses per unit cell, considering only the neighboring and origin cells, there are 29 interactions possible for each mass. To simplify the physical constants dialog, distinct "spring constant" selection is only available for a subset of the possible interactions. This imposes the following respective constraints along the "horizontal," "vertical," "NE diagonal," and "NW diagonal" directions connecting the masses across the cells

(Subscript[K, (OverVector[0]),k;(0,1),j])= (Subscript[K, (OverVector[0]),k;(0,-1),j]), (Subscript[K, (OverVector[0]),k;(1,0),j])= (Subscript[K, (OverVector[0]),k;(-1,0),j]), (Subscript[K, (OverVector[0]),k;(1,1),j])= (Subscript[K, (OverVector[0]),k;(-1,-1),j]), (Subscript[K, (OverVector[0]),k;(1,-1),j])= (Subscript[K, (OverVector[0]),k;(-1,1),j]).

These are labeled (Subscript[k, 1]), (Subscript[k, 2]), (Subscript[k, 3]), and (Subscript[k, 4]) in the parameters tab, respectively. For example, the first identity above is the imposition of an equality constraint on the coupling constants in the "left" and "right" directions separating masses in neighboring "horizontal" cells (when the lattice is square). Those are the interactions directed primarily along the lattice vector directions ±(OverVector[a]) (primarily since these directions are also adjusted for the positions of the masses within the respective cells if different). Similarly, the same coupling constants are used for each of the pairs of directions that are directed primarily along the ±(OverVector[b]), ±((OverVector[a])+(OverVector[b])), ±((OverVector[a])-(OverVector[b])) directions separating the lattice cells.

When there is more than mass per unit cell, distinct coupling constants for the intra-cell interactions between the masses within the origin unit cell may be selected (i.e. (Subscript[K, (OverVector[0]),k;(OverVector[0]),j])). These are labeled (Subscript[k, 5]) in the parameters tab.

General theory describing oscillations around lattice equilibrium points can be found in [1].

Reference

[1] N. W. Ashcroft and N. D. Mermin, Solid State Physics, New York: Holt, Rinehart and Winston, 1976, chapters 21 and 22.

Contributed by: Peeter Joot