# Set up lattice
vector_a = [7.166, 0.0, 0.0]*Angstrom
vector_b = [0.0, 7.166, 0.0]*Angstrom
vector_c = [0.0, 0.0, 7.166]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Silicon, Silicon, Oxygen, Oxygen, Oxygen, Oxygen, Silicon, Silicon,
            Oxygen, Oxygen, Oxygen, Oxygen, Silicon, Silicon, Oxygen, Oxygen,
            Oxygen, Oxygen, Silicon, Silicon, Oxygen, Oxygen, Oxygen, Oxygen]

# Define coordinates
fractional_coordinates = [[ 0.125,  0.125,  0.125],
                          [ 0.875,  0.875,  0.875],
                          [ 0.   ,  0.   ,  0.   ],
                          [ 0.25 ,  0.   ,  0.25 ],
                          [ 0.   ,  0.25 ,  0.25 ],
                          [ 0.25 ,  0.25 ,  0.   ],
                          [ 0.625,  0.625,  0.125],
                          [ 0.375,  0.375,  0.875],
                          [ 0.5  ,  0.5  ,  0.   ],
                          [ 0.75 ,  0.5  ,  0.25 ],
                          [ 0.5  ,  0.75 ,  0.25 ],
                          [ 0.75 ,  0.75 ,  0.   ],
                          [ 0.625,  0.125,  0.625],
                          [ 0.375,  0.875,  0.375],
                          [ 0.5  ,  0.   ,  0.5  ],
                          [ 0.75 ,  0.   ,  0.75 ],
                          [ 0.5  ,  0.25 ,  0.75 ],
                          [ 0.75 ,  0.25 ,  0.5  ],
                          [ 0.125,  0.625,  0.625],
                          [ 0.875,  0.375,  0.375],
                          [ 0.   ,  0.5  ,  0.5  ],
                          [ 0.25 ,  0.5  ,  0.75 ],
                          [ 0.   ,  0.75 ,  0.75 ],
                          [ 0.25 ,  0.75 ,  0.5  ]]

# Set up configuration
sio2_cristobalite = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    ).repeat(3, 3, 3)

# Set up the q-values at which the diffraction intensity should be plotted.
n_points = 100
q_values = numpy.linspace(1.0, 6.0, n_points) / Ang

xray_intensities = DiffractionPeak.diffractionPlotData(sio2_cristobalite, q_values, cutoff_radius=12.0*Ang)
plot_model = Plot.PlotModel(Ang**-1)
line = Plot.Line(q_values, xray_intensities)
plot_model.addItem(line)
plot_model.setLimits()

plot_model.xAxis().setLabel('q')
plot_model.yAxis().setLabel('Intensity')

Plot.show(plot_model)