myBZ[] Brillouin Zone Plotter and myRed2Ext[] scripts are written by Csaba Géza Péterfalvi.
Sent to the Wolfram Library Archive on 18th June 2014.
http://peterfalvi.web.elte.hu

Options[myBZ]={origin->{0,0,0},join->{},print->True,kPoints->True,opacity->-1,resolution->0.01,hight->350,colors->"Rainbow",normLimit->10$MachineEpsilon};
myBZ[g_,csn_,csx_,csl_,OptionsPattern[]]:=Module[{or=OptionValue[origin],join=OptionValue[join],pr=OptionValue[print],kpp=OptionValue[kPoints],
op=OptionValue[opacity],res=OptionValue[resolution],hght=OptionValue[hight],clrs=OptionValue[colors],nl=OptionValue[normLimit],(* nl is supposed to be smaller than any number and distance that matters in this calculation (<<res), but large compared to any possible numerical inaccuracy so that these can be avoided.
It is used at various places independently, and its value could be adjusted separately too, if needed. *)
mn=ToString[$ModuleNumber-1],csn0=Normalize[csn],csx0,csy0,cshh,startY,org,gx,gPoints,kPoints={},gLines,bzi,bzid,bzio=0,bzios=0,bzix,bzin,xy,dir,bx,iy,bzFig,figure,clrBar,string},
csx0=Normalize[csx-csn0(csn0.csx)];csy0=csn0\[Cross]csx0;(* Unit vectors for the xy-coordinates (Just in case, csx is projected to the plane first.) *)
org=(or=={0,0,0});(* Is Γ the origin? *)
gx=Ceiling[2(Norm[csl]+Norm[or])Max[Norm[g[[#[[1]]]]\[Cross]g[[#[[2]]]]]&/@{{1,2},{1,3},{2,3}}]/Det[g]];
(* Here, we determine how many gPoints do we need to find all the BZ borders within the region around the origin: 
(Norm[csl]+Norm[or]) is the maximum possible distance of the furthest separator, and twice this much is that of the furthest gPoint. 
Let us divide this with the smallest thickness of the parallelepiped spanned by the 'g' vectors. This leads us to the maximum surface area divided by the volume. *)
gPoints=Flatten[Table[N[g\[Transpose]].{i,j,k},{i,-gx,gx},{j,-gx,gx},{k,-gx,gx}],2];
gPoints=Select[gPoints,(Norm[#]≤2(Norm[csl]+Norm[or]))&];(* We don't need gPoints that are too far away. *)
If[kpp,kPoints=Select[(#-or)&/@gPoints,(Norm[csn0.#]<nl)&];];(* kPoints are gPoints in The Plane (already relative to the origin. *)
gPoints=Complement[gPoints,{{0.,0.,0.}}];(* We don't need the Γ-point here either. *)
If[!org,(* If the origin is not the Γ point, we'll need bzio, bzios and maybe shift some gPoints... *)
gPoints=Table[If[(Norm[csn0\[Cross]gp]>nl)&&(Abs[gp.(or-gp/2)]<nl),gp+Normalize[gp]res/2,gp],{gp,gPoints}];
(* The gPoints whose separator plane contains the origin (but not The Plain itself) are shifted away from Γ by res/2 to lift the degeneracy. *)
bzio=Total[Table[If[gp.(or-gp/2)>nl,1,0],{gp,gPoints}]];
(* The BZ Index of the origin. We assume that if The Plane is also a BZ border that will not appear in this sum because of the requirement >nl. *)
If[Norm[Round[#]-#]&[2(csn0.or)(csn0.Inverse[g])]<nl,bzios=-Sign[csn0.or];If[bzios==-1,bzio++;];];];(* bzios=-1/0/+1 => The BZ-index right under The Plane (that is also a BZ border now if bzios≠0!) is one less/the same/one more. *)
gPoints=Select[gPoints,(Norm[csn0\[Cross]#]>nl)&];(* We don't need the parallel planes. *)
gPoints=Table[{(Norm[gp]/2-or.#)(#-(csn0.#)csn0)/(1-(csn0.#)^2),#}&[Normalize[gp]],{gp,gPoints}];(* The intersection of The Plane and the separator plane of Γ and the gPoint has its closest point to 'or' at these new coordinates (relative to 'or').
Here we also store the normal vector to the gPoint. *)
gPoints=Table[{gp[[1]],Sign[gp[[1]].gp[[2]]]},{gp,gPoints}];
(* Here we add the sign of change of the BZ Index on crossing the border: + > regular, - > decreasing. *)
gPoints=Select[Table[{csx0.gp[[1]],csy0.gp[[1]],gp[[2]]},{gp,gPoints}],(Min[Flatten[Table[Norm[#[[1;;2]]-{i csl[[1]],j csl[[2]]}/2],{i,{-1,1}},{j,{-1,1}}]]]<Norm[csl]/2)&];(* These are the in-plane coordinates of the separator points already in the xy-system within the Thales' circles determined by the 'csl' intervals. *)
If[kpp,kPoints=Select[Table[{csx0.kp,csy0.kp},{kp,kPoints}],((Abs[#[[1]]]≤csl[[1]])&&(Abs[#[[2]]]≤csl[[2]]))&];];
(* The kPoints' coordinates are transformed here to the system of The Plane and those are selected that are within the given region. *)
gLines=Table[Piecewise[{{{{-csl[[1]],gp[[2]]},{csl[[1]],gp[[2]]}},gp[[1]]==0.},{{{gp[[1]],-csl[[2]]},{gp[[1]],csl[[2]]}},gp[[2]]==0.}},Sort[{{-csl[[1]],-((-csl[[1]] gp[[1]]-((gp[[1]])^2)-((gp[[2]])^2))/(gp[[2]]))},{csl[[1]],-((csl[[1]] gp[[1]]-((gp[[1]])^2)-((gp[[2]])^2))/(gp[[2]]))},{-((-((gp[[1]])^2)-csl[[2]] gp[[2]]-((gp[[2]])^2))/(gp[[1]])),-csl[[2]]},{-((-((gp[[1]])^2)+csl[[2]] gp[[2]]-((gp[[2]])^2))/(gp[[1]])),csl[[2]]}},(#1[[1]]>#2[[1]])&][[2;;3]]],{gp,gPoints}];(* Here, we generate the end-points of the separator lines. *)
(* Summing up the BZI *)
bzid=2 Floor[csl/res]+{1,1};(* The size of the discrete grid for BZI. *)
cshh=(bzid[[2]]+1)/2;startY=If[org,1,cshh];
bzi=If[org,ConstantArray[1,{bzid[[1]],cshh}],ConstantArray[1,bzid]];(* Initializing BZI *)
Do[xy=Sort[gpl[[4;;5]],(#1[[2]]<#2[[2]])&];(* Sorting the end-points of the separator lines *)If[(gpl[[4]]≠gpl[[5]])&&((!org)||(xy[[2,2]]>0)),(* (In case the separator line just touches the corner of the rectangle,
which should not happen btw because we used < instead of ≤ at the Thales' circles.) && (If Γ is the origin, for symmetry reasons, we can only deal with lines that are at least partly above the x-axis.) *)
dir=Sign[gpl];(* Which quarter are we in? *)
If[Abs[gpl[[1]]]<nl,(* We have a horizontal line... (BZI increases if gpl[[3]]=1 or decreases if-1 on crossing the border arriving from the origin.) *)
If[dir[[2]]==1,bzi[[All,startY+Ceiling[xy[[1,2]]/res];;-1]]+=gpl[[3]];,(* in the top half *)
bzi[[All,1;;Floor[(csl[[2]]+xy[[1,2]])/res]+1]]+=gpl[[3]];];,(* in the bottom half *)
If[Abs[gpl[[2]]]<nl,(* or do we have a vertical line? *)If[dir[[1]]==1,bzi[[Ceiling[(csl[[1]]+xy[[1,1]])/res]+1;;-1,All]]+=gpl[[3]];,(* in the right half *)bzi[[1;;Ceiling[(csl[[1]]+xy[[1,1]])/res],All]]+=gpl[[3]];];,(* in the left half *)
(* or are we left with the general situation, with a separator at coordinates given by the function bx[iy]? *)bx[iy_]:=(-xy[[1,2]] xy[[2,1]]+xy[[1,1]] xy[[2,2]]-xy[[1,1]] (iy-startY) res+xy[[2,1]] (iy-startY) res)/(xy[[2,2]]-xy[[1,2]]);
If[gpl[[2]]>0,For[iy=Max[1,startY+Ceiling[xy[[1,2]]/res]],iy<startY+cshh,iy++,(* With a separator point in the *)If[dir[[1]]==1,bzi[[Max[1,Ceiling[(csl[[1]]+bx[iy])/res]+1];;-1,iy]]+=gpl[[3]];,(* top-right *)bzi[[1;;Min[Ceiling[(csl[[1]]+bx[iy])/res],bzid[[1]]],iy]]+=gpl[[3]];];];,(* top-left quarter, *)For[iy=1,iy≤Floor[xy[[2,2]]/res]+startY,iy++,(* and with a separator point in the *)If[dir[[1]]==1,bzi[[Max[1,Ceiling[(csl[[1]]+bx[iy])/res]+1];;-1,iy]]+=gpl[[3]];,(* bottom-right, *)bzi[[1;;Min[Ceiling[(csl[[1]]+bx[iy])/res],bzid[[1]]],iy]]+=gpl[[3]];];];];];];];,(* bottom-left quarter. *)
{gpl,Join[gPoints\[Transpose],gLines\[Transpose]]\[Transpose]}];
bzi=If[org,Join[Reverse[bzi\[Transpose],{1,2}][[;;-2]],bzi\[Transpose]],bzi\[Transpose]+bzio];(* Applying symmetry/Adding BZI of the origin *)
{bzin,bzix}={Min[bzi],Max[bzi]};
(* Constructing the output figure *)
gLines={Pick[gLines,gPoints[[All,3]],1],Pick[gLines,gPoints[[All,3]],-1]};
If[kpp,kPoints=Point/@kPoints];
bzFig=Show[{ListContourPlot[bzi,DataRange->{{-csl[[1]],csl[[1]]},{-csl[[2]],csl[[2]]}},PlotRange->{Full,Full,Full},
Contours->Function[{min,max},Range[min,max-1]+0.5],PerformanceGoal->"Quality",ContourStyle->None,ColorFunction->(Directive[Opacity[If[op==-1,If[Length[join]==0,1,0.5],op]],ColorData[clrs][#]]&),AspectRatio->Automatic,ImageSize->{Automatic,hght}],Graphics[{Line/@gLines[[1]],{Red,Line/@gLines[[2]]},If[bzios≠0,{Thick,Switch[bzios,-1,Blue,1,Red],
Line[{{-csl[[1]],-csl[[2]]},{-csl[[1]],csl[[2]]},{csl[[1]],csl[[2]]},{csl[[1]],-csl[[2]]},{-csl[[1]],-csl[[2]]}}]},{}],{Green,kPoints}}],join}];
(* Colorbar *)
clrBar=ContourPlot[y,{x,0,1},{y,bzin-0.5,bzix+0.5},ContourStyle->None,AspectRatio->25,PlotRange->All,PlotRangePadding->0,
Contours->(Range[bzin,bzix-1]+0.5),FrameTicks->{{Range[bzin,bzix],None},{None,None}},ImageSize->{Automatic,hght},
ColorFunctionScaling->False,ColorFunction->(Directive[Opacity[If[op==-1,If[Length[join]==0,1,0.5],op]],
ColorData[clrs][Rescale[Piecewise[{{bzin,#<bzin},{bzix,#>bzix}},#],{bzin,bzix},{0,1}]]]&)];
(* Tooltips show the BZI of the region under the mouse pointer. *)
figure=Row[{Tooltip[bzFig,Dynamic[bzi[[Apply[Sequence,If[NumericQ[#[[1]]],If[(1≤#[[2]]≤bzid[[1]])&&(1≤#[[1]]≤bzid[[2]]),#,(Reverse[bzid]+{1,1})/2],
(Reverse[bzid]+{1,1})/2]&[Reverse[Round[(MousePosition["Graphics"]+csl)/res]]+{1,1}]]]]]],
(* The red highlighting above suggests that 'Order of evaluation conflicts'. These are all constants, so why the conflict?!
Anyway, these are not saved within the cell together with all the rest of the data for the figure.
This means that the Tooltips are not wotking in a new session, unless we explicitely take care of them.
(Is this in any way related to the 'Temporary' attribute of these variables? But even if so, it does not make any difference if I clear this attribute already within the module.) *)
(* Loading Tooltip data when clicked on the colorbar
(EventHandler[] is better here than Button[]: it is free from fancy effects. Note that practically both have the 'HoldRest' attribute, even though Attributes[EventHandler] tells it otherwise of course. The MouseClicked rule also gets evaluated first when you click!) *)
EventHandler[clrBar,{"MouseClicked":>If[!ValueQ[bzi],SelectionMove[EvaluationNotebook[],All,EvaluationCell];
string="=Uncompress[CurrentValue[NotebookSelection[],{TaggingRules,\"";
ToExpression[ToString[bzi]<>string<>ToString[bzi]<>"\"}]];"];
ToExpression[ToString[bzid]<>string<>ToString[bzid]<>"\"}]];"];
ToExpression[ToString[res]<>string<>ToString[res]<>"\"}]];"];
SelectionMove[EvaluationNotebook[],After,Cell];]}]},Spacer[5]];
(* Writing the output cell and saving the Tooltip data which is lost otherwise when the session is closed *)
NotebookWrite[EvaluationNotebook[],Cell[BoxData@ToBoxes[If[pr,Grid[{{TableForm[If[bzios≠0,(* Vectors, etc to print if required *)
Join[{"The chosen plane is a BZ border. Right under this, the BZ index is one "<>If[bzios==1,"higher.","lower."]},#],#]&
[{{"(Subscript[e, x]) = "<>ToString@N[csx0]<>"   (Subscript[e, y]) = "<>ToString@N[csy0]<>"   n = "<>ToString@N[csn0]<>"   origin = "<>ToString@N[or]},
{"(Subscript[g, 1]) = "<>ToString@N[g[[1]]]<>"   (Subscript[g, 2]) = "<>ToString@N[g[[2]]]<>"   (Subscript[g, 3]) = "<>ToString@N[g[[3]]]},{}}]]},{figure}}],figure]],"Output",TaggingRules->{("bzi$"<>mn)->Compress[bzi],("bzid$"<>mn)->Compress[bzid],("res$"<>mn)->Compress[res]}]]]

(* Try this for example *)
g={{1, -1, 0}, {0, 1, 1.84127}, {-1, 0, 1.84127}}; (* The three reciproc lattice vectors (for white Tin) *)
csn = {0, 1, 0}; (* Normal vector to the cross-section *)
csx = {1,0, 0}; (* Axis x of the cross-section *)
csl = {1.5,1}; (* Cross-section limits in x and y: {{-csl[[1]],csl[[1]]},{-csl[[2]],csl[[2]]}} *)
neworigin = {1.5,1,0}; (* A new origin of the plot, instead of Γ *)

(* Input parameters and optional arguments with their default value (except for the origin). *)
myBZ[g, csn, csx, csl, origin->neworigin, join->{}, print->True,kPoints->True, opacity->-1, resolution->0.01, hight->350, colors->"Rainbow", normLimit->10$MachineEpsilon]

Brillouin Zone Plotter

Optional arguments:
   origin -> can be any point in k-space, but it is the Γ point by default.
   join -> some other graphics to overlay the BZs. (If not specified otherwise, the opacity is 0.5 with a joined figure and 1 without that.)
   print -> the important parameters of the plot or not.
   kPoints -> plots the equivalent, in-plane Γ-points in the figure with green.
   opacity -> of the plot. (If not specified otherwise, the opacity is 0.5 with a joined figure and 1 without that.)
   resolution -> of the k-point sampling in both directions within the plane.
   hight -> is the hight of the figure in pixels.
   colors -> define the colorscheme of the plot. For colorschemes, see the Mathematica manual.
   normLimit -> is supposed to be smaller than any number and distance that matters in the calculation (<< res), but large compared to any possible numerical inaccuracy so that these can be avoided. It is used at various places independently, and its value could be adjusted separately too, if needed.

If the plane of the plot coincides with a BZ border, the plot is framed with blue or red lines, meaning that the BZ index right under the surface is one less or one more respectively. (The Tooltip numbers (and those on the colorbar) give the BZI right above the plane.)

The out-of-plane BZ borders appear as lines in the plot. In some cases of high symmetry, several lines may coincide. If there is one or more among them, crossing which (when coming from the origin) would make the BZI decrease by one, then you see a red line. Otherwise the lines are black.

For technical reasons, Tooltips don’t work by default in new sessions of Mathematica. Some data has to be saved within the cell of the plot and re-loaded in the next session. This is taken care of by the script, all you need to do is to click on the colorbar whenever you start a new session.

Sometimes the plots need lots of storage, consider saving them in separate notebook files.

Options[myRed2Ext]={maxRecursion->3};
myRed2Ext[k_,bzi_,g_,OptionsPattern[]]:=Module[{gx0,gx0max=OptionValue[maxRecursion],kSet={},gx,gPoints},
For[gx0=1,gx0≤gx0max,gx0++,(* First trying to find the solution close to the Γ point and later, we progressively increase the scope of the search. *)
kSet=Complement[Flatten[Table[k+N[g\[Transpose]].{i,j,l},{i,-gx0,gx0},{j,-gx0,gx0},{l,-gx0,gx0}],2],kSet];(* The new k-vectors to try *)
gx=Ceiling[2Max[Norm/@kSet]Max[Norm[g[[#[[1]]]]\[Cross]g[[#[[2]]]]]&/@{{1,2},{1,3},{2,3}}]/Det[g]];(* This many lattice points are needed in each directions to calculate the BZ index for all the vectors in kSet. *)
gPoints=Complement[Flatten[Table[N[g\[Transpose]].{i,j,l},{i,-gx,gx},{j,-gx,gx},{l,-gx,gx}],2],{{0.,0.,0.}}];(* These are the lattice points that determine the BZ borders. *)
(Switch[Length[#],
0,If[gx0==gx0max,Print[Text[
Style["No solution! Check the input or increase 'maxRecursion'. (Maybe we are right on a BZ border? Try to shift 'k' by $MachineEpsilon in some well chosen direction!)",Red,14]]];];,
1,Return[#[[1,1]]];,
_,Print[Text[Style["More than one solution!",Red,14]],"\n",#];Return[Null];])&[Select[Table[{kp,Total[(UnitStep[#.(kp-#/2)])&/@gPoints]},{kp,kSet}],(#[[2]]==bzi-1)&]]]];
(* If one of the vectors in kSet has the right BZ index, then we are done and return with the new vector, otherwise we go on with the search, or throw an error message. *)

myRed2Ext[k, bzi, g, maxRecursion -> 3] transformes the coordinates of the vector ‘k’ to the Extended Zone Scheme from the Reduced Zone Scheme. ‘bzi’ is the Brillouin Zone Index of the point ‘k’ and ‘g’ is a matrix composed by the basis vectors of the crystal lattice: g={r1,r2,r3}. If you increase the optional parameter ‘maxRecursion’, you can find the solution even for larger BZ indeces, but at the same time, you get an error message later when entering a wrong input.