2
$\begingroup$

I am trying to set up Brillouin zone and representing the lines that are formed by the intersection of two planes is part of my step. I have tried the sample sections from the graphics guide book, but I got some wrong output. Here is the code that I copied from the guidebook:

(* to save memory *) $HistoryLength = 0;

The function latticePointListSC generates all lattice points that we want to take into account. It returns all lattice points inside a sphere of radius d around the origin.

latticePointListSC[d_] :=
  With[{n = Ceiling[d]}, 
    Select[
      DeleteCases[
        Flatten[Table[{i, j, k}, {i, -n, n}, {j, -n, n}, {k, -n, n}], 2],
        (* without origin *){0, 0, 0}],
      #. # <= d^2&]]

maxDist = 4;Length[latticePoints = latticePointListSC[maxDist]]

Show[Graphics3D[Cuboid[# + 0.1 {1, 1, 1},# - 0.1 {1, 1, 1}]& /@ latticePoints],
  Axes -> True];

Next, we construct the perpendicular planes in the middle of the lines origin–latticePoint. We denote a plane in the form Plane[onePointOfThePlane, listOfTwoOrthogonalDirections].

The function toPlane generates a plane (head Plane) of the bisector plane formed by the point latticePoint.

toPlane[p:latticePoint_] :=
  Plane[latticePoint/2,
    Which[(* lattice point is on a coordinate axis *)
      p[1] == 0 && p[2] == 0, {{1, 0, 0}, {0, 1, 0}},
      p[1] == 0 && p[3] == 0, {{1, 0, 0}, {0, 0, 1}},
      p[2] == 0 && p[3] == 0, {{0, 1, 0}, {0, 0, 1}},
      (* lattice point is on a coordinate plane *)
      p[1] != 0 && p[2] != 0, {#, Cross[#, p]}&[{p[2], -p[1], 0}],
      p[1] != 0 && p[3] != 0, {#, Cross[#, p]}&[{p[3], 0, -p[1]}], 
      p[2] != 0 && [P3] != 0, {#, Cross[#, p]}&[{0, p[3], -p[2]}],
      (* lattice point is in generic position *)
      True, {#, Cross[#, p]}&[{p[2], -p[1], 0}]]]

We are now adding some planes that guarantee that polygons are divided along the symmetry planes. A “symmetry unit” (for brevity, just called unit in the following) is given by the following domain x>=0, y>= 0, z>= 0, z >= x, x >=y. This is 1/48 of the whole space. The list symmetrySlicingPlanes contains the planes that bound the unit.

symmetrySlicingPlanes =
  Plane[{0, 0, 0}, #]& /@
    {{{1, 0, 0}, {0, 1, 0}}, {{1, 0, 0}, {0, 0, 1}},{{0, 1, 0}, {0, 0, 1}},
     {{0, 0, 1}, {-1, 1, 0}}, {{0, 0, 1}, {1, 1, 0}}, {{1, 0, 0}, {0, -1, 1}}, 
     {{1, 0, 0}, {0, 1, 1}}, {{0, 1, 0}, {-1, 0, 1}},{{0, 1, 0}, {1, 0, 1}}};

Here are the symmetrySlicingPlanes. The red, slightly sticking out polygons mark one symmetry unit. We will in the following concentrate on this unit and only later in the visualization part generate all other 47 units by reflection and rotation.

With[{ε = 0.05},
  Show[
    Graphics3D[{
      (* make polygons *)
      Polygon[{#[1] - #[2, 1] - #[2, 2], 
               #[1] + #[2, 1] - #[2, 2],
               #[1] + #[2, 1] + #[2, 2], 
               #[1] - #[2, 1] + #[2, 2]}]& /@ symmetrySlicingPlanes,
      {SurfaceColor[Hue[0]],
        (* lift a bit up *)
       Map[# + {ε, ε, 2ε}&, 
         (* boundary of the unit cone *)
         {Polygon[{{0, 0, 0}, {0, 0, 1}, {1, 0, 1}}],
          Polygon[{{0, 0, 0}, {0, 0, 1}, {1, 1, 1}}],
          Polygon[{{0, 0, 0}, {1, 0, 1}, {1, 1, 1}}]},
         {-2}]}}], 
    Boxed->False]]

planes = 
  Join[toPlane /@ latticePoints, symmetrySlicingPlanes];
Take[planes, 4]

Next, we construct representations of the lines that are formed by the intersection of two planes. The presentation of the lines will be in the form Line[onePointOfTheLine, lineDirection]. We use here a two-argument version of Line, in distinction to the built-in Line, which takes one argument. The function lineOnPlane[plane1, plane2] calculates the intersection line which is located on plane1 , induced by its intersection with the plane plane2 .

(* three equations cannot be solved for four variables *)
Off[Solve::"svars"];
lineOnPlane[Plane[p_, {dir1_, dir2_}], Plane[q_, {d1_, d2_}]] := 
  Module[{eqs, sol, line, var, aux, P1, P2},
    If[
      (* are the planes parallel? *)
      Length[DeleteCases[RowReduce[{dir1, dir2, d1, d2}], {0, 0, 0}, {1}]] == 2, 
      {},
      (* calculate direction of the intersecting line *)
      eqs = Thread[p + s dir1 + t dir2 == q + u d1 + v d2];
    sol = Solve[eqs, {s, t, u, v}];
    aux = p + s dir1 + t dir2 /. sol[[1]];
    (* two points on the line *)
    {P1, P2} = {aux /. {u -> 0, v -> 0}, aux /. {u -> 1, v -> 1}};Line[P1, P2 - P1]]]

lineOnPlane[planes[[1]], planes[[5]]]

I should have gotten 

Line[{-2, 1, -1}, {0, -20, 10}]

but i got 

Line[{-2, s, -(1/2) - s/2}, {0, 0, 0}]

instead

The lines after this all got messed up beacuse of the wrong out[16] Could anyone please help? Is there a better way I could do compute the lines that are formed by the intersection of two planes?

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ This question looks very similar to the one you last posted (and that went through a lot of edits, last one by @OleksandrR.)? It seems wasteful to have people put effort into edits of multiple, nearly identical questions. $\endgroup$
    – Yves Klett
    Jun 18, 2014 at 5:38

2 Answers 2

Reset to default
8
$\begingroup$

From a first look, late at night.

In the third panel of code, aren't all of the terms really parts of a list p, and should therefore be written as p[[1]], p[[2]], etc., rather than p[1], p[2], etc. The paired single square brackets would surround the argument(s) of a function, e.g. Sin[x] or Sqrt[x]; paired double square brackets designate an element of a list. There is also a typo: what appears as [p3] is almost certainly intended to be p[[3]].

In the fifth panel of code, below the comment (*make polygons*), shouldn't all of the terms be double square brackets , such as #[[1]], #[[2]], etc. rather than #[1], #[2]? When mapping the Polygon onto symmetrySlicingPlanes you are extracting parts from each of the elements of the symmetrySlicingPlanes list, not making function evaluations with paired single square brackets.

Once I changed a lot of [ ] cases to [[ ]] cases and corrected a typo, I got the following for the Graphics3D. Is this what you expect?

Result for the Graphics3D

[I got the same final result as you got for lineOnPlane[planes[[1]], planes[[5]]], but I think printing some intermediate results leading up to that answer, such as

eqs = {-2==-(3/2)-u-6 v, s==-(1/2)+3 u-2 v, t==-1+10 v}

sol = {{t->-(1/2)-s/2, u->1/5+(3 s)/10, v->1/20-s/20}}

aux = {-2,s,-(1/2)-s/2}

Line[{-2, s, -(1/2) - s/2}, {0, 0, 0}]

should provide clues as to what is wrong and how this needs to be fixed. Too late at night now.

P.S. The printing did supply some clues, and I think the code below now gives the answer you desire, but you really need to check this yourself. I'm not familiar enough with this particular type of calculation.

latticePointListSC[d_] :=

  With[{n = Ceiling[d]},

    Select[

      DeleteCases[

        Flatten[

          Table[{i, j, k}, {i, -n, n}, {j, -n, n}, {k, -n, n}],

            2],

       (*without origin*){0, 0, 0}],

     #.# <= d^2 &]]

maxDist = 4;

Length[latticePoints = latticePointListSC[maxDist]]

(* 256 *)

g3d1 = Show[

  Graphics3D[

    Cuboid[# + 0.1 {1, 1, 1}, # - 0.1 {1, 1, 1}] & /@ latticePoints

           ],

    Axes -> True];

Raster[Image[g3d1, ImageSize -> 4*72]]  (*No way to include this graphic in answer.*)


(*Notice the many double square brackets that are in this toPlane function now*)

(*when compared with the original.*)

toPlane[p : latticePoint_] := Plane[latticePoint/2,

  Which[

    (*lattice point is on a coordinate axis*)

    p[[1]] == 0 && p[[2]] == 0,

    {{1, 0, 0}, {0, 1, 0}},

    p[[1]] == 0 && p[[3]] == 0,

    {{1, 0, 0}, {0, 0, 1}},

    p[[2]] == 0 && p[[3]] == 0,

    {{0, 1, 0}, {0, 0, 1}},

    (*lattice point is on a coordinate plane*)

    (*WATCH OUT--IN THE NEXT THREE CONDITIONALS*)

    (*I ADDED A THIRD && CONDITION THAT WAS NOT IN THE ORIGINAL*)

    (*I COULD BE WRONG--I"M NO EXPERT ON THIS*)

    p[[1]] != 0 && p[[2]] != 0 && p[[3]] == 0,

    {#, Cross[#, p]} &[{p[[2]], -p[[1]], 0}],

    p[[1]] != 0 && p[[3]] != 0 && p[[2]] == 0,

    {#, Cross[#, p]} &[{p[[3]], 0, -p[[1]]}],

    p[[2]] != 0 && p[[3]] != 0 && p[[1]] == 0,

    {#, Cross[#, p]} &[{0, p[[3]], -p[[2]]}],

    (*lattice point is in generic position*)

    True,

    {#, Cross[#, p]} &[{p[[2]], -p[[1]], 0}]
]]

symmetrySlicingPlanes =

Plane[{0, 0, 0}, #] & /@

{

{{1, 0, 0}, {0, 1, 0}}, {{1, 0, 0}, {0, 0, 1}}, {{0, 1, 0}, {0, 0,1}},

{{0, 0, 1}, {-1, 1, 0}}, {{0, 0, 1}, {1, 1, 0}}, {{1, 0, 0}, {0, -1, 1}},

{{1, 0, 0}, {0, 1, 1}}, {{0, 1, 0}, {-1, 0, 1}}, {{0, 1, 0}, {1, 0, 1}}

};

g3d2 = With[{\[CurlyEpsilon] = 0.05},

  Show[

  Graphics3D[

    {(*make polygons*)

    Polygon[{#[[1]] - #[[2, 1]] - #[[2, 2]], #[[1]] + #[[2, 1]] - #[[ 2, 2]], 

      #[[1]] + #[[2, 1]] + #[[2, 2]], #[[1]] - #[[2, 1]] + #[[2, 2]]}]

      & /@ symmetrySlicingPlanes,

 {

  SurfaceColor[Hue[0]],(*lift a bit up*)

  Map[# + {\[CurlyEpsilon], \[CurlyEpsilon],

      2 \[CurlyEpsilon]} &,

    (*boundary of the unit cone*)

    {Polygon[{{0, 0, 0}, {0, 0, 1}, {1, 0, 1}}],

    Polygon[{{0, 0, 0}, {0, 0, 1}, {1, 1, 1}}],

    Polygon[{{0, 0, 0}, {1, 0, 1}, {1, 1, 1}}]}, {-2}]

  }

 }

],

   Boxed -> False]

]

Raster[Image[g3d2, ImageSize -> 4*72]]

(* This made the 2D image I uploaded earlier *)

planes = Join[toPlane /@ latticePoints, symmetrySlicingPlanes];

Take[planes, 4];

Off[Solve::"svars"];

lineOnPlane[Plane[p_, {dir1_, dir2_}], Plane[q_, {d1_, d2_}]] :=

  Module[

      {eqs, sol, line, var, aux, P1, P2},

    If[

      (*are the planes parallel?*)

      Length[

        DeleteCases[RowReduce[{dir1, dir2, d1, d2}], {0, 0, 0}, {1}]

            ] == 2,

           {},

           (*calculate direction of the intersecting line*)

           eqs = Thread[p + s dir1 + t dir2 == q + u d1 + v d2];

           Print["eqs = ", eqs];

           (*BIG CHANGE IN SOLVE ON NEXT LINE!*)
           (*Now solving for {s,t,u}, NOT {s, t, u, v}, as in original*)

           sol = Solve[eqs, {s, t, u}];

           Print["sol = ", sol];

           aux = ((p + s dir1 + t dir2) /. sol[[1]]);

           Print["aux = ", aux];

          (*two points on the line*)

          {P1, P2} = {aux /. {u -> 0, v -> 0}, aux /. {u -> 1, v -> 1}};

          Line[P1, P2 - P1]

  ]

] 

lineOnPlane[planes[[1]], planes[[5]]]

(* Prints:  eqs = {-2==-(3/2)-u-6 v, s==-(1/2)+3 u-2 v, t==-1+10 v} *)

(* Prints:  sol = {{s->1-20 v, t->-1+10 v, u->1/2 (1-12 v)}} *)

(* Prints:  aux = {-2,1-20 v,-1+10 v} *)

(* Line[{-2, 1, -1}, {0, -20, 10}] *)

I cannot promise that these changes are correct--you must check them. You should remove the three Print[] functions; I used them for the diagnostic work I was doing, but they are not necessary.

New errors have probably been introduced in the course of all of the copy and paste work necessary to move my Mathematica notebook code to this answer window. This is not an efficient way to answer a question with so much intricate code. Would have been better to e-mail you my .nb, but I think this site is not intended for exchanging .nb files.

I hope this helps you somewhat though.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you for pointing out the typo, it's just that copying onto this website accidentally dropped the bracket. Could you elaborate on how to fix it? Sorry i am a beginner... $\endgroup$
    – Jenny
    Jun 18, 2014 at 16:16
  • $\begingroup$ By the way the figure looked right. $\endgroup$
    – Jenny
    Jun 18, 2014 at 16:18
  • $\begingroup$ Did the code above work for you? $\endgroup$
    – user15996
    Jun 20, 2014 at 2:02
  • $\begingroup$ ,oh my god, you are awesome! The code worked! I really appreciate your help! Thank you so much!!! $\endgroup$
    – Jenny
    Jun 20, 2014 at 13:26
3
$\begingroup$

This problem is now easier to resolve in the current version. I will also take the opportunity to simplify some of the code used above.

Thus:

latticePointListSC[d_] := With[{n = Ceiling[d]}, 
       Select[DeleteCases[Flatten[CoordinateBoundsArray[ConstantArray[{-n, n}, 3]], 
                                  2], {0, 0, 0}], #.# <= d^2 &]]

With[{maxDist = 4}, 
     Graphics3D[Cuboid[# + {1, 1, 1}/10, # - {1, 1, 1}/10] & /@
                (latticePoints = latticePointListSC[maxDist]), Axes -> True]]

cuboids on lattice points

toPlane[p_?VectorQ] := Module[{c, dir, pos, x, y, z},
  {x, y, z} = p;
  If[FreeQ[p, 0], (* normal case *)
     dir = {#, Cross[#, p]} &[{y, -x, 0}],
     c = Count[p, 0];
     If[c == 2 (* on an axis *),
        dir = Pick[IdentityMatrix[3], p, 0],
        (* on a coordinate plane *)
        pos = Position[p, 0];
        dir = {#, Cross[#, p]} &[Insert[Cross[Delete[p, pos]], 0, pos]]]];
  InfinitePlane[p/2, dir]]

symmetrySlicingPlanes = InfinitePlane[{0, 0, 0}, #] & /@
{{{1, 0, 0}, {0, 1, 0}}, {{1, 0, 0}, {0, 0, 1}}, {{0, 1, 0}, {0, 0, 1}},
 {{0, 0, 1}, {-1, 1, 0}}, {{0, 0, 1}, {1, 1, 0}}, {{1, 0, 0}, {0, -1, 1}},
 {{1, 0, 0}, {0, 1, 1}}, {{0, 1, 0}, {-1, 0, 1}}, {{0, 1, 0}, {1, 0, 1}}};

With[{ε = 0.05},
     Graphics3D[{symmetrySlicingPlanes, {Red, 
                 Map[# + {ε, ε, 2 ε} &,
                     {Polygon[{{0, 0, 0}, {0, 0, 1}, {1, 0, 1}}],
                      Polygon[{{0, 0, 0}, {0, 0, 1}, {1, 1, 1}}], 
                      Polygon[{{0, 0, 0}, {1, 0, 1}, {1, 1, 1}}]}, {-2}]}}, 
                Boxed -> False, Method -> {"RelieveDPZFighting" -> True},
                PlotRange -> 1]]

cutting planes

planes = Join[toPlane /@ latticePoints, symmetrySlicingPlanes];

RegionIntersection[planes[[1]], planes[[5]]]
   InfiniteLine[{-2, 1, -1}, {0, -20, 10}]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.