Repository: brmather/pycurious (language: python)

Issue:
Update paper.md
Completes author affiliations, and adds comma after "e.g."
[REVIEW] Notebook is raising an error
On [Ex5-Mapping-Curie-depth-EMAG2.ipynb](https://github.com/brmather/pycurious/blob/master/pycurious/Examples/Notebooks/Bouligand/Ex5-Mapping-Curie-depth-EMAG2.ipynb) an error is raised when running the following cell:
```python
grid = pycurious.CurieOptimise(mag_grid, xmin, xmax, ymin, ymax)

window_size = 400e3

# centroid spacing of 50 km
xc_list, yc_list = grid.create_centroid_list(window_size, spacingX=50e3, spacingY=50e3)
print("number of centroids = {}".format(len(xc_list)))


beta, zt, dz, C = grid.optimise_routine(window_size, xc_list, yc_list)


>>> ValueError: node spacing should be identical (1000.5265929436546, 1000.5885815185403)
```

Seems that when creating the grid, the nodes are not equally spaced between two axis.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 7166ecd..d86c9c1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,7 @@ pycurious-src/pycurious/Examples/Notebooks/Bouligand/.ipynb_checkpoints/*
 pycurious-src/pycurious/Examples/Notebooks/Tanaka/.ipynb_checkpoints/*
 pycurious-src/pycurious/Examples/.ipynb_checkpoints/*
 *.pyc
+.eggs
+build
+*.egg-info
+*.so
diff --git a/.travis.yml b/.travis.yml
index ddaa2c7..e0bf433 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -21,10 +21,6 @@ matrix:
       osx_image: xcode10.2
       language: generic
 
-before_install:
-  - if [ "$TRAVIS_OS_NAME" = "osx" ]; then brew update ; fi
-  - if [ "$TRAVIS_OS_NAME" = "osx" ]; then brew install python3 ; fi
-
 install:
   - python3 -m pip install Cython numpy scipy pytest git+https://github.com/lmoresi/pdoc.git
   - python3 -m pip install -e .
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index ded1122..61e1dcf 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -11,7 +11,7 @@ been working on something similar already.
 
 ## How to create a Pull Request (PR)
 
-Create your own fork of the project on github: log into your github account, navigate to the [`pycurious` repository](https://github.com/brmather/pycurious/tree/dev)
+Create your own fork of the project on GitHub: log into your GitHub account, navigate to the [`pycurious` repository](https://github.com/brmather/pycurious/tree/dev)
 and click on the *fork* button at the top right corner of this repository. You can find detailed instructions and a discussion of how to fork a repository, clone it locally and work on the changes
 [in the GitHub guides](https://guides.github.com/activities/forking/). 
 
@@ -21,6 +21,14 @@ It is much easier to merge in small changes to the code than extensive ones that
 
 It is a good idea to commit frequently and with commit messages that refer to the changes that have been make and why they have been made. The commit message is usually browsed without seeing the actual changes themselves and which sections of the code have been altered so a helpful message provides relevant details. 
 
+## We want your help!
+
+We welcome all contributions to `pycurious`, but if you would like to assist us on our development milestones, then we would appreciate all pull requests pertaining to the following areas:
+
+- Implementing different taper functions in addition to built-in numpy functions (e.g. `numpy.hanning`, `numpy.hanning`)
+- Enhancing I/O capabilities (e.g. exporting to netCDF, HDF5 formats)
+- Expanding the Jupyter notebook examples for specific use cases and study areas.
+
 ## Tests
 
 We use `pytest` for the unit testing framework in `pycurious`. In the source directory this means running:
diff --git a/MANIFEST.in b/MANIFEST.in
index 3d647bf..4cbb4df 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -5,7 +5,7 @@ include COPYING
 include COPYING.LESSER
 include setup.cfg
 include src/*
-include __init__.py
-recursive-include pycurious/*.py
-recursive-include pycurious/*.ipynb
+recursive-include pycurious/ *.py
+recursive-include pycurious/ *.ipynb
+include pycurious/Examples/data/test_mag_data.txt
 include tests/*.py
diff --git a/README.md b/README.md
index 4792bee..da7e01e 100644
--- a/README.md
+++ b/README.md
@@ -1,13 +1,13 @@
-# PyCurious
+![PyCurious](https://github.com/brmather/pycurious/blob/master/pycurious/Examples/Images/pycurious-logo.png?raw=true)
 
 [![Docker Cloud Automated build](https://img.shields.io/docker/cloud/automated/brmather/pycurious.svg)](https://hub.docker.com/r/brmather/pycurious)
 [![PyPI](https://img.shields.io/pypi/v/pycurious.svg)](https://pypi.org/project/pycurious/)
 [![DOI](https://zenodo.org/badge/123281222.svg)](https://zenodo.org/badge/latestdoi/123281222)
 [![Build Status](https://travis-ci.org/brmather/pycurious.svg?branch=master)](https://travis-ci.org/brmather/pycurious)
 
-Magnetic data is one of the most common geophysics datasets available on the surface of the Earth. Curie depth is the depth at which rocks lose their magnetism. The most prevalent magnetic mineral is magnetite, which has a Curie point of 580C, thus the Curie depth is often interpreted as the 580C isotherm.
+Magnetic data is one of the most common geophysics datasets available on the surface of the Earth. Curie depth is the depth at which rocks lose their magnetism. The most prevalent magnetic mineral is magnetite, which has a Curie point of 580°C, thus the Curie depth is often interpreted as the 580°C isotherm.
 
-Current methods to derive Curie depth first compute the (fast) Fourier transform over a square window of a magnetic anomaly that has been reduced to the pole. The depth and thickness of magnetic sources is estimated from the slope of the radial power spectrum. `pycurious` implements the Tanaka *et al.* (1999) and Bouligand *et al.* (2009) methods for computing the thickness of a buried magnetic source, which are covered within Jupyter notebooks.
+Current methods to derive Curie depth first compute the (fast) Fourier transform over a square window of a magnetic anomaly that has been reduced to the pole. The depth and thickness of magnetic sources is estimated from the slope of the radial power spectrum. `pycurious` implements the Tanaka *et al.* (1999) and Bouligand *et al.* (2009) methods for computing the thickness of a buried magnetic source. `pycurious` ingests maps of the magnetic anomaly and distributes the computation of Curie depth across multiple CPUs. Common computational workflows and geospatial manipulation of magnetic data are covered in the Jupyter notebooks bundled with this package.
 
 #### Binder
 
@@ -15,6 +15,12 @@ Launch the demonstration at [mybinder.org](https://mybinder.org/v2/gh/brmather/p
 
 [![badge](https://img.shields.io/badge/launch-pycurious-E66581.svg?logo=data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFkAAABZCAMAAABi1XidAAAB8lBMVEX///9XmsrmZYH1olJXmsr1olJXmsrmZYH1olJXmsr1olJXmsrmZYH1olL1olJXmsr1olJXmsrmZYH1olL1olJXmsrmZYH1olJXmsr1olL1olJXmsrmZYH1olL1olJXmsrmZYH1olL1olL0nFf1olJXmsrmZYH1olJXmsq8dZb1olJXmsrmZYH1olJXmspXmspXmsr1olL1olJXmsrmZYH1olJXmsr1olL1olJXmsrmZYH1olL1olLeaIVXmsrmZYH1olL1olL1olJXmsrmZYH1olLna31Xmsr1olJXmsr1olJXmsrmZYH1olLqoVr1olJXmsr1olJXmsrmZYH1olL1olKkfaPobXvviGabgadXmsqThKuofKHmZ4Dobnr1olJXmsr1olJXmspXmsr1olJXmsrfZ4TuhWn1olL1olJXmsqBi7X1olJXmspZmslbmMhbmsdemsVfl8ZgmsNim8Jpk8F0m7R4m7F5nLB6jbh7jbiDirOEibOGnKaMhq+PnaCVg6qWg6qegKaff6WhnpKofKGtnomxeZy3noG6dZi+n3vCcpPDcpPGn3bLb4/Mb47UbIrVa4rYoGjdaIbeaIXhoWHmZYHobXvpcHjqdHXreHLroVrsfG/uhGnuh2bwj2Hxk17yl1vzmljzm1j0nlX1olL3AJXWAAAAbXRSTlMAEBAQHx8gICAuLjAwMDw9PUBAQEpQUFBXV1hgYGBkcHBwcXl8gICAgoiIkJCQlJicnJ2goKCmqK+wsLC4usDAwMjP0NDQ1NbW3Nzg4ODi5+3v8PDw8/T09PX29vb39/f5+fr7+/z8/Pz9/v7+zczCxgAABC5JREFUeAHN1ul3k0UUBvCb1CTVpmpaitAGSLSpSuKCLWpbTKNJFGlcSMAFF63iUmRccNG6gLbuxkXU66JAUef/9LSpmXnyLr3T5AO/rzl5zj137p136BISy44fKJXuGN/d19PUfYeO67Znqtf2KH33Id1psXoFdW30sPZ1sMvs2D060AHqws4FHeJojLZqnw53cmfvg+XR8mC0OEjuxrXEkX5ydeVJLVIlV0e10PXk5k7dYeHu7Cj1j+49uKg7uLU61tGLw1lq27ugQYlclHC4bgv7VQ+TAyj5Zc/UjsPvs1sd5cWryWObtvWT2EPa4rtnWW3JkpjggEpbOsPr7F7EyNewtpBIslA7p43HCsnwooXTEc3UmPmCNn5lrqTJxy6nRmcavGZVt/3Da2pD5NHvsOHJCrdc1G2r3DITpU7yic7w/7Rxnjc0kt5GC4djiv2Sz3Fb2iEZg41/ddsFDoyuYrIkmFehz0HR2thPgQqMyQYb2OtB0WxsZ3BeG3+wpRb1vzl2UYBog8FfGhttFKjtAclnZYrRo9ryG9uG/FZQU4AEg8ZE9LjGMzTmqKXPLnlWVnIlQQTvxJf8ip7VgjZjyVPrjw1te5otM7RmP7xm+sK2Gv9I8Gi++BRbEkR9EBw8zRUcKxwp73xkaLiqQb+kGduJTNHG72zcW9LoJgqQxpP3/Tj//c3yB0tqzaml05/+orHLksVO+95kX7/7qgJvnjlrfr2Ggsyx0eoy9uPzN5SPd86aXggOsEKW2Prz7du3VID3/tzs/sSRs2w7ovVHKtjrX2pd7ZMlTxAYfBAL9jiDwfLkq55Tm7ifhMlTGPyCAs7RFRhn47JnlcB9RM5T97ASuZXIcVNuUDIndpDbdsfrqsOppeXl5Y+XVKdjFCTh+zGaVuj0d9zy05PPK3QzBamxdwtTCrzyg/2Rvf2EstUjordGwa/kx9mSJLr8mLLtCW8HHGJc2R5hS219IiF6PnTusOqcMl57gm0Z8kanKMAQg0qSyuZfn7zItsbGyO9QlnxY0eCuD1XL2ys/MsrQhltE7Ug0uFOzufJFE2PxBo/YAx8XPPdDwWN0MrDRYIZF0mSMKCNHgaIVFoBbNoLJ7tEQDKxGF0kcLQimojCZopv0OkNOyWCCg9XMVAi7ARJzQdM2QUh0gmBozjc3Skg6dSBRqDGYSUOu66Zg+I2fNZs/M3/f/Grl/XnyF1Gw3VKCez0PN5IUfFLqvgUN4C0qNqYs5YhPL+aVZYDE4IpUk57oSFnJm4FyCqqOE0jhY2SMyLFoo56zyo6becOS5UVDdj7Vih0zp+tcMhwRpBeLyqtIjlJKAIZSbI8SGSF3k0pA3mR5tHuwPFoa7N7reoq2bqCsAk1HqCu5uvI1n6JuRXI+S1Mco54YmYTwcn6Aeic+kssXi8XpXC4V3t7/ADuTNKaQJdScAAAAAElFTkSuQmCC)](https://mybinder.org/v2/gh/brmather/pycurious/binder?filepath=Notebooks%2F0-StartHere.ipynb)
 
+#### Citation
+
+[![DOI](http://joss.theoj.org/papers/10.21105/joss.01544/status.svg)](https://doi.org/10.21105/joss.01544)
+
+Mather, B. and Delhaye, R. (2019). PyCurious: A Python module for computing the Curie depth from the magnetic anomaly. _Journal of Open Source Software_, 4(39), 1544, https://doi.org/10.21105/joss.01544
+
 ## Navigation / Notebooks
 
 There are two matching sets of Jupyter notebooks - one set for the [Tanaka](#Tanaka) and one for [Bouligand](#Bouligand) implementations. The Bouligand set of noteboks are a natural choice for Bayesian inference applications.
@@ -28,17 +34,17 @@ pycurious.install_documentation(path="Notebooks")
 
 ### Tanaka
 
-- [Ex1-Plot-power-spectrum.ipynb](pycurious-src/pycurious/Examples/Notebooks/Tanaka/Ex1-Plot-power-spectrum.ipynb)
-- [Ex2-Compute-Curie-depth.ipynb](pycurious-src/pycurious/Examples/Notebooks/Tanaka/Ex2-Compute-Curie-depth.ipynb)
-- [Ex3-Parameter-exploration.ipynb](pycurious-src/pycurious/Examples/Notebooks/Tanaka/Ex3-Parameter-exploration.ipynb)
+- [Ex1-Plot-power-spectrum.ipynb](pycurious/Examples/Notebooks/Tanaka/Ex1-Plot-power-spectrum.ipynb)
+- [Ex2-Compute-Curie-depth.ipynb](pycurious/Examples/Notebooks/Tanaka/Ex2-Compute-Curie-depth.ipynb)
+- [Ex3-Parameter-exploration.ipynb](pycurious/Examples/Notebooks/Tanaka/Ex3-Parameter-exploration.ipynb)
 
 ### Bouligand
 
-- [Ex1-Plot-power-spectrum.ipynb](pycurious-src/pycurious/Examples/Notebooks/Bouligand/Ex1-Plot-power-spectrum.ipynb)
-- [Ex2-Compute-Curie-depth.ipynb](pycurious-src/pycurious/Examples/Notebooks/Bouligand/Ex2-Compute-Curie-depth.ipynb)
-- [Ex3-Posing-the-inverse-problem.ipynb](pycurious-src/pycurious/Examples/Notebooks/Bouligand/Ex3-Posing-the-inverse-problem.ipynb)
-- [Ex4-Spatial-variation-of-Curie-depth.ipynb](pycurious-src/pycurious/Examples/Notebooks/Bouligand/Ex4-Spatial-variation-of-Curie-depth.ipynb)
-- [Ex5-Mapping-Curie-depth-EMAG2.ipynb](pycurious-src/pycurious/Examples/Notebooks/Bouligand/Ex5-Mapping-Curie-depth-EMAG2.ipynb)
+- [Ex1-Plot-power-spectrum.ipynb](pycurious/Examples/Notebooks/Bouligand/Ex1-Plot-power-spectrum.ipynb)
+- [Ex2-Compute-Curie-depth.ipynb](pycurious/Examples/Notebooks/Bouligand/Ex2-Compute-Curie-depth.ipynb)
+- [Ex3-Posing-the-inverse-problem.ipynb](pycurious/Examples/Notebooks/Bouligand/Ex3-Posing-the-inverse-problem.ipynb)
+- [Ex4-Spatial-variation-of-Curie-depth.ipynb](pycurious/Examples/Notebooks/Bouligand/Ex4-Spatial-variation-of-Curie-depth.ipynb)
+- [Ex5-Mapping-Curie-depth-EMAG2.ipynb](pycurious/Examples/Notebooks/Bouligand/Ex5-Mapping-Curie-depth-EMAG2.ipynb)
 
 
 ## Installation
@@ -54,6 +60,7 @@ Also, the following packages are required:
 
 __Optional dependencies__ for mapping module and running the Notebooks:
 
+- [`jupyter`](https://jupyter.org/)
 - [`matplotlib`](https://matplotlib.org/)
 - [`pyproj`](https://github.com/jswhit/pyproj)
 - [`cartopy`](https://scitools.org.uk/cartopy/docs/latest/)
@@ -69,14 +76,74 @@ python3 -m pip install pycurious
 ```
 All the dependencies will be automatically installed by `pip`.
 
+### Installing with conda
+
+You can install `pycurious` using the [conda package manager](https://conda.io).
+Its required dependencies can be easily installed with:
+
+```bash
+conda install numpy scipy cython
+```
+
+And the full set of dependencies with:
+
+```bash
+conda install numpy scipy cython matplotlib pyproj cartopy
+```
+
+Then `pycurious` can be installed with `pip`:
+
+```bash
+pip install pycurious
+```
+
+#### Conda environment
+
+Alternatively, you can create a custom
+[conda environment](https://conda.io/docs/user-guide/tasks/manage-environments.html)
+where `pycurious` can be installed along with its dependencies.
+
+Clone the repository:
+```bash
+git clone https://github.com/brmather/pycurious
+cd pycurious
+```
+
+Create the environment from the `environment.yml` file:
+```bash
+conda env create -f environment.yml
+```
+
+Activate the newly created environment:
+```bash
+conda activate pycurious
+```
+
+And install `pycurious` with `pip`:
+```bash
+pip install pycurious
+```
+
+#### Issue with gcc
+
+If the `pycurious` installation fails due to [an issue with `gcc` and
+Anaconda](https://github.com/Anaconda-Platform/anaconda-project/issues/183), you just
+need to install `gxx_linux-64` with conda:
+
+```bash
+conda install gxx_linux-64
+```
+
+And then install `pycurious` normally.
+
+
 ### Installing using Docker
 
 A more straightforward installation for `pycurious` and all of its dependencies may be deployed with [Docker](https://www.docker.com).
 To install the docker image and start the Jupyter notebook examples:
 
 ```bash
-docker pull brmather/pycurious:latest
-docker run --name pycurious -p 8888:8888 brmather/pycurious:latest
+docker run --name pycurious -p 127.0.0.1:8888:8888 brmather/pycurious:latest
 ```
 
 ## Usage
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 0000000..13a75b9
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,16 @@
+name: pycurious
+channels:
+    - default
+    - conda-forge
+dependencies:
+    - python=3.7
+    - pip
+    - numpy
+    - scipy
+    - cython
+    - pytest
+    # Optional dependencies
+    - jupyter
+    - matplotlib
+    - pyproj
+    - cartopy
diff --git a/paper/citation.bib b/paper/citation.bib
new file mode 100644
index 0000000..4d83960
--- /dev/null
+++ b/paper/citation.bib
@@ -0,0 +1,18 @@
+@article{Mather2019,
+	journal = {Journal of Open Source Software},
+	doi = {10.21105/joss.01544},
+	issn = {2475-9066},
+	number = {39},
+	publisher = {The Open Journal},
+	title = {PyCurious: A Python module for computing the Curie depth from the magnetic anomaly.},
+	url = {http://dx.doi.org/10.21105/joss.01544},
+	volume = {4},
+	author = {Mather, Ben and Delhaye, Robert},
+	pages = {1544},
+	date = {2019-07-26},
+	year = {2019},
+	month = {7},
+	day = {26},
+}
+
+
diff --git a/paper/figure.png b/paper/figure.png
new file mode 100644
index 0000000..d8b4cc1
Binary files /dev/null and b/paper/figure.png differ
diff --git a/paper/paper.bib b/paper/paper.bib
new file mode 100644
index 0000000..ce54825
--- /dev/null
+++ b/paper/paper.bib
@@ -0,0 +1,49 @@
+@article{Bouligand2009,
+author = {Bouligand, Claire and Glen, Jonathan M. G. and Blakely, Richard J.},
+doi = {10.1029/2009JB006494},
+issn = {0148-0227},
+journal = {Journal of Geophysical Research},
+month = {nov},
+number = {B11104},
+pages = {1--25},
+title = {{Mapping Curie temperature depth in the western United States with a fractal model for crustal magnetization}},
+url = {http://doi.wiley.com/10.1029/2009JB006494},
+volume = {114},
+year = {2009}
+}
+
+@article{Tanaka1999,
+author = {Tanaka, A and Okubo, Y and Matsubayashi, O},
+doi = {10.1016/S0040-1951(99)00072-4},
+issn = {00401951},
+journal = {Tectonophysics},
+month = {June},
+number = {3-4},
+pages = {461--470},
+title = {{Curie point depth based on spectrum analysis of the magnetic anomaly data in East and Southeast Asia}},
+url = {https://dx.doi.org/10.1016/S0040-1951(99)00072-4},
+volume = {306},
+year = {1999}
+}
+
+@article{Meyer2017,
+author = {Meyer, Brian and Saltus, Richard and Chulliat, Arnaud},
+journal = {National Centers for Environmental Information, NOAA.},
+doi = {10.7289/V5H70CVX},
+title = {{EMAG2: Earth Magnetic Anomaly Grid (2-arc-minute resolution) Version 3}},
+url = {https://dx.doi.org/10.7289/V5H70CVX},
+urldate = {2018-09-01},
+year = {2017}
+}
+
+@article{Mather2019,
+author = {Mather, Ben and Fullea, Javier},
+doi = {10.5194/se-10-839-2019},
+journal = {Solid Earth},
+number = {3},
+pages = {839--850},
+title = {{Constraining the geotherm beneath the British Isles from Bayesian inversion of Curie depth: integrated modelling of magnetic, geothermal, and seismic data}},
+url = {https://doi.org/10.5194/se-10-839-2019},
+volume = {10},
+year = {2019}
+}
diff --git a/paper/paper.md b/paper/paper.md
new file mode 100644
index 0000000..4cbd888
--- /dev/null
+++ b/paper/paper.md
@@ -0,0 +1,83 @@
+---
+title: 'PyCurious: A Python module for computing the Curie depth from the magnetic anomaly.'
+tags:
+  - Python
+  - Curie depth
+  - magnetism
+  - magnetic anomaly
+  - Bayesian inference
+authors:
+  - name: Ben Mather
+    orcid: 0000-0003-3566-1557
+    affiliation: "1,2" # (Multiple affiliations must be quoted)
+  - name: Robert Delhaye
+    orcid: 0000-0003-2128-4295
+    affiliation: "2"
+affiliations:
+ - name: School of Geoscience, The University of Sydney, Sydney, Australia
+   index: 1
+ - name: School of Cosmic Physics, Dublin Institute for Advanced Studies,  Dublin, Ireland
+   index: 2
+date: 30 April 2019
+bibliography: paper.bib
+---
+
+# Summary
+
+Multiple geophysical methods have been proposed to resolve the thermal structure of the Earth's lithosphere with varying degrees of precision.
+Geotherms may be constructed from heat flow or xenolith data, but the spatial coverage of these are often limited.
+Seismic velocity has proven effective to infer upper-mantle temperature, but its application relies on building a compositional model suitable for the geological context and estimating attenuation from grainsize and water content.
+In contrast, magnetic data are among the most widespread geophysical datasets available on the surface of the Earth and can offer useful insight into its thermal structure. 
+
+The magnetic anomaly is the observed magnetic field less the contribution from the core and external fields.
+Geologic features enhance or depress the local magnetic field; this response is primarily controlled by the concentration of magnetite in the Earth's crust.
+Magnetite is the most prevalent magnetic mineral, in terms of both its magnetic susceptibility and concentration in the crust, and has a Curie point of 580°C.
+The depth at which rocks lose their permanent magnetic properties is referred to as the Curie depth.
+Due to the prevalence of magnetite in the crust, the Curie depth is often interpreted as the 580°C isotherm.
+Current methods to derive the Curie depth first compute the (fast) Fourier transform over a square window of a magnetic anomaly that has been reduced to the pole.
+The depth and thickness of magnetic sources is estimated from the slope of the radial power spectrum.
+
+Our Python package, `pycurious`, ingests a map of the Earth's magnetic anomaly and distributes the computation of Curie depth across multiple CPUs.
+`pycurious` implements the @Tanaka1999 and @Bouligand2009 methods for computing the thickness of a buried magnetic source.
+The former selects portions of the radial power spectrum in the low and high frequency domain to compute the depth of magnetic sources, while the latter fits an analytical solution to the entire power spectrum.
+We cast the @Bouligand2009 method within a Bayesian framework to estimate the uncertainty of Curie depth calculations [@Mather2019].
+Common computational workflows and geospatial manipulation of magnetic data are covered in the Jupyter notebooks bundled with this package.
+The `mapping` module includes a set of functions that help to wrangle maps of the magnetic anomaly into a useful form for `pycurious`.
+Such an approach is commonly encountered for transforming global compilations of the magnetic anomaly, e.g., EMAG2 [@Meyer2017], from latitudinal/longitudinal coordinates to a local projection in eastings/northings.
+
+![Radial power spectra (right) computed from different sized windows of a synthetic magnetic anomaly (left) using `pycurious`. Larger window sizes increase the number of points in the low wavenumber range of the radial power spectrum, which is important for estimating the thickness of magnetic sources in the crust.](figure.png)
+
+`pycurious` includes the following functionality:
+
+- Importing and exporting GeoTiff files.
+- Converting between geospatial coordinate reference systems (CRS).
+- Decomposition of magnetic grids into square windows over which to compute Curie depth.
+- Parallel distribution of computation across computing resources.
+- Bayesian inversion framework and sensitivity analysis.
+
+
+## Documentation and examples
+
+`pycurious` is bundled with a linked collection of Jupyter notebooks that can act as a user guide and an introduction to the package.
+The notebooks are split into matching sets for frequentist and Bayesian estimates of Curie depth.
+The notebooks cover:
+
+- Plotting the radial and azimuthal power spectrum.
+- Computing Curie depth from a synthetic magnetic anomaly.
+- Exploring parameter sensitivity.
+- Posing the inverse problem and objective function.
+- Mapping Curie depth using the EMAG2 magnetic anomaly dataset [@Meyer2017].
+
+All documentation can be accessed from within the module via a Python function that installs the notebooks at a filesystem location specified by the user at run time.
+The API documentation is kept up-to-date on [GitHub pages](https://brmather.github.io/pycurious/).
+`pycurious` and all Python dependencies can be installed through the pypi.org `pip` package manager, however, several of the dependencies for `cartopy` may cause problems for inexperienced users.
+We therefore provide a Docker image built with all required and optional dependencies and a deployment of the documentation / examples on [mybinder.org](https://mybinder.org/v2/gh/brmather/pycurious/binder?filepath=Notebooks%2F0-StartHere.ipynb).
+
+
+# Acknowledgements
+
+This work was made possible by the G.O.THERM.3D project, supported by an Irish Research Council Research for Policy & Society grant (RfPS/2016/50) co-funded by Geological Survey Ireland and by Sustainable Energy Authority Of Ireland.
+Development of `pycurious` was inspired from [`pycpd`](https://github.com/groupeLIAMG/pycpd) and [`fatiando a terra`](https://github.com/fatiando/fatiando).
+
+
+# References
diff --git a/pycurious/Examples/0-StartHere.ipynb b/pycurious/Examples/0-StartHere.ipynb
index d1df166..888bab7 100644
--- a/pycurious/Examples/0-StartHere.ipynb
+++ b/pycurious/Examples/0-StartHere.ipynb
@@ -4,9 +4,11 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# PyCurious\n",
+    "![PyCurious](https://github.com/brmather/pycurious/blob/master/pycurious/Examples/Images/pycurious-logo.png?raw=true)\n",
     "\n",
-    "Magnetic data is one of the most common geophysics datasets available on the surface of the Earth. Curie depth is the depth at which rocks lose their magnetism. The most prevalent magnetic mineral is magnetite, which has a Curie point of 580C, thus the Curie depth is often interpreted as the 580C isotherm.\n",
+    "---\n",
+    "\n",
+    "Magnetic data is one of the most common geophysics datasets available on the surface of the Earth. Curie depth is the depth at which rocks lose their magnetism. The most prevalent magnetic mineral is magnetite, which has a Curie point of 580°C, thus the Curie depth is often interpreted as the 580°C isotherm.\n",
     "\n",
     "Current methods to derive Curie depth first compute the (fast) Fourier transform over a square window of a magnetic anomaly that has been reduced to the pole. The depth and thickness of magnetic sources is estimated from the slope of the radial power spectrum. `pycurious` implements the Tanaka *et al.* (1999) and Bouligand *et al.* (2009) methods for computing the thickness of a buried magnetic source, which are covered within Jupyter notebooks.\n",
     "\n",
@@ -14,9 +16,18 @@
     "\n",
     "Launch the demonstration at [mybinder.org](https://mybinder.org/v2/gh/brmather/pycurious/binder?filepath=Notebooks%2F0-StartHere.ipynb)\n",
     "\n",
+    "\n",
     "## Navigation / Notebooks\n",
     "\n",
-    "There are two matching sets of `pycurious` notebooks - one set for the [Tanaka](#Tanaka) and one for [Bouligand](#Bouligand) implementations. The Bouligand set of noteboks are a natural choice for Bayesian inference applications.\n",
+    "There are two matching sets of Jupyter notebooks - one set for the [Tanaka](#Tanaka) and one for [Bouligand](#Bouligand) implementations. The Bouligand set of noteboks are a natural choice for Bayesian inference applications.\n",
+    "\n",
+    "Note, these examples can be installed from the package itself by running:\n",
+    "\n",
+    "```python\n",
+    "import pycurious\n",
+    "pycurious.install_documentation(path=\"Notebooks\")\n",
+    "```\n",
+    "\n",
     "\n",
     "### Tanaka\n",
     "\n",
@@ -32,13 +43,42 @@
     "- [Ex4-Spatial-variation-of-Curie-depth.ipynb](./Notebooks/Bouligand/Ex4-Spatial-variation-of-Curie-depth.ipynb)\n",
     "- [Ex5-Mapping-Curie-depth-EMAG2.ipynb](./Notebooks/Bouligand/Ex5-Mapping-Curie-depth-EMAG2.ipynb)\n",
     "\n",
-    "### Examples\n",
+    "## Installation\n",
     "\n",
-    "Note, these examples can be installed from the package itself by running:\n",
+    "### Dependencies\n",
     "\n",
-    "```python\n",
-    "import pycurious\n",
-    "pycurious.install_documentation(path=\"Notebooks\")\n",
+    "You will need **Python 2.7 or 3.5+**.\n",
+    "Also, the following packages are required:\n",
+    "\n",
+    "- [`numpy`](http://numpy.org)\n",
+    "- [`scipy`](https://scipy.org)\n",
+    "- [`cython`](https://cython.org/)\n",
+    "\n",
+    "__Optional dependencies__ for mapping module and running the Notebooks:\n",
+    "\n",
+    "- [`matplotlib`](https://matplotlib.org/)\n",
+    "- [`pyproj`](https://github.com/jswhit/pyproj)\n",
+    "- [`cartopy`](https://scitools.org.uk/cartopy/docs/latest/)\n",
+    "\n",
+    "### Installing using pip\n",
+    "\n",
+    "You can install `pycurious` using the\n",
+    "[`pip package manager`](https://pypi.org/project/pip/) with either version of Python:\n",
+    "\n",
+    "```bash\n",
+    "python2 -m pip install pycurious\n",
+    "python3 -m pip install pycurious\n",
+    "```\n",
+    "All the dependencies will be automatically installed by `pip`.\n",
+    "\n",
+    "### Installing using Docker\n",
+    "\n",
+    "A more straightforward installation for `pycurious` and all of its dependencies may be deployed with [Docker](https://www.docker.com).\n",
+    "To install the docker image and start the Jupyter notebook examples:\n",
+    "\n",
+    "```bash\n",
+    "docker pull brmather/pycurious:latest\n",
+    "docker run --name pycurious -p 8888:8888 brmather/pycurious:latest\n",
     "```\n",
     "\n",
     "## Usage\n",
@@ -65,49 +105,32 @@
     "k, Phi, sigma_Phi = grid.radial_spectrum(subgrid)\n",
     "```\n",
     "\n",
-    "## Installation\n",
-    "\n",
-    "To install with **setuptools**:\n",
+    "A series of tests are located in the *tests* subdirectory.\n",
+    "In order to perform these tests, clone the repository and run [`pytest`](https://pypi.org/project/pytest/):\n",
     "\n",
     "```bash\n",
-    "python setup.py install --user\n",
+    "git checkout https://github.com/brmather/pycurious.git\n",
+    "cd pycurious\n",
+    "pytest -v\n",
     "```\n",
     "\n",
-    "This will compile all C source code and install them to the user directory (omit the `--user` flag to install to the system directory). Remember to delete the build folder if you are upgrading this package.\n",
-    "\n",
-    "Alternatively, install using **pip**:\n",
-    "\n",
-    "```bash\n",
-    "pip install pycurious --user\n",
-    "```\n",
+    "### API Documentation\n",
     "\n",
-    "### Dependencies\n",
+    "The API for all functions and classes in `pycurious` can be accessed from [https://brmather.github.io/pycurious/](https://brmather.github.io/pycurious/).\n",
     "\n",
-    "- Python 2.7+ and 3.6+\n",
-    "- Numpy 1.9+\n",
-    "- Scipy 0.14+\n",
-    "- Cython\n",
-    "\n",
-    "Optional dependencies for mapping module:\n",
-    "\n",
-    "- Matplotlib\n",
-    "- [pyproj](https://github.com/jswhit/pyproj)\n",
-    "- [Cartopy](http://scitools.org.uk/cartopy/docs/latest/index.html)\n",
-    "\n",
-    "### Docker\n",
-    "\n",
-    "A more straightforward installation for `pycurious` and all of its dependencies may be deployed with [Docker](https://www.docker.com). To install the docker image and test it is working:\n",
-    "\n",
-    "```bash\n",
-    "docker pull brmather/pycurious:latest\n",
-    "docker run --rm brmather/pycurious:latest help\n",
-    "```\n",
     "\n",
     "## References\n",
     "\n",
     "1. Bouligand, C., Glen, J. M. G., & Blakely, R. J. (2009). Mapping Curie temperature depth in the western United States with a fractal model for crustal magnetization. Journal of Geophysical Research, 114(B11104), 1–25. https://doi.org/10.1029/2009JB006494\n",
     "2. Tanaka, A., Okubo, Y., & Matsubayashi, O. (1999). Curie point depth based on spectrum analysis of the magnetic anomaly data in East and Southeast Asia. Tectonophysics, 306(3–4), 461–470. https://doi.org/10.1016/S0040-1951(99)00072-4\n"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
diff --git a/pycurious/Examples/Images/pycurious-logo.pdf b/pycurious/Examples/Images/pycurious-logo.pdf
new file mode 100644
index 0000000..52b5a3f
Binary files /dev/null and b/pycurious/Examples/Images/pycurious-logo.pdf differ
diff --git a/pycurious/Examples/Images/pycurious-logo.png b/pycurious/Examples/Images/pycurious-logo.png
new file mode 100644
index 0000000..2c50377
Binary files /dev/null and b/pycurious/Examples/Images/pycurious-logo.png differ
diff --git a/pycurious/Examples/Notebooks/Bouligand/Ex1-Plot-power-spectrum.ipynb b/pycurious/Examples/Notebooks/Bouligand/Ex1-Plot-power-spectrum.ipynb
index 671aea3..50ae868 100644
--- a/pycurious/Examples/Notebooks/Bouligand/Ex1-Plot-power-spectrum.ipynb
+++ b/pycurious/Examples/Notebooks/Bouligand/Ex1-Plot-power-spectrum.ipynb
@@ -237,13 +237,13 @@
    "source": [
     "## Azimuthal power spectrum\n",
     "\n",
-    "The azimuthal spectrum computes the FFT on a square window that is polarised at a range of radii. This is useful to explore linear trends in the magnetic anomaly that may align with a particular foliation or strike orientation.\n",
+    "The azimuthal spectrum computes the FFT on a square window that is polarised at a range of radii. Subdividing the transforms into bins by azimuth is useful to explore linear trends in the magnetic anomaly that may align with a particular foliation or strike orientation.\n",
     "\n",
     "```python\n",
-    "azimuthal_spectrum(subgrid, theta=5.0)\n",
+    "azimuthal_spectrum(subgrid, power=2.0, theta=5.0)\n",
     "```\n",
     "\n",
-    "`theta` controls the thickness of each azimuth (in degrees)."
+    "`theta` controls the bin size of each azimuth (in degrees). The FFT of the magnetic anomaly is raised to the `power=2` (default) which is compatible with Bouligand _et al._ (2009) computation of Curie depth."
    ]
   },
   {
@@ -282,6 +282,13 @@
     "ax1.legend(bbox_to_anchor=(1,1))\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -300,7 +307,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.6.9"
   }
  },
  "nbformat": 4,
diff --git a/pycurious/Examples/Notebooks/Bouligand/Ex2-Compute-Curie-depth.ipynb b/pycurious/Examples/Notebooks/Bouligand/Ex2-Compute-Curie-depth.ipynb
index d3386be..ca8607d 100644
--- a/pycurious/Examples/Notebooks/Bouligand/Ex2-Compute-Curie-depth.ipynb
+++ b/pycurious/Examples/Notebooks/Bouligand/Ex2-Compute-Curie-depth.ipynb
@@ -61,7 +61,15 @@
    "source": [
     "## Plot radial power spectrum\n",
     "\n",
-    "The radial power spectrum is computed from a square window of the magnetic anomaly. Methods to select window sizes and compute the Fast Fourier Transform (FFT) belong to the `CurieGrid` object. We apply the default `np.hanning` taper to the power spectrum as in [Ex1-Plot-power-spectrum](#./Ex1-Plot-power-spectrum.ipynb)."
+    "The radial power spectrum is computed from a square window of the magnetic anomaly. Methods to select window sizes and compute the Fast Fourier Transform (FFT) belong to the `CurieGrid` object. We apply the default `np.hanning` taper to the power spectrum as in [Ex1-Plot-power-spectrum](#./Ex1-Plot-power-spectrum.ipynb).\n",
+    "\n",
+    "By default the FFT of the magnetic anomaly is raised to the power 2:\n",
+    "\n",
+    "```python\n",
+    "grid.radial_spectrum(subgrid, taper=None, power=2.0, **kwargs)\n",
+    "```\n",
+    "\n",
+    "which is compatible with Bouligand _et al._ (2009) computation of Curie depth."
    ]
   },
   {
@@ -249,7 +257,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.6.9"
   }
  },
  "nbformat": 4,
diff --git a/pycurious/Examples/Notebooks/Bouligand/Ex3-Posing-the-inverse-problem.ipynb b/pycurious/Examples/Notebooks/Bouligand/Ex3-Posing-the-inverse-problem.ipynb
index b11dcca..1709620 100644
--- a/pycurious/Examples/Notebooks/Bouligand/Ex3-Posing-the-inverse-problem.ipynb
+++ b/pycurious/Examples/Notebooks/Bouligand/Ex3-Posing-the-inverse-problem.ipynb
@@ -342,7 +342,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.6.9"
   }
  },
  "nbformat": 4,
diff --git a/pycurious/Examples/Notebooks/Tanaka/Ex1-Plot-power-spectrum.ipynb b/pycurious/Examples/Notebooks/Tanaka/Ex1-Plot-power-spectrum.ipynb
index 8676c09..d1728c5 100644
--- a/pycurious/Examples/Notebooks/Tanaka/Ex1-Plot-power-spectrum.ipynb
+++ b/pycurious/Examples/Notebooks/Tanaka/Ex1-Plot-power-spectrum.ipynb
@@ -130,7 +130,7 @@
    "outputs": [],
    "source": [
     "\n",
-    "k1, Phi1, sigma_Phi2 = grid.radial_spectrum(subgrid, taper=None)\n",
+    "k1, Phi1, sigma_Phi1 = grid.radial_spectrum(subgrid, taper=None)\n",
     "k2, Phi2, sigma_Phi2 = grid.radial_spectrum(subgrid, taper=np.hanning)\n",
     "k3, Phi3, sigma_Phi2 = grid.radial_spectrum(subgrid, taper=np.hamming)\n",
     "\n",
@@ -150,13 +150,13 @@
    "source": [
     "## Azimuthal power spectrum\n",
     "\n",
-    "The azimuthal spectrum computes the FFT on a square window that is polarised at a range of radii. This is useful to explore linear trends in the magnetic anomaly that may align with a particular foliation or strike orientation.\n",
+    "The azimuthal spectrum computes the FFT on a square window that is polarised at a range of radii. Subdividing the transforms into bins by azimuth is useful to explore linear trends in the magnetic anomaly that may align with a particular foliation or strike orientation.\n",
     "\n",
     "```python\n",
-    "azimuthal_spectrum(subgrid, theta=5.0)\n",
+    "azimuthal_spectrum(subgrid, power=2.0, theta=5.0)\n",
     "```\n",
     "\n",
-    "`theta` controls the thickness of each azimuth (in degrees)."
+    "`theta` controls the bin size of each azimuth (in degrees). The FFT of the magnetic anomaly is raised to the `power=2`(default) which is compatible with Bouligand _et al._ (2009) computation of Curie depth. For Tanaka *et al.* (1999), however, we need to take the square root which is set with `power=0.5`."
    ]
   },
   {
@@ -165,7 +165,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "k, Phi, theta = grid.azimuthal_spectrum(subgrid, theta=20.0)\n",
+    "k, Phi, theta = grid.azimuthal_spectrum(subgrid, power=0.5, theta=20.0)\n",
     "\n",
     "fig = plt.figure()\n",
     "ax1 = fig.add_subplot(111, ylabel=r'azimuth $\\theta$', xlabel='wavenumber (rad/km)')\n",
@@ -195,6 +195,13 @@
     "ax1.legend(bbox_to_anchor=(1,1))\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -213,7 +220,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.6.9"
   }
  },
  "nbformat": 4,
diff --git a/pycurious/Examples/Notebooks/Tanaka/Ex2-Compute-Curie-depth.ipynb b/pycurious/Examples/Notebooks/Tanaka/Ex2-Compute-Curie-depth.ipynb
index e2df8b5..be69111 100644
--- a/pycurious/Examples/Notebooks/Tanaka/Ex2-Compute-Curie-depth.ipynb
+++ b/pycurious/Examples/Notebooks/Tanaka/Ex2-Compute-Curie-depth.ipynb
@@ -76,7 +76,15 @@
    "metadata": {},
    "source": [
     "## Plot radial power spectrum\n",
-    "The radial power spectrum is computed from a square window of the magnetic anomaly. Methods to select window sizes and compute the Fast Fourier Transform (FFT) belong to the `CurieGrid` object. We apply the default `np.hanning` taper to the power spectrum as in [Ex1-Plot-power-spectrum](#./Ex1-Plot-power-spectrum.ipynb)."
+    "The radial power spectrum is computed from a square window of the magnetic anomaly. Methods to select window sizes and compute the Fast Fourier Transform (FFT) belong to the `CurieGrid` object. We apply the default `np.hanning` taper to the power spectrum as in [Ex1-Plot-power-spectrum](#./Ex1-Plot-power-spectrum.ipynb).\n",
+    "\n",
+    "By default the FFT of the magnetic anomaly is raised to the power 2:\n",
+    "\n",
+    "```python\n",
+    "grid.radial_spectrum(subgrid, taper=None, power=2.0, **kwargs)\n",
+    "```\n",
+    "\n",
+    "which is compatible with Bouligand _et al._ (2009) computation of Curie depth. For Tanaka *et al.* (1999), however, we need to take the square root which is set with `power=0.5`."
    ]
   },
   {
@@ -103,13 +111,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "k, Phi, sigma_Phi = grid.radial_spectrum_log(subgrid, taper=None)\n",
+    "k, Phi, sigma_Phi = grid.radial_spectrum(subgrid, taper=None, power=0.5)\n",
     "\n",
     "# Plot of power spectrum, as function of spatial frequency\n",
     "plt.plot(k/(2*np.pi), Phi/(2*np.pi), '-o')\n",
     "plt.title('Power Spectrum')\n",
     "plt.xlabel(r'$|k|/(2\\pi)$ [km$^{-1}$]')\n",
-    "plt.ylabel(r'$\\ln(\\Phi_{\\Delta T}(|k|)^{1/2})/(2\\pi)$')"
+    "plt.ylabel(r'$\\ln(\\Phi_{\\Delta T}(|k|)^{1/2})/(2\\pi)$')\n",
+    "plt.show()"
    ]
   },
   {
@@ -122,7 +131,8 @@
     "plt.plot(k/(2*np.pi), Phi_n/(2*np.pi), '-o')\n",
     "plt.title('Power Spectrum, weighted by k')\n",
     "plt.xlabel(r'$|k|/(2\\pi)$ [km$^{-1}$]')\n",
-    "plt.ylabel(r'$\\ln(\\Phi_{\\Delta T}(|k|)^{1/2}/|k|)/(2\\pi)$')"
+    "plt.ylabel(r'$\\ln(\\Phi_{\\Delta T}(|k|)^{1/2}/|k|)/(2\\pi)$')\n",
+    "plt.show()"
    ]
   },
   {
@@ -167,8 +177,16 @@
     "\n",
     "ax2.plot(k/(2*np.pi), Phi_n/(2*np.pi), '-o',k/(2*np.pi),(Zor*k/(2*np.pi)+bor)/(2*np.pi),'r-')\n",
     "ax2.set_xlabel(r'$|k|/(2\\pi)$ [km$^{-1}$]')\n",
-    "ax2.set_ylabel(r'$\\ln(\\Phi_{\\Delta T}(|k|)^{1/2}/|k|)/(2\\pi)$')\n"
+    "ax2.set_ylabel(r'$\\ln(\\Phi_{\\Delta T}(|k|)^{1/2}/|k|)/(2\\pi)$')\n",
+    "plt.show()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -187,7 +205,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.6.9"
   }
  },
  "nbformat": 4,
diff --git a/pycurious/Examples/Notebooks/Tanaka/Ex3-Parameter-exploration.ipynb b/pycurious/Examples/Notebooks/Tanaka/Ex3-Parameter-exploration.ipynb
index 717a8f2..7399a49 100644
--- a/pycurious/Examples/Notebooks/Tanaka/Ex3-Parameter-exploration.ipynb
+++ b/pycurious/Examples/Notebooks/Tanaka/Ex3-Parameter-exploration.ipynb
@@ -78,7 +78,7 @@
     "# 1) Where in k-space do we locate window?\n",
     "window_size = 304e3\n",
     "subgrid = grid.subgrid(window_size, xpt, ypt)\n",
-    "k, Phi, sigma_Phi = grid.radial_spectrum_log(subgrid, taper=None)\n",
+    "k, Phi, sigma_Phi = grid.radial_spectrum(subgrid, taper=None, power=0.5)\n",
     "\n",
     "# Number of bins to divide k-range into\n",
     "nbins = 10\n",
@@ -145,7 +145,7 @@
    "source": [
     "# 2) How big a bandwidth in k-space do we want?\n",
     "\n",
-    "k, Phi, sigma_Phi = grid.radial_spectrum_log(subgrid, taper=None)\n",
+    "k, Phi, sigma_Phi = grid.radial_spectrum(subgrid, taper=None, power=0.5)\n",
     "\n",
     "# We are testing this by taking one central spatial frequency (k=1), and estimating Zb with increasing k-bandwidths.\n",
     "# Recall k=2*pi*nu\n",
@@ -210,7 +210,7 @@
     "for i in range(0, nwin):\n",
     "    ws = baseW + i*inc\n",
     "    subgrid = grid.subgrid(ws, xpt, ypt)\n",
-    "    k, Phi, sigma_Phi = grid.radial_spectrum_log(subgrid, taper=None)\n",
+    "    k, Phi, sigma_Phi = grid.radial_spectrum(subgrid, taper=None, power=0.5)\n",
     "    (Ztr,btr,dZtr), (Zor, bor, dZor) = pycurious.tanaka1999(k, Phi, sigma_Phi, (0.11, 0.21), (0.11, 0.21))\n",
     "    Zb, eZb = pycurious.ComputeTanaka(Ztr, dZtr, Zor, dZor)\n",
     "    \n",
@@ -256,7 +256,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.6.9"
   }
  },
  "nbformat": 4,
diff --git a/pycurious/__init__.py b/pycurious/__init__.py
index 51783ff..87fb5dc 100644
--- a/pycurious/__init__.py
+++ b/pycurious/__init__.py
@@ -1,7 +1,26 @@
+# Copyright 2018-2019 Ben Mather, Robert Delhaye
+#
+# This file is part of PyCurious.
+#
+# PyCurious is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any later version.
+#
+# PyCurious is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+
 """
-Copyright 2018 Ben Mather, Robert Delhaye
 
-This file is part of PyCurious.
+> PyCurious is a Python package for computing the Curie depth from the magnetic anomaly.
+
+## Licence
+
+Copyright 2018-2019 Ben Mather, Robert Delhaye
 
 PyCurious is free software: you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
@@ -14,6 +33,112 @@
 
 You should have received a copy of the GNU Lesser General Public License
 along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+
+## Installation
+
+### Dependencies
+
+You will need **Python 2.7 or 3.5+**.
+Also, the following packages are required:
+
+- [`numpy`](http://numpy.org)
+- [`scipy`](https://scipy.org)
+- [`cython`](https://cython.org/)
+
+__Optional dependencies__ for mapping module and running the Notebooks:
+
+- [`matplotlib`](https://matplotlib.org/)
+- [`pyproj`](https://github.com/jswhit/pyproj)
+- [`cartopy`](https://scitools.org.uk/cartopy/docs/latest/)
+
+### Installing using pip
+
+You can install `pycurious` using the
+[`pip package manager`](https://pypi.org/project/pip/) with either
+version of Python:
+
+>>> python2 -m pip install pycurious
+>>> python3 -m pip install pycurious
+
+All the dependencies will be automatically installed by `pip`.
+
+### Installing using Docker
+
+A more straightforward installation for `pycurious` and all of its
+dependencies may be deployed with [Docker](https://www.docker.com).
+To install the docker image and start the Jupyter notebook examples:
+
+>>> docker pull brmather/pycurious:latest
+>>> docker run --name pycurious -p 8888:8888 brmather/pycurious:latest
+
+## Documenation / Notebooks
+
+Jupyter notebooks that demonstrate the functionality of PyCurious and
+some common workflows may be installed to a local directory with
+
+```python
+import pycurious
+pycurious.install_documentation(path="Notebooks")
+```
+
+## Contributing to pycurious
+
+We welcome contributions to `pycurious`, large or small [\\*](#footnote1).
+That can be in the form of new code, improvements to the documentation, 
+helping with a missing test, or it may even just be pointing out a bug or
+potential improvement to the team. 
+
+For bugs and suggestions, the most effective way to reach the team is by
+raising an issue on the github issue tracker. Github allows you to classify
+your issues so that we know if it is a bug report, feature request or
+feedback to the authors.
+
+If you wish to contribute some changes to the code then you should submit a
+*pull request*. On GitHub we can review the code that you are contributing
+and discuss it before the changes are merged into the development version of
+the code. Before embarking on changes to the code, please first take a look
+at the [`dev` branch](https://github.com/brmather/pycurious/tree/dev) which
+is where unreleased changes are staged: we may have been working on something
+similar already. 
+
+### How to create a Pull Request (PR)
+
+Create your own fork of the project on github: log into your github account,
+navigate to the [`pycurious` repository](https://github.com/brmather/pycurious/tree/dev)
+and click on the *fork* button at the top right corner of this repository.
+You can find detailed instructions and a discussion of how to fork a
+repository, clone it locally and work on the changes
+[in the GitHub guides](https://guides.github.com/activities/forking/). 
+
+### When to make your pull request
+
+It is much easier to merge in small changes to the code than extensive ones
+that touch multiple files. If a small incremental change is possible, please
+issue a request for that change rather than saving everything up. 
+
+It is a good idea to commit frequently and with commit messages that refer
+to the changes that have been make and why they have been made. The commit
+message is usually browsed without seeing the actual changes themselves and
+which sections of the code have been altered so a helpful message provides
+relevant details. 
+
+### Tests
+
+We use [`pytest`](https://pypi.org/project/pytest/) for the unit testing
+framework in `pycurious`. In the source directory this means running:
+
+>>> python setup.py test
+
+The existing tests should be passing before you start coding (help us out with
+an issue if that is not the case !) and when you have finished. Any new
+functionality should also have tests that we can use to verify the code.
+It is important that you make it clear if the original tests have had to
+change to accomodate new code / functionality.
+
+<a name="footnote1">*</a> _Small changes are our favorites as they are much
+easier to quickly merge into the existing code. Proofreading typos, bug fixes,
+... the little stuff is especially welcome._
+
 """
 
 # -*- coding: utf-8 -*-
@@ -21,4 +146,4 @@
 from .grid import CurieGrid, bouligand2009, tanaka1999, maus1995, ComputeTanaka
 from .optimise import CurieOptimise
 from . import mapping
-from . import download
\ No newline at end of file
+from . import download
diff --git a/pycurious/documentation.py b/pycurious/documentation.py
index ff346ed..8415e62 100644
--- a/pycurious/documentation.py
+++ b/pycurious/documentation.py
@@ -1,19 +1,23 @@
-"""
-Copyright 2018 Ben Mather, Robert Delhaye
-
-This file is part of PyCurious.
-
-PyCurious is free software: you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation, either version 3 of the License, or any later version.
+# Copyright 2018-2019 Ben Mather, Robert Delhaye
+#
+# This file is part of PyCurious.
+#
+# PyCurious is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any later version.
+#
+# PyCurious is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
 
-PyCurious is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU Lesser General Public License for more details.
+"""
+Use the `pycurious.documentation.install_documentation` function to copy all
+Jupyter Notebooks and example data to a local directory.
 
-You should have received a copy of the GNU Lesser General Public License
-along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 import pkg_resources as _pkg_resources
@@ -33,19 +37,30 @@ def install_documentation(path="./PyCurious-Examples"):
 
     Some dependencies exist for the notebooks to be useful:
 
-       - matplotlib: for some diagrams
-       - cartopy: for mapping and visualisation
-       - pyepsg: for converting between map projections
+       - `matplotlib`: for some diagrams
+       - `cartopy`: for mapping and visualisation
+       - `pyepsg`: for converting between map projections
 
     PyCurious dependencies may be explicitly imported in the notebooks including:
 
-       - numpy
-       - scipy
+       - `numpy`
+       - `scipy`
 
     """
 
-    Notebooks_Path = _pkg_resources.resource_filename("pycurious", os.path.join("Examples"))
-
-    ct = _dir_util.copy_tree(Notebooks_Path, path, preserve_mode=1, preserve_times=1, preserve_symlinks=1, update=0, verbose=1, dry_run=0)
+    Notebooks_Path = _pkg_resources.resource_filename(
+        "pycurious", os.path.join("Examples")
+    )
+
+    ct = _dir_util.copy_tree(
+        Notebooks_Path,
+        path,
+        preserve_mode=1,
+        preserve_times=1,
+        preserve_symlinks=1,
+        update=0,
+        verbose=1,
+        dry_run=0,
+    )
 
     return
diff --git a/pycurious/download.py b/pycurious/download.py
index bec09d5..a860e6e 100644
--- a/pycurious/download.py
+++ b/pycurious/download.py
@@ -1,22 +1,48 @@
+# Copyright 2018-2019 Ben Mather, Robert Delhaye
+#
+# This file is part of PyCurious.
+#
+# PyCurious is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any later version.
+#
+# PyCurious is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+
 """
-Copyright 2018 Ben Mather, Robert Delhaye
+This module provides some simple functions to download or install files.
+
+It is good practise to store a checksum for each file. A checksum is a
+unique signature used to identify a file, and ensure that the data
+contained within the file is legitimate. Here is a sample workflow for
+downloading EMAG2 (v3) and evaluating its checksum:
 
-This file is part of PyCurious.
+```python
+resource = {
+     "local_file":"../../data/EMAG2_V3_20170530.npz",
+     "md5":'c0898b6a91efb3f13783873a8b67380c',
+     "url":"https://zenodo.org/record/3245551/files/EMAG2_V3_20170530.npz?download=1",
+     "expected_size":"500Mb"
+    }
 
-PyCurious is free software: you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation, either version 3 of the License, or any later version.
+download_cached_file(resource["url"], resource["local_file"], resource["md5"], resource["expected_size"])
+```
 
-PyCurious is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU Lesser General Public License for more details.
+The file will commence downloading if it does not already exist in the local directory
+or if the checksum (md5) does not match. If you do not know the checksum for a file, run
 
-You should have received a copy of the GNU Lesser General Public License
-along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+```python
+md5sum(filename)
+```
+
+to return its unique identifier.
 """
 
-## Here are some functions to download or install files. 
 
 def download_file(url, local_filename, expected_size="Unknown"):
     """
@@ -35,9 +61,9 @@ def download_file(url, local_filename, expected_size="Unknown"):
 
     # We might want to bundle some files if they are small / compressed or not readily available for download
 
-    if (os.path.isfile(url)):
+    if os.path.isfile(url):
         shutil.copy(url, local_filename)
-    
+
     else:
         r = requests.get(url, stream=True)
 
@@ -45,17 +71,23 @@ def download_file(url, local_filename, expected_size="Unknown"):
         last_time = start_time
         datasize = 0
 
-        with open(local_filename, 'wb') as f:
-            for chunk in r.iter_content(chunk_size=10000000): 
-                if chunk: # filter out keep-alive new chunks
+        with open(local_filename, "wb") as f:
+            for chunk in r.iter_content(chunk_size=10000000):
+                if chunk:  # filter out keep-alive new chunks
                     f.write(chunk)
                     f.flush()
                     datasize += len(chunk)
                     if (time.time() - last_time) > 2.5:
-                        print ("{:.1f} Mb in {:.1f}s / {}".format(datasize / 1.0e6, time.time() - start_time, expected_size))
+                        print(
+                            "{:.1f} Mb in {:.1f}s / {}".format(
+                                datasize / 1.0e6,
+                                time.time() - start_time,
+                                expected_size,
+                            )
+                        )
                         last_time = time.time()
 
-    return 
+    return
 
 
 def md5sum(filename):
@@ -65,14 +97,16 @@ def md5sum(filename):
     import hashlib
     from functools import partial
 
-    with open(filename, mode='rb') as f:
+    with open(filename, mode="rb") as f:
         d = hashlib.md5()
-        for buf in iter(partial(f.read, 4096), b''):
+        for buf in iter(partial(f.read, 4096), b""):
             d.update(buf)
     return d.hexdigest()
 
 
-def download_cached_file(location_url, local_file, expected_md5, expected_size="Unknown"):
+def download_cached_file(
+    location_url, local_file, expected_md5, expected_size="Unknown"
+):
     """
     Download files from a URL to a local filename.
 
@@ -93,21 +127,21 @@ def download_cached_file(location_url, local_file, expected_md5, expected_size="
 
     try:
         assert md5sum(local_file) == expected_md5
-        print ("Using cached file - {}".format(local_file))
-        return(2)
-    
+        print("Using cached file - {}".format(local_file))
+        return 2
+
     except (IOError, AssertionError) as error_info:
         # No file or the wrong file ... best go download it
         # print "Assertion failed - {}".format(sys.exc_info())
-        
+
         try:
             data_file = download_file(location_url, local_file, expected_size)
-            print ("Downloaded from {}".format( location_url ))
-            return(1)
+            print("Downloaded from {}".format(location_url))
+            return 1
 
         except:
-            print ("Unable to download {} [{}] ".format( location_url, sys.exc_info() ))
-            return(0)
+            print("Unable to download {} [{}] ".format(location_url, sys.exc_info()))
+            return 0
 
 
 def report_cached_file(local_file, expected_md5):
@@ -122,12 +156,15 @@ def report_cached_file(local_file, expected_md5):
     """
     import os
     import os.path
-    
+
     if not os.path.isfile(local_file):
-        print ("Local file {} does not exist".format(local_file))
+        print("Local file {} does not exist".format(local_file))
     else:
         if len(expected_md5) == 0 or md5sum(local_file) == expected_md5:
-            print ("Cached file {} is valid".format(local_file))
+            print("Cached file {} is valid".format(local_file))
         else:
-            print ("Cached file {} failed, checksum {}".format(local_file, md5sum(local_file)))
-
+            print(
+                "Cached file {} failed, checksum {}".format(
+                    local_file, md5sum(local_file)
+                )
+            )
diff --git a/pycurious/grid.py b/pycurious/grid.py
index 72bd171..3978894 100644
--- a/pycurious/grid.py
+++ b/pycurious/grid.py
@@ -1,30 +1,49 @@
+# Copyright 2018-2019 Ben Mather, Robert Delhaye
+#
+# This file is part of PyCurious.
+#
+# PyCurious is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any later version.
+#
+# PyCurious is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+
 """
-Copyright 2018 Ben Mather, Robert Delhaye
+This PyCurious module contains the `pycurious.grid.CurieGrid` class,
+which can be initialised with a magnetic grid of equal spacing in the x and y direction.
+It contains methods for the following functionality:
 
-This file is part of PyCurious.
+- Decomposition of subgrids for processing square windows of the magnetic anomaly
+- Removing linear trends from the magnetic anomaly
+- Upward continuation
+- Reduction to the pole
 
-PyCurious is free software: you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation, either version 3 of the License, or any later version.
+Other functions within this module are useful to compute analytical solutions
+of the radial power spectrum, \\( \\Phi \\) according to Bouligand *et al.* (2009),
+Maus and Dimri (1995), and the decomposition of \\( \\Phi \\) from the magnetic
+anomaly according to Tanaka *et al.* (1999):
 
-PyCurious is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU Lesser General Public License for more details.
+- `bouligand2009`: analytic solution used in `pycurious.optimise.CurieOptimise`
+- `maus1995`: simplified version of `bouligand2009` without higher order integration.
+- `tanaka1999`: to be used in conjunction with `ComputeTanaka`
 
-You should have received a copy of the GNU Lesser General Public License
-along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
 """
 
 # -*- coding: utf-8 -*-
 import numpy as np
-from scipy.special import gamma, kv, lambertw
-from scipy.optimize import minimize
-from scipy.signal import tukey
+from scipy.special import gamma, kv
 import warnings
 
-try: range = xrange
-except: pass
+try:
+    range = xrange
+except:
+    pass
 
 
 class CurieGrid(object):
@@ -32,17 +51,52 @@ class CurieGrid(object):
     Accepts a 2D array and Cartesian coordinates specifying the
     bounding box of the array
 
-    Grid must be projected in metres
+    Grid must be projected in metres.
 
     Args:
-     grid : 2D array of magnetic data
-     xmin : minimum x bound in metres
-     xmax : maximum x bound in metres
-     ymin : minimum y bound in metres
-     ymax : maximum y bound in metres
+        grid : 2D numpy array
+            2D array of magnetic data
+        xmin : float
+            minimum x bound in metres
+        xmax : float
+            maximum x bound in metres
+        ymin : float
+            minimum y bound in metres
+        ymax : float
+            maximum y bound in metres
+
+    Attributes:
+        grid : 2D numpy array
+            2D array of magnetic data
+        xmin : float
+            minimum x bound in metres
+        xmax : float
+            maximum x bound in metres
+        ymin : float
+            minimum y bound in metres
+        ymax : float
+            maximum y bound in metres
+        dx : float
+            grid spacing in the x-direction in metres
+        dy : float
+            grid spacing in the y-direction in metres
+        nx : int
+            number of nodes in the x-direction
+        ny : int
+            number of nodes in the y-direction
+        xcoords : 1D numpy array
+            1D numpy array of coordinates in the x-direction
+        ycoords : 1D numpy array
+            1D numpy array of coordinates in the y-direction
+
+    Notes:
+        In all instances `x` indicates eastings in metres and `y` indicates northings.
+        Using a grid of longitude / latitudinal coordinates (degrees) will result
+        in incorrect Curie depth calculations.
     """
+
     def __init__(self, grid, xmin, xmax, ymin, ymax):
-        
+
         self.data = np.array(grid)
         ny, nx = self.data.shape
         self.xmin, self.xmax = xmin, xmax
@@ -52,15 +106,14 @@ def __init__(self, grid, xmin, xmax, ymin, ymax):
         self.nx, self.ny = nx, ny
         self.dx, self.dy = dx, dy
 
-        if abs(dx - dy) > 1.0:
-            warnings.warn("node spacing should be identical {}".format((dx,dy)), RuntimeWarning)
-
+        if not np.allclose(dx, dy, 1.0):
+            raise ValueError("node spacing should be identical {}".format((dx, dy)))
 
     def subgrid(self, window, xc, yc):
         """
         Extract a subgrid from the data at a window around
         the point (xc,yc)
-        
+
         Args:
             xc : float
                 x coordinate
@@ -74,16 +127,16 @@ def subgrid(self, window, xc, yc):
                 subgrid encompassing window size
         """
 
-        # check in coordinate in grid
+        # check whether coordinate is inside grid
         if xc < self.xmin or xc > self.xmax or yc < self.ymin or yc > self.ymax:
-            raise ValueError("Point ({},{}) outside data range".format(xc,yc))
+            raise ValueError("Point {} outside data range".format((xc, yc)))
 
         # find nearest index to xc,yc
         ix = np.abs(self.xcoords - xc).argmin()
         iy = np.abs(self.ycoords - yc).argmin()
 
-        nw = int(round(window/self.dx))
-        n2w = nw//2
+        nw = int(round(window / self.dx))
+        n2w = nw // 2
 
         # extract a square window from the data
         imin = ix - n2w
@@ -91,17 +144,18 @@ def subgrid(self, window, xc, yc):
         jmin = iy - n2w
         jmax = iy + n2w + 1
 
-        # safeguard if window size is larger than grid
-        imin = max(imin, 0)
-        imax = min(imax, self.nx)
-        jmin = max(jmin, 0)
-        jmax = min(jmax, self.ny)
-        
+        # check whether window fits inside grid
+        if imin < 0 or imax > self.nx or jmin < 0 or jmax > self.ny:
+            raise ValueError(
+                "Window size {} at centroid {} exceeds the data range".format(
+                    window, (xc, yc)
+                )
+            )
+
         data = self.data[jmin:jmax, imin:imax]
 
         return data
 
-
     def create_centroid_list(self, window, spacingX=None, spacingY=None):
         """
         Create a list of xc,yc values to extract subgrids.
@@ -122,28 +176,25 @@ def create_centroid_list(self, window, spacingX=None, spacingY=None):
             yc_list : 1D array
                 array of y coordinates
         """
-        
-        nx, ny = self.nx, self.ny
         xcoords = self.xcoords
         ycoords = self.ycoords
-        
-        nw = int(round(window/self.dx))
-        n2w = nw//2
-        
+
+        nw = int(round(window / self.dx))
+        n2w = nw // 2
+
         # this is the densest spacing possible given the data
         xc = xcoords[n2w:-n2w]
         yc = ycoords[n2w:-n2w]
-        
+
         # but we can alter it if required
-        if type(spacingX) != type(None):
+        if spacingX is not None:
             xc = np.arange(xc.min(), xc.max(), spacingX)
-        if type(spacingY) != type(None):
+        if spacingY is not None:
             yc = np.arange(yc.min(), yc.max(), spacingY)
 
         xq, yq = np.meshgrid(xc, yc)
-            
-        return xq.ravel(), yq.ravel()
 
+        return xq.ravel(), yq.ravel()
 
     def remove_trend_linear(self, data):
         """
@@ -156,165 +207,150 @@ def remove_trend_linear(self, data):
             data : 2D numpy array
 
         Returns:
-         data : 2D numpy array
+            data : 2D numpy array
         """
-        data = data.copy()
         nr, nc = data.shape
-        xq, yq = np.meshgrid(np.arange(0, nc), np.arange(0, nr))
-
+        yq, xq = np.mgrid[0:nc, 0:nr]
         A = np.c_[xq.ravel(), yq.ravel(), np.ones(xq.size)]
         c, resid, rank, sigma = np.linalg.lstsq(A, data.ravel(), rcond=None)
+        return data - np.dot(A, c).reshape(data.shape)
 
-        data.flat[:] -= np.dot(A, c)
-        return data
-
-    def radial_spectrum(self, subgrid, taper=np.hanning, scale=0.001, **kwargs):
+    def _taper_spectrum(self, subgrid, taper=np.hanning, scale=0.001, **kwargs):
         """
-        Compute the radial spectrum for a square grid.
-
-        Args:
-            subgrid : 2D array
-                window of the original data (see subgrid method)
-            taper : function (default=np.hanning)
-                taper function, set to None for no taper function
-            scale : float (default: 0.001)
-                scaling factor to get k into rad/km
-            kwargs : keyword arguments
-                keyword arguments to pass to `taper`
+        Template for tapering the power spectrum used in:
 
-        Returns:
-            k : 1D array shape (n,)
-                wavenumber in rad/km
-            Phi : 1D array shape (n,)
-                Radial power spectrum
-            sigma_Phi : 1D array shape (n,)
-                Standard deviation of Phi
+        - `radial_spectrum`
+        - `radial_spectrum_log`
+        - `azimuthal_spectrum`
         """
-
         data = subgrid
         nr, nc = data.shape
-        nw = nr
 
         if nr != nc:
-            warnings.warn("subgrid is not square {}".format((nr,nc)), RuntimeWarning)
-
+            warnings.warn("subgrid is not square {}".format((nr, nc)), RuntimeWarning)
 
         # control taper
-        if type(taper) == type(None):
+        if taper is None:
             vtaper = 1.0
         else:
             rt = taper(nr, **kwargs)
             ct = taper(nc, **kwargs)
             xq, yq = np.meshgrid(ct, rt)
-            vtaper = xq*yq
+            vtaper = xq * yq
 
-        dx_scale = self.dx*scale
-        dk = 2.0*np.pi/(nw - 1)/dx_scale
+        # scaling factor to transform wavenumber into units of rad/km
+        dx_scale = self.dx * scale
+        dk = 2.0 * np.pi / (nr - 1) / dx_scale
 
-        # fast Fourier transform and shift
-        FT = np.abs(np.fft.fft2(data*vtaper))
-        FT = np.fft.fftshift(FT)
+        kbins = np.arange(dk, dk * nr / 2, dk)
+        return vtaper, dk, kbins
+
+    def _FFT_spectrum(self, subgrid, vtaper, dk, kbins, const):
+        """
+        Template for computing the (fast) Fourier transform used in:
+
+        - `radial_spectrum`
+        - `radial_spectrum_log`
+
+        A constant `const` should be applied to the FFT of the magnetic anomaly
+        to convert `S` and `sigma` to specific units for further analysis.
+
+        It is useful to remember that:
+
+        ```python
+        2*log(FFT) == log(FFT**2)
+        ```
+        """
+        data = subgrid
+        nr, nc = data.shape
 
-        kbins = np.arange(dk, dk*nw/2, dk)
         nbins = kbins.size - 1
 
+        # fast Fourier transform and shift
+        FT = np.abs(np.fft.fft2(data * vtaper))
+        FT = np.fft.fftshift(FT)
+
         S = np.empty(nbins)
         k = np.empty(nbins)
         sigma = np.empty(nbins)
 
-        i0 = int((nw - 1)//2)
-        ix, iy = np.mgrid[0:nw,0:nw]
-        kk = np.hypot((ix - i0)*dk, (iy - i0)*dk)
+        i0 = int((nr - 1) // 2)
+        ix, iy = np.mgrid[0:nr, 0:nr]
+        kk = np.hypot((ix - i0) * dk, (iy - i0) * dk)
 
-        for i in range(0, nbins):
-            mask = np.logical_and(kk >= kbins[i], kk <= kbins[i+1])
-            rr = 2.0*np.log(FT[mask])
+        for i in range(nbins):
+            mask = np.logical_and(kk >= kbins[i], kk <= kbins[i + 1])
+            rr = const * np.log(FT[mask])
             S[i] = rr.mean()
             k[i] = kk[mask].mean()
             sigma[i] = np.std(rr)
 
         return k, S, sigma
 
-
-    def radial_spectrum_log(self, subgrid, taper=np.hanning, scale=0.001, **kwargs):
+    def radial_spectrum(self, subgrid, taper=np.hanning, power=2.0, **kwargs):
         """
-        Compute the log of the radial spectrum for a square grid.
+        Compute the radial spectrum for a square grid.
+
+        > Wavenumber is returned in values of __rad/km__
 
         Args:
-         subgrid : 2D array
-            window of the original data (see subgrid method)
-         taper : function (default=np.hanning)
-            taper function, set to None for no taper function
-         scale : float (detault=0.001)
-            scaling factor to get k into rad/km
-         kwargs : keyword arguments
-            keyword arguments to pass to `taper`
+            subgrid : 2D array
+                window of the original data (see subgrid method)
+            taper : function (default=np.hanning)
+                taper function, set to None for no taper function
+            power : float
+                raise the FFT of the magnetic anomaly to the power.
+                - 2.0 for Bouligand _et al._ (2009) use cases
+                - 0.5 for Tanaka _et al.__ (1999) use cases
+            kwargs : keyword arguments
+                keyword arguments to pass to `taper`
 
         Returns:
             k : 1D array shape (n,)
                 wavenumber in rad/km
-            lnPhi : 1D array
-                log of the radial power spectrum in ln(sqrt(S))
-            lnsigma_Phi : 1D array
-                standard deviation of lnPhi
-        """
-
-        data = subgrid
-        nr, nc = data.shape
-        nw = nr
-
-        if nr != nc:
-            warnings.warn("subgrid is not square {}".format((nr,nc)), RuntimeWarning)
-
-
-        # control taper
-        if type(taper) == type(None):
-            vtaper = 1.0
-        else:
-            rt = taper(nr, **kwargs)
-            ct = taper(nc, **kwargs)
-            xq, yq = np.meshgrid(ct, rt)
-            vtaper = xq*yq
-
-        dx_scale = self.dx*scale
-        dk = 2.0*np.pi/(nw - 1)/dx_scale
-
-        # fast Fourier transform and shift
-        FT = np.abs(np.fft.fft2(data*vtaper))
-        FT = np.fft.fftshift(FT)
+            Phi : 1D array shape (n,)
+                Radial power spectrum
+            sigma_Phi : 1D array shape (n,)
+                Standard deviation of Phi
 
-        kbins = np.arange(dk, dk*nw/2, dk)
-        nbins = kbins.size - 1
+        Notes:
+            While `subgrid` is projected in eastings / northings (in metres),
+            the wavenumber, \\( k \\), is returned in units of rad/km.
+            This is because both Bouligand *et al.* (2009) and Tanaka *et al.*
+            (1999) require the computation of Curie depth in these units.
 
-        S = np.empty(nbins)
-        k = np.empty(nbins)
-        sigma = np.empty(nbins)
+        References:
+            Bouligand, C., J. M. G. Glen, and R. J. Blakely (2009), Mapping Curie
+            temperature depth in the western United States with a fractal model for
+            crustal magnetization, J. Geophys. Res., 114, B11104,
+            doi:10.1029/2009JB006494
+
+            Tanaka, A., Okubo, Y., & Matsubayashi, O. (1999). Curie point depth
+            based on spectrum analysis of the magnetic anomaly data in East and
+            Southeast Asia. Tectonophysics, 306(3–4), 461–470.
+            doi:10.1016/S0040-1951(99)00072-4
+        """
 
-        i0 = int((nw - 1)//2)
-        ix, iy = np.mgrid[0:nw,0:nw]
-        kk = np.hypot((ix - i0)*dk, (iy - i0)*dk)
+        # bin the spectrum and compute the taper
+        vtaper, dk, kbins = self._taper_spectrum(subgrid, taper, **kwargs)
 
-        for i in range(0, nbins):
-            mask = np.logical_and(kk >= kbins[i], kk <= kbins[i+1])
-            rr = np.log(np.sqrt(FT[mask]))
-            S[i] = rr.mean()
-            k[i] = kk[mask].mean()
-            sigma[i] = np.std(rr)
+        # calculate the Fourier transform and apply scaling constant to retrieve
+        # values compatible with Bouligand or Tanaka analysis
+        return self._FFT_spectrum(subgrid, vtaper, dk, kbins, power)
 
-        return k, S, sigma
-        
-        
-    def azimuthal_spectrum(self, subgrid, taper=np.hanning, scale=0.001, theta=5.0, **kwargs):
+    def azimuthal_spectrum(
+        self, subgrid, taper=np.hanning, power=2.0, theta=5.0, **kwargs
+    ):
         """
         Compute azimuthal spectrum for a square grid.
 
+        > Wavenumber is returned in values of __rad/km__
+
         Args:
             subgrid : 2D array
                 window of the original data (see subgrid method)
             taper : function (default=np.hanning)
                 taper function, set to None for no taper function
-            scale : float (default=0.001)
-                scaling factor to get k into rad/km
             theta : float
                 angle increment in degrees
             args : arguments
@@ -327,24 +363,31 @@ def azimuthal_spectrum(self, subgrid, taper=np.hanning, scale=0.001, theta=5.0,
                 Radial power spectrum
             sigma_Phi : 1D array shape (n,)
                 Standard deviation of Phi
-        """
-        from pycurious import radon
-        data = subgrid
-        nr, nc = data.shape
-        nw = nr
 
-        if nr != nc:
-            raise RuntimeWarning("subgrid is not square {}".format((nr,nc)))
+        Notes:
+            While `subgrid` is projected in eastings / northings (in metres),
+            the wavenumber, \\( k \\), is returned in units of rad/km.
+            This is because both Bouligand *et al.* (2009) and Tanaka *et al.*
+            (1999) require the computation of Curie depth in these units.
 
-        dx_scale = self.dx*scale
-        dk = 2.0*np.pi/(nw - 1)/dx_scale
+        References:
+            Bouligand, C., J. M. G. Glen, and R. J. Blakely (2009), Mapping Curie
+            temperature depth in the western United States with a fractal model for
+            crustal magnetization, J. Geophys. Res., 114, B11104,
+            doi:10.1029/2009JB006494
+
+            Tanaka, A., Okubo, Y., & Matsubayashi, O. (1999). Curie point depth
+            based on spectrum analysis of the magnetic anomaly data in East and
+            Southeast Asia. Tectonophysics, 306(3–4), 461–470.
+            doi:10.1016/S0040-1951(99)00072-4
+        """
+        from pycurious import radon
 
-        kbins = np.arange(dk, dk*nw/2, dk)
+        vtaper, dk, kbins = self._taper_spectrum(subgrid, taper, **kwargs)
 
         dtheta = np.arange(0.0, 180.0, theta)
-        sinogram = radon.radon2d(data, np.pi*dtheta/180.0)
+        sinogram = radon.radon2d(subgrid, np.pi * dtheta / 180.0)
         S = np.zeros((dtheta.size, kbins.size))
-        sigma2 = np.zeros((dtheta.size, kbins.size))
 
         # control taper
         if taper is None:
@@ -352,16 +395,15 @@ def azimuthal_spectrum(self, subgrid, taper=np.hanning, scale=0.001, theta=5.0,
         else:
             vtaper = taper(sinogram.shape[0], **kwargs)
 
-        nk = 1 + 2*kbins.size
+        nk = 1 + 2 * kbins.size
         for i in range(0, dtheta.size):
-            PSD = np.abs(np.fft.fft(vtaper*sinogram[:,i], n=nk))
-            S[i,:] = np.log( np.sqrt(PSD[1:kbins.size+1] ))
-        
-        return kbins, S, dtheta
+            PSD = np.abs(np.fft.fft(vtaper * sinogram[:, i], n=nk))
+            S[i, :] = power * np.log(PSD[1 : kbins.size + 1])
 
+        return kbins, S, dtheta
 
     def reduce_to_pole(self, data, inc, dec, sinc=None, sdec=None):
-        r"""
+        """
         Reduce total field magnetic anomaly data to the pole.
 
         The reduction to the pole if a phase transformation that can be
@@ -387,7 +429,7 @@ def reduce_to_pole(self, data, inc, dec, sinc=None, sdec=None):
         Returns:
             rtp : 2D array
                 the data reduced to the pole.
-        
+
         References:
             Blakely, R. J. (1996), Potential Theory in Gravity and Magnetic
             Applications, Cambridge University Press.
@@ -395,26 +437,35 @@ def reduce_to_pole(self, data, inc, dec, sinc=None, sdec=None):
         Notes:
             This functions performs the reduction in the frequency domain
             (using the FFT). The transform filter is (in the freq domain):
-                RTP(k_x, k_y) = \frac{|k|}{
-                    a_1 k_x^2 + a_2 k_y^2 + a_3 k_x k_y +
-                    i|k|(b_1 k_x + b_2 k_y)}
-            in which k_x and k_y are the wave-numbers in the x and y
+
+            \\( RTP(k_x, k_y) = \\frac{|k|}{
+                a_1 k_x^2 + a_2 k_y^2 + a_3 k_x k_y +
+                i|k|(b_1 k_x + b_2 k_y)}    \\)
+
+            in which \\( k_x, k_y \\) are the wave-numbers in the x and y
             directions and
-                |k| = \sqrt{k_x^2 + k_y^2} \\
-                a_1 = m_z f_z - m_x f_x \\
-                a_2 = m_z f_z - m_y f_y \\
-                a_3 = -m_y f_x - m_x f_y \\
-                b_1 = m_x f_z + m_z f_x \\
-                b_2 = m_y f_z + m_z f_y
-            \mathbf{m} = (m_x, m_y, m_z)` is the unit-vector of the total
+
+            \\( |k| = \\sqrt{k_x^2 + k_y^2} \\)
+
+            \\( a_1 = m_z f_z - m_x f_x     \\)
+
+            \\( a_2 = m_z f_z - m_y f_y     \\)
+
+            \\( a_3 = -m_y f_x - m_x f_y    \\)
+
+            \\( b_1 = m_x f_z + m_z f_x     \\)
+
+            \\( b_2 = m_y f_z + m_z f_y     \\)
+
+            \\( \\mathbf{m} = (m_x, m_y, m_z) \\) is the unit-vector of the total
             magnetization of the source and
-            \mathbf{f} = (f_x, f_y, f_z)` is the unit-vector of the
+            \\( \\mathbf{f} = (f_x, f_y, f_z) \\) is the unit-vector of the
             Geomagnetic field.
         """
         nr, nc = data.shape
 
         if nr != nc:
-            warnings.warn("subgrid is not square {}".format((nr,nc)), RuntimeWarning)
+            warnings.warn("subgrid is not square {}".format((nr, nc)), RuntimeWarning)
 
         fx, fy, fz = ang2vec(1.0, inc, dec)
         if sinc is None or sdec is None:
@@ -425,52 +476,50 @@ def reduce_to_pole(self, data, inc, dec, sinc=None, sdec=None):
         kx, ky = [k for k in _fftfreqs(self.dx, self.dy, data.shape)]
         kz = np.hypot(kx, ky)
 
-        a1 = mz*fz - mx*fx
-        a2 = mz*fz - my*fy
-        a3 = -my*fx - mx*fy
-        b1 = mx*fz + mz*fx
-        b2 = my*fz + mz*fy
+        a1 = mz * fz - mx * fx
+        a2 = mz * fz - my * fy
+        a3 = -my * fx - mx * fy
+        b1 = mx * fz + mz * fx
+        b2 = my * fz + mz * fy
 
         # The division gives a RuntimeWarning because of the zero frequency term.
         # This suppresses the warning.
-        with np.errstate(divide='ignore', invalid='ignore'):
-            rtp = (kz)/(a1*kx**2 + a2*ky**2 + a3*kx*ky + \
-                            1j*np.sqrt(kz)*(b1*kx + b2*ky))
-        
+        with np.errstate(divide="ignore", invalid="ignore"):
+            rtp = (kz) / (
+                a1 * kx ** 2
+                + a2 * ky ** 2
+                + a3 * kx * ky
+                + 1j * np.sqrt(kz) * (b1 * kx + b2 * ky)
+            )
+
         rtp[0, 0] = 0
-        ft_pole = rtp*np.fft.fft2(data)
+        ft_pole = rtp * np.fft.fft2(data)
         return np.real(np.fft.ifft2(ft_pole))
 
-
     def upward_continuation(self, data, height):
-        r"""
+        """
         Upward continuation of potential field data.
 
         Calculates the continuation through the Fast Fourier Transform in
         the wavenumber domain (Blakely, 1996):
 
-            F\{h_{up}\} = F\{h\} e^{-\Delta z |k|}
+        \\( F\\{h_{up}\\} = F\\{h\\} e^{-\\Delta z |k|} \\)
 
-        and then transformed back to the space domain. :math:`h_{up}` is the
-        upward continue data, \Delta z is the height increase, F denotes
-        the Fourier Transform, |k| is the wavenumber modulus.
+        and then transformed back to the space domain. \\( h_{up} \\) is the
+        upward continue data, \\( \\Delta z \\) is the height increase,
+        \\( F \\) denotes the Fourier Transform,
+        \\( |k| \\) is the wavenumber modulus.
 
         Args:
-            data   : 2D array
+            data : 2D array
                 potential field at the grid points
             height : float
                 height increase (delta z) in meters.
-        
+
         Returns:
             cont : array
                 upward continued data
 
-        Notes:
-            It is not possible to get the FFT of a masked grid. The default
-            :func:`fatiando.gridder.interp` call using minimum curvature will
-            not be suitable.  Use extrapolate=True or algorithm='nearest' to
-            get an unmasked grid.
-
         References:
             Blakely, R. J. (1996), Potential Theory in Gravity and Magnetic
             Applications, Cambridge University Press.
@@ -478,25 +527,28 @@ def upward_continuation(self, data, height):
         nr, nc = data.shape
 
         if nr != nc:
-            warnings.warn("subgrid is not square {}".format((nr,nc)), RuntimeWarning)
+            warnings.warn("subgrid is not square {}".format((nr, nc)), RuntimeWarning)
 
         if height <= 0:
-            warnings.warn("Using 'height' <= 0 means downward continuation, " +
-                          "which is known to be unstable.")
+            warnings.warn(
+                "Using 'height' <= 0 means downward continuation, "
+                + "which is known to be unstable."
+            )
 
-        fx = 2.0*np.pi*np.fft.fftfreq(nr, self.dx)
-        fy = 2.0*np.pi*np.fft.fftfreq(nc, self.dy)
+        fx = 2.0 * np.pi * np.fft.fftfreq(nr, self.dx)
+        fy = 2.0 * np.pi * np.fft.fftfreq(nc, self.dy)
 
         kx, ky = np.meshgrid(fy, fx)[::-1]
         kz = np.hypot(kx, ky)
 
-        upcont_ft = np.fft.fft2(data)*np.exp(-height*kz)
+        upcont_ft = np.fft.fft2(data) * np.exp(-height * kz)
         cont = np.real(np.fft.ifft2(upcont_ft))
         return cont
 
 
 # Helper functions to calculate Curie depth
 
+
 def bouligand2009(kh, beta, zt, dz, C):
     """
     Calculate the synthetic radial power spectrum of
@@ -531,13 +583,19 @@ def bouligand2009(kh, beta, zt, dz, C):
         129, 163-168, doi:10.1111/j.1365-246X.1997.tb00945.x
     """
     # from scipy.special import kv
-    khdz = kh*dz
+    khdz = kh * dz
     coshkhdz = np.cosh(khdz)
 
-    Phi1d = C - 2.0*kh*zt - (beta-1.0)*np.log(kh) - khdz
-    A = np.sqrt(np.pi)/gamma(1.0+0.5*beta) * \
-        (0.5*coshkhdz*gamma(0.5*(1.0+beta)) - \
-        kv((-0.5*(1.0+beta)), khdz) * np.power(0.5*khdz,(0.5*(1.0+beta)) ))
+    Phi1d = C - 2.0 * kh * zt - (beta - 1.0) * np.log(kh) - khdz
+    A = (
+        np.sqrt(np.pi)
+        / gamma(1.0 + 0.5 * beta)
+        * (
+            0.5 * coshkhdz * gamma(0.5 * (1.0 + beta))
+            - kv((-0.5 * (1.0 + beta)), khdz)
+            * np.power(0.5 * khdz, (0.5 * (1.0 + beta)))
+        )
+    )
     Phi1d += np.log(A)
     return Phi1d
 
@@ -554,7 +612,7 @@ def tanaka1999(k, lnPhi, sigma_lnPhi, kmin_range=(0.05, 0.2), kmax_range=(0.05,
             expected in ln(sqrt(S)) form
         sigma_lnPhi : standard deviation of lnPhi
         kmin_range : tuple (default:(0.05, 0.2))
-            minimum and maximum range of spatial frequencies to fit for the 
+            minimum and maximum range of spatial frequencies to fit for the
             top of magnetic sources - ideally low frequency, straight line
         kmax_range : tuple (default:(0.05, 0.2))
             minimum and maximum range of spatial frequencies to fit for the
@@ -569,56 +627,53 @@ def tanaka1999(k, lnPhi, sigma_lnPhi, kmin_range=(0.05, 0.2), kmax_range=(0.05,
     """
     # for now...
     S = lnPhi
-    sigma2 = sigma_lnPhi**2
+    sigma2 = sigma_lnPhi ** 2
 
     def compute_coefficients(X, Y, E):
-        X2 = X**2
-        Y2 = Y**2
-        E2 = E**2
+        X2 = X ** 2
+        Y2 = Y ** 2
+        E2 = E ** 2
 
         XY = np.multiply(X, Y)
-        XE2sum = np.sum(X/E2)
-        YE2sum = np.sum(Y/E2)
-        rE2sum = np.sum(1.0/E2)
-        X2E2sum = np.sum(X2/E2)
-        Y2E2sum = np.sum(Y2/E2)
+        XE2sum = np.sum(X / E2)
+        YE2sum = np.sum(Y / E2)
+        rE2sum = np.sum(1.0 / E2)
+        X2E2sum = np.sum(X2 / E2)
 
-        #TL = XE2sum*YE2sum - np.sum(XY/E2*rE2sum)
+        # TL = XE2sum*YE2sum - np.sum(XY/E2*rE2sum)
         # I think summation in second TL term needed to be split
-        TL = XE2sum*YE2sum - np.sum(XY/E2)*rE2sum
-        BL = XE2sum**2 - X2E2sum*rE2sum
+        TL = XE2sum * YE2sum - np.sum(XY / E2) * rE2sum
+        BL = XE2sum ** 2 - X2E2sum * rE2sum
 
-        Z  = TL/BL
-        b  = (np.sum(XY/E2) - Z*X2E2sum)/XE2sum
-        #dZ = np.sqrt( rE2sum/(X2E2sum*rE2sum - XE2sum) )
+        Z = TL / BL
+        b = (np.sum(XY / E2) - Z * X2E2sum) / XE2sum
+        # dZ = np.sqrt( rE2sum/(X2E2sum*rE2sum - XE2sum) )
         ## There was a missing **2 term at end of error term.
-        dZ = np.sqrt( rE2sum/(X2E2sum*rE2sum - XE2sum**2) )
+        dZ = np.sqrt(rE2sum / (X2E2sum * rE2sum - XE2sum ** 2))
         return Z, b, dZ
 
-
-    sf=k/(2.0*np.pi)
-    S2=np.log(np.exp(S)/sf)
+    sf = k / (2.0 * np.pi)
 
     # mask low wavenumbers
     kmin, kmax = kmin_range
-    mask1 = np.logical_and(sf >=kmin, sf <=kmax)
+    mask1 = np.logical_and(sf >= kmin, sf <= kmax)
     X1 = sf[mask1]
     Y1 = S[mask1]
     E1 = sigma2[mask1]
-    
+
     # mask high wavenumbers
     kmin, kmax = kmax_range
-    mask2 = np.logical_and(sf >=kmin, sf <=kmax)
+    mask2 = np.logical_and(sf >= kmin, sf <= kmax)
     X2 = sf[mask2]
-    Y2 = np.log(np.exp(S[mask2])/(X2*2*np.pi))
-    E2 = np.log(np.exp(sigma2[mask2])/(X2*2*np.pi))
+    Y2 = np.log(np.exp(S[mask2]) / (X2 * 2 * np.pi))
+    E2 = np.log(np.exp(sigma2[mask2]) / (X2 * 2 * np.pi))
 
     # compute top and bottom of magnetic layer
     Ztr, btr, dZtr = compute_coefficients(X1, Y1, E1)
     Zor, bor, dZor = compute_coefficients(X2, Y2, E2)
-    return (Ztr,btr,dZtr), (Zor, bor, dZor)
-    
-    
+    return (Ztr, btr, dZtr), (Zor, bor, dZor)
+
+
 def ComputeTanaka(Ztr, dZtr, Zor, dZor):
     """
     Compute the Curie depth from the results of tanaka1999
@@ -628,19 +683,19 @@ def ComputeTanaka(Ztr, dZtr, Zor, dZor):
             top of the magnetic source
         dZtr : float / 1D array
             error of Ztr
-        Zor : float / 1D array 
-            bottom of the magnetic source
+        Zor : float / 1D array
+            centroid depth of the magnetic source
         dZor : float / 1D array
             error of Zor
 
     Returns:
         Zb : float / 1D array
-            estimated Curie point depth
+            estimated Curie point depth at bottom of magnetic source
         eZb : float / 1D array
             error of `Zb`
     """
-    Zb = 2.0*Zor - Ztr
-    dZb  = 2.0*dZor - dZtr
+    Zb = 2.0 * Zor - Ztr
+    dZb = 2.0 * dZor + dZtr
     return abs(Zb), dZb
 
 
@@ -677,17 +732,18 @@ def maus1995(beta, zt, kh, C=0.0):
         estimation using a self-similar magnetization model, Geophys. J. Int.,
         129, 163-168, doi:10.1111/j.1365-246X.1997.tb00945.x
     """
-    return C - 2.0*kh*zt - (beta-1.0)*np.log(kh)
+    return C - 2.0 * kh * zt - (beta - 1.0) * np.log(kh)
 
 
 def _fftfreqs(dx, dy, shape):
     """
     Get two 2D-arrays with the wave numbers in the x and y directions.
     """
-    fx = 2.0*np.pi*np.fft.fftfreq(shape[0], dx)
-    fy = 2.0*np.pi*np.fft.fftfreq(shape[1], dy)
+    fx = 2.0 * np.pi * np.fft.fftfreq(shape[0], dx)
+    fy = 2.0 * np.pi * np.fft.fftfreq(shape[1], dy)
     return np.meshgrid(fy, fx)[::-1]
 
+
 def ang2vec(intensity, inc, dec):
     """
     Convert intensity, inclination and  declination to a 3-component vector
@@ -743,8 +799,10 @@ def dircos(inc, dec):
         Inclination is positive down and declination is measured with respect
         to x (North).
     """
-    d2r = np.pi/180.0
-    vect = [np.cos(d2r * inc) * np.cos(d2r * dec), \
-            np.cos(d2r * inc) * np.sin(d2r * dec), \
-            np.sin(d2r * inc)]
+    d2r = np.pi / 180.0
+    vect = [
+        np.cos(d2r * inc) * np.cos(d2r * dec),
+        np.cos(d2r * inc) * np.sin(d2r * dec),
+        np.sin(d2r * inc),
+    ]
     return vect
diff --git a/pycurious/mapping.py b/pycurious/mapping.py
index 6513265..84bc7e6 100644
--- a/pycurious/mapping.py
+++ b/pycurious/mapping.py
@@ -1,36 +1,85 @@
+# Copyright 2018-2019 Ben Mather, Robert Delhaye
+#
+# This file is part of PyCurious.
+#
+# PyCurious is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any later version.
+#
+# PyCurious is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+
 """
-Copyright 2018 Ben Mather, Robert Delhaye
+The `pycurious.mapping` module of PyCurious contains various functions to help
+manipulate geospatial data into common formats. It handles commonly encounted
+operations, such as:
+
+- Gridding scattered data points
+- Converting between coordinate reference systems (CRS)
+- Importing and exporting GeoTiff files
+
+It requires some **additional dependencies**:
 
-This file is part of PyCurious.
+- [`matplotlib`](https://matplotlib.org/) - for plotting
+- [`pyproj`](https://github.com/jswhit/pyproj) - for transforming between different CRS
+- [`cartopy`](https://scitools.org.uk/cartopy/docs/latest/) - for generating maps
 
-PyCurious is free software: you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation, either version 3 of the License, or any later version.
+Beware that most global data are georeferenced in WGS84 (EPSG: 4326).
+The radial power spectrum must be in rad/km, which requires a transformation
+from longitude / latitude to a local projection in eastings / northings.
 
-PyCurious is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU Lesser General Public License for more details.
+For example, EMAG2 is a global compilation of the magnetic anomaly
+georeferenced in WGS84 longitude / latitude. This will need to be projected
+in a local CRS to use with PyCurious. If, for example, we are interested in a
+region across Ireland we could use the IRENET95 local CRS (EPSG: 2157),
 
-You should have received a copy of the GNU Lesser General Public License
-along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+```python
+transform_coordinates(lons, lats, epsg_in=4326, epsg_out=2157)
+```
+
+which would return a list of eastings and northings in IRENET95 projection.
 """
 
 # -*- coding: utf-8 -*-
 import numpy as np
 
-try: range=xrange
-except: pass
+try:
+    range = xrange
+except:
+    pass
+
 
 def transform_coordinates(x, y, epsg_in, epsg_out):
     """
     Transform between any coordinate system.
 
-    Requires `pyproj`
+    **Requires `pyproj`** - install using pip.
+
+    Args:
+        x : float / 1D array
+            x coordinates (may be in degrees or metres/eastings)
+        y : float / 1D array
+            y coordinates (may be in degrees or metres/northings)
+        epsg_in : int
+            CRS of x and y coordinates
+        epsg_out : int
+            CRS of output
+
+    Returns:
+        x_out : float / list of floats
+            x coordinates projected in `epsg_out`
+        y_out : float / list of floats
+            y coordinates projected in `epsg_out`
     """
     import pyproj
-    proj_in  = pyproj.Proj("+init=EPSG:"+str(epsg_in))
-    proj_out = pyproj.Proj("+init=EPSG:"+str(epsg_out))
+
+    proj_in = pyproj.Proj("+init=EPSG:" + str(epsg_in))
+    proj_out = pyproj.Proj("+init=EPSG:" + str(epsg_out))
     return pyproj.transform(proj_in, proj_out, x, y)
 
 
@@ -62,7 +111,7 @@ def trim(coords, data, extent, buffer_amount=0.0):
     """
     Trim a smaller section of a large dataset taking into
     consideration transformations into various coordinate
-    reference systems (CRS)
+    reference systems (CRS).
     
     Args:
         coords : array shape (n,2)
@@ -83,21 +132,18 @@ def trim(coords, data, extent, buffer_amount=0.0):
     xmin, xmax, ymin, ymax = extent
 
     # Extract only the data within the extent
-    data_mask = np.ones(data.shape[0], dtype=bool)
+    data_mask = np.zeros(data.shape[0], dtype=bool)
 
     # Add a 1 percent buffer zone
-    x_buffer = buffer_amount*(xmax - xmin)
-    y_buffer = buffer_amount*(ymax - ymin)
-
-    mask_e = coords[:,0] < xmin - x_buffer
-    mask_w = coords[:,0] > xmax + x_buffer
-    mask_n = coords[:,1] < ymin - y_buffer
-    mask_s = coords[:,1] > ymax + y_buffer
-    data_mask[mask_n] = False
-    data_mask[mask_s] = False
-    data_mask[mask_e] = False
-    data_mask[mask_w] = False
-    
+    x_buffer = buffer_amount * (xmax - xmin)
+    y_buffer = buffer_amount * (ymax - ymin)
+
+    data_mask += coords[:, 0] < xmin - x_buffer
+    data_mask += coords[:, 0] > xmax + x_buffer
+    data_mask += coords[:, 1] < ymin - y_buffer
+    data_mask += coords[:, 1] > ymax + y_buffer
+    data_mask = np.invert(data_mask)
+
     data_trim = data[data_mask]
     coords_trim = coords[data_mask]
 
@@ -108,7 +154,9 @@ def grid(coords, data, extent, shape=None, epsg_in=None, epsg_out=None, **kwargs
     """
     Grid a smaller section of a large dataset taking into
     consideration transformations into various coordinate
-    reference systems (CRS)
+    reference systems (CRS).
+
+    **Requires `scipy.interpolate.griddata`**
     
     Args:
         coords : array shape (n,2)
@@ -116,7 +164,7 @@ def grid(coords, data, extent, shape=None, epsg_in=None, epsg_out=None, **kwargs
         data : array shape (n,) 
             values corresponding to coordinates
         extent : tuple
-           box contained within the data in espg_out coordinates
+           bounding box in espg_out coordinates
         shape : tuple (nrows,ncols)
            size of the box, if None, shape is estimated from coords spacing
         epsg_in : int
@@ -132,12 +180,16 @@ def grid(coords, data, extent, shape=None, epsg_in=None, epsg_out=None, **kwargs
             rectangular section of data bounded by extent
     """
     from scipy.interpolate import griddata
+
     xmin, xmax, ymin, ymax = extent
-    
-    if type(epsg_in) != type(None):
-        xt, yt = transform_coordinates(np.array([xmin, xmin, xmax, xmax]),\
-                                       np.array([ymin, ymax, ymin, ymax]),\
-                                       epsg_out, epsg_in)
+
+    if epsg_in is not None:
+        xt, yt = transform_coordinates(
+            np.array([xmin, xmin, xmax, xmax]),
+            np.array([ymin, ymax, ymin, ymax]),
+            epsg_out,
+            epsg_in,
+        )
         # find the coordinates that will completely
         # engulf the extent
         xtmin, xtmax = min(xt), max(xt)
@@ -151,24 +203,24 @@ def grid(coords, data, extent, shape=None, epsg_in=None, epsg_out=None, **kwargs
     # trim data - buffer = 5%
     coords_trim, data_trim = trim(coords, data, xtextent, 0.05)
 
-
-    if type(epsg_in) != type(None):
+    if epsg_in is not None:
         # convert back to output CRS
-        xtrim, ytrim = transform_coordinates(coords_trim[:,0],\
-                                             coords_trim[:,1],\
-                                             epsg_in, epsg_out)
+        xtrim, ytrim = transform_coordinates(
+            coords_trim[:, 0], coords_trim[:, 1], epsg_in, epsg_out
+        )
         coords_trim = np.column_stack([xtrim, ytrim])
 
-
     if shape == None:
         # estimate based on the data spacing
-        xunique = np.unique(coords_trim[:,0])
-        yunique = np.unique(coords_trim[:,1])
+        xunique = np.unique(coords_trim[:, 0])
+        yunique = np.unique(coords_trim[:, 1])
         dx = np.diff(xunique).mean()
         dy = np.diff(yunique).mean()
-        nc = int((xtmax - xtmin)/dx)
-        nr = int((ytmax - ytmin)/dy)
-        print("using nrows={}, ncols={} with cell spacing of {}".format(nr,nc,(dy,dx)))
+        nc = int((xtmax - xtmin) / dx)
+        nr = int((ytmax - ytmin) / dy)
+        print(
+            "using nrows={}, ncols={} with cell spacing of {}".format(nr, nc, (dy, dx))
+        )
     else:
         nr, nc = shape
 
@@ -182,10 +234,52 @@ def grid(coords, data, extent, shape=None, epsg_in=None, epsg_out=None, **kwargs
     return vq
 
 
+def ungrid(grid, extent, coordinates, **kwargs):
+    """
+    Maps the value at unstructured coordinates from a regular 2D grid.
+    Uses the `scipy.ndimage.map_coordinates` function with cubic interpolation
+    
+    Args:
+        grid : array shape (ny,nx)
+            regularly spaced grid 
+        extent : tuple
+            bounding box of `grid` e.g. [xmin,xmax,ymin,ymax]
+        coordinates : array shape (n,2)
+            coordinates (x,y) to interpolate `grid`
+        kwargs : keyword arguments
+            keyword arguments to pass to `map_coordinates`
+     
+    Returns:
+        grid_interp : array shape (n,)
+            interpolated values at coordinates
+    """
+    from scipy.ndimage import map_coordinates
+    
+    ny, nx = grid.shape
+    xmin, xmax, ymin, ymax = extent
+    
+    # normalise coordinates within extent
+    icoords = coordinates.copy()
+    icoords[:,0] -= xmin
+    icoords[:,1] -= ymin
+    icoords[:,0] /= (xmax - xmin)
+    icoords[:,1] /= (ymax - ymin)
+    
+    # icoords now somewhere within range [0, 1]
+    # project coordinates to the number of grid indices
+    icoords[:,0] *= nx - 1
+    icoords[:,1] *= ny - 1
+    
+    # now interpolate
+    return map_coordinates(grid, icoords.T, **kwargs)
+
+
 def import_geotiff(file_path):
     """
     Import a GeoTIFF to a numpy array and prints
-    information of the Coordinate Reference System (CRS)
+    information of the Coordinate Reference System (CRS).
+
+    **Requires `osgeo`.**
 
     Args:
         file_path : str
@@ -194,11 +288,11 @@ def import_geotiff(file_path):
     Returns:
         data : 2D array
         extent : tuple
-            extent in the projection of the GeoTIFF
+            bounding box in the projection of the GeoTIFF
             e.g. [xmin, xmax, ymin, ymax]
     """
     from osgeo import gdal, osr
-    
+
     gtiff = gdal.Open(file_path)
     data = gtiff.ReadAsArray()
     gt = gtiff.GetGeoTransform()
@@ -207,8 +301,12 @@ def import_geotiff(file_path):
     inproj = osr.SpatialReference()
     inproj.ImportFromWkt(gtproj)
 
-    gtextent = (gt[0], gt[0] + gtiff.RasterXSize*gt[1],\
-                gt[3], gt[3] + gtiff.RasterYSize*gt[5])
+    gtextent = (
+        gt[0],
+        gt[0] + gtiff.RasterXSize * gt[1],
+        gt[3],
+        gt[3] + gtiff.RasterYSize * gt[5],
+    )
 
     # print projection information
     print(inproj)
@@ -222,7 +320,9 @@ def import_geotiff(file_path):
 def export_geotiff(file_path, array, extent, epsg):
     """
     Export a GeoTIFF from a numpy array projected in a
-    predefined Coordinate Reference System (CRS)
+    predefined Coordinate Reference System (CRS).
+
+    **Requires `osgeo`.**
 
     Args:
         file_path : str
@@ -230,7 +330,7 @@ def export_geotiff(file_path, array, extent, epsg):
         array: 2D array
             array to save to GeoTiff
         extent : tuple
-            extent in the projection of the GeoTIFF
+            bounding box in the projection of the GeoTIFF
             e.g. [xmin, xmax, ymin, ymax]
         epsg : int
             CRS of the GeoTIFF
@@ -238,16 +338,17 @@ def export_geotiff(file_path, array, extent, epsg):
 
     """
     from osgeo import gdal, osr
+
     # import ogr, gdal, osr, os
 
     cols = array.shape[1]
     rows = array.shape[0]
 
     xmin, xmax, ymin, ymax = extent
-    spacingX = (xmax - xmin)/cols
-    spacingY = (ymax - ymin)/rows
+    spacingX = (xmax - xmin) / cols
+    spacingY = (ymax - ymin) / rows
 
-    driver = gdal.GetDriverByName('GTiff')
+    driver = gdal.GetDriverByName("GTiff")
     outRaster = driver.Create(file_path, cols, rows, 1, gdal.GDT_Float64)
     outRaster.SetGeoTransform((xmin, spacingX, 0, ymin, 0, spacingY))
     outband = outRaster.GetRasterBand(1)
@@ -256,4 +357,86 @@ def export_geotiff(file_path, array, extent, epsg):
     outRasterSRS.ImportFromEPSG(epsg)
     outRaster.SetProjection(outRasterSRS.ExportToWkt())
     outband.FlushCache()
-    return
\ No newline at end of file
+    return
+
+
+def export_netcdf4(file_path, array, extent):
+    """
+    Export a netCDF4 file from a numpy array
+    over a user-defined extent.
+
+    **Requires `netcdf4`.**
+    `pip install netcdf4`
+
+    Args:
+        file_path : str
+            path to write the netCDF4
+        array: 2D array
+            array to save to netCDF4
+        extent : tuple
+            bounding box in the projection of the netCDF4 file
+            e.g. [xmin, xmax, ymin, ymax]
+    """
+    import netCDF4
+
+    ny, nx = array.shape
+    xmin, xmax, ymin, ymax = extent
+
+    with netCDF4.Dataset(str(file_path), 'w') as cdf:
+        cdf.createDimension('x', nx)
+        cdf.createDimension('y', ny)
+        cdf_x = cdf.createVariable('x', np.float, ('x',), zlib=True)
+        cdf_y = cdf.createVariable('y', np.float, ('y',), zlib=True)
+        cdf_x[:] = np.linspace(xmin, xmax, nx)
+        cdf_y[:] = np.linspace(ymin, ymax, ny)
+
+        cdf_data = cdf.createVariable('z', np.float, ('y','x'), zlib=True)
+        cdf_data[:,:]  = array
+
+
+def import_netcdf4(file_path):
+    """
+    Import a netCDF4 file from a numpy array
+    over a user-defined extent.
+
+    **Requires `netcdf4`.**
+    `pip install netcdf4`
+
+    Args:
+        file_path : str
+            path to write the netCDF4
+
+    Returns:
+        array : ndarray
+            numpy array containing the gridded netCDF4 data
+        extent : tuple
+            bounding box in the projection of the netCDF4 file
+            e.g. [xmin, xmax, ymin, ymax]
+
+    Note:
+        The layout of the netCDF4 file must follow the standard
+        naming convention such as 'lon', 'lat', 'z', 'data', etc 
+        otherwise this import function will not work as intended.
+    """
+    import netCDF4
+
+    with netCDF4.Dataset(str(file_path), 'r') as cdf:
+        try:
+            lons = cdf['lon']
+            lats = cdf['lat']
+        except:
+            lons = cdf['x']
+            lats = cdf['y']
+        else:
+            raise ValueError(print(cdf))
+
+        try:
+            data = cdf['z'][:]
+        except:
+            data = cdf['data'][:]
+        else:
+            raise ValueError(print(cdf))
+
+        extent = [lons.min(), lons.max(), lats.min(), lats.max()]
+
+    return data, extent
\ No newline at end of file
diff --git a/pycurious/optimise.py b/pycurious/optimise.py
index 6de0f10..7439968 100644
--- a/pycurious/optimise.py
+++ b/pycurious/optimise.py
@@ -1,19 +1,35 @@
+# Copyright 2018-2019 Ben Mather, Robert Delhaye
+#
+# This file is part of PyCurious.
+#
+# PyCurious is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any later version.
+#
+# PyCurious is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+
 """
-Copyright 2018 Ben Mather, Robert Delhaye
+This PyCurious module contains the `pycurious.optimise.CurieOptimise` class,
+which inherits the `pycurious.grid.CurieGrid` class with added functionality for:
 
-This file is part of PyCurious.
+- Fitting the synthetic power spectrum \\( \\Phi \\) computed with `pycurious.grid.bouligand2009`
+- Defining an objective function with a flexible interface for adding *a priori* and likelihood functions
+- Parallel decomposition of routines to compute the optimal radial power spectrum, and thus Curie depth
+- Metropolis-Hastings algorithm and sensitivity analysis to estimate the uncertainty of the posterior
 
-PyCurious is free software: you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation, either version 3 of the License, or any later version.
+The posterior is defined as
 
-PyCurious is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU Lesser General Public License for more details.
+\\( P(\\mathbf{m}|\\mathbf{d}) = P(\\beta, z_t, \\Delta z, C|\\Phi_d) \\)
 
-You should have received a copy of the GNU Lesser General Public License
-along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
+where \\( \\beta, z_t, \\Delta z, C \\) are input parameters to `pycurious.grid.bouligand2009` and
+\\( \\Phi_d \\) is the radial power spectrum computed from a FFT over square windows of the magnetic
+anomaly in `pycurious.grid.CurieGrid.radial_spectrum`.
 """
 
 # -*- coding: utf-8 -*-
@@ -25,14 +41,17 @@
 from scipy import stats
 from multiprocessing import Pool, Process, Queue, cpu_count
 
-try: range=xrange
-except: pass
+try:
+    range = xrange
+except:
+    pass
 
 
 class CurieOptimise(CurieGrid):
     """
-    Extends the CurieGrid class to include optimisation routines
-    see scipy.optimize.minimize for a description of the algorithm
+    Extends the `pycurious.grid.CurieGrid` class to include
+    optimisation routines see `scipy.optimize.minimize` for
+    a description of the algorithm.
 
     Args:
         grid : 2D numpy array
@@ -48,10 +67,38 @@ class CurieOptimise(CurieGrid):
     
     Attributes:
         bounds : list of tuples
-            lower and upper bounds for beta, zt, dz, C
-        prior  : dict
-            dictionary of priors for beta, zt, dz, C
+            lower and upper bounds for \\( \\beta, z_t, \\Delta z, C \\)
+        prior : dict
+            dictionary of priors for \\( \\beta, z_t, \\Delta z, C \\)
+        grid : 2D numpy array
+            2D array of magnetic data
+        xmin : float
+            minimum x bound in metres
+        xmax : float
+            maximum x bound in metres
+        ymin : float
+            minimum y bound in metres
+        ymax : float
+            maximum y bound in metres
+        dx : float
+            grid spacing in the x-direction in metres
+        dy : float
+            grid spacing in the y-direction in metres
+        nx : int
+            number of nodes in the x-direction
+        ny : int
+            number of nodes in the y-direction
+        xcoords : 1D numpy array
+            1D numpy array of coordinates in the x-direction
+        ycoords : 1D numpy array
+            1D numpy array of coordinates in the y-direction
+
+    Notes:
+        In all instances `x` indicates eastings in metres and `y` indicates northings.
+        Using a grid of longitude / latitudinal coordinates (degrees) will result
+        in incorrect Curie depth calculations.
     """
+
     def __init__(self, grid, xmin, xmax, ymin, ymax, **kwargs):
 
         super(CurieOptimise, self).__init__(grid, xmin, xmax, ymin, ymax)
@@ -62,8 +109,8 @@ def __init__(self, grid, xmin, xmax, ymin, ymax, **kwargs):
         # lower / upper bounds
         # [beta, zt, dz, C]
         lb = [0.0, 0.0, 0.0, None]
-        ub = [None]*len(lb)
-        bounds = list(zip(lb,ub))
+        ub = [None] * len(lb)
+        bounds = list(zip(lb, ub))
         self.bounds = bounds
 
         self.max_processors = kwargs.pop("max_processors", cpu_count())
@@ -73,10 +120,10 @@ def __init__(self, grid, xmin, xmax, ymin, ymax, **kwargs):
     def add_prior(self, **kwargs):
         """
         Add a prior to the dictionary (tuple)
-        Available priors are beta, zt, dz, C
+        Available priors are \\( \\beta, z_t, \\Delta z, C \\)
 
         Assumes a normal distribution or
-        define another distribution from scipy.stats
+        define another distribution from `scipy.stats`
 
         Usage:
             >>> add_prior(beta=(p, sigma_p))
@@ -102,14 +149,12 @@ def add_prior(self, **kwargs):
             else:
                 raise ValueError("prior must be one of {}".format(self.prior.keys()))
 
-
     def reset_priors(self):
         """
         Reset priors to uniform distribution
         """
-        self.prior = {'beta':None, 'zt':None, 'dz':None, 'C':None}
-        self.prior_pdf = {'beta':None, 'zt':None, 'dz':None, 'C':None}
-
+        self.prior = {"beta": None, "zt": None, "dz": None, "C": None}
+        self.prior_pdf = {"beta": None, "zt": None, "dz": None, "C": None}
 
     def objective_routine(self, **kwargs):
         """
@@ -120,7 +165,7 @@ def objective_routine(self, **kwargs):
             >>> objective_routine(beta=2.5)
 
         Returns:
-            misfit  : float
+            misfit : float
                 misfit integrated over all observations and priors
         """
         c = 0.0
@@ -144,10 +189,9 @@ def objective_function(self, x, x0, sigma_x0, *args):
             sigma_x0 : float, ndarray
 
         Returns:
-            misfit   : float
+            misfit : float
         """
-        return 0.5*np.sum((x - x0)**2/sigma_x0**2)
-
+        return 0.5 * np.sum((x - x0) ** 2 / sigma_x0 ** 2)
 
     def min_func(self, x, kh, Phi, sigma_Phi):
         """
@@ -155,20 +199,20 @@ def min_func(self, x, kh, Phi, sigma_Phi):
 
         Args:
             x : array shape (n,)
-                array of variables [beta, zt, dz, C]
+                array of variables \\( \\beta, z_t, \\Delta z, C \\)
             kh : array shape (n,)
                 wavenumbers (rad/km)
             Phi : array shape (n,)
-                radial power spectrum
+                radial power spectrum \\( \\Phi \\)
             sigma_Phi : array shape (n,)
-                standard deviation of Phi
+                standard deviation of Phi, \\( \\sigma_{\\Phi} \\)
 
         Returns:
-            misfit    : float
-                sum of misfit (float)
+            misfit : float
+                sum of misfit (scalar)
 
         Notes:
-            We purposely ignore all warnings raised by the bouligand2009
+            We purposely ignore all warnings raised by the `pycurious.grid.bouligand2009`
             function because some combinations of input parameters will
             trigger an out-of-range warning that will crash the minimiser.
             Instead, the misfit is set to a very large number when this occurs.
@@ -185,14 +229,25 @@ def min_func(self, x, kh, Phi, sigma_Phi):
             misfit += self.objective_routine(beta=beta, zt=zt, dz=dz, C=C)
         return misfit
 
-
-    def optimise(self, window, xc, yc, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=np.hanning, process_subgrid=None, **kwargs):
+    def optimise(
+        self,
+        window,
+        xc,
+        yc,
+        beta=3.0,
+        zt=1.0,
+        dz=10.0,
+        C=5.0,
+        taper=np.hanning,
+        process_subgrid=None,
+        **kwargs
+    ):
         """
-        Find the optimal parameters of beta, zt, dz, C for a given
-        centroid (xc,yc) and window size.
+        Find the optimal parameters of \\( \\beta, z_t, \\Delta z, C \\)
+        for a given centroid (xc,yc) and window size.
 
         Args:
-            window  : float
+            window : float
                 size of window in metres
             xc : float
                 centroid x values
@@ -206,7 +261,7 @@ def optimise(self, window, xc, yc, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=np.ha
                 thickness of magnetic layer (starting value)
             C : float
                 field constant (starting value)
-            taper : taper (default=np.hanning)
+            taper : taper (default=`numpy.hanning`)
                 taper function, set to None for no taper function
             process_subgrids : function
                 a custom function to process the subgrid
@@ -224,7 +279,7 @@ def optimise(self, window, xc, yc, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=np.ha
                 field constant
         """
 
-        if type(process_subgrid) == type(None):
+        if process_subgrid is None:
             # dummy function
             def process_subgrid(subgrid):
                 return subgrid
@@ -245,7 +300,6 @@ def process_subgrid(subgrid):
         res = minimize(self.min_func, x0, args=(k, Phi, sigma_Phi), bounds=self.bounds)
         return res.x
 
-
     def _func_queue(self, func, q_in, q_out, window, *args, **kwargs):
         """ Retrive processes from the queue """
         while True:
@@ -260,7 +314,6 @@ def _func_queue(self, func, q_in, q_out, window, *args, **kwargs):
             q_out.put((pos, res))
         return
 
-
     def parallelise_routine(self, window, xc_list, yc_list, func, *args, **kwargs):
         """
         Implements shared memory multiprocessing to split multiple
@@ -270,22 +323,22 @@ def parallelise_routine(self, window, xc_list, yc_list, func, *args, **kwargs):
         along with any additional arguments or keyword arguments.
 
         Args:
-         window : float
-            size of window in metres
-         xc_list : array shape (l,)
-            centroid x values
-         yc_list : array shape (l,)
-            centroid y values
-         func : function
-            Python function to evaluate in parallel
-         args : arguments
-            additional arguments to pass to func
-         kwargs : keyword arguments
-            additional keyword arguments to pass to func
+            window : float
+                size of window in metres
+            xc_list : array shape (l,)
+                centroid x values
+            yc_list : array shape (l,)
+                centroid y values
+            func : function
+                Python function to evaluate in parallel
+            args : arguments
+                additional arguments to pass to `func`
+            kwargs : keyword arguments
+                additional keyword arguments to pass to `func`
 
         Returns:
-            out : list
-                (depends on output of func - see notes)
+            out : list of lists
+                (depends on output of `func` - see notes)
 
         Usage:
             An obvious use case is to compute the Curie depth for many
@@ -295,16 +348,33 @@ def parallelise_routine(self, window, xc_list, yc_list, func, *args, **kwargs):
         
             Each centroid is assigned a new process and sent to a free processor
             to compute. In this case, the output is separate lists of shape(l,)
-            for beta, zt, dz, and C.
+            for \\( \\beta, z_t, \\Delta z, C \\). If `len(xc_list)=2` then,
+
+            >>> self.parallelise_routine(window, [x1,x2], [y1, y2], self.optimise)
+            [[beta1  beta2], [zt1  zt2], [dz1  dz2], [C1  C2]]
 
             Another example is to parallelise the sensitivity analysis:
 
             >>> self.parallelise_routine(window, xc_list, yc_list, self.sensitivity, nsim)
 
-            This time the output will be a list of lists for beta, zt, dz, and C
-            i.e. if nc=2 is the number of centroids and nsim=4 is the number of
-            simulations then separatee lists [[k1, k2, k3, k4], [k1, k2, k3, k4]]
-            will be returned for beta, zt, dz, and C.
+            This time the output will be a list of lists for \\( \\beta, z_t, \\Delta z, C \\)
+            i.e. if `len(xc_list)=2` is the number of centroids and `nsim=4` is the number of
+            simulations then separatee lists will be returned for \\( \\beta, z_t, \\Delta z, C \\).
+
+            >>> self.parallelise_routine(window, [x1,x2], [y1,y2], self.sensitivity, 4)
+
+            which would return:
+
+            ```python
+            [[[ beta1a , beta1b , beta1c , beta1d ],   # centroid 1 (x1,y1)
+              [ beta2a , beta2b , beta2c , beta2d ]],  # centroid 2 (x2,y2)
+             [[   zt1a ,   zt1b ,   zt1c ,   zt1d ],   # centroid 1 (x1,y1)
+              [   zt2a ,   zt2b ,   zt2c ,   zt2d ]],  # centroid 2 (x2,y2)
+             [[   dz1a ,   dz1b ,   dz1c ,   dz1d ],   # centroid 1 (x1,y1)
+              [   dz2a ,   dz2b ,   dz2c ,   dz2d ]]   # centroid 2 (x2,y2)
+             [[    C1a ,    C1b ,    C1c ,    C1d ],   # centroid 1 (x1,y1)
+              [    C2a ,    C2b ,    C2c ,    C2d ]]]  # centroid 2 (x2,y2)
+            ```
         """
 
         n = len(xc_list)
@@ -318,32 +388,43 @@ def parallelise_routine(self, window, xc_list, yc_list, func, *args, **kwargs):
 
         nprocs = self.max_processors
 
-        for i in range(nprocs):
-            pass_args = [func, q_in, q_out, window]
-            pass_args.extend(args)
+        if nprocs == 1:
+            # skip all the OpenMP cruft
+            for i in range(n):
+                xc = xc_list[i]
+                yc = yc_list[i]
 
-            p = Process(target=self._func_queue,\
-                        args=tuple(pass_args),\
-                        kwargs=kwargs)
+                res = func(window, xc, yc, *args, **kwargs)
+                xOpt[i] = res
 
-            processes.append(p)
+        elif nprocs > 1:
+            # more than one processor
+            for i in range(nprocs):
+                pass_args = [func, q_in, q_out, window]
+                pass_args.extend(args)
 
-        for p in processes:
-            p.daemon = True
-            p.start()
+                p = Process(target=self._func_queue, args=tuple(pass_args), kwargs=kwargs)
 
-        # put items in the queue
-        sent = [q_in.put((i, xc_list[i], yc_list[i])) for i in range(n)]
-        [q_in.put((None, None, None)) for _ in range(nprocs)]
+                processes.append(p)
 
-        # get the results
-        for i in range(len(sent)):
-            i, res = q_out.get()
-            xOpt[i] = res
+            for p in processes:
+                p.daemon = True
+                p.start()
 
+            # put items in the queue
+            sent = [q_in.put((i, xc_list[i], yc_list[i])) for i in range(n)]
+            [q_in.put((None, None, None)) for _ in range(nprocs)]
 
-        # wait until each processor has finished
-        [p.join() for p in processes]
+            # get the results
+            for i in range(len(sent)):
+                index, res = q_out.get()
+                xOpt[index] = res
+
+            # wait until each processor has finished
+            [p.join() for p in processes]
+
+        else:
+            raise ValueError("{} processors is invalid, specify a positive integer value".format(nprocs))
 
         # process dimensions of output
         ndim = np.array(res).ndim
@@ -362,12 +443,22 @@ def parallelise_routine(self, window, xc_list, yc_list, func, *args, **kwargs):
         else:
             raise ValueError("Cannot determine shape of output")
 
-
-
-    def optimise_routine(self, window, xc_list, yc_list, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=np.hanning, process_subgrid=None, **kwargs):
+    def optimise_routine(
+        self,
+        window,
+        xc_list,
+        yc_list,
+        beta=3.0,
+        zt=1.0,
+        dz=10.0,
+        C=5.0,
+        taper=np.hanning,
+        process_subgrid=None,
+        **kwargs
+    ):
         """
         Iterate through a list of centroids to compute the optimal values
-        of beta, zt, dz, C for a given window size.
+        of \\( \\beta, z_t, \\Delta z, C \\) for a given window size.
         
         Args:
             window : float
@@ -385,12 +476,12 @@ def optimise_routine(self, window, xc_list, yc_list, beta=3.0, zt=1.0, dz=10.0,
             C : float
                 field constant
             taper : function
-                taper function (default=np.hanning)
+                taper function (default=`numpy.hanning`)
                 set to None for no taper function
-         process_subgrids : func
-            a custom function to process the subgrid
-         kwargs : keyword arguments
-            to pass to radial_spectrum.
+            process_subgrids : func
+                a custom function to process the subgrid
+            kwargs : keyword arguments
+                to pass to radial_spectrum.
 
         Returns:
             beta : ndarray shape (l,)
@@ -403,15 +494,42 @@ def optimise_routine(self, window, xc_list, yc_list, beta=3.0, zt=1.0, dz=10.0,
                 field constant
 
         """
-        return self.parallelise_routine(window, xc_list, yc_list, self.optimise, beta, zt, dz, C, taper, process_subgrid, **kwargs)
-
-
-    def metropolis_hastings(self, window, xc, yc, nsim, burnin, x_scale=None, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=np.hanning, process_subgrid=None, **kwargs):
+        return self.parallelise_routine(
+            window,
+            xc_list,
+            yc_list,
+            self.optimise,
+            beta,
+            zt,
+            dz,
+            C,
+            taper,
+            process_subgrid,
+            **kwargs
+        )
+
+    def metropolis_hastings(
+        self,
+        window,
+        xc,
+        yc,
+        nsim,
+        burnin,
+        x_scale=None,
+        beta=3.0,
+        zt=1.0,
+        dz=10.0,
+        C=5.0,
+        taper=np.hanning,
+        process_subgrid=None,
+        **kwargs
+    ):
         """
         MCMC algorithm using a Metropolis-Hastings sampler.
 
-        Evaluates a Markov-Chain for starting values of beta, zt, dz, C
-        and returns the ensemble of model realisations.
+        Evaluates a Markov-Chain for starting values of
+        \\( \\beta, z_t, \\Delta z, C \\) and returns the
+        ensemble of model realisations.
         
         Args:
             window : float
@@ -426,7 +544,8 @@ def metropolis_hastings(self, window, xc, yc, nsim, burnin, x_scale=None, beta=3
                 number of burn-in simulations before to nsim
             x_scale: float(4) (optional)
                 scaling factor for new proposals
-                default is [1,1,1,1] for [beta, zt, dz, C] - see notes
+                (default=`[1,1,1,1]` for `[beta, zt, dz, C]`)
+                - see notes
             beta : float
                 fractal parameter (starting value)
             zt : float
@@ -447,9 +566,10 @@ def metropolis_hastings(self, window, xc, yc, nsim, burnin, x_scale=None, beta=3
                 field constant
 
         Notes:
-            nsim, burnin, and x_scale should be tweaked for optimal performance
-            Use starting values of beta, zt, dz, C relatively close to the solution
-            - C can easily found from the mean of the radial power spectrum.
+            `nsim`, `burnin`, and `x_scale` should be tweaked for optimal performance
+            Use starting values of \\( \\beta, z_t, \\Delta z, C \\) relatively
+            close to the solution - \\( C \\) can easily found from the mean of the
+            radial power spectrum.
 
             During the burn-in stage we apply tempering to the PDF to iterate closer
             towards the solution. This has the effect of smoothing out the posterior
@@ -457,19 +577,17 @@ def metropolis_hastings(self, window, xc, yc, nsim, burnin, x_scale=None, beta=3
             large portions of the posterior probability are zero.
             see see Sambridge 2013, DOI:10.1093/gji/ggt342 for more information.
         """
-        if type(process_subgrid) == type(None):
+        if process_subgrid is None:
             # dummy function
             def process_subgrid(subgrid):
                 return subgrid
 
-
         samples = np.empty((nsim, 4))
         x0 = np.array([beta, zt, dz, C])
 
-        if type(x_scale) == type(None):
+        if x_scale is None:
             x_scale = np.ones(4)
 
-
         # get subgrid
         subgrid = self.subgrid(window, xc, yc)
         subgrid = process_subgrid(subgrid)
@@ -477,15 +595,15 @@ def process_subgrid(subgrid):
         # compute radial spectrum
         k, Phi, sigma_Phi = self.radial_spectrum(subgrid, taper=taper, **kwargs)
 
-        P0 = np.exp(-self.min_func(x0, k, Phi, sigma_Phi)/1000)
+        P0 = np.exp(-self.min_func(x0, k, Phi, sigma_Phi) / 1000)
 
         # Burn-in phase
         for i in range(burnin):
             # add random perturbation
-            x1 = x0 + np.random.normal(size=4)*x_scale
+            x1 = x0 + np.random.normal(size=4) * x_scale
 
             # evaluate proposal probability + tempering
-            P1 = np.exp(-self.min_func(x1, k, Phi, sigma_Phi)/1000)
+            P1 = np.exp(-self.min_func(x1, k, Phi, sigma_Phi) / 1000)
 
             # iterate towards MAP estimate
             if P1 > P0:
@@ -497,13 +615,13 @@ def process_subgrid(subgrid):
         # Now sample posterior
         for i in range(nsim):
             # add random perturbation
-            x1 = x0 + np.random.normal(size=4)*x_scale
+            x1 = x0 + np.random.normal(size=4) * x_scale
 
             # evaluate proposal probability
             P0 = max(P0, 1e-99)
             P1 = np.exp(-self.min_func(x1, k, Phi, sigma_Phi))
 
-            P = min(P1/P0, 1.0)
+            P = min(P1 / P0, 1.0)
 
             # randomly accept probability
             if np.random.rand() <= P:
@@ -514,13 +632,25 @@ def process_subgrid(subgrid):
 
         return list(samples.T)
 
-
-
-    def sensitivity(self, window, xc, yc, nsim, beta=3.0, zt=1.0, dz=10.0, C=5.0, taper=np.hanning, process_subgrid=None, **kwargs):
+    def sensitivity(
+        self,
+        window,
+        xc,
+        yc,
+        nsim,
+        beta=3.0,
+        zt=1.0,
+        dz=10.0,
+        C=5.0,
+        taper=np.hanning,
+        process_subgrid=None,
+        **kwargs
+    ):
         """
         Iterate through a list of centroids to compute the mean and
-        standard deviation of beta, zt, dz, C by perturbing their
-        prior distributions (if provided by the user - see add_prior).
+        standard deviation of \\( \\beta, z_t, \\Delta z, C \\) by
+        perturbing their prior distributions
+        (if provided by the user - see add_prior).
         
         Args:
             nsim : int
@@ -532,7 +662,7 @@ def sensitivity(self, window, xc, yc, nsim, beta=3.0, zt=1.0, dz=10.0, C=5.0, ta
             yc : float
                 centroid y values
             nsim : int
-                - number of simulations
+                number of simulations
             beta : float
                 starting fractal parameter 
             zt : float
@@ -553,21 +683,20 @@ def sensitivity(self, window, xc, yc, nsim, beta=3.0, zt=1.0, dz=10.0, C=5.0, ta
             C : ndarray shape (nsim,)
                 field constant
         """
-        if type(process_subgrid) == type(None):
+        if process_subgrid is None:
             # dummy function
             def process_subgrid(subgrid):
                 return subgrid
 
-
         samples = np.empty((nsim, 4))
         x0 = np.array([beta, zt, dz, C])
-        
+
         use_keys = []
         for key in self.prior_pdf:
             prior_pdf = self.prior_pdf[key]
             if prior_pdf is not None:
                 use_keys.append(key)
-        
+
         # get subgrid
         subgrid = self.subgrid(window, xc, yc)
         subgrid = process_subgrid(subgrid)
@@ -583,10 +712,11 @@ def process_subgrid(subgrid):
 
             # minimise function
             rPhi = np.random.normal(Phi, sigma_Phi)
-            res = minimize(self.min_func, x0, args=(k, rPhi, sigma_Phi), bounds=self.bounds)
+            res = minimize(
+                self.min_func, x0, args=(k, rPhi, sigma_Phi), bounds=self.bounds
+            )
             samples[sim] = res.x
 
-
         # restore priors
         for key in use_keys:
             prior_pdf = self.prior_pdf[key]
diff --git a/pycurious/parallel.py b/pycurious/parallel.py
deleted file mode 100644
index cb75af7..0000000
--- a/pycurious/parallel.py
+++ /dev/null
@@ -1,204 +0,0 @@
-"""
-Copyright 2018 Ben Mather, Robert Delhaye
-
-This file is part of PyCurious.
-
-PyCurious is free software: you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation, either version 3 of the License, or any later version.
-
-PyCurious is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU Lesser General Public License for more details.
-
-You should have received a copy of the GNU Lesser General Public License
-along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
-"""
-
-# -*- coding: utf-8 -*-
-import numpy as np
-import warnings
-from .optimise import CurieOptimise
-
-try: range=xrange
-except: pass
-
-
-class CurieParallel(CurieOptimise):
-    """
-    A parallel implementation built on the PETSc DMDA mesh structure
-    Almost a trivially parallel problem: each processor works locally
-    on their portion of the grid. A stipulation being that the window
-    size must not exceed the thickness of shadow zones or else the
-    subgrid will be truncated.
-
-    This inherits from CurieOptimise and CurieGrid classes.
-
-
-    Parameters
-    ----------
-     grid       : 2D array of magnetic data
-     xmin       : minimum x bound in metres
-     xmax       : maximum x bound in metres
-     ymin       : minimum y bound in metres
-     ymax       : maximum y bound in metres
-     max_window : maximum size of the windows
-                : this will be enforced for all methods
-    
-    Attributes
-    ----------
-     bounds     : lower and upper bounds for beta, zt, dz, C
-     prior      : dictionary of priors for beta, zt, dz, C
-     coords     : local xy coordinates
-    """
-
-    def __init__(self, grid, xmin, xmax, ymin, ymax, max_window):
-
-        from petsc4py import PETSc
-        from mpi4py import MPI
-        comm = MPI.COMM_WORLD
-        self.comm = comm
-        self.MPI = MPI
-
-        # super(CurieParallel, self).__init__(grid, xmin, xmax, ymin, ymax)
-
-        # determine stencil width (should be window size/2)
-        ny, nx = grid.shape
-        dx = (xmax - xmin)/nx
-        dy = (ymax - ymin)/ny
-        nw = max_window/dx
-        n2w = int(nw/2) + 1 # add some buffer
-
-        if abs(dx - dy) > 1.0:
-            warnings.warn("node spacing should be identical {}".format((dx,dy)), RuntimeWarning)
-
-
-        dm = PETSc.DMDA().create(dim=2, sizes=(nx,ny), stencil_width=n2w, comm=comm)
-        dm.setUniformCoordinates(xmin, xmax, ymin, ymax)
-
-        self.dm = dm
-        self.lgmap = dm.getLGMap()
-        self.lvec = dm.createLocalVector()
-        self.gvec = dm.createGlobalVector()
-
-        coords = dm.getCoordinatesLocal().array.reshape(-1, 2)
-        xmin, ymin = coords.min(axis=0)
-        xmax, ymax = coords.max(axis=0)
-        self.coords = coords
-        self.max_window = max_window
-
-        (imin, imax), (jmin, jmax) = dm.getGhostRanges()
-
-        # now decompose grid for each processor
-
-        grid_chunk = grid[jmin:jmax, imin:imax]
-
-        super(CurieParallel, self).__init__(grid_chunk, xmin, xmax, ymin, ymax)
-
-        reduce_methods = {'sum': MPI.SUM, 'max': MPI.MAX, 'min': MPI.MIN, 'mean': MPI.SUM}
-        self._reduce_methods = reduce_methods
-
-
-    def subgrid(self, xc, yc, window):
-        """
-        Extract a subgrid from the data at a window around
-        the point (xc,yc)
-        
-        Parameters
-        ----------
-         xc      : x coordinate
-         yc      : y coordinate
-         window  : size of window in metres
-
-        Returns
-        -------
-         data    : subgrid
-        """
-        if window > self.max_window:
-            raise ValueError("Max window size is {}".format(self.max_window))
-
-        return super(CurieParallel, self).subgrid(window, xc, yc)
-
-    
-    def sync(self, vector):
-        """
-        Synchronise a local vector across all processors
-        """
-        self.lvec.setArray(vector)
-        self.dm.localToGlobal(self.lvec, self.gvec)
-        self.dm.globalToLocal(self.gvec, self.lvec)
-        return self.lvec.array.copy()
-
-
-    def add_parallel_prior(self, **kwargs):
-        """
-        Add a prior to the dictionary (tuple)
-        This version broadcasts just the values from the root processor
-
-        Available priors are beta, zt, dz, C
-
-        Usage
-        -----
-         add_parallel_prior(beta=(p, sigma_p))
-        """
-
-        comm = self.comm
-        MPI = self.MPI
-
-        for key in kwargs:
-            if key in self.prior:
-                prior = kwargs[key]
-                p  = np.array(prior[0]) # prior
-                dp = np.array(prior[1]) # uncertainty
-
-                comm.Bcast([p, MPI.DOUBLE], root=0)
-                comm.Bcast([dp,MPI.DOUBLE], root=0)
-
-                self.prior[key] = (float(p), float(dp))
-            else:
-                raise ValueError("prior must be one of {}".format(self.prior.keys()))
-
-
-    def distribute_prior(self, method, **kwargs):
-        """
-        Distribute priors across all processors using a specific method.
-
-        Parameters
-        ----------
-         method  : operation to reduce the local priors to a single value
-                 : choose from one of 'sum', 'mean', 'min', 'max'
-         kwargs  : (prior, sigma_prior) tuple
-
-        Notes
-        -----
-         add_prior will broadcast the prior on the root processor (rank=0)
-         distribute_prior enacts a MPI Allreduce operation.
-        """
-
-        comm = self.comm
-        MPI = self.MPI
-
-        if method in self._reduce_methods:
-            op = self._reduce_methods[method]
-        else:
-            raise ValueError("choose one of the following methods {}".format(self._reduce_methods.keys()))
-
-        for key in kwargs:
-            if key in self.prior:
-                prior = kwargs[key]
-                local_p  = np.array(prior[0]) # prior
-                local_dp = np.array(prior[1]) # uncertainty
-
-                global_p = np.array(0.0)
-                global_dp = np.array(0.0)
-
-                comm.Allreduce([local_p, MPI.DOUBLE], [global_p, MPI.DOUBLE], op=op)
-                comm.Allreduce([local_dp, MPI.DOUBLE], [global_dp, MPI.DOUBLE], op=op)
-                if method == 'mean':
-                    global_p  /= comm.size
-                    global_dp /= comm.size
-
-                self.prior[key] = (float(global_p), float(global_dp))
-            else:
-                raise ValueError("prior must be one of {}".format(self.prior.keys()))
\ No newline at end of file
diff --git a/pycurious/therm.py b/pycurious/therm.py
deleted file mode 100644
index ceba4d3..0000000
--- a/pycurious/therm.py
+++ /dev/null
@@ -1,377 +0,0 @@
-"""
-Copyright 2018 Ben Mather, Robert Delhaye
-
-This file is part of PyCurious.
-
-PyCurious is free software: you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation, either version 3 of the License, or any later version.
-
-PyCurious is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU Lesser General Public License for more details.
-
-You should have received a copy of the GNU Lesser General Public License
-along with PyCurious.  If not, see <http://www.gnu.org/licenses/>.
-"""
-
-import numpy as np
-from scipy import sparse
-from scipy.sparse.linalg import spsolve
-from multiprocessing import Pool, Process, Queue, cpu_count
-
-try: range=xrange
-except: pass
-
-
-class CurieTherm(object):
-    """
-    Implicit 1d solver for the steady-state heat equation
-    over a structured grid using sparse matrices.
-
-    Parameters
-    ----------
-     minZ : float, minimum Cartesian coordinate
-     maxZ : float, maximum Cartesian coordinate
-     resZ : int, resolution in the z direction
-     kwargs : dict, keyword arguments to pass to control
-       optional functionality e.g.
-       - stencil_width=int : number of node neighbours to
-           include in each matrix row
-
-    Attributes
-    ----------
-     diffusivity  : float shape(n,) thermal conductivity field
-     heat_sources : float shape(n,) heat source field
-     temperature  : float shape(n,) temperature field
-
-     npoints      : int, number of nodes in the mesh
-     sizes        : (int,int) dimensions of the sparse matrix
-     
-    """
-
-    def __init__(self, minZ, maxZ, resZ, **kwargs):
-
-        self.npoints = resZ
-        self.extent = (minZ, maxZ)
-        self.coords = np.linspace(minZ, maxZ, resZ)
-        self.nodes = np.arange(0, resZ, dtype=np.int)
-        self.sizes = (resZ, resZ)
-
-        width = kwargs.pop('stencil_width', 1)
-        self.stencil_width = 2*width + 1
-
-        closure = [(0,-2), (2,0), (1,-1)]
-        self.closure = self._create_closure_object(closure)
-
-
-        self._initialise_COO_vectors()
-        self._initialise_boundary_dictionary()
-        self.diffusivity  = np.zeros(resZ)
-        self.heat_sources = np.zeros(resZ)
-        self.temperature  = np.zeros(resZ)
-
-
-    def _initialise_COO_vectors(self):
-
-        nn = self.npoints
-
-        index = np.empty(nn + 2, dtype=np.int)
-        index.fill(-1)
-        self.index = index
-
-        self.rows = np.empty((self.stencil_width, nn), dtype=np.int)
-        self.cols = np.empty((self.stencil_width, nn), dtype=np.int)
-        self.vals = np.empty((self.stencil_width, nn))
-
-
-    def _initialise_boundary_dictionary(self):
-
-        nn = self.npoints
-        coords = self.coords
-        minZ, maxZ = self.extent
-
-        # Setup boundary dictionary
-        bc = dict()
-
-        m0 = coords == minZ
-        m1 = coords == maxZ
-        d0 = abs(coords[1] - coords[0])
-        d1 = abs(coords[-1] - coords[-2])
-
-        bc["minZ"] = {"val": 0.0, "delta": d0, "flux": True, "mask": m0}
-        bc["maxZ"] = {"val": 0.0, "delta": d1, "flux": True, "mask": m1}
-
-        self.bc = bc
-        self.dirichlet_mask = np.zeros(nn, dtype=bool)
-
-
-    def _create_closure_object(self, closure):
-
-        nn = self.npoints
-        obj = [[0] for i in range(self.stencil_width)]
-
-        for i in range(0, self.stencil_width):
-            # construct slicing object
-            start, end = closure[i-j]
-            obj[i] = slice(start, nn+end+2)
-
-        return obj
-
-
-    def update_properties(self, diffusivity, heat_sources):
-        """
-        Update diffusivity and heat sources
-        """
-        self.diffusivity[:] = diffusivity
-        self.heat_sources[:] = heat_sources
-        return
-
-
-    def boundary_condition(self, wall, val, flux=False):
-        """
-        Set the boundary conditions on each wall of the domain.
-        By default each wall is a Dirichlet condition.
-
-        Parameters
-        ----------
-         wall : str, wall to assign bc - 'minZ' or 'maxZ'
-         val  : float or array(n,) value(s) to assign to wall
-         flux : bool, toggle type of boundary condition
-           True = Neumann flux boundary condition
-           False = Dirichlet boundary condition (default)
-
-        Notes
-        -----
-         If flux=True, positive val indicates a flux vector towards the centre
-         of the domain.
-
-         val can be a vector with the same number of elements as the wall
-        """
-        wall = str(wall)
-
-        if wall in self.bc:
-            self.bc[wall]["val"]  = np.array(val, copy=True)
-            self.bc[wall]["flux"] = flux
-            d = self.bc[wall]
-
-            mask = d['mask']
-
-            if flux:
-                self.dirichlet_mask[mask] = False
-                self.bc[wall]["val"] /= -d['delta']
-            else:
-                self.dirichlet_mask[mask] = True
-
-        else:
-            raise ValueError("Wall should be one of {}".format(self.bc.keys()))
-
-
-    def construct_matrix(self):
-        """
-        Construct a sparse coefficient matrix
-        i.e. matrix A in AT = b
-
-        Notes
-        -----
-         We vectorise the 7-point stencil for fast matrix insertion.
-         An extra border of dummy values around the domain allows for automatic
-         Neumann (flux) boundary creation.
-         These are stomped on if there are any Dirichlet conditions.
-        """
-
-        nodes = self.nodes
-        nn = self.npoints
-
-        index = self.index
-
-        rows = self.rows
-        cols = self.cols
-        vals = self.vals
-
-        dirichlet_mask = self.dirichlet_mask
-        coords = self.coords.reshape(-1,1)
-
-        u = self.diffusivity
-
-        k = np.zeros(nn + 2)
-        k[1:-1] = u
-
-        for i in range(0, self.stencil_width):
-            obj = self.closure[i]
-
-            rows[i] = nodes
-            cols[i] = index[obj].ravel()
-
-            distance = np.linalg.norm(coords[cols[i]] - coords, axis=1)
-            distance[distance==0] = 1e-12 # protect against dividing by zero
-            delta = 1.0/(2.0*distance**2)
-
-            vals[i] = delta*(k[obj] + u).ravel()
-
-
-        # Dirichlet boundary conditions (duplicates are summed)
-        cols[:,dirichlet_mask] = nodes[dirichlet_mask]
-        vals[:,dirichlet_mask] = 0.0
-
-        # zero off-grid coordinates
-        vals[cols < 0] = 0.0
-
-        # centre point
-        vals[-1] = 0.0
-        vals[-1][dirichlet_mask] = -1.0
-
-
-        row = rows.ravel()
-        col = cols.ravel()
-        val = vals.ravel()
-
-
-        # mask off-grid entries and sum duplicates
-        mask = col >= 0
-        row = row[mask]
-        col = col[mask]
-        val = val[mask]
-
-        mat = sparse.coo_matrix((val, (row, col)), shape=self.sizes).tocsr()
-        mat.sum_duplicates()
-        diag = np.ravel(mat.sum(axis=1))
-        diag *= -1
-        mat.setdiag(diag)
-
-        return mat
-
-
-    def construct_rhs(self):
-        """
-        Construct the right-hand-side vector
-        i.e. vector b in AT = b
-
-        Notes
-        -----
-         Boundary conditions are grabbed from the dictionary and
-         summed to the rhs.
-         Be careful of duplicate entries on the corners!!
-        """
-        
-        rhs = -1.0*self.heat_sources
-
-        for wall in self.bc:
-            val  = self.bc[wall]['val']
-            flux = self.bc[wall]['flux']
-            mask = self.bc[wall]['mask']
-            if flux:
-                rhs[mask] += val
-            else:
-                rhs[mask] = val
-
-        return rhs
-
-
-    def solve(self, matrix=None, rhs=None):
-        """
-        Construct the matrix A and vector b in AT = b and solve for T
-        (i.e. temperature field)
-
-        Arguments
-        ---------
-         matrix : (optional) scipy sparse matrix object
-                 build using construct_matrix()
-         rhs    : (optional) numpy right-hand-side vector
-                 build using construct_rhs()
-
-        Returns
-        -------
-         sol    : 1d array shape(n,) temperature solution
-
-        Notes
-        -----
-         The solution to the system of linear equations AT = b is solved
-         with spsolve in the scipy.sparse.linalg module
-        """
-        if matrix is None:
-            matrix = self.construct_matrix()
-        if rhs is None:
-            rhs = self.construct_rhs()
-        # res = self.temperature
-
-        T = spsolve(matrix, rhs)
-        self.temperature[:] = T
-
-        return T
-
-
-    def gradient(self, vector, **kwargs):
-        """
-        Calculate gradient of a vector
-        
-        Arguments
-        ---------
-         vector : 1d array of shape(n,)
-
-        Returns
-        -------
-         dvdz   : 1d array of shape(n,)
-                 derivative of vector in z direction
-        """
-        return np.gradient(vector, self.coords, **kwargs)
-
-
-    def heatflux(self):
-        """
-        Calculate the heat flux from the conductivity field
-        and temperature solution
-
-        Returns
-        -------
-         qz  : 1d array shape(n,) heat flux vector
-        """
-
-        T = self.temperature
-        k = self.diffusivity
-        divT = self.gradient(T)
-        return -k*divT
-
-
-    def isosurface(self, vector, isoval, interp='linear'):
-        """
-        Calculate an isosurface along the z axis
-
-        Parameters
-        ----------
-         vector : array, the same size as the mesh (n,)
-         isoval : float, isosurface value
-         interp : str, interpolation method can be either
-            'nearest' - nearest neighbour interpolation
-            'linear'  - linear interpolation
-        
-        Returns
-        -------
-         z_interp : isosurface
-        """
-        coords = self.coords
-        sort_idx = np.abs(vector - isoval).argsort()    
-        i0 = sort_idx[0]
-        z0 = coords[i0]
-
-        if interp == 'linear':
-            v0 = vector[i0]
-            
-            # identify next nearest node
-            i1 = sort_idx[1]
-            z1 = coords[i1]
-            v1 = vector[i1]
-
-            vmin = min(v0,v1)
-            vmax = max(v0,v1)
-
-            ratio = np.array([isoval, vmin, vmax])
-            ratio -= ratio.min()
-            ratio /= ratio.max()
-            z_interp = ratio[0]*z1 + (1.0 - ratio[0])*z0
-            return z_interp
-        elif interp == 'nearest':
-            return z0
-        else:
-            raise ValueError("enter a valid interp method: 'linear' or 'nearest'")
diff --git a/setup.cfg b/setup.cfg
index 52985f1..7295b60 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,7 +1,12 @@
+[aliases]
+test=pytest
+
 [bdist_wheel]
+
 # This flag says that the code is written to work on both Python 2 and Python
 # 3. If at all possible, it is good practice to do this. If you cannot, you
 # will need to generate wheels for each Python version that you support.
+
 universal=1
 [metadata]
 description-file = README.md
diff --git a/setup.py b/setup.py
index 854f975..a5e3123 100644
--- a/setup.py
+++ b/setup.py
@@ -19,27 +19,27 @@
 import io
 
 ## in development set version to none and ...
-PYPI_VERSION = "0.3"
+PYPI_VERSION = "1.1"
 
 # Return the git revision as a string (from numpy)
 def git_version():
     def _minimal_ext_cmd(cmd):
         # construct minimal environment
         env = {}
-        for k in ['SYSTEMROOT', 'PATH']:
+        for k in ["SYSTEMROOT", "PATH"]:
             v = os.environ.get(k)
             if v is not None:
                 env[k] = v
         # LANGUAGE is used on win32
-        env['LANGUAGE'] = 'C'
-        env['LANG'] = 'C'
-        env['LC_ALL'] = 'C'
-        out = subprocess.Popen(cmd, stdout = subprocess.PIPE, env=env).communicate()[0]
+        env["LANGUAGE"] = "C"
+        env["LANG"] = "C"
+        env["LC_ALL"] = "C"
+        out = subprocess.Popen(cmd, stdout=subprocess.PIPE, env=env).communicate()[0]
         return out
 
     try:
-        out = _minimal_ext_cmd(['git', 'rev-parse', '--short', 'HEAD'])
-        GIT_REVISION = out.strip().decode('ascii')
+        out = _minimal_ext_cmd(["git", "rev-parse", "--short", "HEAD"])
+        GIT_REVISION = out.strip().decode("ascii")
     except OSError:
         GIT_REVISION = "Unknown"
 
@@ -50,40 +50,49 @@ def _minimal_ext_cmd(cmd):
     PYPI_VERSION = git_version()
 
 this_directory = path.abspath(path.dirname(__file__))
-with open(path.join(this_directory, 'README.md')) as f:
+with open(path.join(this_directory, "README.md")) as f:
     long_description = f.read()
 
 
-ext = Extension(name               = 'pycurious.radon',
-                sources            = ['src/radon.pyx','src/cradon.c'],
-                include_dirs       = [np.get_include(),'src'],
-                library_dirs       = ['pycurious'],
-                extra_compile_args = ['-std=c99'])
+ext = Extension(
+    name="pycurious.radon",
+    sources=["src/radon.pyx", "src/cradon.c"],
+    include_dirs=[np.get_include(), "src"],
+    library_dirs=["pycurious"],
+    extra_compile_args=["-std=c99"],
+)
 
 if __name__ == "__main__":
-    setup(name              = 'pycurious',
-          long_description  = long_description,
-          long_description_content_type = 'text/markdown',
-          author            = "Ben Mather",
-          author_email      = "brmather1@gmail.com",
-          url               = "https://github.com/brmather/pycurious",
-          version           = PYPI_VERSION,
-          description       = "Python tool for computing the Curie depth from magnetic data",
-          ext_modules       = cythonize([ext]),
-          install_requires  = ['numpy', 'scipy>=0.15.0', 'Cython>=0.25.2'],
-          python_requires   = '>=2.7, >=3.5',
-          setup_requires    = ['pytest-runner', 'webdav'],
-          tests_require     = ['pytest', 'webdav'],
-          packages          = ['pycurious'],
-          package_data      = {'pycurious': ['Examples/data/*',
-                                             'Examples/*.ipynb',
-                                             'Examples/Notebooks/Tanaka/*.ipynb',
-                                             'Examples/Notebooks/Bouligand/*.ipynb',
-                                             'Examples/Scripts/*.py']},
-          classifiers       = ['Programming Language :: Python :: 2',
-                               'Programming Language :: Python :: 2.7',
-                               'Programming Language :: Python :: 3',
-                               'Programming Language :: Python :: 3.5',
-                               'Programming Language :: Python :: 3.6',
-                               'Programming Language :: Python :: 3.7']
-          )
+    setup(
+        name="pycurious",
+        long_description=long_description,
+        long_description_content_type="text/markdown",
+        author="Ben Mather",
+        author_email="brmather1@gmail.com",
+        url="https://github.com/brmather/pycurious",
+        version=PYPI_VERSION,
+        description="Python tool for computing the Curie depth from magnetic data",
+        ext_modules=cythonize([ext]),
+        install_requires=["numpy", "scipy>=0.15.0", "Cython>=0.25.2"],
+        python_requires=">=2.7, >=3.5",
+        setup_requires=["pytest-runner", "webdav"],
+        tests_require=["pytest", "webdav"],
+        packages=["pycurious"],
+        package_data={
+            "pycurious": [
+                "Examples/data/test_mag_data.txt",
+                "Examples/*.ipynb",
+                "Examples/Notebooks/Tanaka/*.ipynb",
+                "Examples/Notebooks/Bouligand/*.ipynb",
+                "Examples/Scripts/*.py",
+            ]
+        },
+        classifiers=[
+            "Programming Language :: Python :: 2",
+            "Programming Language :: Python :: 2.7",
+            "Programming Language :: Python :: 3",
+            "Programming Language :: Python :: 3.5",
+            "Programming Language :: Python :: 3.6",
+            "Programming Language :: Python :: 3.7",
+        ],
+    )


Test patch:
diff --git a/tests/Bouligand_forward.py b/tests/Bouligand_forward.py
index d78b039..f398b12 100644
--- a/tests/Bouligand_forward.py
+++ b/tests/Bouligand_forward.py
@@ -1,156 +1,162 @@
 ## Code to generate test magnetic data for PyCurious examples
-## Converted from MATLAB code by R Delhaye. 
-#Original code:
-#Claire Bouligand
-#03-19-2007
-#05-09-2014
-#Modified to use less memory
+## Converted from MATLAB code by R Delhaye.
+# Original code:
+# Claire Bouligand
+# 03-19-2007
+# 05-09-2014
+# Modified to use less memory
 
 import numpy as np
 import pycurious
 import time
 
-start_time=time.time()
-cm=1e-7
+start_time = time.time()
+cm = 1e-7
 
 ##########################################################
-#A VERIFIER: multiplier facteur cm par 4*pi
+# A VERIFIER: multiplier facteur cm par 4*pi
 ## TO VERIFY: multiplication factor "cm for 4*pi"
 ###########################################################
 ## Degree of fractal
-b=-3.0
+b = -3.0
 ## Dimension of pixels (km), should be multiples of the smallest
-dlx=1.0
-dly=1.0
-dlz=1.0
-dl=min([dlx,dly,dlz])
-dx=dlx/dl
-dy=dly/dl
-dz=dlz/dl
-
-#nombre de pixels
+dlx = 1.0
+dly = 1.0
+dlz = 1.0
+dl = min([dlx, dly, dlz])
+dx = dlx / dl
+dy = dly / dl
+dz = dlz / dl
+
+# nombre de pixels
 ## Number of pixels.
-nx=305
-ny=305
-nz=10
-nmax=int(max([nx*dx,ny*dy,nz*dz]))
+nx = 305
+ny = 305
+nz = 10
+nmax = int(max([nx * dx, ny * dy, nz * dz]))
 
-#texte="dimensions x : %d / y : %d / z : %d" % (nx*dlx,ny*dly,nz*dlz)
-#print(texte)
-#disp(texte)
+# texte="dimensions x : %d / y : %d / z : %d" % (nx*dlx,ny*dly,nz*dlz)
+# print(texte)
+# disp(texte)
 
-#Parametres de l'aimantation
+# Parametres de l'aimantation
 ## Parameters of magnetisation
-#M1=1         # aimantation moyenne = 1 A/m ## Average magnetisation
-#M2=10^(0.25) # std de l'aimantation = 10^(0.25) A/m ## Standard dev of
-#magnetisation
-M1=0
-M2=0.2
-#prof sommet zt # Depth Zt
-zt=0.305
-#prof base zb ## Depth Zb
-zb=zt+nz*dlz
+# M1=1         # aimantation moyenne = 1 A/m ## Average magnetisation
+# M2=10^(0.25) # std de l'aimantation = 10^(0.25) A/m ## Standard dev of
+# magnetisation
+M1 = 0
+M2 = 0.2
+# prof sommet zt # Depth Zt
+zt = 0.305
+# prof base zb ## Depth Zb
+zb = zt + nz * dlz
 
 ############################################################
-#GENERATION DE LA DISTRIB 3D D'AIMANTATION FRACTAL 
-#A PARTIR VOLUME DE DIM ID. (nmax,dl) SELON X, Y ET Z
+# GENERATION DE LA DISTRIB 3D D'AIMANTATION FRACTAL
+# A PARTIR VOLUME DE DIM ID. (nmax,dl) SELON X, Y ET Z
 
 if __name__ == "__main__":
 
-    M=np.zeros((nmax,nmax,nmax))
-    l=[0,(nmax-1)*dl,dl]
-
-    #x=np.arange(0,dlx*(nx-1),dlx)
-    #y=np.arange(0,dly*(nx-1),dly)
-    #z=np.arange(0,dlz*(nx-1),dlz)
-
-    df=1.0/(nmax*dl)
-    dfx=1.0/(nx*dlx)
-    dfy=1.0/(ny*dly)
-    dfz=1.0/(nz*dlz)
-
-    #Center indices of matrix = Nyquist
-    c=(nmax+1)/2
-    f=np.zeros((int(np.floor((nmax+1))),1))
-    fx=np.zeros((int(np.floor((nmax+1))),1))
-    fy=np.zeros((int(np.floor((nmax+1))),1))
-    fz=np.zeros((int(np.floor((nmax+1))),1))
-    for i in range(0,int(np.floor((nmax+1)/2))-1):
-        f[i]=(i)*df
-    for i in range(int(np.floor((nmax+1)/2))-1,int(nmax)):
-        f[i]=-1.0*(nmax-i)*df
-    for i in range(0,int(np.floor((nx+1)/2))-1):
-        fx[i]=(i)*dfx
-    for i in range(int(np.floor((nx+1)/2))-1,int(nx)):
-        fx[i]=-1.0*(nx-i)*dfx
-    for i in range(1,int(np.floor((ny+1)/2))-1):
-        fy[i]=(i)*dfy
-    for i in range(int(np.floor((ny+1)/2))-1,int(ny)):
-        fy[i]=-1.0*(ny-i)*dfy
-    for i in range(1,int(np.floor((nz+1)/2))-1):
-        fz[i]=(i)*dfz
-    for i in range(int(np.floor((nz+1)/2))-1,int(nz)):
-        fz[i]=-1.0*(nz-i)*dfz
+    M = np.zeros((nmax, nmax, nmax))
+    l = [0, (nmax - 1) * dl, dl]
+
+    # x=np.arange(0,dlx*(nx-1),dlx)
+    # y=np.arange(0,dly*(nx-1),dly)
+    # z=np.arange(0,dlz*(nx-1),dlz)
+
+    df = 1.0 / (nmax * dl)
+    dfx = 1.0 / (nx * dlx)
+    dfy = 1.0 / (ny * dly)
+    dfz = 1.0 / (nz * dlz)
+
+    # Center indices of matrix = Nyquist
+    c = (nmax + 1) / 2
+    f = np.zeros((int(np.floor((nmax + 1))), 1))
+    fx = np.zeros((int(np.floor((nmax + 1))), 1))
+    fy = np.zeros((int(np.floor((nmax + 1))), 1))
+    fz = np.zeros((int(np.floor((nmax + 1))), 1))
+    for i in range(0, int(np.floor((nmax + 1) / 2)) - 1):
+        f[i] = (i) * df
+    for i in range(int(np.floor((nmax + 1) / 2)) - 1, int(nmax)):
+        f[i] = -1.0 * (nmax - i) * df
+    for i in range(0, int(np.floor((nx + 1) / 2)) - 1):
+        fx[i] = (i) * dfx
+    for i in range(int(np.floor((nx + 1) / 2)) - 1, int(nx)):
+        fx[i] = -1.0 * (nx - i) * dfx
+    for i in range(1, int(np.floor((ny + 1) / 2)) - 1):
+        fy[i] = (i) * dfy
+    for i in range(int(np.floor((ny + 1) / 2)) - 1, int(ny)):
+        fy[i] = -1.0 * (ny - i) * dfy
+    for i in range(1, int(np.floor((nz + 1) / 2)) - 1):
+        fz[i] = (i) * dfz
+    for i in range(int(np.floor((nz + 1) / 2)) - 1, int(nz)):
+        fz[i] = -1.0 * (nz - i) * dfz
 
     # Normal distribution
-    M=M2*np.random.randn(nmax,nmax,nmax)+M1
-
-    Mf=np.fft.fftn(M)
-    for i in range(0,nmax-1):
-        for k in range(0,nmax-1):
-            for l in range(0,nmax-1):
-                Mf[i,k,l]=Mf[i,k,l]*((f[i])**2+(f[k])**2+(f[l])**2+0.0000001)**(b/4)
-    ## DC correction - not needed in python, req. in MATLAB        
+    M = M2 * np.random.randn(nmax, nmax, nmax) + M1
+
+    Mf = np.fft.fftn(M)
+    for i in range(0, nmax - 1):
+        for k in range(0, nmax - 1):
+            for l in range(0, nmax - 1):
+                Mf[i, k, l] = Mf[i, k, l] * (
+                    (f[i]) ** 2 + (f[k]) ** 2 + (f[l]) ** 2 + 0.0000001
+                ) ** (b / 4)
+    ## DC correction - not needed in python, req. in MATLAB
     ##Mf(1,1,1)=mean([Mf(2,1,1) Mf(1,2,1) Mf(1,1,2)])
 
-    Mi=np.fft.ifftn(Mf)
-    TFANO=np.zeros((nx,ny), dtype=np.complex)
-    ANO=np.zeros((nx,ny), dtype=np.complex)
-    #input("Press Enter to continue.")
+    Mi = np.fft.ifftn(Mf)
+    TFANO = np.zeros((nx, ny), dtype=np.complex)
+    ANO = np.zeros((nx, ny), dtype=np.complex)
+    # input("Press Enter to continue.")
     ## executes fine up to here.
-    for k in range(0,nz):
-        TFHM=np.fft.fft2(Mi[:,:,k])
-        z1=zt+k*dlz
-        z2=z1+dlz
-    #    disp(z1)
-    #    disp(z2)
-        for i in range(0,nx):
-            for j in range(0,ny):
-                fH=float(np.sqrt((fx[i])**2+(fy[j])**2))
-                TFANO[i,j]=TFHM[i,j]*2*np.pi*cm*(np.exp(-2*np.pi*fH*z1)-np.exp(-2*np.pi*fH*z2))
-        ANO=ANO+np.fft.ifft2(TFANO)
+    for k in range(0, nz):
+        TFHM = np.fft.fft2(Mi[:, :, k])
+        z1 = zt + k * dlz
+        z2 = z1 + dlz
+        #    disp(z1)
+        #    disp(z2)
+        for i in range(0, nx):
+            for j in range(0, ny):
+                fH = float(np.sqrt((fx[i]) ** 2 + (fy[j]) ** 2))
+                TFANO[i, j] = (
+                    TFHM[i, j]
+                    * 2
+                    * np.pi
+                    * cm
+                    * (np.exp(-2 * np.pi * fH * z1) - np.exp(-2 * np.pi * fH * z2))
+                )
+        ANO = ANO + np.fft.ifft2(TFANO)
 
     print("--- %s seconds ---" % (time.time() - start_time))
 
-    x=np.arange(0.5,0.5+dlx*(nx),dlx)*1000
-    y=np.arange(0.5,0.5+dly*(nx),dly)*1000
-    fid=open("test_mag_data.txt","w+")
-    for i in range(0,nx):
-        for j in range(0,ny):
-            fid.write(' '.join([str(x[i]),str(y[j]),str(np.real(ANO[i,j])),"\n"]))
-            
+    x = np.arange(0.5, 0.5 + dlx * (nx), dlx) * 1000
+    y = np.arange(0.5, 0.5 + dly * (nx), dly) * 1000
+    fid = open("test_mag_data.txt", "w+")
+    for i in range(0, nx):
+        for j in range(0, ny):
+            fid.write(" ".join([str(x[i]), str(y[j]), str(np.real(ANO[i, j])), "\n"]))
+
     fid.close()
-    #plt.pcolor(np.real(ANO))
-    #plt.show()
-    #clear TFHM TFANO
+    # plt.pcolor(np.real(ANO))
+    # plt.show()
+    # clear TFHM TFANO
 
-    #figure (7)
+    # figure (7)
     ## RD - PROBLEM! "Data inputs must be real". Trying abs(), or real()
     ##pcolor(ANO)
-    #pcolor(log10(abs(ANO)))
-    #shading flat
-    #colorbar
+    # pcolor(log10(abs(ANO)))
+    # shading flat
+    # colorbar
 
     ##Sauvegarde carte :
-    #ncol=nx  #Number of columns
-    #nrow=ny  #Number of rows
-    #xleft=x(1) #x coordinate for left hand corner of grid
-    #dx=dlx    #delta x
-    #yleft=y(1) #y coordinate for left hand corner of grid
-    #dy=dly    #delta y
-    #gridC=ANO
+    # ncol=nx  #Number of columns
+    # nrow=ny  #Number of rows
+    # xleft=x(1) #x coordinate for left hand corner of grid
+    # dx=dlx    #delta x
+    # yleft=y(1) #y coordinate for left hand corner of grid
+    # dy=dly    #delta y
+    # gridC=ANO
 
     ##save carte_synth_format.mat ncol nrow xleft yleft dx dy gridC
-    ##save carte_synth.mat b M1 M2 dlx dly dlz zt zb MM2 ANO 
-
-
+    ##save carte_synth.mat b M1 M2 dlx dly dlz zt zb MM2 ANO
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..c818161
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,37 @@
+import pytest
+import pycurious
+import numpy as np
+
+
+@pytest.fixture(scope="module")
+def load_magnetic_anomaly():
+    # load magnetic anomaly - i.e. random fractal noise
+    try:
+        mag_data = np.loadtxt("tests/test_mag_data.txt")
+    except:
+        mag_data = np.loadtxt("test_mag_data.txt")
+
+    nx, ny = 305, 305
+
+    x = mag_data[:, 0]
+    y = mag_data[:, 1]
+    d = mag_data[:, 2].reshape(ny, nx)
+
+    max_window = 300e3
+
+    xmin, xmax = x.min(), x.max()
+    ymin, ymax = y.min(), y.max()
+
+    xc = 0.5 * (xmin + xmax)
+    yc = 0.5 * (ymin + ymax)
+
+    # store inside dictionary
+    mag_dict = {
+        "extent": [xmin, xmax, ymin, ymax],
+        "mag_data": d,
+        "xc": xc,
+        "yc": yc,
+        "max_window": max_window,
+    }
+
+    return mag_dict
diff --git a/tests/test_0_imports.py b/tests/test_0_imports.py
index bdd6081..39423ed 100644
--- a/tests/test_0_imports.py
+++ b/tests/test_0_imports.py
@@ -2,17 +2,24 @@
 
 ## ===================
 
+
 def test_numpy_import():
     import numpy
+
     return
 
+
 def test_scipy_import():
     import scipy
+
     print("\t\t You have scipy version {}".format(scipy.__version__))
 
+
 def test_cython_import():
-	import Cython
-	return
+    import Cython
+
+    return
+
 
 def test_pycurious_modules():
     import pycurious
@@ -40,4 +47,4 @@ def test_pycurious_modules():
     test_numpy_import()
     test_scipy_import()
     test_cython_import()
-    test_pycurious_modules()
\ No newline at end of file
+    test_pycurious_modules()
diff --git a/tests/test_grid.py b/tests/test_grid.py
new file mode 100644
index 0000000..4e7b93e
--- /dev/null
+++ b/tests/test_grid.py
@@ -0,0 +1,106 @@
+import pytest
+import pycurious
+import numpy as np
+
+from conftest import load_magnetic_anomaly
+
+
+def test_subgrid(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+
+    grid = pycurious.CurieGrid(d, xmin, xmax, ymin, ymax)
+
+    window_size = 100e3
+    subgrid = grid.subgrid(window_size, xc, yc)
+
+    error_msg = "FAILED! Subgrid is of shape {} and domain is of shape {}".format(
+        subgrid.shape, grid.data.shape
+    )
+    assert subgrid.shape[0] < grid.data.shape[0], error_msg
+    assert subgrid.shape[1] < grid.data.shape[1], error_msg
+
+
+def test_FFT(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+
+    grid = pycurious.CurieGrid(d, xmin, xmax, ymin, ymax)
+
+    # Take Fourier transform
+    k, Phi, sigma_Phi = grid.radial_spectrum(grid.data, taper=None)
+
+    # radial power spectrum should decrease with wavenumber
+    # divide Phi into three sections
+    i3 = len(k) // 3
+    Phi1 = Phi[:i3]
+    Phi2 = Phi[i3 : 2 * i3]
+    Phi3 = Phi[2 * i3 :]
+
+    # also sigma_Phi should decrease with wavenumber
+    sigma_Phi1 = sigma_Phi[:i3]
+    sigma_Phi2 = sigma_Phi[i3 : 2 * i3]
+    sigma_Phi3 = sigma_Phi[2 * i3 :]
+
+    error_msg = "FAILED! Fast Fourier Transform did not produce a valid power spectrum"
+    assert Phi1.mean() > Phi2.mean() > Phi3.mean(), error_msg
+    assert sigma_Phi1.mean() > sigma_Phi2.mean() > sigma_Phi3.mean(), error_msg
+
+
+def test_taper_functions(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+
+    grid = pycurious.CurieGrid(d, xmin, xmax, ymin, ymax)
+
+    # Take Fourier transform using different taper functions
+    k, Phi1, sigma_Phi1 = grid.radial_spectrum(grid.data, taper=None)
+    k, Phi2, sigma_Phi2 = grid.radial_spectrum(grid.data, taper=np.hanning)
+    k, Phi3, sigma_Phi3 = grid.radial_spectrum(grid.data, taper=np.hamming)
+
+    grad_Phi1 = np.gradient(Phi1, k)
+    grad_Phi2 = np.gradient(Phi2, k)
+    grad_Phi3 = np.gradient(Phi3, k)
+
+    assert (
+        Phi1.mean() > Phi2.mean()
+    ), "FAILED! 'taper=np.hanning' not significantly different from 'taper=None'"
+    assert (
+        Phi1.mean() > Phi3.mean()
+    ), "FAILED! 'taper=np.hamming' not significantly different from 'taper=None'"
+    assert (
+        grad_Phi1.mean() > grad_Phi2.mean()
+    ), "FAILED! 'taper=np.hanning' has steeper gradient from 'taper=None'"
+    assert (
+        grad_Phi1.mean() > grad_Phi3.mean()
+    ), "FAILED! 'taper=np.hamming' has steeper gradient from 'taper=None'"
+
+
+def test_Tanaka(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+
+    grid = pycurious.CurieGrid(d, xmin, xmax, ymin, ymax)
+
+    # wavenumber bands for Z0 and Zt, respectively
+    kwin_Z0 = (0.005, 0.03)
+    kwin_Zt = (0.03, 0.7)
+
+    k, Phi, sigma_Phi = grid.radial_spectrum(grid.data, taper=np.hanning, power=0.5)
+    (Ztr, btr, dZtr), (Zor, bor, dZor) = pycurious.tanaka1999(
+        k, Phi, sigma_Phi, kwin_Z0, kwin_Zt
+    )
+    Zb, eZb = pycurious.ComputeTanaka(Ztr, dZtr, Zor, dZor)
+
+    error_msg = "FAILED! Tanaka CPD is {:.4f} different from expected, uncertainty is {:.4f}".format(
+        Zb - 10.0, eZb
+    )
+    assert np.abs(Zb - 10.0) < 2.0 and eZb < Zb, error_msg
diff --git a/tests/test_mag_data.txt b/tests/test_mag_data.txt
new file mode 120000
index 0000000..13d9149
--- /dev/null
+++ b/tests/test_mag_data.txt
@@ -0,0 +1,1 @@
+../pycurious/Examples/data/test_mag_data.txt
\ No newline at end of file
diff --git a/tests/test_optimise.py b/tests/test_optimise.py
index 5139e73..712816a 100644
--- a/tests/test_optimise.py
+++ b/tests/test_optimise.py
@@ -1,94 +1,208 @@
-
 import pytest
 import pycurious
 import numpy as np
+import numpy.testing as npt
 from scipy.optimize import minimize
 
-# global variables
-max_window = 100e3
+from conftest import load_magnetic_anomaly
+
+
+def test_optimisation(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+    max_window = load_magnetic_anomaly["max_window"]
+
+    grid = pycurious.CurieOptimise(d, xmin, xmax, ymin, ymax)
+    beta, zt, dz, C = grid.optimise(max_window, xc, yc, taper=np.hanning)
 
-# square grid 100km x 100km
-xmin = 0.0
-xmax = 100e3
-ymin = 0.0
-ymax = 100e3
+    x_opt = np.array([beta, zt, dz])
 
-nx = 100
-ny = 100
+    # hard-coded parameters used to generate the magnetic anomaly
+    zt0 = 0.305
+    dz0 = 10.0 + zt0
+    beta0 = 3.0
+
+    x0 = np.array([beta0, zt0, dz0])
 
-# random noise
-mag_grid = np.random.randn(ny,nx)
+    # compare if they are close or not
+    # some parameters should be more similar than others
+    tol = np.array([0.3, 0.1, 2.0])
 
-xc = 50e3
-yc = 50e3
+    parameters = ["beta", "zt", "dz"]
+    err_msg = "FAILED! {} = {:.4f} is not within an acceptable tolerance of {}"
 
+    for i in range(x0.size):
+        npt.assert_allclose(
+            x_opt[i],
+            x0[i],
+            atol=tol[i],
+            err_msg=err_msg.format(parameters[i], x_opt[i], tol[i]),
+        )
 
-def test_optimisation():
-    grid = pycurious.CurieOptimise(mag_grid, xmin, xmax, ymin, ymax)
-    beta, zt, dz, C = grid.optimise(max_window, xc, yc)
 
-def test_priors():
-    grid = pycurious.CurieOptimise(mag_grid, xmin, xmax, ymin, ymax)
+def test_priors(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+    max_window = load_magnetic_anomaly["max_window"]
+
+    grid = pycurious.CurieOptimise(d, xmin, xmax, ymin, ymax)
     beta0, zt0, dz0, C0 = grid.optimise(max_window, xc, yc)
 
     grid.add_prior(beta=(1.0, 0.1))
     beta1, zt1, dz1, C1 = grid.optimise(max_window, xc, yc)
 
-    if abs(beta1 - 1.0) < abs(beta0 - 1.0):
-        print("PASSED! Optimisation is sensitive to prior information")
-    else:
-        assert False, "FAILED! Optimisation with priors failed"
-
-def test_valid_numbers():
-    S = np.array([ 22.16409774,  19.95258494,  18.27873722,  17.10575637,\
-                   16.53959747,  16.31539575,  15.69619005,  15.29953307,\
-                   14.83475976,  14.54031396,  14.33361716,  13.81764026,\
-                   13.5176055 ,  13.27386563,  13.03493328,  12.88581369,\
-                   12.61998358,  12.48616749,  12.11261083,  12.13079154,\
-                   11.85440661,  11.79244826,  11.66823202,  11.40231744,\
-                   11.32521296,  11.13634007,  11.10650999,  10.94822598,\
-                   10.78032794,  10.66593304,  10.55815845,  10.56805594,\
-                   10.33514462,  10.22026537,  10.22945756,  10.09275259,\
-                   10.11562101,   9.85061009,   9.87165772,   9.85976847,\
-                    9.73954992,   9.72021054,   9.52959744,   9.59582531,\
-                    9.50927273,   9.44691364,   9.39293966,   9.33097387,   9.33191784])
-
-    k = np.array([ 0.09237068,  0.15443902,  0.21486191,  0.282031  ,  0.33452175,\
-                   0.4020083 ,  0.46381582,  0.5290359 ,  0.59611689,  0.65875047,\
-                   0.71442694,  0.77929408,  0.84320843,  0.90939865,  0.96863721,\
-                   1.03155349,  1.09130906,  1.15699532,  1.22338685,  1.28109263,\
-                   1.34606572,  1.40861721,  1.46920824,  1.53421046,  1.59675833,\
-                   1.66109235,  1.7233503 ,  1.78492941,  1.84835837,  1.90933735,\
-                   1.973275  ,  2.03650008,  2.0983359 ,  2.16386249,  2.22416603,\
-                   2.2876692 ,  2.34841931,  2.41225232,  2.47774604,  2.5379226 ,\
-                   2.59980401,  2.66346692,  2.72735801,  2.79146076,  2.85220715,\
-                   2.91443307,  2.97643404,  3.04008782,  3.10369815])
+    assert abs(beta1 - 1.0) < abs(
+        beta0 - 1.0
+    ), "FAILED! Optimisation with priors failed"
+
+
+def test_valid_numbers(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+    max_window = load_magnetic_anomaly["max_window"]
+
+    # create phoney power spectrum
+    S = np.array(
+        [
+            22.16409774,
+            19.95258494,
+            18.27873722,
+            17.10575637,
+            16.53959747,
+            16.31539575,
+            15.69619005,
+            15.29953307,
+            14.83475976,
+            14.54031396,
+            14.33361716,
+            13.81764026,
+            13.5176055,
+            13.27386563,
+            13.03493328,
+            12.88581369,
+            12.61998358,
+            12.48616749,
+            12.11261083,
+            12.13079154,
+            11.85440661,
+            11.79244826,
+            11.66823202,
+            11.40231744,
+            11.32521296,
+            11.13634007,
+            11.10650999,
+            10.94822598,
+            10.78032794,
+            10.66593304,
+            10.55815845,
+            10.56805594,
+            10.33514462,
+            10.22026537,
+            10.22945756,
+            10.09275259,
+            10.11562101,
+            9.85061009,
+            9.87165772,
+            9.85976847,
+            9.73954992,
+            9.72021054,
+            9.52959744,
+            9.59582531,
+            9.50927273,
+            9.44691364,
+            9.39293966,
+            9.33097387,
+            9.33191784,
+        ]
+    )
+
+    k = np.array(
+        [
+            0.09237068,
+            0.15443902,
+            0.21486191,
+            0.282031,
+            0.33452175,
+            0.4020083,
+            0.46381582,
+            0.5290359,
+            0.59611689,
+            0.65875047,
+            0.71442694,
+            0.77929408,
+            0.84320843,
+            0.90939865,
+            0.96863721,
+            1.03155349,
+            1.09130906,
+            1.15699532,
+            1.22338685,
+            1.28109263,
+            1.34606572,
+            1.40861721,
+            1.46920824,
+            1.53421046,
+            1.59675833,
+            1.66109235,
+            1.7233503,
+            1.78492941,
+            1.84835837,
+            1.90933735,
+            1.973275,
+            2.03650008,
+            2.0983359,
+            2.16386249,
+            2.22416603,
+            2.2876692,
+            2.34841931,
+            2.41225232,
+            2.47774604,
+            2.5379226,
+            2.59980401,
+            2.66346692,
+            2.72735801,
+            2.79146076,
+            2.85220715,
+            2.91443307,
+            2.97643404,
+            3.04008782,
+            3.10369815,
+        ]
+    )
 
     sigma_S = np.ones_like(S)
 
-    grid = pycurious.CurieOptimise(mag_grid, xmin, xmax, ymin, ymax)
+    grid = pycurious.CurieOptimise(d, xmin, xmax, ymin, ymax)
 
     beta0 = 3.0
     zt0 = 1.0
     dz0 = 40.0
     C0 = 5.0
     x0 = np.array([beta0, zt0, dz0, C0])
-    
+
     lower_bound = np.zeros_like(x0)
-    upper_bound = [None]*len(lower_bound)
+    upper_bound = [None] * len(lower_bound)
 
     # xi = func([5, 0., 45., 6.], S, k)
     xi = grid.min_func(x0, k, S, sigma_S)
 
-    res = minimize(grid.min_func, x0, args=(k, S, sigma_S), method='TNC', bounds=list(zip(lower_bound,upper_bound)))
+    res = minimize(
+        grid.min_func,
+        x0,
+        args=(k, S, sigma_S),
+        method="TNC",
+        bounds=list(zip(lower_bound, upper_bound)),
+    )
     print("beta={:.2f}, zt={:.2f}, dz={:.2f}, C={:.2f}".format(*res.x))
 
-    if np.isfinite(res.x).all():
-        print("PASSED! All numbers are finite")
-    else:
-        assert False, "FAILED! Not all numbers are finite"
+    parameters = ["beta", "zt", "dz", "C"]
+    err_msg = "FAILED! {} = {} is not a finite number"
 
-if __name__ == "__main__":
-    test_optimisation()
-    test_priors()
-    test_valid_numbers()
\ No newline at end of file
+    for i in range(res.x.size):
+        assert np.isfinite(res.x[i]), err_msg.format(parameters[i], res.x[i])
diff --git a/tests/test_routines.py b/tests/test_routines.py
index 4f2479e..3622045 100644
--- a/tests/test_routines.py
+++ b/tests/test_routines.py
@@ -3,62 +3,57 @@
 import numpy as np
 from time import time
 
-# global variables
-max_window = 100e3
-
-# square grid 100km x 100km
-xmin = 0.0
-xmax = 100e3
-ymin = 0.0
-ymax = 100e3
-
-nx = 100
-ny = 100
-
-# random noise
-mag_grid = np.random.randn(ny,nx)
+from conftest import load_magnetic_anomaly
 
 
 def time_routine(routine, *args, **kwargs):
     t = time()
     routine(*args, **kwargs)
     str_fmt = "- {:30} completed in  {:1.6f} seconds"
-    print(str_fmt.format(str(routine.__name__), time()-t))
+    print(str_fmt.format(str(routine.__name__), time() - t))
+    return True
 
 
 ##
 ## Test PyCurious routines
 ##
 
-def test_CurieGrid_routines():
-    cpd = pycurious.CurieGrid(mag_grid, xmin, xmax, ymin, ymax)
-    xc = 50e3
-    yc = 50e3
-
-    time_routine(cpd.__init__, mag_grid, xmin, xmax, ymin, ymax)
-    time_routine(cpd.subgrid, max_window, xc, yc)
-    time_routine(cpd.create_centroid_list, 0.1*max_window)
-    time_routine(cpd.remove_trend_linear, cpd.data)
-    time_routine(cpd.radial_spectrum, cpd.data)
-    time_routine(cpd.radial_spectrum_log, cpd.data)
-    time_routine(cpd.azimuthal_spectrum, cpd.data)
-    time_routine(cpd.reduce_to_pole, cpd.data, 30., 30.)
-    time_routine(cpd.upward_continuation, cpd.data, 10e3)
 
-
-def test_CurieOptimise_routines():
-    cpd = pycurious.CurieOptimise(mag_grid, xmin, xmax, ymin, ymax)
-    xc_list, yc_list = cpd.create_centroid_list(0.5*max_window, spacingX=5e3, spacingY=5e3)
+def test_CurieGrid_routines(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+    max_window = load_magnetic_anomaly["max_window"]
+
+    cpd = pycurious.CurieGrid(d, xmin, xmax, ymin, ymax)
+
+    assert time_routine(cpd.__init__, d, xmin, xmax, ymin, ymax)
+    assert time_routine(cpd.subgrid, max_window, xc, yc)
+    assert time_routine(cpd.create_centroid_list, 0.1 * max_window)
+    assert time_routine(cpd.remove_trend_linear, cpd.data)
+    assert time_routine(cpd.radial_spectrum, cpd.data)
+    assert time_routine(cpd.azimuthal_spectrum, cpd.data)
+    assert time_routine(cpd.reduce_to_pole, cpd.data, 30.0, 30.0)
+    assert time_routine(cpd.upward_continuation, cpd.data, 10e3)
+
+
+def test_CurieOptimise_routines(load_magnetic_anomaly):
+    d = load_magnetic_anomaly["mag_data"]
+    xc = load_magnetic_anomaly["xc"]
+    yc = load_magnetic_anomaly["yc"]
+    xmin, xmax, ymin, ymax = load_magnetic_anomaly["extent"]
+    max_window = load_magnetic_anomaly["max_window"]
+
+    cpd = pycurious.CurieOptimise(d, xmin, xmax, ymin, ymax)
+    xc_list, yc_list = cpd.create_centroid_list(
+        0.5 * max_window, spacingX=5e3, spacingY=5e3
+    )
     xc = xc_list[0]
     yc = yc_list[0]
 
-    time_routine(cpd.add_prior, beta=(3.0, 0.5))
-    time_routine(cpd.optimise, 0.5*max_window, xc, yc)
-    time_routine(cpd.optimise_routine, 0.5*max_window, xc_list, yc_list)
-    time_routine(cpd.metropolis_hastings, 0.5*max_window, xc, yc, 100, 10)
-    time_routine(cpd.sensitivity, 0.5*max_window, xc, yc, 100)
-
-
-if __name__ == "__main__":
-    test_CurieGrid_routines()
-    test_CurieOptimise_routines()
\ No newline at end of file
+    assert time_routine(cpd.add_prior, beta=(3.0, 0.5))
+    assert time_routine(cpd.optimise, 0.5 * max_window, xc, yc)
+    assert time_routine(cpd.optimise_routine, 0.5 * max_window, xc_list, yc_list)
+    assert time_routine(cpd.metropolis_hastings, 0.5 * max_window, xc, yc, 100, 10)
+    assert time_routine(cpd.sensitivity, 0.5 * max_window, xc, yc, 100)