Repository: ornlneutronimaging/braggedgemodeling

Issue:
clean up and add doc

Discussion:

Patch:
diff --git a/.travis.yml b/.travis.yml
index 569a751..b238c14 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -24,6 +24,9 @@ before_install:
   - chmod +x miniconda.sh
   - ./miniconda.sh -b -p /home/travis/mc
   - export PATH=/home/travis/mc/bin:$PATH
+  - "export DISPLAY=:99.0" # configure a headless display to test plot generation
+  - "sh -e /etc/init.d/xvfb start"
+  - sleep 3 # give xvfb some time to start
 
 install:
   - export GIT_FULL_HASH=`git rev-parse HEAD`
diff --git a/README.md b/README.md
index 6425bab..83a3cc0 100644
--- a/README.md
+++ b/README.md
@@ -2,14 +2,43 @@
 [![Coverage Status](https://coveralls.io/repos/github/ornlneutronimaging/braggedgemodeling/badge.svg?branch=master)](https://coveralls.io/github/ornlneutronimaging/braggedgemodeling?branch=master)
 # Bragg Edge Modeling
 
+This python package provides tools to model and help analyze neutron Bragg Edge imaging data.
+
+**Main functionality:** given lattice structure of a material and optionally a texture model and
+an instrument beam model,
+calculate neutron Bragg Edge spectrum as a function of neutron wavelength.
+
+## Features
 * Cross section calculation from crystal structure specification
-  * Calculate diffraction peaks data (d-spacing, form factors, etc.) according to sample crystal structure
-  * Estimate inelastic scattering using incoherent approximation
-* Support of texture model: March Dollase
-* Peak profile modeling: Jorgensen model
-* Flexible design to allow future extension to the texture and peak profile models
-* Allow easy fitting to measured Bragg Edge data
+  - Calculate diffraction peaks data (d-spacing, form factors, etc.) according to sample crystal structure
+  - Estimate inelastic scattering using incoherent approximation
+* Modeling of texture:
+  - March Dollase
+* Modeling of peak profile:
+  - Jorgensen model
+* Flexible design to allow future extension to texture and peak profile models
+* Enable easy fitting to measured Bragg Edge data
+
+## Documentation
+
+Please refer to https://ornlneutronimaging.github.io/braggedgemodeling for documentation
+for installation, usage, and API.
+
+## Community guidelines
+
+**How to contribute**
+
+Please clone the repository, make changes and make a pull request.
+
+**How to report issues**
+
+Please use [the github issues](https://github.com/ornlneutronimaging/braggedgemodeling/issues) to report issues or bug reports.
+
+**Support**
+
+Please either use [the github issues](https://github.com/ornlneutronimaging/braggedgemodeling/issues) to ask for support, or contact the authors directly using email.
 
 
-Known problems:
-* Debye temperature is a table. It is missing data for some elements.
+## Known problems
+* Debye temperatures are listed in a table, which is missing data for some elements.
+  However, users can provide their own table in a configuration file.
diff --git a/bem/__init__.py b/bem/__init__.py
index 51f935d..0061cfb 100644
--- a/bem/__init__.py
+++ b/bem/__init__.py
@@ -1,6 +1,8 @@
 # -*- Python -*-
 
 """
+Bragg edge modeling
+
 deps:
   * diffpy.Structure
   * periodictable
diff --git a/bem/matter.py b/bem/matter.py
index 0709a4a..ac48555 100644
--- a/bem/matter.py
+++ b/bem/matter.py
@@ -9,7 +9,20 @@ else:
     
 
 def Structure(*args, **kwds):
-    """a wrapper for Structure method that injects "sg" data member"""
+    """a wrapper for Structure method that injects "sg" data member
+
+    Structure(atoms, lattice, sgid=)
+    
+    Parameters
+    ----------
+    atoms : list
+        list of atoms
+
+    lattice : Lattice
+    
+    sgid : int
+        space group id. For example 225 for FCC
+    """
     if 'sgid' in kwds:
         sgid = kwds.pop('sgid')
     else:
@@ -21,6 +34,8 @@ def Structure(*args, **kwds):
 
 
 def loadCif(path):
+    """load CIF file from given path to create a Structure instance
+    """
     if sys.version_info[0] < 3:
         from diffpy.Structure.Parsers import getParser
     else:
diff --git a/bem/peak_profile.py b/bem/peak_profile.py
index 8affcbc..74124f0 100644
--- a/bem/peak_profile.py
+++ b/bem/peak_profile.py
@@ -30,7 +30,23 @@ class DeltaFunction(AbstractPeakProfile):
 
 class Jorgensen(AbstractPeakProfile):
 
+    """Jorgensen peak profile
+    """
+
     def __init__(self, alpha=None, beta=None, sigma=None):
+        """Jorgensen peak profile constructor
+        
+        Parameters
+        ----------
+        alpha :
+          2-tuple of floats
+
+        beta :
+          2-tuple of floats
+
+        sigma :
+          3-tuple of floats        
+        """
         self.alpha = alpha if alpha is not None else [1., 0.]
         self.beta = beta if beta is not None else [1., 0.]
         self.sigma = sigma if sigma is not None else [0., 1., 0.]
diff --git a/bem/xscalc.py b/bem/xscalc.py
index 30de03a..be1bc16 100644
--- a/bem/xscalc.py
+++ b/bem/xscalc.py
@@ -8,24 +8,34 @@ from .xtaloriprobmodel import IsotropicXOPM
 class XSCalculator:
 
     """Cross section calculator
+    """
 
-    xc = XSCalculator(structure, T, xopm=, size=, max_diffraction_index=)
-    xc.xs(wavelen)                         # compute total cross section
-    xc.diffpeaks                           # show diffraction peaks
-    xc.xopm(peak=, wavelen=)               # compute factor due to orientation distribution (texture)
-    xc.extinction_factor(wavelen=, peak=)  # compute factor due to dynamical diffraction
+    example_code = """
+    >>> xc = XSCalculator(structure, T, xopm=, size=, max_diffraction_index=)
+    >>> xc.xs(wavelen)                         # compute total cross section
+    >>> xc.diffpeaks                           # show diffraction peaks
+    >>> xc.xopm(peak=, wavelen=)               # compute factor due to orientation distribution (texture)
     """
 
     def __init__(self, structure, T, xopm=None, max_diffraction_index=5, size=0):
         """constructor
-  Required args
-    - structure: lattice structure
-    - T: temperature
-
-  Optional args
-    - max_diffraction_index
-    - xopm: xtal orientation probability model
-    - size: size of crystallites along beam (for extinction effect calculation)
+
+        Parameters
+        ----------
+        structure : diffpy.structure
+            lattice structure
+
+        T : float
+            temperature (Kelvin)
+
+        max_diffraction_index : int, optional
+            max integer for h,k, and l
+
+        xopm : xtaloriprobmodel.XtalOriProbModel, optional
+            crystal orientation probability model
+
+        size : int, optional
+            size of crystallites along beam (for extinction effect calculation)
         """
         self.name = structure.title
         occs = np.array([atom.occupancy for atom in structure])
@@ -52,7 +62,10 @@ class XSCalculator:
     def xs(self, wavelen):
         """calculate total cross section in barns
 
-    - wavelen: a single float or a list of floats. unit: angstom
+        Parameters
+        ----------
+        wavelen : float, floats
+            wavelength or wavelengths. unit: angstrom
         """
         abs = self.xs_abs(wavelen)
         coh_el = self.xs_coh_el(wavelen)
@@ -61,18 +74,46 @@ class XSCalculator:
         return abs+coh_el+inc_el+inel
 
     def xs_inel(self, wavelen):
+        """calculate inelastic cross section in barns
+
+        Parameters
+        ----------
+        wavelen : float, floats
+            wavelength or wavelengths. unit: angstrom
+        """
         S = self._S_inel(wavelen)
         return (self.coh_xs + self.inc_xs)*S
 
     def xs_coh_inel(self, wavelen):
+        """calculate coherent inelastic cross section in barns
+
+        Parameters
+        ----------
+        wavelen : float, floats
+            wavelength or wavelengths. unit: angstrom
+        """
         S = self._S_inel(wavelen)
         return self.coh_xs*S
 
     def xs_inc_inel(self, wavelen):
+        """calculate incoherent inelastic cross section in barns
+
+        Parameters
+        ----------
+        wavelen : float, floats
+            wavelength or wavelengths. unit: angstrom
+        """
         S = self._S_inel(wavelen)
         return self.inc_xs*S
 
     def xs_inc_el(self, wavelen):
+        """calculate incoherent elastic cross section in barns
+
+        Parameters
+        ----------
+        wavelen : float, floats
+            wavelength or wavelengths. unit: angstrom
+        """
         sctts = self.sctts
         Sarr = np.array([sc.S_el_inc(wavelen, self.T)*sc.occupancy for sc in sctts])
         if len(Sarr.shape) == 1:
@@ -82,6 +123,13 @@ class XSCalculator:
         return self.inc_xs * S
 
     def xs_coh_el(self, wavelen):
+        """calculate coherent elastic cross section in barns
+
+        Parameters
+        ----------
+        wavelen : float, floats
+            wavelength or wavelengths. unit: angstrom
+        """
         vs = [np.abs(p.F)**2*p.d*p.mult*self.xopm(p, wavelen)*self.extinction_factor(wavelen, p)*(p.d*2>wavelen) for p in self.diffpeaks]
         vs = np.array(vs) # unit fm^2
         # print wavelen, vs * wavelen*wavelen/(2*self.uc_vol)
@@ -90,17 +138,36 @@ class XSCalculator:
         return np.sum(vs, axis=0)/100 * wavelen*wavelen/(2*self.uc_vol) # unit: barn
 
     def xs_coh_el__peak(self, wavelen, peak):
+        """calculate coherent elastic cross section for one peak in barns
+
+        Parameters
+        ----------
+        wavelen : float, floats
+            wavelength or wavelengths. unit: angstrom
+
+        peak : diffraction.DiffrPeak
+            diffraction peak instance
+        """
         v = np.abs(peak.F)**2*peak.d*peak.mult*self.xopm(peak, wavelen)*self.extinction_factor(wavelen, peak)*(peak.d*2>wavelen)
         return v/100 * wavelen*wavelen/(2*self.uc_vol) # unit: barn
 
     def xs_abs(self, wavelen):
+        """calculate absorption cross section in barns
+
+        Parameters
+        ----------
+        wavelen : float, floats
+            wavelength or wavelengths. unit: angstrom
+        """
         Q = 2*pi/wavelen
         from .constants import K2V
         v = K2V*Q
         return self.abs_xs_at2200/v*2200
 
     def extinction_factor(self, wavelen, pk):
-        "compute extinction factor for given wavelength and diffraction peak"
+        """
+        compute extinction factor for given wavelength and diffraction peak
+        """
         size = self.size
         if size == 0:
             return 1.
@@ -121,4 +188,31 @@ class XSCalculator:
         return S
 
 
+    def plotAll(self, wavelen):
+        """plot all cross sections w.r.t the given wavelength array
+
+        Parameters
+        ----------
+        wavelen : floats
+            wavelengths. unit: angstrom
+        """
+        coh_el_xs = self.xs_coh_el(wavelen)
+        inc_el_xs = self.xs_inc_el(wavelen)
+        abs_xs = self.xs_abs(wavelen)
+        coh_inel_xs = self.xs_coh_inel(wavelen)
+        inc_inel_xs = self.xs_inc_inel(wavelen)
+        # and the total cross section
+        total = self.xs(wavelen)
+        # plot
+        from matplotlib import pyplot as plt
+        plt.plot(wavelen, coh_el_xs, label='coh el')
+        plt.plot(wavelen, inc_el_xs, label='inc el')
+        plt.plot(wavelen, coh_inel_xs, label='coh inel')
+        plt.plot(wavelen, inc_inel_xs, label='inc inel')
+        plt.plot(wavelen, abs_xs, label='abs')
+        plt.plot(wavelen, total, label='total')
+        plt.legend()
+        return
+
+
 # End of file
diff --git a/bem/xtaloriprobmodel.py b/bem/xtaloriprobmodel.py
index 109d3da..2e1eb4e 100644
--- a/bem/xtaloriprobmodel.py
+++ b/bem/xtaloriprobmodel.py
@@ -18,8 +18,32 @@ class IsotropicXOPM(XtalOriProbModel):
 
 class MarchDollase(XtalOriProbModel):
 
+    """March Dollase model for texture
+
+    Usage
+    -----
+
+    Create an instance
+
+    >>> md = MarchDollase()
+    
+    Now md.r and md.beta are parameters.
+    The default value for r is 1.
+    The default value for beta is 0.
+    
+    To change r and beta value for a specific hkl, do
+
+    >>> md.r[hkl] = new_r_value
+    >>> md.beta[hkl] = new_beta_value
+    """
+
     class Texture(dict):
 
+        def __init__(self, defaultvalue):
+            self.defaultvalue = defaultvalue
+            dict.__init__(self)
+            return
+
         def __call__(self, hkl):
             if hkl in self: return self[hkl]
             for hkl1 in self:
@@ -28,12 +52,12 @@ class MarchDollase(XtalOriProbModel):
                 if np.allclose(hkl, np.array(hkl1)*div):
                     return self[hkl1]
                 continue
-            return 1.
+            return self.defaultvalue
     
 
     def __init__(self, r=None, beta=None):
-        self.r = r or self.Texture()
-        self.beta = beta or self.Texture()
+        self.r = r or self.Texture(1.)
+        self.beta = beta or self.Texture(0.)
         return
 
 
diff --git a/docs/Makefile b/docs/Makefile
new file mode 100644
index 0000000..893168b
--- /dev/null
+++ b/docs/Makefile
@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = sphinx-build
+SPHINXPROJ    = bem
+SOURCEDIR     = .
+BUILDDIR      = _build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
diff --git a/docs/api.rst b/docs/api.rst
new file mode 100644
index 0000000..febaeb3
--- /dev/null
+++ b/docs/api.rst
@@ -0,0 +1,30 @@
+API
+===
+
+The braggedgemodeling package provides tools to calculate neutron Bragg Edge
+spectrum for a material.
+
+Atomic structure
+----------------
+
+.. autofunction:: bem.matter.Structure
+.. autofunction:: bem.matter.loadCif
+
+Cross section calculator
+------------------------
+
+.. autoclass:: bem.xscalc.XSCalculator
+   :members:
+   :special-members: __init__
+		  
+Texture
+-------
+
+.. autoclass:: bem.xtaloriprobmodel.MarchDollase
+
+
+Peak profile
+------------
+.. autoclass:: bem.peak_profile.Jorgensen
+   :members:
+   :special-members: __init__
diff --git a/docs/conf.py b/docs/conf.py
new file mode 100644
index 0000000..c3b3717
--- /dev/null
+++ b/docs/conf.py
@@ -0,0 +1,163 @@
+# -*- coding: utf-8 -*-
+#
+# Configuration file for the Sphinx documentation builder.
+#
+# This file does only contain a selection of the most common options. For a
+# full list see the documentation:
+# http://www.sphinx-doc.org/en/master/config
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+# import os
+# import sys
+# sys.path.insert(0, os.path.abspath('.'))
+
+
+# -- Project information -----------------------------------------------------
+
+project = 'bem'
+copyright = '2017-2018, ORNL Neutron Data Sciences Group'
+author = 'ORNL Neutron Data Sciences Group'
+
+# The short X.Y version
+version = ''
+# The full version, including alpha/beta/rc tags
+release = '0.1.1'
+
+
+# -- General configuration ---------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#
+# needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+    'sphinx.ext.autodoc',
+    'sphinx.ext.mathjax',
+    'sphinx.ext.napoleon',
+    'sphinx.ext.viewcode',
+    'sphinx.ext.githubpages',
+]
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+#
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# The master toctree document.
+master_doc = 'index'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#
+# This is also used if you do content translation via gettext catalogs.
+# Usually you set "language" from the command line for these cases.
+language = None
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path .
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_theme = 'sphinx_rtd_theme'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further.  For a list of options available for each theme, see the
+# documentation.
+#
+# html_theme_options = {}
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# Custom sidebar templates, must be a dictionary that maps document names
+# to template names.
+#
+# The default sidebars (for documents that don't match any pattern) are
+# defined by theme itself.  Builtin themes are using these templates by
+# default: ``['localtoc.html', 'relations.html', 'sourcelink.html',
+# 'searchbox.html']``.
+#
+# html_sidebars = {}
+
+
+# -- Options for HTMLHelp output ---------------------------------------------
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'bemdoc'
+
+
+# -- Options for LaTeX output ------------------------------------------------
+
+latex_elements = {
+    # The paper size ('letterpaper' or 'a4paper').
+    #
+    # 'papersize': 'letterpaper',
+
+    # The font size ('10pt', '11pt' or '12pt').
+    #
+    # 'pointsize': '10pt',
+
+    # Additional stuff for the LaTeX preamble.
+    #
+    # 'preamble': '',
+
+    # Latex figure (float) alignment
+    #
+    # 'figure_align': 'htbp',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title,
+#  author, documentclass [howto, manual, or own class]).
+latex_documents = [
+    (master_doc, 'bem.tex', 'bem Documentation',
+     'ORNL Neutron Data Sciences Group', 'manual'),
+]
+
+
+# -- Options for manual page output ------------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+    (master_doc, 'bem', 'bem Documentation',
+     [author], 1)
+]
+
+
+# -- Options for Texinfo output ----------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+#  dir menu entry, description, category)
+texinfo_documents = [
+    (master_doc, 'bem', 'bem Documentation',
+     author, 'bem', 'One line description of project.',
+     'Miscellaneous'),
+]
+
+
+# -- Extension configuration -------------------------------------------------
diff --git a/docs/Bragg_edge_elastic_coherent_scattering_calculation.xlsx b/docs/excel/Bragg_edge_elastic_coherent_scattering_calculation.xlsx
similarity index 100%
rename from docs/Bragg_edge_elastic_coherent_scattering_calculation.xlsx
rename to docs/excel/Bragg_edge_elastic_coherent_scattering_calculation.xlsx
diff --git a/docs/Structure_factor_calculation_with_Debye_Waller_factor_Ni.xlsx b/docs/excel/Structure_factor_calculation_with_Debye_Waller_factor_Ni.xlsx
similarity index 100%
rename from docs/Structure_factor_calculation_with_Debye_Waller_factor_Ni.xlsx
rename to docs/excel/Structure_factor_calculation_with_Debye_Waller_factor_Ni.xlsx
diff --git a/docs/TDS_and_Absorption_scattering_calculations_Ni.xlsx b/docs/excel/TDS_and_Absorption_scattering_calculations_Ni.xlsx
similarity index 100%
rename from docs/TDS_and_Absorption_scattering_calculations_Ni.xlsx
rename to docs/excel/TDS_and_Absorption_scattering_calculations_Ni.xlsx
diff --git a/docs/index.rst b/docs/index.rst
new file mode 100644
index 0000000..b08f760
--- /dev/null
+++ b/docs/index.rst
@@ -0,0 +1,29 @@
+.. bem documentation master file, created by
+   sphinx-quickstart on Tue Aug 21 21:37:10 2018.
+   You can adapt this file completely to your liking, but it should at least
+   contain the root `toctree` directive.
+
+Welcome to braggedgemodeling!
+=============================
+
+.. toctree::
+   :maxdepth: 1
+   :caption: Contents:
+
+   installation
+   tutorial
+   api
+
+
+Links
+=====
+
+* `source repo <http://github.com/ornlneutronimaging/braggedgemodeling>`_
+
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
diff --git a/docs/installation.rst b/docs/installation.rst
new file mode 100644
index 0000000..4bdfad0
--- /dev/null
+++ b/docs/installation.rst
@@ -0,0 +1,31 @@
+.. _installation:
+
+Installation
+============
+
+The braggedgemodeling package depends on python version 3 (3.5 and newer).
+
+Using conda
+-----------
+
+The braggedgemodeling package can be installed using conda::
+
+      $ conda config --add channels conda-forge
+      $ conda install braggedgemodeling
+
+Information on dependencies of this code can be found at `the conda recipe <https://github.com/conda-forge/braggedgemodeling-feedstock/blob/master/recipe/meta.yaml>`_.
+
+      
+Using setup.py
+--------------
+
+This package is a pure python package and can be installed simply by 
+
+* Checking out the code::
+
+  $ git clone https://github.com/ornlneutronimaging/braggedgemodeling.git
+
+* Run setup.py::
+
+  $ cd braggedgemodeling && python setup.py install
+
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
new file mode 100644
index 0000000..e7a25e8
--- /dev/null
+++ b/docs/tutorial.rst
@@ -0,0 +1,95 @@
+.. _tutorial:
+
+Tutorial
+========
+
+Please see :ref:`installation` before start here.
+
+Atomic structure
+----------------
+
+Let us start by creating an atomic structure::
+
+  from bem.matter import Atom, Lattice, Structure
+  atoms = [Atom('Ni', (0,0,0)), Atom('Ni', (0.5, 0.5, 0)),
+           Atom('Ni', (0.5,0,0.5)), Atom('Ni', (0, 0.5, 0.5))]
+  a=3.5238
+  alpha = 90.
+  lattice = Lattice(a=a, b=a, c=a, alpha=alpha, beta=alpha, gamma=alpha)
+  fccNi = Structure(atoms, lattice, sgid=225)
+
+.. note:: You can also use :code:`bem.matter.loadCif(path)` to load an atomic structure
+   from a CIF file.
+
+Cross section
+-------------
+
+Then we can perform a simple Bragg Edge neutron cross section calculation and plot them::
+
+  # define wavelength axis
+  import numpy as np
+  wavelengths = np.arange(0.05, 5.5, 0.005)
+  T = 300
+  # create calculator
+  from bem import xscalc
+  xscalculator = xscalc.XSCalculator(fccNi, T, max_diffraction_index=7)
+  # compute various contributions
+  # In neutron Bragg Edge data analysis, it may not be necessary to calculate all these
+  # contributions, but it is useful to see them when exploring.
+  coh_el_xs = xscalculator.xs_coh_el(wavelengths)
+  inc_el_xs = xscalculator.xs_inc_el(wavelengths)
+  abs_xs = xscalculator.xs_abs(wavelengths)
+  coh_inel_xs = xscalculator.xs_coh_inel(wavelengths)
+  inc_inel_xs = xscalculator.xs_inc_inel(wavelengths)
+  # and the total cross section
+  total = xscalculator.xs(wavelengths)
+  # plot
+  from matplotlib import pyplot as plt
+  plt.plot(wavelengths, coh_el_xs, label='coh el')
+  plt.plot(wavelengths, inc_el_xs, label='inc el')
+  plt.plot(wavelengths, coh_inel_xs, label='coh inel')
+  plt.plot(wavelengths, inc_inel_xs, label='inc inel')
+  plt.plot(wavelengths, abs_xs, label='abs')
+  plt.plot(wavelengths, total, label='total')
+  plt.ylim(-0.2, None)
+  plt.xlim(0,7)
+  plt.legend()
+  plt.show()
+
+
+Texture
+-------
+
+To introduce texture into the sample, we can use a texture model::
+
+  from bem import xtaloriprobmodel as xopm
+  texture_model = xopm.MarchDollase()
+  texture_model.r[(0,0,1)] = 2
+
+Now we recreate the calculator using this texture model::
+  
+  xscalculator = xscalc.XSCalculator(fccNi, T, texture_model)
+
+And replot::
+    
+  xscalculator.plotAll(wavelengths)
+  plt.show()
+
+The "plotAll" method simplifies plotting.
+
+
+Peak profile
+------------
+
+To take instrument broadening into account::
+  
+  from bem import peak_profile as pp, calc
+  jorgensen = pp.Jorgensen(alpha=[50, 0.], beta=[10, 0], sigma=[0, .003, 0])
+  spectrum_calculator = calc.BraggEdgeSpectrumCalculator(xscalculator, jorgensen)
+  # calculate total cross section convolved with peak profile
+  spectrum = spectrum_calculator('total', wavelengths)
+  # plot it
+  plt.plot(wavelengths, spectrum)
+  # also plot the cross sections
+  xscalculator.plotAll(wavelengths)
+  plt.show()
diff --git a/docs/tutorial_Ni.py b/docs/tutorial_Ni.py
new file mode 100644
index 0000000..7fd7168
--- /dev/null
+++ b/docs/tutorial_Ni.py
@@ -0,0 +1,55 @@
+from bem.matter import Atom, Lattice, Structure
+atoms = [Atom('Ni', (0,0,0)), Atom('Ni', (0.5, 0.5, 0)),
+         Atom('Ni', (0.5,0,0.5)), Atom('Ni', (0, 0.5, 0.5))]
+a=3.5238
+alpha = 90.
+lattice = Lattice(a=a, b=a, c=a, alpha=alpha, beta=alpha, gamma=alpha)
+fccNi = Structure(atoms, lattice, sgid=225)
+
+# define wavelength axis
+import numpy as np
+wavelengths = np.arange(0.05, 5.5, 0.005)
+T = 300
+# create calculator
+from bem import xscalc
+xscalculator = xscalc.XSCalculator(fccNi, T, max_diffraction_index=7)
+# compute various contributions
+# In neutron Bragg Edge data analysis, it may not be necessary to calculate all these
+# contributions, but it is useful to see them when exploring.
+coh_el_xs = xscalculator.xs_coh_el(wavelengths)
+inc_el_xs = xscalculator.xs_inc_el(wavelengths)
+abs_xs = xscalculator.xs_abs(wavelengths)
+coh_inel_xs = xscalculator.xs_coh_inel(wavelengths)
+inc_inel_xs = xscalculator.xs_inc_inel(wavelengths)
+# and the total cross section
+total = xscalculator.xs(wavelengths)
+# plot
+from matplotlib import pyplot as plt
+plt.plot(wavelengths, coh_el_xs, label='coh el')
+plt.plot(wavelengths, inc_el_xs, label='inc el')
+plt.plot(wavelengths, coh_inel_xs, label='coh inel')
+plt.plot(wavelengths, inc_inel_xs, label='inc inel')
+plt.plot(wavelengths, abs_xs, label='abs')
+plt.plot(wavelengths, total, label='total')
+plt.ylim(-0.2, None)
+plt.xlim(0,7)
+plt.legend()
+plt.show()
+
+# texture
+from bem import xtaloriprobmodel as xopm
+texture_model = xopm.MarchDollase()
+texture_model.r[(0,0,1)] = 2
+xscalculator = xscalc.XSCalculator(fccNi, T, texture_model)
+xscalculator.plotAll(wavelengths)
+plt.show()
+
+
+# peak profile
+from bem import peak_profile as pp, calc
+jorgensen = pp.Jorgensen(alpha=[50, 0.], beta=[10, 0], sigma=[0, .003, 0])
+spectrum_calculator = calc.BraggEdgeSpectrumCalculator(xscalculator, jorgensen)  
+spectrum = spectrum_calculator('total', wavelengths)
+plt.plot(wavelengths, spectrum)
+xscalculator.plotAll(wavelengths)
+plt.show()
diff --git a/setup.py b/setup.py
index 6e879ac..b01a70b 100644
--- a/setup.py
+++ b/setup.py
@@ -17,6 +17,8 @@ setup(
     data_files = [],
     test_suite = 'tests',
     install_requires = [
+        'pyyaml', 'numpy', 'scipy', 'matplotlib',
+        'diffpy.structure', 'periodictable'
     ],
     dependency_links = [
     ],


Test patch:
diff --git a/tests/test_texture.py b/tests/test_texture.py
index eed96ef..a9d0135 100755
--- a/tests/test_texture.py
+++ b/tests/test_texture.py
@@ -12,6 +12,7 @@ def test_fccNi():
     T = 300
     texture_model = xopm.MarchDollase()
     texture_model.r[(0,0,1)] = 2
+    texture_model.beta[(0,0,1)] = 1.
     calc = xscalc.XSCalculator(fccNi, T, texture_model)
     coh_el_xs = calc.xs_coh_el(lambdas)
     # coh_el_xs = [calc.xs_coh_el(l) for l in lambdas]
diff --git a/tests/test_xs_Ni.py b/tests/test_xs_Ni.py
index dda67d8..e1efc60 100755
--- a/tests/test_xs_Ni.py
+++ b/tests/test_xs_Ni.py
@@ -39,21 +39,12 @@ def test_fccNi():
     expected = np.load(os.path.join(thisdir, 'expected', 'fccNi-coh-el-xs.npy'))
     assert np.isclose(data, expected).all()
 
-    inc_el_xs = calc.xs_inc_el(lambdas)
-    inel_xs = calc.xs_inel(lambdas)
-    abs_xs = calc.xs_abs(lambdas)
-    coh_inel_xs = calc.xs_coh_inel(lambdas)
-    inc_inel_xs = calc.xs_inc_inel(lambdas)
+    if not interactive:
+        import matplotlib as mpl
+        mpl.use('Agg')
+    from matplotlib import pyplot as plt        
+    calc.plotAll(lambdas)
     if interactive:
-        from matplotlib import pyplot as plt
-        plt.plot(lambdas, coh_el_xs, label='coh el')
-        plt.plot(lambdas, inc_el_xs, label='inc el')
-        plt.plot(lambdas, coh_inel_xs, label='inc inel')
-        plt.plot(lambdas, inc_inel_xs, label='inc inel')
-        plt.plot(lambdas, inel_xs, label='inel')
-        plt.plot(lambdas, abs_xs, label='abs')
-        plt.plot(lambdas, abs_xs+coh_el_xs+inc_el_xs+inel_xs)
-        plt.legend()
         plt.show()
     return