It is convenient to index a crystal using a 'surface coordinate' system. In this system we chose the basis vectors a_s and b_s to be within the surface plane, and the c_s vector is defined perpendicular to the surface plane (parallel to the surface normal).
A non-rational termination is one in which there is no lattice point that falls on the line perpendicular to a surface. In that case the unit cell cannot be simply repeated in the surface normal direction - rather the indexing needs to include an additional vector that defines the lateral shift of 'slabs' defined with respect to the surface plane
The below examples show using tdl for indexing some non-rational surfaces. (you can follow the same examples for a rational surface)
Also see tdl/Modules/XtalCalcs for more info on the lattice transforms

Olivine (101)

The surface_vectors function determine the possible set of in-plane vectors that can be used to define the surface mesh for the olivine (101) plane:

from xtal import Lattice, LatticeTransform
from xtal.surf_indexing import surface_vectors, calc_surf_cell

#  This is the cell defined for pnma setting
cell = [10.223, 5.992, 4.761,90,90,90]
lat = Lattice(*cell)

# calculate inplane vectors to choose Va and Vb
hkl = [1.,0.,1.]
Vs,Vr = surface_vectors(hkl,lat)

Vs is the list of real space lattice vectors that are contained within the surface plane. If you list the values of Vs they will shown [x,y,z,length], sorted by length.
You need to chose two of these vectors to form your surface mesh (ie a_s and b_s vectors). You should chose the two smallest non-collinear vectors that will form a right handed system (surface normal defining 'up').
```
Va = [1.,0.,-1.]
Vb = [0.,1.,0.]
```
Now we make Vc parallel to the surface normal. In this case we chose its length to be five times the d-spacing
```
# calculate vector d in bulk real space basis
d = lat.dvec(hkl)
Vc = 5*d
```

Now use Lattice Transform to get the surface unit cell

trans = LatticeTransform(lat,Va=Va,Vb=Vb,Vc=Vc)
slat = trans.plat()

For a non-rational termination we need to chose a repeat vector. Look at the output of Vr, this shows all the lattice vectors terminating at planes below the surface, their lengths, which plane they terminate on and the angle with respect to the surface normal. The format is [x,y,z,length,n,angle]
Note that these vectors are sorted by angle - you'll see there is not a zero degree angle, meaning we have a non-rational termination. You'll also see that the smallest angle vector (about 3degrees) terminates on the 5th plane below the surface. This is why we used Vc = 5*d above.
```
Vrpt = [-1.,0.,-4.]

# find repeat vector in surface basis
Vrpt_s = trans.xp(Vrpt)
del_1 = -1*Vrpt_s[0]
del_2 = -1*Vrpt_s[1]

# check repeat angle 
repeat_angle = slat.angle(Vrpt_s, [0.,0.,-1.])
```
The del values can be used directly to figure out how far off integer L's the bulk bragg peaks will be in the surface setting: ie Ls = n -del_1*Hs -del_2*Ks

Note that once the tansform is defined we can compute the indicies of reflections defined in the original indexing in the surface indexing and visa versa. For example:

# given the (101) in the bulk indexing, compute it in the surface basis (ie the primed basis)
h = [1,0,1]
hs = trans.hp(h)  # this should be the (005)

# take the (001) in the surface indexing and compute the bulk
hs = [0,0,10]
h = trans.h(hs)   # this should be the (202)

Hematite (1-12)

The surface_vectors function determine the possible set of in-plane vectors that can be used to define the surface mesh for the hematite (1-12) plane:

from xtal import Lattice, LatticeTransform
from xtal.surf_indexing import surface_vectors, calc_surf_cell

#  This is bulk hexagonal cell
cell = [5.035, 5.035, 13.747, 90, 90, 120.]
lat = Lattice(*cell)

# work on the (1-12) surface
hkl = [1.,-1,2.]

# calculate inplane vectors to choose Va and Vb
Vs,Vr = surface_vectors(hkl,lat)

Vs is the list of real space lattice vectors that are contained within the surface plane. If you list the values of Vs they will shown [x,y,z,length], sorted by length.
You need to chose two of these vectors to form your surface mesh (ie a_s and b_s vectors). You should chose the two smallest non-collinear vectors that will form a right handed system (surface normal defining 'up').
```
Va = [1., 1., 0.]
Vb = [-1.,1.,1.]
```
Note for hematite the following is also a lattice vector (ie in the primitive rhombohedral setting). Since its shorter its a better chioce for Vb (note they are parallel).
```
Vb = [-1/3., 1/3., 1/3.]
```

Now we make Vc parallel to the surface normal. In this case we chose its length to be twice the d-spacing

# calculate d spacing in Angstroms
d_space = lat.d(hkl)

# calculate vector d in bulk real space basis and define Vc
d = lat.dvec(hkl)
Vc = 2*d

Now use Lattice Transform to get the surface unit cell

trans = LatticeTransform(lat,Va=Va,Vb=Vb,Vc=Vc)
slat = trans.plat()

If you look at the output of Vr you'll see lots of large repeat vectors - again there is a simpler repeat vector in the rhombohedral setting

#Vrpt is a lattice vector that terminates two planes below the surface (note chose the same plane as used to define Vc
Vrpt = [-2/3., 2/3., -1/3.]
Vrpt_s = trans.xp(Vrpt)
del_1 = -1*Vrpt_s[0]
del_2 = -1*Vrpt_s[1]

# check repeat angle 
repeat_angle = slat.angle(Vrpt_s, [0.,0.,-1.])

The del values can be used directly to figure out how far off integer L's the bulk bragg peaks will be in the surface setting: ie Ls = n -del_1*Hs -del_2*Ks

Note should expand on the hematite example by first transforming to the primitive rhombohedral then finding the surface cell... This would be a nice addition to the routine - if you pass it the cell centering then the primitive unit cell could always be calculated prior to searching for the in-plane vectors.