# User guide: indexing reflections¶

New in version 1.3.0.

The topic of indexing reflections is a large one. The indexing functionality of
`crystals`

is intended to be general enough to be useable with x-ray,
electron, and neutron diffraction data. As such, all indexing routines presented
here expect that you have a list of peaks, already oriented in reciprocal space.
This way, the details of the scattering experiments are not important to `crystals`

.

## General purpose direct indexing with DirAx¶

The first and simplest indexing routine, `index_dirax()`

, is based on the DirAx algorithm presented in:

For this example, we will need a list of reflections oriented in reciprocal space. Let’s start with a very simple example using the simplest crystal structure, cubic polonium:

```
>>> from crystals import Crystal
>>> import numpy as np
>>>
>>> indices = [ (0,0,0), (1,0,0), (0,1,0), (0,0,1) ]
>>>
>>> cryst = Crystal.from_database('Pu-epsilon') # cubic polonium
>>> peaks = [cryst.scattering_vector(hkl) for hkl in indices]
```

The list of `peaks`

are peak positions in three-dimensional reciprocal space. To index:

```
>>> from crystals import index_dirax
>>> lattice, hkls = index_dirax(peaks)
>>> lattice
< Lattice object with parameters 3.638Å, 3.638Å, 3.638Å, 90.00°, 90.00°, 90.00° >
>>> hkls.astype(int)
array([[0, 0, 0],
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
```

Here, `hkls`

are the indices of the reflections in `peaks`

according to the indexing lattice, `lattice`

.
As expected, the Miller indices in `hkls`

match the ones we created in `indices`

.

Indexing a cubic lattice with perfect peak placement is pretty easy: we only needed four reflections.
The DirAx algorithm is robust against missing reflections, but also *alien* reflections that do not belong there.
In the next example, we will introduce 10% of non-fitting reflections:

```
>>> from crystals import Crystal, index_dirax
>>> import numpy as np
>>>
>>> cryst = Crystal.from_database('Pu-epsilon') # cubic polonium
>>> indices = list(cryst.bounded_reflections(3))
>>> num_aliens = len(indices) // 10 # 10% alien reflections
>>> aliens = [ np.random.random(size=(3,)) for _ in range(num_aliens) ]
>>>
>>> peaks = [cryst.scattering_vector(hkl) for hkl in indices + aliens]
>>> lattice, hkls = index_dirax(peaks)
>>> lattice
< Lattice object with parameters 3.638Å, 3.638Å, 3.638Å, 90.00°, 90.00°, 90.00° >
>>> hkls.round(decimals=3)
array([[-1. , -1. , -1. ],
[-1. , -1. , 0. ],
[-1. , -1. , 1. ],
[-1. , 0. , -1. ],
[-1. , 0. , 0. ],
[-1. , -0. , 1. ],
[-1. , 1. , -1. ],
[-1. , 1. , 0. ],
[-1. , 1. , 1. ],
[ 0. , -1. , -1. ],
[ 0. , -1. , 0. ],
[ 0. , -1. , 1. ],
[ 0. , 0. , -1. ],
[ 0. , 0. , 0. ],
[-0. , -0. , 1. ],
[-0. , 1. , -1. ],
[-0. , 1. , 0. ],
[-0. , 1. , 1. ],
[ 1. , -1. , -1. ],
[ 1. , -1. , 0. ],
[ 1. , -1. , 1. ],
[ 1. , 0. , -1. ],
[ 1. , 0. , 0. ],
[ 1. , -0. , 1. ],
[ 1. , 1. , -1. ],
[ 1. , 1. , 0. ],
[ 1. , 1. , 1. ],
[ 0.549, 0.715, 0.603],
[ 0.545, 0.424, 0.646]])
```

As you can see, the last indexed reflections in `hkls`

don’t have integer
Miller indices; those are the non-fitting (alien) reflections!

Of course, the DirAx algorithm is also resistant to noise. Here is an example where we add noise to our perfect peaks:

```
>>> from crystals import Crystal, Lattice, index_dirax
>>> import numpy as np
>>>
>>> cryst = Crystal.from_database("Pu-epsilon")
>>> indices = list(cryst.bounded_reflections(2))
>>>
>>> peaks = [
... cryst.scattering_vector(hkl) + np.random.normal(0, scale=0.01, size=(3,))
... for hkl in indices
... ]
>>> lattice, hkls = index_dirax(peaks)
>>> hkls.round(decimals=2)
array([[-0.98, -0. , -0. ],
[ 0. , -1. , 0.01],
[ 0. , -0. , -1. ],
[ 0. , 0.01, -0. ],
[ 0. , -0. , 0.99],
[ 0. , 1.01, -0. ],
[ 1. , -0.01, 0. ]])
```

### Dealing with a restricted set of reflections¶

`index_dirax()`

can also deal with a restricted set of reflections. For example, in
electron diffraction measurements, it might occur that the only reflections measured are
of the form \(\{ (hk0) \}\). This happens, for example, when measuring diffraction from layered compounds
like graphite.

Let’s try to index **without** prior knowledge first:

```
>>> from crystals import Crystal, index_dirax
>>> import numpy as np
>>> from itertools import product
>>>
>>> cryst = Crystal.from_database('C') # graphite
>>> indices = product(
... [-2,-1,0,1,2], #-3 < h < 3
... [-2,-1,0,1,2], #-3 < k < 3
... [0] # l = 0
... )
>>> peaks = [cryst.scattering_vector(hkl) for hkl in indices]
>>> lattice, hkls = index_dirax(peaks)
Traceback (most recent call last):
crystals.indexing.common.IndexingError: No candidate lattice vectors could be determined.
```

As you can see, there not enough information to index. Let’s now use an initial guess:

```
>>> lattice, hkls = index_dirax(peaks, initial=cryst) # Lattice or Crystal works here
>>> lattice
< Lattice object with parameters 2.464Å, 2.464Å, 6.711Å, 90.00°, 90.00°, 120.00° >
```

Great!

## Other indexing methods¶

More specialized and performant indexing methods are planned to be included in `crystals`

.
If there’s a particular method you would like to see included, please do not hesitate to
raise an issue!