I am writing a code to generate the Wigner-Seitz cell of the reciprocal lattice for a given set of lattice translation vectors. For example, consider the Body Centered Cubic (BCC) grid, whose base translation vectors are given by

```
a1 = {-1, 1, 1} / 2;
a2 = {1, -1, 1} / 2;
a3 = {1, 1, -1} / 2;
```

The reciprocal basis vectors are then defined according to

```
d = 2 Pi;
v = a1. (a2 [Cross]a3);
b1 = d / v (a2 [Cross]a3);
b2 = d / v (a3 [Cross]a1);
b3 = d / v (a1 [Cross]a2);
```

The reciprocal lattice is then defined by the set of reciprocal lattice vectors, the set of all linear combinations of integer multiples of reciprocal base vectors, i.

$$ vec {G} = n_1 vec {b} _1 + n_2 vec {b} _2 + n_3 vec {b} _3, qquad n_i in mathbb {Z} $$

The Wigner-Seitz cell (in this case the first Brillouin zone) is defined as the region containing the origin delimited by the vertical bisector levels of the reciprocal lattice vectors. We can generally accomplish this by considering only the first, second, and perhaps third, closest mutual grid points to the origin. For example, in the case of BCC, the following vectors suffice:

```
recipvecs =
Choose[Flatten[
Table[n1 b1 + n2 b2 + n3 b3, {n1, -1, 1}, {n2, -1, 1}, {n3, -1, 1}], 2].
standard[#] <= 2 d &];
```

**Question:** How can I construct the Wigner-Seitz cell with these vectors?

For example, one possibility is to construct the equations for all levels

```
plan = ({x, y, z} - (# / 2)). # == 0 & / @ reciplattice
```

(Note that there is a redundancy for the origin that only exists `True`

can be removed). Now it is about rewriting each of these equations as inequalities, so that the half space defined by the inequality contains the origin. I do not think that would be too difficult, but not every one of the equations can be solved for one of the coordinates, eg. we can not solve every equation for $ z $, to like

```
To solve[#, z] & / @ Planes
```

Some of the equations have to be solved $ x $ or $ y $ before they are transformed into inequalities. I think I could find a brute force solution, but I hope there is something more elegant.

Finally, I want to get the inequalities that define the region so that I can imagine it `RegionPlot3D`

and use it `Choose`

Points from a network.

**To update:**

So far I have done a simple one `If`

construct to solve the equations for both $ z $. $ x $, or $ y $.

```
sols = flattening[
If[Sz=release[Sz=Solve[sz=Löse[sz=Solve[#, z]; sz! = {}, sz,
If[sx = Solve[#, x]; sx! = {}, sx, solve[#, y]]]& / @ Airplanes]
```

This prints out the following list of rules

```
{z -> -1, z -> -1 - y, z -> -1 + x, y -> -1, x -> 1 + y, x -> 1, z -> -1 - x, z -> -1 + y, x -> -1 - y, x -> 1 - y, z -> 1 + y, z -> 1 - x, x -> -1, x -> -1 + y, y -> 1, z -> 1 + x, z -> 1 - y, z -> 1}
```

Now I just have to convert each of these into $ LHS <RHS $ and $ LHS> RHS $ and test, which contains the origin, but I do not yet know how to achieve that.