Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
2  
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. –  Yves Klett Jun 18 at 5:38

1 Answer 1

up vote 5 down vote accepted

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.

share|improve this answer
    
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... –  Jenny Jun 18 at 16:16
    
By the way the figure looked right. –  Jenny Jun 18 at 16:18
    
Did the code above work for you? –  user15996 Jun 20 at 2:02
    
,oh my god, you are awesome! The code worked! I really appreciate your help! Thank you so much!!! –  Jenny Jun 20 at 13:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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