Repository: SHTOOLS/SHTOOLS (language: python)

Issue:
setup.py instead of Makefile for Python package
The main reason is the possibility of using true pythonic and recommended  way to distribute (pip), install, remove and upgrade package since all routines are fully wrapped now.



Discussion / hints:
Agree, this needs to be done! A main objective of the next update is to make everything more 'pythonic'.

Great! Does it also means support of Python 3?

That's the other thing that needs to be done! :)

I've already done (or I think so) this with the current 'master' branch and it is not hard with 2to3 tools if there is no need for Python 2. In other case it is also possible, of course. I can help with this, if required.

This would be great! It should go to the develop branch though which already has a few other additions for Python (a class interface). I don't know what the best way is to organize this, maybe you make your own fork of the project and then somehow push to the main repository?

For info, here is our planned dev schedule.

https://github.com/SHTOOLS/SHTOOLS/wiki/Development-schedule

I will be writing doc and man pages for the next month, so any help would be appreciated as I'm not going to get anything else done... I'd suggest using the develop branch if you want to contribute.

Yes, I think to fork is the best way. The main question is: do we need Python 2? 

I think we need to support 2 and 3, as almost everyone used version 2 now. What if we just add some flag to the setup file for --with-python3 ?

In the Makefile it can be flag, yes. In the setup.py file there is no need for this, the version will depend on user.

Here is a good introduction about the stuff that we need to add:

https://packaging.python.org/en/latest/distributing/#configuring-your-project


Fix patch:
diff --git a/AUTHORS.md b/AUTHORS.md
index 166b73784..8c8012211 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -15,4 +15,7 @@
 * Akihisa Hattori / [HattoriAkihisa](https://github.com/HattoriAkihisa)
 * Aaryaman Vasishta / [jammm](https://github.com/jammm)
 * Alex Kalinin / [alkalinin](https://github.com/alkalinin)
-* Armin Corbin / [rainbowsend](https://github.com/rainbowsend)
\ No newline at end of file
+* Armin Corbin / [rainbowsend](https://github.com/rainbowsend)
+* Tristan Brugere
+* Erik Schnetter / [eschnett](https://github.com/eschnett)
+* Martin Reinecke / [mreineck](https://github.com/mreineck)
diff --git a/docs/Gemfile.lock b/docs/Gemfile.lock
index 0c6b4ba85..8082fce5f 100644
--- a/docs/Gemfile.lock
+++ b/docs/Gemfile.lock
@@ -1,13 +1,13 @@
 GEM
   remote: https://rubygems.org/
   specs:
-    activesupport (6.0.3.4)
+    activesupport (6.0.4)
       concurrent-ruby (~> 1.0, >= 1.0.2)
       i18n (>= 0.7, < 2)
       minitest (~> 5.1)
       tzinfo (~> 1.1)
       zeitwerk (~> 2.2, >= 2.2.2)
-    addressable (2.7.0)
+    addressable (2.8.0)
       public_suffix (>= 2.0.2, < 5.0)
     coffee-script (2.4.1)
       coffee-script-source
@@ -16,26 +16,40 @@ GEM
     colorator (1.1.0)
     commonmarker (0.17.13)
       ruby-enum (~> 0.5)
-    concurrent-ruby (1.1.8)
-    dnsruby (1.61.5)
+    concurrent-ruby (1.1.9)
+    dnsruby (1.61.7)
       simpleidn (~> 0.1)
     em-websocket (0.5.2)
       eventmachine (>= 0.12.9)
       http_parser.rb (~> 0.6.0)
-    ethon (0.12.0)
-      ffi (>= 1.3.0)
+    ethon (0.14.0)
+      ffi (>= 1.15.0)
     eventmachine (1.2.7)
-    execjs (2.7.0)
-    faraday (1.3.0)
+    execjs (2.8.1)
+    faraday (1.7.0)
+      faraday-em_http (~> 1.0)
+      faraday-em_synchrony (~> 1.0)
+      faraday-excon (~> 1.1)
+      faraday-httpclient (~> 1.0.1)
       faraday-net_http (~> 1.0)
+      faraday-net_http_persistent (~> 1.1)
+      faraday-patron (~> 1.0)
+      faraday-rack (~> 1.0)
       multipart-post (>= 1.2, < 3)
-      ruby2_keywords
+      ruby2_keywords (>= 0.0.4)
+    faraday-em_http (1.0.0)
+    faraday-em_synchrony (1.0.0)
+    faraday-excon (1.1.0)
+    faraday-httpclient (1.0.1)
     faraday-net_http (1.0.1)
-    ffi (1.14.2)
+    faraday-net_http_persistent (1.2.0)
+    faraday-patron (1.0.0)
+    faraday-rack (1.0.0)
+    ffi (1.15.3)
     forwardable-extended (2.6.0)
     gemoji (3.0.1)
-    github-pages (209)
-      github-pages-health-check (= 1.16.1)
+    github-pages (218)
+      github-pages-health-check (= 1.17.2)
       jekyll (= 3.9.0)
       jekyll-avatar (= 0.7.0)
       jekyll-coffeescript (= 1.1.1)
@@ -50,39 +64,39 @@ GEM
       jekyll-readme-index (= 0.3.0)
       jekyll-redirect-from (= 0.16.0)
       jekyll-relative-links (= 0.6.1)
-      jekyll-remote-theme (= 0.4.2)
+      jekyll-remote-theme (= 0.4.3)
       jekyll-sass-converter (= 1.5.2)
-      jekyll-seo-tag (= 2.6.1)
+      jekyll-seo-tag (= 2.7.1)
       jekyll-sitemap (= 1.4.0)
       jekyll-swiss (= 1.0.0)
-      jekyll-theme-architect (= 0.1.1)
-      jekyll-theme-cayman (= 0.1.1)
-      jekyll-theme-dinky (= 0.1.1)
-      jekyll-theme-hacker (= 0.1.2)
-      jekyll-theme-leap-day (= 0.1.1)
-      jekyll-theme-merlot (= 0.1.1)
-      jekyll-theme-midnight (= 0.1.1)
-      jekyll-theme-minimal (= 0.1.1)
-      jekyll-theme-modernist (= 0.1.1)
-      jekyll-theme-primer (= 0.5.4)
-      jekyll-theme-slate (= 0.1.1)
-      jekyll-theme-tactile (= 0.1.1)
-      jekyll-theme-time-machine (= 0.1.1)
+      jekyll-theme-architect (= 0.2.0)
+      jekyll-theme-cayman (= 0.2.0)
+      jekyll-theme-dinky (= 0.2.0)
+      jekyll-theme-hacker (= 0.2.0)
+      jekyll-theme-leap-day (= 0.2.0)
+      jekyll-theme-merlot (= 0.2.0)
+      jekyll-theme-midnight (= 0.2.0)
+      jekyll-theme-minimal (= 0.2.0)
+      jekyll-theme-modernist (= 0.2.0)
+      jekyll-theme-primer (= 0.6.0)
+      jekyll-theme-slate (= 0.2.0)
+      jekyll-theme-tactile (= 0.2.0)
+      jekyll-theme-time-machine (= 0.2.0)
       jekyll-titles-from-headings (= 0.5.3)
       jemoji (= 0.12.0)
-      kramdown (= 2.3.0)
+      kramdown (= 2.3.1)
       kramdown-parser-gfm (= 1.1.0)
       liquid (= 4.0.3)
       mercenary (~> 0.3)
       minima (= 2.5.1)
       nokogiri (>= 1.10.4, < 2.0)
-      rouge (= 3.23.0)
+      rouge (= 3.26.0)
       terminal-table (~> 1.4)
-    github-pages-health-check (1.16.1)
+    github-pages-health-check (1.17.2)
       addressable (~> 2.3)
       dnsruby (~> 1.60)
       octokit (~> 4.0)
-      public_suffix (~> 3.0)
+      public_suffix (>= 2.0.2, < 5.0)
       typhoeus (~> 1.3)
     html-pipeline (2.14.0)
       activesupport (>= 2)
@@ -136,57 +150,57 @@ GEM
       jekyll (>= 3.3, < 5.0)
     jekyll-relative-links (0.6.1)
       jekyll (>= 3.3, < 5.0)
-    jekyll-remote-theme (0.4.2)
+    jekyll-remote-theme (0.4.3)
       addressable (~> 2.0)
       jekyll (>= 3.5, < 5.0)
       jekyll-sass-converter (>= 1.0, <= 3.0.0, != 2.0.0)
       rubyzip (>= 1.3.0, < 3.0)
     jekyll-sass-converter (1.5.2)
       sass (~> 3.4)
-    jekyll-seo-tag (2.6.1)
-      jekyll (>= 3.3, < 5.0)
+    jekyll-seo-tag (2.7.1)
+      jekyll (>= 3.8, < 5.0)
     jekyll-sitemap (1.4.0)
       jekyll (>= 3.7, < 5.0)
     jekyll-swiss (1.0.0)
-    jekyll-theme-architect (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-architect (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-cayman (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-cayman (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-dinky (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-dinky (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-hacker (0.1.2)
+    jekyll-theme-hacker (0.2.0)
       jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-leap-day (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-leap-day (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-merlot (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-merlot (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-midnight (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-midnight (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-minimal (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-minimal (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-modernist (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-modernist (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-primer (0.5.4)
+    jekyll-theme-primer (0.6.0)
       jekyll (> 3.5, < 5.0)
       jekyll-github-metadata (~> 2.9)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-slate (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-slate (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-tactile (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-tactile (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
-    jekyll-theme-time-machine (0.1.1)
-      jekyll (~> 3.5)
+    jekyll-theme-time-machine (0.2.0)
+      jekyll (> 3.5, < 5.0)
       jekyll-seo-tag (~> 2.0)
     jekyll-titles-from-headings (0.5.3)
       jekyll (>= 3.3, < 5.0)
@@ -196,12 +210,12 @@ GEM
       gemoji (~> 3.0)
       html-pipeline (~> 2.2)
       jekyll (>= 3.0, < 5.0)
-    kramdown (2.3.0)
+    kramdown (2.3.1)
       rexml
     kramdown-parser-gfm (1.1.0)
       kramdown (~> 2.0)
     liquid (4.0.3)
-    listen (3.4.1)
+    listen (3.6.0)
       rb-fsevent (~> 0.10, >= 0.10.3)
       rb-inotify (~> 0.9, >= 0.9.10)
     mercenary (0.3.6)
@@ -209,26 +223,26 @@ GEM
       jekyll (>= 3.5, < 5.0)
       jekyll-feed (~> 0.9)
       jekyll-seo-tag (~> 2.1)
-    minitest (5.14.3)
+    minitest (5.14.4)
     multipart-post (2.1.1)
-    nokogiri (1.11.1-x86_64-darwin)
+    nokogiri (1.12.3-x86_64-darwin)
       racc (~> 1.4)
-    octokit (4.20.0)
+    octokit (4.21.0)
       faraday (>= 0.9)
       sawyer (~> 0.8.0, >= 0.5.3)
     pathutil (0.16.2)
       forwardable-extended (~> 2.6)
-    public_suffix (3.1.1)
+    public_suffix (4.0.6)
     racc (1.5.2)
-    rb-fsevent (0.10.4)
+    rb-fsevent (0.11.0)
     rb-inotify (0.10.1)
       ffi (~> 1.0)
-    rexml (3.2.4)
-    rouge (3.23.0)
-    ruby-enum (0.8.0)
+    rexml (3.2.5)
+    rouge (3.26.0)
+    ruby-enum (0.9.0)
       i18n
-    ruby2_keywords (0.0.4)
-    rubyzip (2.3.0)
+    ruby2_keywords (0.0.5)
+    rubyzip (2.3.2)
     safe_yaml (1.0.5)
     sass (3.7.4)
       sass-listen (~> 4.0.0)
@@ -254,10 +268,10 @@ GEM
     zeitwerk (2.4.2)
 
 PLATFORMS
-  x86_64-darwin-19
+  x86_64-darwin-20
 
 DEPENDENCIES
   github-pages
 
 BUNDLED WITH
-   2.2.5
+   2.2.26
diff --git a/docs/_data/sidebars/mydoc_sidebar.yml b/docs/_data/sidebars/mydoc_sidebar.yml
index 765a22b12..817ae98d3 100755
--- a/docs/_data/sidebars/mydoc_sidebar.yml
+++ b/docs/_data/sidebars/mydoc_sidebar.yml
@@ -64,8 +64,8 @@ entries:
       url: /file-formats.html
       output: web, pdf
 
-    - title: Implementation details
-      url: /implementation-details.html
+    - title: Backends
+      url: /backends.html
       output: web, pdf
 
   - title: Python Reference
@@ -104,6 +104,10 @@ entries:
       url: /python-utilities.html
       output: web, pdf
 
+    - title: pysh.backends
+      url: /python-backends.html
+      output: web, pdf
+
   - title: Support
     output: web, pdf
 
diff --git a/docs/index.md b/docs/index.md
index 7fcb3a7c5..62961aafb 100755
--- a/docs/index.md
+++ b/docs/index.md
@@ -17,11 +17,13 @@ toc: false
 
 * Use of both regularly sampled geographic grids and grids appropriate for Gauss-Legendre quadrature.
 
-* Spherical harmonic transforms proven to be accurate up to about degree 2800.
+* Spherical harmonic transforms proven to be fast and accurate.
+
+* Supports multiple backends for performing basic mathematical operations, including [shtools](index-fortran.html) and the [Distinctly Useful Code Collection (DUCC)](https://mtr.pages.mpcdf.de/ducc/).
 
 * Perform localized multitaper spectral analyses, or expand functions in terms of localized Slepian bases.
 
-* Support for standard data and file formats, including *xarray* and *netcdf*.
+* Supports standard data and file formats, including *xarray* and *netcdf*.
 
 * Import research-grade gravity, topography, and magnetic field datasets with a single command.
 
diff --git a/docs/pages/fortran/fortran-contributors.md b/docs/pages/fortran/fortran-contributors.md
index 1004ea88a..d759adcb6 100644
--- a/docs/pages/fortran/fortran-contributors.md
+++ b/docs/pages/fortran/fortran-contributors.md
@@ -24,3 +24,6 @@ folder: fortran
 * Aaryaman Vasishta / [jammm](https://github.com/jammm)
 * Alex Kalinin / [alkalinin](https://github.com/alkalinin)
 * Armin Corbin / [rainbowsend](https://github.com/rainbowsend)
+* Tristan Brugere
+* Erik Schnetter / [eschnett](https://github.com/eschnett)
+* Martin Reinecke / [mreineck](https://github.com/mreineck)
\ No newline at end of file
diff --git a/docs/pages/mydoc/backends.md b/docs/pages/mydoc/backends.md
new file mode 100644
index 000000000..a9c870c78
--- /dev/null
+++ b/docs/pages/mydoc/backends.md
@@ -0,0 +1,71 @@
+---
+title: "Backends"
+keywords: spherical harmonics software package, backends, shtools, ducc
+sidebar: mydoc_sidebar
+permalink: backends.html
+summary: pyshtools supports multiple backends for performing the most computationaly intensive mathematical operations.
+toc: true
+folder: mydoc
+---
+
+<style>
+table:nth-of-type(n) {
+    display:table;
+    width:100%;
+}
+</style>
+
+
+## pyshtools backends
+
+*pyshtools* is a collection of utilities that facilitates working with data expressed in spherical harmonics. Though *pyshtools* is accessed through python, many of the mathematical operations that are called are computationally demanding and are instead perfomed by external libraries. In order to access the functions of these external libraries, an intermediary python wrapper function is used that converts the python call and parameters into the form expected by the library. In principle, any pre-existing numerical software package could be used as a numerical backend to *pyshtools* by the creation of such wrapper functions.
+
+The following software packages are supported as numerical backends to *pysthools*:
+
+| Class name | Description |
+| ---------- | ----------- |
+| [shtools](index-fortran.html) (default) | Spherical Harmonic Tools (Fortran 95, installed as part of pyshtools) |
+| [ducc](https://mtr.pages.mpcdf.de/ducc/) | Distinctly Useful Code Collection (C++17) |
+
+The *shtools* package (written in Fortran 95) suppports all the functionality of *pyshtools*. Other packages only implement a subset of the *shtools* routines, and whenever an implementation of a function does not exist, the *shtools* function is used as a fallback.
+
+There are two ways to control which backend is used. The first method is to set the preferred backend that will be used in all subsequent calls. As an example, here we select the *ducc* backend using the `select_preferred_backend()` function, and also set the number of threads to use to 4:
+```python
+In [1]: pysh.backends.select_preferred_backend(backend='ducc', nthreads=4)
+```
+The current preferred backend can be inspected by using the `preferred_backend()` function:
+```python
+In [2]: pysh.backends.preferred_backend()
+Out[2]: 'ducc'
+```
+
+Alternatively, if one has an initialized *pyshtools* class instance (such as from `SHGrid` or `SHCoeffs`), the backend can be specified for an single function call by use of the optional `backend` parameter. In this example, `grid` is expanded into spherical harmonic coefficients using the same backend parameters as above:
+```python
+In [3]: clm = grid.expand(backend='ducc', nthreads=4)
+```
+Following such a call, the preferred backend remains unchanged.
+
+
+## Supported backend key characteristics
+
+### shtools
+
+* Supports all *pyshtools* functionality, including localized spectral analyses, Slepian analyses, and operations on gravity and magnetic field data.
+* Spherical harmonic transforms and reconstructions are accurate up to degree 2800.
+* Spherical harmonic rotations are accurate to about degree 1200.
+
+### ducc
+* Spherical harmonic transforms and rotations are more than 10 times faster than the default *shtools* routines.
+* Supports the use of multiple threads to speed up computations.
+* Spherical harmonic transforms are accurate beyond degree 25,000.
+* Does not implement functions involving localized spectral analyses, Slepian analyses, nor gravity and magnetic field data.
+* The native routines make use orthonormalized spherical harmonics that exclude the Condon-Shortley phase factor: Other normalizations are supported, but require preprocessing in python.
+
+
+## Speed comparisons
+
+The speeds of the spherical harmonic transforms using the above backends were tested for both real and complex data using equidistant 'DH2' grids. The amount of time in seconds required to perform individually the forward and inverse operations is plotted in the figure below as a function of the spherical harmonic bandwidth of the function. These calculations were performed on a 2018 MacBook Pro with a 2.7 GHz Quad-Core Intel i7 processor and 16 GB of memory.
+
+{% include image.html file="backend-timing-tests.png" alt="Timing tests" caption="Figure 1. Time to perform the reconstruction of a function from its spherical harmonic coefficients (solid lines) and the spherical harmonic transform of the gridded function (dashed lines). The spherical harmonic transforms made use of orthonormalized spherical harmonics that exclude the Condon-Shortley phase factor (which is the native implementation in the 'ducc' backend)." %}
+
+For the 'shtools' backend, the transform time is on the order of one second for degrees close to 700 and just under a minute for degree 2800. The 'ducc' backend has transform speeds that are about an order of magnitude faster than the 'shtools' backend. Transforms are about 1 second for degrees near 2000, and between 1 and 3 minutes for degrees near 10,000. In general, the complex routines are slower than then real routines by a factor close to 2.
diff --git a/docs/pages/mydoc/contributors.md b/docs/pages/mydoc/contributors.md
index cd653776c..9f5373259 100644
--- a/docs/pages/mydoc/contributors.md
+++ b/docs/pages/mydoc/contributors.md
@@ -23,4 +23,7 @@ folder: mydoc
 * Akihisa Hattori / [HattoriAkihisa](https://github.com/HattoriAkihisa)
 * Aaryaman Vasishta / [jammm](https://github.com/jammm)
 * Alex Kalinin / [alkalinin](https://github.com/alkalinin)
-* Armin Corbin / [rainbowsend](https://github.com/rainbowsend)
\ No newline at end of file
+* Armin Corbin / [rainbowsend](https://github.com/rainbowsend)
+* Tristan Brugere
+* Erik Schnetter / [eschnett](https://github.com/eschnett)
+* Martin Reinecke / [mreineck](https://github.com/mreineck)
\ No newline at end of file
diff --git a/docs/pages/mydoc/implementation-details.md b/docs/pages/mydoc/implementation-details.md
deleted file mode 100644
index 56f37362a..000000000
--- a/docs/pages/mydoc/implementation-details.md
+++ /dev/null
@@ -1,117 +0,0 @@
----
-title: "Implementation details"
-keywords: spherical harmonics software package, spherical harmonic transform, legendre functions, multitaper spectral analysis, fortran, Python, gravity, magnetic field
-sidebar: mydoc_sidebar
-permalink: implementation-details.html
-summary: The spherical harmonic transforms in pyshtools make use of integrations over longitude that involve fast Fourier transforms and integrations over latitude that utilize either Gauss-Legendre quadrature or exact quadrature rules for regularly spaced grids. The transforms and reconstructions are accurate to about degree 2800.
-toc: true
-folder: mydoc
----
-
-<style>
-table:nth-of-type(n) {
-    display:table;
-    width:100%;
-}
-</style>
-
-## Numerical computations
-
-### Spherical harmonic reconstructions
-
-If the spherical harmonic coefficients of a function $$f$$ were known up to a maximum degree $$L$$, the function could be evaluated at arbitrary colatitude $$\theta$$ and longitude $$\phi$$ using the equation
-
-\begin{equation}
-f\left(\theta,\phi\right) = \sum_{l=0}^{L} \sum_{m=-l}^l f_{lm} \, Y_{lm}\left(\theta,\phi \right).
-\end{equation}
-
-Using the separate variables $$C_{lm}$$ and $$S_{lm}$$ for the cosine and sine spherical harmonic coefficients, respectively, and after interchanging the order of summations over $$l$$ and $$m$$, the function can be expressed equivalently as
-
-$$ \begin{equation}
-f\left(\theta,\phi\right) = \sum_{m=0}^L \sum_{l=m}^{L} \left( C_{lm} \, \bar{P}_{lm}(\cos \theta) \, \cos m\phi + S_{lm} \, \bar{P}_{lm}(\cos \theta) \, \sin m\phi\right).
-\end{equation} $$
-
-Defining the two component vector
-
-$$ \begin{equation}
-\left[a_m(\theta), b_m(\theta)\right] = \sum_{l=m}^{L} \left[C_{lm}, S_{lm}\right] \, \bar{P}_{lm}(\cos \theta),
-\end{equation} $$
-
-The function $$f$$ can be written more simply as
-
-$$ \begin{equation}
-f\left(\theta,\phi\right) = \sum_{m=0}^L \left( a_m(\theta) \, \cos m\phi + b_m(\theta) \, \sin m\phi\right).
-\end{equation} $$
-
-For a given latitude band, the function $$f$$ can thus be evaluated on a series of grid nodes all at the same time using an inverse fast Fourier transform. For this operation, SHTOOLS makes use of the highly optimized software package [FFTW](http://www.fftw.org) [Frigo and Johnson 2005] that supports grids of arbitrary length and both real and complex data. To adequately sample the function in longitude, a minimum of $$2L+1$$ data points should be employed.
-
-The coefficients $$a_m$$ and $$b_m$$ depend upon the associated Legendre functions, which are calculated efficiently using standard three-term recursion relations over adjacent spherical harmonic degrees. For a given colatitude, the sectoral term $$\bar{P}_{mm}$$ is first calculated using an analytic equation, and then $$\bar{P}_{lm}$$ is calculated for all values of $$l>m$$. To circumvent numerical underflows near the poles for large values of $$m$$, the associated Legendre functions are calculated using the approach of Holmes and Featherstone [2002], where the sectoral $$\bar{P}_{mm}$$ terms are multiplied by $$10^{280} \sin^m \theta$$ prior to performing the recursions, and then appropriately unscaled at the end of the recursion. This ensures that the Legendre functions are accurate to about degree 2800.
-
-### Spherical harmonic transforms
-
-The spherical harmonic coefficients of a function can be calculated from the relation
-
-$$ \begin{equation}
-\left[C_{lm}, S_{lm}\right] = \frac{1}{4\pi} \int_{0}^{2\pi} \int_0^{\pi} f(\theta,\phi) \, \bar{P}_{lm}(\theta) \left[\cos m\phi, \sin m\phi\right]\, \sin\theta \, d\theta \, d\phi.
-\end{equation} $$
-
-Defining the two intermediary variables $$c$$ and $$s$$
-
-$$ \begin{equation}
-\left[c_{lm}^{(i)}, s_{lm}^{(i)}\right] = \int_0^{2\pi} f(\theta_i,\phi) \left[ \cos m\phi, \sin m\phi\right] d\phi,
-\end{equation} $$
-
-it is seen that for a given colatitude $$\theta_i$$ and degree $$l$$, all of the
-angular orders can be calculated at once by making use of a fast Fourier
-transform of the function $$f$$. Replacing the integral over latitude with a numerical quadrature rule, the spherical harmonic
-coefficients can be calculated as
-
-$$ \begin{equation}
-\left[C_{lm}, S_{lm}\right] = \frac{1}{4 \pi} \sum_{i=1}^N w_i \, \bar{P}_{lm}(\cos \theta_i) \left[c_{lm}^{(i)}, s_{lm}^{(i)}\right].
-\label{eq:quadrature}
-\end{equation} $$
-
-Here, $$w$$ is the latitudinal weight, and $$N$$ is the number of latitudinal
-points over which the integration is performed. It is possible to choose the weights $$w_i$$ and the locations of the latitudinal sampling points $$\theta_i$$ such that the quadrature in equation \eqref{eq:quadrature} is exact.
-
-As is described in more detail [here](grid-formats.html), SHTOOLS makes use of grid formats that accommodate exact quadrature.These include regularly spaced and sampled grids that satisfy the *Driscoll and Healy* (1994) sampling theorem, and grids for exact quadrature using Gauss-Legendre quadrature [e.g., Press et al. 1992]. For all grids, the quadrature is exact only when the function being integrated is a terminating polynomial. The functions $$c_{lm}$$ and $$s_{lm}$$ can be approximated as polynomials of maximum degree $$L$$, and when multiplied by the associated Legendre function, the integrand is approximately a polynomial of maximum degree $$2L.$$
-
-## Accuracy
-
-To test the accuracy of the spherical harmonic transforms and reconstructions,
-two sets of synthetic spherical harmonic coefficients were created. Each
-coefficient was chosen to be a random Gaussian distributed number with unit
-variance, and the coefficients were then scaled such that the power spectrum
-was proportional to either $$l^2$$ or $$l^{-2}$$. For a given spherical harmonic bandwidth $$L$$, the function was reconstructed in the space domain using the Gauss-Legendre quadrature implementation and then re-expanded into spherical harmonics. The maximum and root-mean square (rms) relative errors between the initial and final set of coefficients were then computed as a function of spherical harmonic degree, as plotted in the figure below.
-
-{% include image.html file="relative-errors.png" alt="Relative errors of the spherical harmonic coefficients" caption="Figure 1. Maximum and rms relative errors of the spherical harmonic coefficients as a function of spherical harmonic bandwidth following a reconstruction of the function in the space domain and a subsequent spherical harmonic transform." %}
-
-The errors associated with the transform and inverse pair increase in an quasi-exponential manner, with the maximum relative error being approximately 1 part in a billion for degrees close to 400, and about 1 part in
-a million for degrees close to 2600. Though the errors are negligible to about
-degree 2600, they then grow somewhat between degrees 2700 and 2800. The rms
-errors for the coefficients as a function of the bandwidth $$L$$ are typically
-three orders of magnitude smaller than the maximum relative errors.
-
-The relative errors are nearly the same for the two tested power spectra, implying that the form of the data do not strongly affect the accuracy of the routines. Relative errors using the *Driscoll and Healy* [1994] sampled grids are nearly identical to those using Gauss-Legendre quadrature. The errors associated with the routines for complex data are lower by a few orders of magnitude.
-
-## Timing tests
-
-The speed of the spherical harmonic reconstructions and transforms were tested for both real and complex data using the Gauss-Legendre and
-*Driscoll and Healy* [1994] quadrature implementations. The amount of time in seconds required to perform these operations is plotted in the figure below as
-a function of the spherical harmonic bandwidth of the function. These
-calculations were performed on a modern Mac Pro 2.7 GHz 12-Core Intel Xeon E5
-using 64 bit executables and level 3 optimizations.
-
-{% include image.html file="timing-tests.png" alt="Timing tests" caption="Figure 2. Time to perform the reconstruction of a function from its spherical harmonic coefficients (solid lines) and the spherical harmonic transform of the function (dashed lines)." %}
-
-For the real Gauss-Legendre quadrature routines, the transform time is on the order of one second for degrees close to 800 and about 30 seconds for degree 2600. For the real *Driscoll and Healy* [1994] sampled grids, the transform time is close to a second for degree 600 and about 1 minute for degrees close to 2600. The complex routines are slower by a factor of about 1.4.
-
-## References
-
-* Driscoll, J. R. and D. M. Healy, Computing Fourier transforms and convolutions on the 2-sphere, Adv. Appl. Math., 15, 202-250, doi:[10.1006/aama.1994.1008](https://doi.org/10.1006/aama.1994.1008), 1994.
-
-* Frigo, M., and S. G. Johnson, The design and implementation of FFTW3, Proc. IEEE, 93, 216-231, doi:[10.1109/JPROC.2004.840301](https://doi.org/10.1109/JPROC.2004.840301), 2005.
-
-* Holmes, S. A., and W. E. Featherstone, A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions, J. Geodesy, 76, 279-299, doi:[10.1007/s00190-002-0216-2](https://doi.org/10.1007/s00190-002-0216-2), 2002.
-
-* Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, "Numerical Recipes in FORTRAN: The Art of Scientific Computing," 2nd ed., Cambridge Univ. Press, Cambridge, UK, 1992.
diff --git a/docs/pages/mydoc/pydoc/backend_module.md b/docs/pages/mydoc/pydoc/backend_module.md
new file mode 100644
index 000000000..1bbbf8c49
--- /dev/null
+++ b/docs/pages/mydoc/pydoc/backend_module.md
@@ -0,0 +1,30 @@
+---
+title: backend_module()
+keywords: spherical harmonics software package, spherical harmonic transform, legendre functions, multitaper spectral analysis, fortran, Python, gravity, magnetic field
+sidebar: mydoc_sidebar
+permalink: backend_module.html
+summary:
+tags: [python]
+toc: false
+editdoc: pydoc
+---
+
+Return the specified backend module used for the spherical harmonic
+transforms in pyshtools.
+
+## Usage
+
+```python
+module = backend_module([backend, nthreads])
+```
+
+## Parameters
+
+**backend : str, optional, default = preferred_backend()**
+:   Name of the preferred backend, either 'shtools' or 'ducc'.
+
+**nthreads : int, optional, default = 1**
+:   Number of threads to use for the 'ducc' backend. Setting this parameter
+        to 0 will use as many threads as there are hardware threads on the
+        system.
+    
\ No newline at end of file
diff --git a/docs/pages/mydoc/pydoc/preferred_backend.md b/docs/pages/mydoc/pydoc/preferred_backend.md
new file mode 100644
index 000000000..7001d9f93
--- /dev/null
+++ b/docs/pages/mydoc/pydoc/preferred_backend.md
@@ -0,0 +1,20 @@
+---
+title: preferred_backend()
+keywords: spherical harmonics software package, spherical harmonic transform, legendre functions, multitaper spectral analysis, fortran, Python, gravity, magnetic field
+sidebar: mydoc_sidebar
+permalink: preferred_backend.html
+summary:
+tags: [python]
+toc: false
+editdoc: pydoc
+---
+
+Return the name of the current preferred backend used for the spherical
+harmonic transforms in pyshtools as a string.
+
+## Usage
+
+```python
+preferred_backend()
+```
+    
\ No newline at end of file
diff --git a/docs/pages/mydoc/pydoc/preferred_backend_module.md b/docs/pages/mydoc/pydoc/preferred_backend_module.md
new file mode 100644
index 000000000..4c1cb7e95
--- /dev/null
+++ b/docs/pages/mydoc/pydoc/preferred_backend_module.md
@@ -0,0 +1,20 @@
+---
+title: preferred_backend_module()
+keywords: spherical harmonics software package, spherical harmonic transform, legendre functions, multitaper spectral analysis, fortran, Python, gravity, magnetic field
+sidebar: mydoc_sidebar
+permalink: preferred_backend_module.html
+summary:
+tags: [python]
+toc: false
+editdoc: pydoc
+---
+
+Return a reference to the current preferred backend module used for the
+spherical harmonic transforms in pyshtools.
+
+## Usage
+
+```python
+module = preferred_backend_module()
+```
+    
\ No newline at end of file
diff --git a/docs/pages/mydoc/pydoc/pymakegradientdh.md b/docs/pages/mydoc/pydoc/pymakegradientdh.md
index 43f2c9760..1f851b77b 100644
--- a/docs/pages/mydoc/pydoc/pymakegradientdh.md
+++ b/docs/pages/mydoc/pydoc/pymakegradientdh.md
@@ -14,7 +14,7 @@ Compute the gradient of a scalar function and return grids of the two horizontal
 ## Usage
 
 ```python
-`theta`, `phi` = MakeGradientDH (`cilm`, [`lmax`, `sampling`, `lmax_calc`, `extend`])
+`theta`, `phi` = MakeGradientDH (`cilm`, [`lmax`, `sampling`, `lmax_calc`, `extend`, `radius`])
 ```
 
 ## Returns
@@ -39,16 +39,20 @@ Compute the gradient of a scalar function and return grids of the two horizontal
 `lmax_calc` : optional, integer, default = `lmax`
 :   The maximum spherical harmonic degree used in evaluating the functions. This must be less than or equal to `lmax`.
 
-`extend` : input, optional, bool, default = False
+`extend` : optional, bool, default = False
 :   If True, compute the longitudinal band for 360 E and the latitudinal band for 90 S. This increases each of the dimensions of `griddh` by 1.
 
+`radius` : optional, float, default = 1.0
+:   The radius of the sphere used when computing the gradient of the function.
+
+
 ## Description
 
 `MakeGradientDH` will compute the horizontal gradient of a scalar function on a sphere defined by the spherical harmonic coefficients `cilm`. The output grids of the theta and phi components of the gradient are either equally sampled (`n` by `n`) or equally spaced (`n` by 2`n`) in latitude and longitude. The gradient is given by the formula
 
 `Grad F = 1/r dF/theta theta-hat + 1/(r sin theta) dF/dphi phi-hat`.
 
-where theta is colatitude and phi is longitude. The radius r is taken from the degree zero coefficient of the input function.
+where theta is colatitude and phi is longitude. The radius r is by default set to 1, but this can be modified by use of the optional parameter `radius`.
 
 The default is to use an input grid that is equally sampled (`n` by `n`), but this can be changed to use an equally spaced grid (`n` by 2`n`) by the optional argument `sampling`. The redundant longitudinal band for 360 E and the latitudinal band for 90 S are excluded by default, but these can be computed by specifying the optional argument `extend`.
 
diff --git a/docs/pages/mydoc/pydoc/pyshread.md b/docs/pages/mydoc/pydoc/pyshread.md
index 312e6074f..cbde74507 100644
--- a/docs/pages/mydoc/pydoc/pyshread.md
+++ b/docs/pages/mydoc/pydoc/pyshread.md
@@ -15,8 +15,9 @@ Read shtools-formatted spherical harmonic coefficients from a text file.
 
 ```python
 coeffs, [errors], lmaxout, [header], [header2] = shread(
-    filename, [error=True, header=True, header2=True, lmax, skip])
+    filename, [error=True, header=True, header2=True, lmax, skip,
 ```
+        encoding])
 
 ## Returns
 
@@ -64,6 +65,9 @@ coeffs, [errors], lmaxout, [header], [header2] = shread(
 **skip : int, optional, default = 0**
 :   The number of lines to skip before parsing the file.
 
+**encoding : str, optional, default = None**
+:   Encoding of the input file. The default is to use the system default.
+
 ## Notes
 
 This function will read spherical harmonic coefficients from an
diff --git a/docs/pages/mydoc/pydoc/read_dov.md b/docs/pages/mydoc/pydoc/read_dov.md
index b750a30d7..82a8bfc7c 100644
--- a/docs/pages/mydoc/pydoc/read_dov.md
+++ b/docs/pages/mydoc/pydoc/read_dov.md
@@ -16,8 +16,9 @@ Read spherical harmonic coefficients from a text file formatted as
 
 ```python
 coeffs, [errors], lmaxout, [header], [header2] = read_dov(
-    filename, [error=True, header=True, header2=True, lmax, skip])
+    filename, [error=True, header=True, header2=True, lmax, skip,
 ```
+        encoding])
 
 ## Returns
 
@@ -65,6 +66,9 @@ coeffs, [errors], lmaxout, [header], [header2] = read_dov(
 **skip : int, optional, default = 0**
 :   The number of lines to skip before parsing the file.
 
+**encoding : str, optional, default = None**
+:   Encoding of the input file. The default is to use the system default.
+
 ## Notes
 
 This function will read spherical harmonic coefficients from a 'dov'-
diff --git a/docs/pages/mydoc/pydoc/read_icgem_gfc.md b/docs/pages/mydoc/pydoc/read_icgem_gfc.md
index d62588d00..2e9389df4 100644
--- a/docs/pages/mydoc/pydoc/read_icgem_gfc.md
+++ b/docs/pages/mydoc/pydoc/read_icgem_gfc.md
@@ -16,8 +16,9 @@ file.
 
 ```python
 cilm, gm, r0, [errors] = read_icgem_gfc(filename,
-    [errors, lmax, epoch, encoding)
+    [errors, lmax, epoch, encoding,
 ```
+                                             quiet)
 
 ## Returns
 
@@ -56,7 +57,40 @@ cilm, gm, r0, [errors] = read_icgem_gfc(filename,
         If the format of the file is 'icgem2.0' then the epoch must be
         specified.
 
-**encoding : str, optional**
-:   Encoding of the input file. Try to use 'iso-8859-1' if the default
-        (UTF-8) fails.
+**encoding : str, optional, default = None**
+:   Encoding of the input file. The default is to use the system default.
+        If the default encoding doesn't work, try 'iso-8859-1'.
+
+**quiet : bool, default = False**
+:   If True, suppress warnings about undefined keywords when reading the
+        file.
+
+## Notes
+
+This routine reads ICGEM formatted files of gravitational potential models
+and outputs arrays of the gravitational potential coefficients, errors, GM,
+and the reference radius. If epoch is specified, the coefficients will make
+use of the time variable terms in order to compute and return the potential
+coefficients for the specified epoch. Otherwise, the coefficients will be
+returned for the reference epoch of the model.
+
+Valid keys in the header section include:
+        modelname (not used)
+        product_type (only 'gravity_field' is allowed)
+        earth_gravity_constant or gravity_constant
+        radius
+        max_degree
+        errors ('unknown', 'formal', 'calibrated' or 'calibrated_and_formal')
+        tide_system (not used)
+        norm (not used)
+        format (either None or 'icgem2.0')
+
+Valid keys in the data section include:
+        gfc
+        gfct
+        trnd or dot
+        asin
+        acos
+
+Data lines starting with an unknown key are ignored.
     
\ No newline at end of file
diff --git a/docs/pages/mydoc/pydoc/read_igrf.md b/docs/pages/mydoc/pydoc/read_igrf.md
index f53594d42..4dcc7a504 100644
--- a/docs/pages/mydoc/pydoc/read_igrf.md
+++ b/docs/pages/mydoc/pydoc/read_igrf.md
@@ -34,6 +34,9 @@ read_igrm(filename, [year])
 **year : float, optional, default = 2020.**
 :   The year to compute the coefficients.
 
+**encoding : str, optional, default = None**
+:   Encoding of the input file. The default is to use the system default.
+
 ## Notes
 
 The current International Geomagnetic Reference Field (IGRF-13) is a
diff --git a/docs/pages/mydoc/pydoc/select_preferred_backend.md b/docs/pages/mydoc/pydoc/select_preferred_backend.md
new file mode 100644
index 000000000..149e62d92
--- /dev/null
+++ b/docs/pages/mydoc/pydoc/select_preferred_backend.md
@@ -0,0 +1,30 @@
+---
+title: select_preferred_backend()
+keywords: spherical harmonics software package, spherical harmonic transform, legendre functions, multitaper spectral analysis, fortran, Python, gravity, magnetic field
+sidebar: mydoc_sidebar
+permalink: select_preferred_backend.html
+summary:
+tags: [python]
+toc: false
+editdoc: pydoc
+---
+
+Select the preferred backend module used for the spherical harmonic
+transforms in pyshtools.
+
+## Usage
+
+```python
+select_preferred_backend_module([backend, nthreads])
+```
+
+## Parameters
+
+**backend : str, optional, default = 'shtools'**
+:   Name of the preferred backend, either 'shtools' or 'ducc'.
+
+**nthreads : int, optional, default = 1**
+:   Number of threads to use for the 'ducc' backend. Setting this parameter
+        to 0 will use as many threads as there are hardware threads on the
+        system.
+    
\ No newline at end of file
diff --git a/docs/pages/mydoc/pydoc/shwrite.md b/docs/pages/mydoc/pydoc/shwrite.md
index 8ef763865..bdf563827 100644
--- a/docs/pages/mydoc/pydoc/shwrite.md
+++ b/docs/pages/mydoc/pydoc/shwrite.md
@@ -14,7 +14,7 @@ Write shtools-formatted spherical harmonic coefficients to a text file.
 ## Usage
 
 ```python
-shwrite(filename, coeffs, [errors, header, header2, lmax])
+shwrite(filename, coeffs, [errors, header, header2, lmax, encoding])
 ```
 
 ## Parameters
@@ -41,6 +41,9 @@ shwrite(filename, coeffs, [errors, header, header2, lmax])
 **lmax : int, optional, default = None**
 :   The maximum spherical harmonic degree to write to the file.
 
+**encoding : str, optional, default = None**
+:   Encoding of the output file. The default is to use the system default.
+
 ## Notes
 
 This function will write spherical harmonic coefficients (and optionally
diff --git a/docs/pages/mydoc/pydoc/write_dov.md b/docs/pages/mydoc/pydoc/write_dov.md
index 7848a2bec..b76f048c0 100644
--- a/docs/pages/mydoc/pydoc/write_dov.md
+++ b/docs/pages/mydoc/pydoc/write_dov.md
@@ -15,7 +15,7 @@ Write spherical harmonic coefficients to a text file formatted as
 ## Usage
 
 ```python
-write_dov(filename, coeffs, [errors, header, header2, lmax])
+write_dov(filename, coeffs, [errors, header, header2, lmax, encoding])
 ```
 
 ## Parameters
@@ -42,6 +42,9 @@ write_dov(filename, coeffs, [errors, header, header2, lmax])
 **lmax : int, optional, default = None**
 :   The maximum spherical harmonic degree to write to the file.
 
+**encoding : str, optional, default = None**
+:   Encoding of the output file. The default is to use the system default.
+
 ## Notes
 
 This function will write spherical harmonic coefficients (and optionally
diff --git a/docs/pages/mydoc/pydoc/write_icgem_gfc.md b/docs/pages/mydoc/pydoc/write_icgem_gfc.md
index beabcb688..7a8286478 100644
--- a/docs/pages/mydoc/pydoc/write_icgem_gfc.md
+++ b/docs/pages/mydoc/pydoc/write_icgem_gfc.md
@@ -16,8 +16,9 @@ file.
 
 ```python
 write_icgem_gfc(filename, coeffs, [errors, header, lmax, modelname, gm, r0,
-    product_type, earth_gm, error_kind, tide_system, normalization, format)
+    product_type, earth_gm, error_kind, tide_system, normalization, format,
 ```
+        encoding)
 
 ## Parameters
 
@@ -55,7 +56,7 @@ write_icgem_gfc(filename, coeffs, [errors, header, lmax, modelname, gm, r0,
 :   Reference radius of the model, in meters.
 
 **error_kind : str, optional, default = None**
-:   Which errors to read. Can be either 'unknown', 'calibrated', or
+:   Which errors to write. Can be either 'unknown', 'calibrated', or
         'formal'.
 
 **tide_system : str, optional, default = 'unknown'**
@@ -67,4 +68,7 @@ write_icgem_gfc(filename, coeffs, [errors, header, lmax, modelname, gm, r0,
 
 **format : str, optional, default = None**
 :   The format of the ICGEM spherical harmonic coefficients.
+
+**encoding : str, optional, default = None**
+:   Encoding of the output file. The default is to use the system default.
     
\ No newline at end of file
diff --git a/docs/pages/mydoc/python-backends.md b/docs/pages/mydoc/python-backends.md
new file mode 100644
index 000000000..80027ab30
--- /dev/null
+++ b/docs/pages/mydoc/python-backends.md
@@ -0,0 +1,26 @@
+---
+title: pysh.backends
+keywords: pyshtools, backends, shtools, ducc
+sidebar: mydoc_sidebar
+permalink: python-backends.html
+summary: This module defines functions for accessing and setting optional parameters of the backends used for the spherical harmonic transforms in pyshtools.
+toc: false
+folder: mydoc
+---
+
+<style>
+table:nth-of-type(n) {
+    display:table;
+    width:100%;
+}
+table:nth-of-type(n) th:nth-of-type(2) {
+    width:75%;
+}
+</style>
+
+| Function name | Description |
+| ------------- | ----------- |
+| [select_preferred_backend](select_preferred_backend.html) | Set the preferred backend module and options. |
+| [preferred_backend](preferred_backend.html) | Return the name of the current preferred backend. |
+| [preferred_backend_module](preferred_backend_module.html) | Return a reference to the preferred backend module. |
+| [backend_module](backend_module.html) | Return a reference to the specified backend module. |
diff --git a/docs/pages/mydoc/python-installing.md b/docs/pages/mydoc/python-installing.md
index dfc54a556..e9959d7f4 100644
--- a/docs/pages/mydoc/python-installing.md
+++ b/docs/pages/mydoc/python-installing.md
@@ -85,10 +85,11 @@ When installing pyshtools using `pip` or `conda`, the following packages should
 When installing pyshtools using `conda`, the following will also be installed automatically:
 
 * [cartopy](https://scitools.org.uk/cartopy/docs/latest/): required for Cartopy map projections. Cartopy requires (see below) *proj*, *geos*, *cython*, *pyshp*, *six*, and *shapely*.
-* [pygmt](https://www.pygmt.org) (>=0.2): required for pygmt map projections. pygmt requires (see below) *gmt (>=6.1.1)*.
+* [pygmt](https://www.pygmt.org) (>=0.3): required for pygmt map projections. pygmt requires (see below) *gmt (>=6.1.1)*.
+* [ducc0](https://mtr.pages.mpcdf.de/ducc/) (>=0.15): required for using the 'ducc' backend for spherical harmonic transforms.
 * [palettable](https://jiffyclub.github.io/palettable/): scientific color maps required by one of the tutorials.
 
-The above three packages will need to be installed separately when installing pyshtools with `pip`, as described in the following subsections.
+The above four packages will need to be installed separately when installing pyshtools with `pip`, as described in the following subsections.
 
 ### How to install Cartopy
 
@@ -108,17 +109,26 @@ pip install cartopy
 See [these instructions](https://scitools.org.uk/cartopy/docs/latest/installing.html#installing) for further details.
 
 ### How to install pygmt
-In order to use the *pygmt* plotting routines, it will be necessary to install both *pygmt (>=0.2)* and the *gmt (>=6.1.1)* library. This is most easily achieved using conda with
+In order to use the *pygmt* plotting routines, it will be necessary to install both *pygmt (>=0.3)* and the *gmt (>=6.1.1)* library. This is most easily achieved using conda with
 ```bash
 conda install -c conda-forge pygmt gmt
 ```
 Alternatively, *pygmt* can be installed using `pip`
 ```bash
-pip install pygmt
+pip install pyshtools[pygmt]  # installs pygmt at the same time as pyshtools
+pip install pygmt  # if pyshtools is already installed
 ```
-The *gmt* library will then need to be installed using other means, such as with brew, macports or apt-get:
+For this case, the *gmt* library will then need to be installed using other means, such as with brew, macports or apt-get:
 ```bash
 brew install gmt  # using brew on macOS
 sudo port install gmt6  # using macports on macOS
 sudo apt-get install gmt  # using apt-get on linux
-```
\ No newline at end of file
+```
+
+### How to install ducc
+To make use of the 'ducc' backend for the spherical harmonic transforms, it will be necessary to install the *ducc0 (>=0.15)* package, preferably via `pip`:
+```bash
+pip install pyshtools[ducc]  # installs ducc at the same time as pyshtools
+pip install ducc0>=0.15  # if pyshtools is already installed
+```
+A conda version of *ducc0* exists, but it is not custom-built for the user's hardware and will therefore be roughly half as fast on most machines.
diff --git a/docs/pages/mydoc/using-with-python.md b/docs/pages/mydoc/using-with-python.md
index 89b319e40..803d94378 100644
--- a/docs/pages/mydoc/using-with-python.md
+++ b/docs/pages/mydoc/using-with-python.md
@@ -46,6 +46,7 @@ This will load the following classes and subpackages into the `pysh` namespace:
 | [rotate](python-spherical-harmonic-rotations.html) | Spherical harmonic rotation routines |
 | [gravmag](python-gravity-magnetics.html) | Gravity and magnetics routines |
 | [utils](python-utilities.html) | Utilities |
+| [backends](python-backends.html) | Functions to control which backend to use for the spherical harmonic transforms. |
 | shtools | All Python wrapped Fortran 95 routines |
 
 
diff --git a/environment.yml b/environment.yml
index bb4968ff6..515573bb8 100644
--- a/environment.yml
+++ b/environment.yml
@@ -19,3 +19,5 @@ dependencies:
     - pypandoc
     - flake8
     - pip
+    - pip:
+      - ducc0>=0.15
diff --git a/pyshtools/__init__.py b/pyshtools/__init__.py
index fd1ed6bdd..af6051cbe 100644
--- a/pyshtools/__init__.py
+++ b/pyshtools/__init__.py
@@ -29,6 +29,8 @@
     rotate            : Spherical harmonic rotation routines.
     gravmag           : Gravity and magnetics routines.
     utils             : pyshtools utilities.
+    backends          : Routines for selecting the spherical harmonic
+                        transform backend.
 
 For further information, consult the web documentation at
 
@@ -52,6 +54,7 @@
 from . import rotate
 from . import gravmag
 from . import utils
+from . import backends
 
 # ---- Import principal classes into pyshtools namespace
 from .shclasses import SHCoeffs
@@ -67,6 +70,6 @@
 
 # ---- Define __all__ for use with: from pyshtools import * ----
 __all__ = ['constants', 'shclasses', 'legendre', 'expand', 'shio', 'shtools',
-           'spectralanalysis', 'rotate', 'gravmag', 'utils', 'SHCoeffs',
-           'SHGrid', 'SHWindow', 'Slepian', 'SHGravCoeffs', 'SHMagCoeffs',
-           'datasets']
+           'spectralanalysis', 'rotate', 'gravmag', 'utils', 'backends',
+           'SHCoeffs', 'SHGrid', 'SHWindow', 'Slepian', 'SHGravCoeffs',
+           'SHMagCoeffs', 'datasets']
diff --git a/pyshtools/backends.py b/pyshtools/backends.py
new file mode 100644
index 000000000..2ab9b8933
--- /dev/null
+++ b/pyshtools/backends.py
@@ -0,0 +1,142 @@
+"""
+backends
+
+This module defines functions for accessing and setting optional parameters of
+the backends used for the spherical harmonic transforms in pyshtools.
+
+Supported backends
+------------------
+shtools (default)            Spherical Harmonic Tools (Fortran 95), installed
+                             as part of pyshtools.
+ducc0                        Distinctly Useful Code Collection (C++17).
+
+Functions
+---------
+select_preferred_backend()   Set the preferred backend module and options.
+preferred_backend()          Return the name of the current preferred backend.
+preferred_backend_module()   Return a reference to the preferred backend
+                             module.
+backend_module()             Return a reference to the specified backend
+                             module.
+"""
+
+# defaults
+from . import shtools as _preferred_backend_module
+_preferred_backend = "shtools"
+
+
+def preferred_backend():
+    """
+    Return the name of the current preferred backend used for the spherical
+    harmonic transforms in pyshtools as a string.
+
+    Usage
+    -----
+    preferred_backend()
+    """
+    return _preferred_backend
+
+
+def preferred_backend_module():
+    """
+    Return a reference to the current preferred backend module used for the
+    spherical harmonic transforms in pyshtools.
+
+    Usage
+    -----
+    module = preferred_backend_module()
+    """
+    return _preferred_backend_module
+
+
+def backend_module(backend=None, nthreads=None):
+    """
+    Return the specified backend module used for the spherical harmonic
+    transforms in pyshtools.
+
+    Usage
+    -----
+    module = backend_module([backend, nthreads])
+
+    Parameters
+    ----------
+    backend : str, optional, default = preferred_backend()
+        Name of the preferred backend, either 'shtools' or 'ducc'.
+    nthreads : int, optional, default = 1
+        Number of threads to use for the 'ducc' backend. Setting this parameter
+        to 0 will use as many threads as there are hardware threads on the
+        system.
+    """
+    backend = backend.lower()
+    if backend == "shtools":
+        from . import shtools
+
+        return shtools
+    elif backend == "ducc":
+        from .wrappers import ducc0_wrapper
+        if not ducc0_wrapper.available():
+            raise ImportError('"ducc" backend requested, but not installed.')
+        if nthreads is not None:
+            ducc0_wrapper.set_nthreads(nthreads)
+        return ducc0_wrapper
+    elif backend is None:
+        return preferred_backend_module()
+    else:
+        print("Unknown backend '{}' requested.".format(backend))
+        raise RuntimeError
+
+
+def select_preferred_backend(backend="shtools", nthreads=None):
+    """
+    Select the preferred backend module used for the spherical harmonic
+    transforms in pyshtools.
+
+    Usage
+    -----
+    select_preferred_backend_module([backend, nthreads])
+
+    Parameters
+    ----------
+    backend : str, optional, default = 'shtools'
+        Name of the preferred backend, either 'shtools' or 'ducc'.
+    nthreads : int, optional, default = 1
+        Number of threads to use for the 'ducc' backend. Setting this parameter
+        to 0 will use as many threads as there are hardware threads on the
+        system.
+    """
+    global _preferred_backend, _preferred_backend_module
+    backend = backend.lower()
+    if backend == "shtools":
+        _preferred_backend = backend
+        from . import shtools
+
+        _preferred_backend_module = shtools
+    elif backend == "ducc":
+        try:
+            import ducc0
+
+            major, minor, patch = ducc0.__version__.split(".")
+            if int(major) < 1 and int(minor) < 15:
+                print(
+                    "ducc0 installation found, but it is too old. "
+                    "Need at least version 0.15"
+                )
+                raise RuntimeError
+        except:
+            print(
+                "DUCC backend requested, but the relevant package cannot be "
+                "imported. Leaving backend unchanged."
+            )
+        _preferred_backend = backend
+        from .wrappers import ducc0_wrapper
+
+        _preferred_backend_module = ducc0_wrapper
+        if nthreads is not None:
+            ducc0_wrapper.set_nthreads(nthreads)
+    else:
+        print("Unknown backend '{}' requested.".format(backend))
+        raise RuntimeError
+
+
+__all__ = ['preferred_backend', 'preferred_backend_module', 'backend_module',
+           'select_preferred_backend']
diff --git a/pyshtools/shclasses/shcoeffs.py b/pyshtools/shclasses/shcoeffs.py
index a55d50b76..8cbdd7737 100644
--- a/pyshtools/shclasses/shcoeffs.py
+++ b/pyshtools/shclasses/shcoeffs.py
@@ -23,6 +23,8 @@
 from ..shio import write_dov as _write_dov
 from ..shio import read_bshc as _read_bshc
 from ..shio import write_bshc as _write_bshc
+from ..backends import backend_module
+from ..backends import preferred_backend
 
 
 class SHCoeffs(object):
@@ -785,7 +787,7 @@ def from_netcdf(self, filename, lmax=None, normalization='4pi', csphase=1,
     @classmethod
     def from_cap(self, theta, lmax, clat=None, clon=None, normalization='4pi',
                  csphase=1, kind='real', name=None, units=None, degrees=True,
-                 copy=True):
+                 copy=True, backend=None, nthreads=None):
         """
         Initialize the class with spherical harmonic coefficients of a
         spherical cap centered at the north pole.
@@ -793,7 +795,8 @@ def from_cap(self, theta, lmax, clat=None, clon=None, normalization='4pi',
         Usage
         -----
         x = SHCoeffs.from_cap(theta, lmax, [clat, clon, normalization, csphase,
-                                            kind, name, units, degrees, copy])
+                                            kind, name, units, degrees, copy,
+                                            backend, nthreads])
 
         Returns
         -------
@@ -826,6 +829,12 @@ def from_cap(self, theta, lmax, clat=None, clon=None, normalization='4pi',
         copy : bool, optional, default = True
             If True, make a copy of array when initializing the class instance.
             If False, initialize the class instance with a reference to array.
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         Notes
         -----
@@ -884,11 +893,15 @@ def from_cap(self, theta, lmax, clat=None, clon=None, normalization='4pi',
                            csphase=csphase, name=name, units=units, copy=copy)
 
         if clat is not None and clon is not None:
+            if backend is None:
+                backend = preferred_backend()
             if degrees is True:
-                temp = temp.rotate(0., -90 + clat, -clon, degrees=True)
+                temp = temp.rotate(0., -90 + clat, -clon, degrees=True,
+                                   backend=backend, nthreads=nthreads)
             else:
                 temp = temp.rotate(0., -_np.pi/2. + clat, -clon,
-                                   degrees=False)
+                                   degrees=False, backend=backend,
+                                   nthreads=nthreads)
 
         return temp
 
@@ -1739,14 +1752,14 @@ def admitcorr(self, hlm, errors=True, lmax=None):
                 admit = _np.column_stack((admit, _np.sqrt(sigma)))
             return admit, corr
 
-    def volume(self, lmax=None):
+    def volume(self, lmax=None, backend=None, nthreads=None):
         """
         If the function is the real shape of an object, calculate the volume
         of the body.
 
         Usage
         -----
-        volume = x.volume([lmax])
+        volume = x.volume([lmax, backend, nthreads])
 
         Returns
         -------
@@ -1758,6 +1771,12 @@ def volume(self, lmax=None):
         lmax : int, optional, default = x.lmax
             The maximum spherical harmonic degree to use when calculating the
             volume.
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         Notes
         -----
@@ -1779,11 +1798,18 @@ def volume(self, lmax=None):
 
         if lmax is None:
             lmax = self.lmax
+        if backend is None:
+            backend = preferred_backend()
+        if backend == 'shtools':
+            lmax = min(3*lmax, 2800)
 
         r0 = self.coeffs[0, 0, 0]
-        grid = self.expand(lmax=min(3*lmax, 2800)) - r0
-        h200 = (grid**2).expand(lmax_calc=0).coeffs[0, 0, 0]
-        h300 = (grid**3).expand(lmax_calc=0).coeffs[0, 0, 0]
+        grid = self.expand(lmax=lmax, backend=backend,
+                           nthreads=nthreads) - r0
+        h200 = (grid**2).expand(lmax_calc=0, backend=backend,
+                                nthreads=nthreads).coeffs[0, 0, 0]
+        h300 = (grid**3).expand(lmax_calc=0, backend=backend,
+                                nthreads=nthreads).coeffs[0, 0, 0]
 
         volume = 4 * _np.pi / 3 * (h300 + 3 * r0 * h200 + r0**3)
         return volume
@@ -1816,7 +1842,7 @@ def centroid(self):
 
     # ---- Operations that return a new SHGravCoeffs class instance ----
     def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
-               body=False, dj_matrix=None):
+               body=False, dj_matrix=None, backend=None, nthreads=None):
         """
         Rotate either the coordinate system used to express the spherical
         harmonic coefficients or the physical body, and return a new class
@@ -1825,7 +1851,7 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
         Usage
         -----
         x_rotated = x.rotate(alpha, beta, gamma, [degrees, convention,
-                             body, dj_matrix])
+                             body, dj_matrix, backend, nthreads])
 
         Returns
         -------
@@ -1845,7 +1871,14 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
         body : bool, optional, default = False
             If true, rotate the physical body and not the coordinate system.
         dj_matrix : ndarray, optional, default = None
-            The djpi2 rotation matrix computed by a call to djpi2.
+            The djpi2 rotation matrix computed by a call to djpi2 (not used if
+            the backend is 'ducc').
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         Notes
         -----
@@ -1914,13 +1947,10 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
         if degrees:
             angles = _np.radians(angles)
 
-        if self.lmax > 1200:
-            _warnings.warn("The rotate() method is accurate only to about" +
-                           " spherical harmonic degree 1200. " +
-                           "lmax = {:d}".format(self.lmax),
-                           category=RuntimeWarning)
+        if backend is None:
+            backend = preferred_backend()
 
-        rot = self._rotate(angles, dj_matrix)
+        rot = self._rotate(angles, dj_matrix, backend, nthreads)
         return rot
 
     def convert(self, normalization=None, csphase=None, lmax=None, kind=None,
@@ -2068,14 +2098,15 @@ def pad(self, lmax, copy=True):
 
     # ---- Expand the coefficients onto a grid ----
     def expand(self, grid='DH2', lat=None, colat=None, lon=None, degrees=True,
-               zeros=None, lmax=None, lmax_calc=None, extend=True):
+               zeros=None, lmax=None, lmax_calc=None, extend=True,
+               backend=None, nthreads=None):
         """
         Evaluate the spherical harmonic coefficients either on a global grid
         or for a list of coordinates.
 
         Usage
         -----
-        f = x.expand([grid, lmax, lmax_calc, zeros])
+        f = x.expand([grid, lmax, lmax_calc, zeros, backend, nthreads])
         g = x.expand(lat=lat, lon=lon, [lmax_calc, degrees])
         g = x.expand(colat=colat, lon=lon, [lmax_calc, degrees])
 
@@ -2110,6 +2141,12 @@ def expand(self, grid='DH2', lat=None, colat=None, lon=None, degrees=True,
         zeros : ndarray, optional, default = None
             The cos(colatitude) nodes used in the Gauss-Legendre Quadrature
             grids.
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         Notes
         -----
@@ -2148,6 +2185,8 @@ def expand(self, grid='DH2', lat=None, colat=None, lon=None, degrees=True,
                 lmax = self.lmax
             if lmax_calc is None:
                 lmax_calc = lmax
+            if backend is None:
+                backend = preferred_backend()
 
             if type(grid) != str:
                 raise ValueError('grid must be a string. Input type is {:s}.'
@@ -2155,13 +2194,16 @@ def expand(self, grid='DH2', lat=None, colat=None, lon=None, degrees=True,
 
             if grid.upper() in ('DH', 'DH1'):
                 gridout = self._expandDH(sampling=1, lmax=lmax,
-                                         lmax_calc=lmax_calc, extend=extend)
+                                         lmax_calc=lmax_calc, extend=extend,
+                                         backend=backend, nthreads=nthreads)
             elif grid.upper() == 'DH2':
                 gridout = self._expandDH(sampling=2, lmax=lmax,
-                                         lmax_calc=lmax_calc, extend=extend)
+                                         lmax_calc=lmax_calc, extend=extend,
+                                         backend=backend, nthreads=nthreads)
             elif grid.upper() == 'GLQ':
                 gridout = self._expandGLQ(zeros=zeros, lmax=lmax,
-                                          lmax_calc=lmax_calc, extend=extend)
+                                          lmax_calc=lmax_calc, extend=extend,
+                                          backend=backend, nthreads=nthreads)
             else:
                 raise ValueError(
                     "grid must be 'DH', 'DH1', 'DH2', or 'GLQ'. " +
@@ -2171,14 +2213,14 @@ def expand(self, grid='DH2', lat=None, colat=None, lon=None, degrees=True,
 
     # ---- Compute the horizontal gradient ----
     def gradient(self, grid='DH2', lmax=None, lmax_calc=None, units=None,
-                 extend=True, radius=None):
+                 extend=True, radius=None, backend=None, nthreads=None):
         """
         Compute the horizontal gradient of the function and return an
         SHGradient class instance.
 
         Usage
         -----
-        g = x.gradient([grid, lmax, lmax_calc, units])
+        g = x.gradient([grid, lmax, lmax_calc, units, backend, nthreads])
 
         Returns
         -------
@@ -2203,12 +2245,20 @@ def gradient(self, grid='DH2', lmax=None, lmax_calc=None, units=None,
         radius : float, optional, default = 1.0
             The radius of the sphere used when computing the gradient of the
             function.
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         """
         if lmax is None:
             lmax = self.lmax
         if lmax_calc is None:
             lmax_calc = lmax
+        if backend is None:
+            backend = preferred_backend()
 
         if type(grid) != str:
             raise ValueError('grid must be a string. Input type is {:s}.'
@@ -2217,11 +2267,13 @@ def gradient(self, grid='DH2', lmax=None, lmax_calc=None, units=None,
         if grid.upper() in ('DH', 'DH1'):
             gradientout = self._gradientDH(sampling=1, lmax=lmax,
                                            lmax_calc=lmax_calc, units=units,
-                                           extend=extend, radius=radius)
+                                           extend=extend, radius=radius,
+                                           backend=backend, nthreads=nthreads)
         elif grid.upper() == 'DH2':
             gradientout = self._gradientDH(sampling=2, lmax=lmax,
                                            lmax_calc=lmax_calc, units=units,
-                                           extend=extend, radius=radius)
+                                           extend=extend, radius=radius,
+                                           backend=backend, nthreads=nthreads)
         elif grid.upper() == 'GLQ':
             raise NotImplementedError('gradient() does not support the use '
                                       'of GLQ grids.')
@@ -3968,13 +4020,20 @@ def _make_complex(self):
                                    csphase=self.csphase, units=self.units,
                                    copy=False)
 
-    def _rotate(self, angles, dj_matrix):
+    def _rotate(self, angles, dj_matrix, backend, nthreads):
         """Rotate the coefficients by the Euler angles alpha, beta, gamma."""
-        if dj_matrix is None:
-            dj_matrix = _shtools.djpi2(self.lmax + 1)
+        if self.lmax > 1200 and backend.lower() == "shtools":
+            _warnings.warn("The rotate() method is accurate only to about" +
+                           " spherical harmonic degree 1200 when using the" +
+                           " shtools backend. " +
+                           "lmax = {:d}".format(self.lmax),
+                           category=RuntimeWarning)
 
         # The coefficients need to be 4pi normalized with csphase = 1
-        coeffs = _shtools.SHRotateRealCoef(
+        if backend == "shtools" and dj_matrix is None:
+            dj_matrix = _shtools.djpi2(self.lmax + 1)
+        coeffs = backend_module(
+            backend=backend, nthreads=nthreads).SHRotateRealCoef(
             self.to_array(normalization='4pi', csphase=1, errors=False),
             angles, dj_matrix)
 
@@ -3991,7 +4050,7 @@ def _rotate(self, angles, dj_matrix):
             return SHCoeffs.from_array(coeffs, errors=self.errors,
                                        units=self.units, copy=False)
 
-    def _expandDH(self, sampling, lmax, lmax_calc, extend):
+    def _expandDH(self, sampling, lmax, lmax_calc, extend, backend, nthreads):
         """Evaluate the coefficients on a Driscoll and Healy (1994) grid."""
         from .shgrid import SHGrid
         if self.normalization == '4pi':
@@ -4008,14 +4067,16 @@ def _expandDH(self, sampling, lmax, lmax_calc, extend):
                 "'unnorm'. Input value is {:s}."
                 .format(repr(self.normalization)))
 
-        data = _shtools.MakeGridDH(self.coeffs, sampling=sampling, norm=norm,
-                                   csphase=self.csphase, lmax=lmax,
-                                   lmax_calc=lmax_calc, extend=extend)
+        data = backend_module(
+            backend=backend, nthreads=nthreads).MakeGridDH(
+                self.coeffs, sampling=sampling, norm=norm,
+                csphase=self.csphase, lmax=lmax,
+                lmax_calc=lmax_calc, extend=extend)
         gridout = SHGrid.from_array(data, grid='DH', units=self.units,
                                     copy=False)
         return gridout
 
-    def _expandGLQ(self, zeros, lmax, lmax_calc, extend):
+    def _expandGLQ(self, zeros, lmax, lmax_calc, extend, backend, nthreads):
         """Evaluate the coefficients on a Gauss Legendre quadrature grid."""
         from .shgrid import SHGrid
         if self.normalization == '4pi':
@@ -4032,12 +4093,13 @@ def _expandGLQ(self, zeros, lmax, lmax_calc, extend):
                 "'unnorm'. Input value is {:s}."
                 .format(repr(self.normalization)))
 
-        if zeros is None:
+        if backend == "shtools" and zeros is None:
             zeros, weights = _shtools.SHGLQ(self.lmax)
-
-        data = _shtools.MakeGridGLQ(self.coeffs, zeros, norm=norm,
-                                    csphase=self.csphase, lmax=lmax,
-                                    lmax_calc=lmax_calc, extend=extend)
+        data = backend_module(
+            backend=backend, nthreads=nthreads).MakeGridGLQ(
+                self.coeffs, zero=zeros, norm=norm,
+                csphase=self.csphase, lmax=lmax,
+                lmax_calc=lmax_calc, extend=extend)
         gridout = SHGrid.from_array(data, grid='GLQ', units=self.units,
                                     copy=False)
         return gridout
@@ -4097,14 +4159,17 @@ def _expand_coord(self, lat, lon, lmax_calc, degrees):
                              'ndarray, or list. Input types are {:s} and {:s}.'
                              .format(repr(type(lat)), repr(type(lon))))
 
-    def _gradientDH(self, sampling, lmax, lmax_calc, units, extend, radius):
+    def _gradientDH(self, sampling, lmax, lmax_calc, units, extend, radius,
+                    backend, nthreads):
         """Evaluate the gradient on a Driscoll and Healy (1994) grid."""
         from .shgradient import SHGradient
 
-        theta, phi = _shtools.MakeGradientDH(
-            self.to_array(normalization='4pi', csphase=1, errors=False),
-            sampling=sampling, lmax=lmax, lmax_calc=lmax_calc, extend=extend,
-            radius=radius)
+        theta, phi = backend_module(
+                backend=backend, nthreads=nthreads).MakeGradientDH(
+                    self.to_array(
+                        normalization='4pi', csphase=1, errors=False),
+                    sampling=sampling, lmax=lmax,
+                    lmax_calc=lmax_calc, extend=extend, radius=radius)
 
         return SHGradient(theta, phi, lmax, lmax_calc, units=units)
 
@@ -4198,18 +4263,42 @@ def _make_real(self, check=True):
                                    normalization=self.normalization,
                                    csphase=self.csphase, units=self.units)
 
-    def _rotate(self, angles, dj_matrix):
+    def _rotate(self, angles, dj_matrix, backend, nthreads):
         """Rotate the coefficients by the Euler angles alpha, beta, gamma."""
+        if backend == 'ducc':
+            coeffs = backend_module(
+                backend=backend, nthreads=nthreads).SHRotateComplexCoef(
+                    self.to_array(normalization='4pi', csphase=1,
+                                  errors=False), angles)
+            # Convert 4pi normalized coefficients to the same normalization
+            # as the unrotated coefficients.
+            if self.normalization != '4pi' or self.csphase != 1:
+                temp = _convert(coeffs, normalization_in='4pi', csphase_in=1,
+                                normalization_out=self.normalization,
+                                csphase_out=self.csphase)
+                return SHCoeffs.from_array(
+                    temp, errors=self.errors, normalization=self.normalization,
+                    csphase=self.csphase, units=self.units, copy=False)
+            else:
+                return SHCoeffs.from_array(coeffs, errors=self.errors,
+                                           units=self.units, copy=False)
+
         # Note that the current method is EXTREMELY inefficient. The complex
         # coefficients are expanded onto real and imaginary grids, each of
         # the two components are rotated separately as real data, the rotated
         # real data are re-expanded on new real and complex grids, they are
         # combined to make a complex grid, and the resultant is expanded
         # in complex spherical harmonics.
+        if self.lmax > 1200 and backend.lower() == "shtools":
+            _warnings.warn("The rotate() method is accurate only to about" +
+                           " spherical harmonic degree 1200 when using the" +
+                           " shtools backend. " +
+                           "lmax = {:d}".format(self.lmax),
+                           category=RuntimeWarning)
         if dj_matrix is None:
             dj_matrix = _shtools.djpi2(self.lmax + 1)
 
-        cgrid = self.expand(grid='DH')
+        cgrid = self.expand(grid='DH', extend=False)
         rgrid, igrid = cgrid.data.real, cgrid.data.imag
         rgridcoeffs = _shtools.SHExpandDH(rgrid, norm=1, sampling=1, csphase=1)
         igridcoeffs = _shtools.SHExpandDH(igrid, norm=1, sampling=1, csphase=1)
@@ -4242,12 +4331,12 @@ def _rotate(self, angles, dj_matrix):
         coeffs_rot = _shtools.SHExpandDHC(grid_rot, norm=norm,
                                           csphase=self.csphase)
 
-        return SHCoeffs.from_array(coeffs_rot, errros=self.errors,
+        return SHCoeffs.from_array(coeffs_rot, errors=self.errors,
                                    normalization=self.normalization,
                                    csphase=self.csphase, units=self.units,
                                    copy=False)
 
-    def _expandDH(self, sampling, lmax, lmax_calc, extend):
+    def _expandDH(self, sampling, lmax, lmax_calc, extend, backend, nthreads):
         """Evaluate the coefficients on a Driscoll and Healy (1994) grid."""
         from .shgrid import SHGrid
         if self.normalization == '4pi':
@@ -4264,14 +4353,16 @@ def _expandDH(self, sampling, lmax, lmax_calc, extend):
                 "'unnorm'. Input value is {:s}."
                 .format(repr(self.normalization)))
 
-        data = _shtools.MakeGridDHC(self.coeffs, sampling=sampling,
-                                    norm=norm, csphase=self.csphase, lmax=lmax,
-                                    lmax_calc=lmax_calc, extend=extend)
+        data = backend_module(
+                backend=backend, nthreads=nthreads).MakeGridDHC(
+                    self.coeffs, sampling=sampling, norm=norm,
+                    csphase=self.csphase, lmax=lmax,
+                    lmax_calc=lmax_calc, extend=extend)
         gridout = SHGrid.from_array(data, grid='DH', units=self.units,
                                     copy=False)
         return gridout
 
-    def _expandGLQ(self, zeros, lmax, lmax_calc, extend):
+    def _expandGLQ(self, zeros, lmax, lmax_calc, extend, backend, nthreads):
         """Evaluate the coefficients on a Gauss-Legendre quadrature grid."""
         from .shgrid import SHGrid
         if self.normalization == '4pi':
@@ -4288,12 +4379,13 @@ def _expandGLQ(self, zeros, lmax, lmax_calc, extend):
                 "'unnorm'. Input value is {:s}."
                 .format(repr(self.normalization)))
 
-        if zeros is None:
+        if backend == "shtools" and zeros is None:
             zeros, weights = _shtools.SHGLQ(self.lmax)
-
-        data = _shtools.MakeGridGLQC(self.coeffs, zeros, norm=norm,
-                                     csphase=self.csphase, lmax=lmax,
-                                     lmax_calc=lmax_calc, extend=extend)
+        data = backend_module(
+                backend=backend, nthreads=nthreads).MakeGridGLQC(
+                    self.coeffs, zero=zeros, norm=norm,
+                    csphase=self.csphase, lmax=lmax,
+                    lmax_calc=lmax_calc, extend=extend)
         gridout = SHGrid.from_array(data, grid='GLQ', units=self.units,
                                     copy=False)
         return gridout
diff --git a/pyshtools/shclasses/shgravcoeffs.py b/pyshtools/shclasses/shgravcoeffs.py
index f8929d914..ba274782a 100644
--- a/pyshtools/shclasses/shgravcoeffs.py
+++ b/pyshtools/shclasses/shgravcoeffs.py
@@ -37,8 +37,9 @@
 from ..shtools import MakeGravGradGridDH as _MakeGravGradGridDH
 from ..shtools import MakeGeoidGridDH as _MakeGeoidGridDH
 from ..shtools import djpi2 as _djpi2
-from ..shtools import SHRotateRealCoef as _SHRotateRealCoef
 from ..shtools import MakeGravGridPoint as _MakeGravGridPoint
+from ..backends import backend_module
+from ..backends import preferred_backend
 
 
 class SHGravCoeffs(object):
@@ -976,7 +977,8 @@ def from_netcdf(self, filename, lmax=None, normalization='4pi', csphase=1,
 
     @classmethod
     def from_shape(self, shape, rho, gm, nmax=7, lmax=None, lmax_grid=None,
-                   lmax_calc=None, omega=None, name=None, epoch=None):
+                   lmax_calc=None, omega=None, name=None, epoch=None,
+                   backend=None, nthreads=None):
         """
         Initialize a class of gravitational potential spherical harmonic
         coefficients by calculuting the gravitational potential associatiated
@@ -986,7 +988,7 @@ def from_shape(self, shape, rho, gm, nmax=7, lmax=None, lmax_grid=None,
         -----
         x = SHGravCoeffs.from_shape(shape, rho, gm, [nmax, lmax, lmax_grid,
                                                      lmax_calc, omega, name,
-                                                     epoch])
+                                                     epoch, backend, nthreads])
 
         Returns
         -------
@@ -1028,6 +1030,12 @@ def from_shape(self, shape, rho, gm, nmax=7, lmax=None, lmax_grid=None,
         epoch : str or float, optional, default = None
             The epoch time of the spherical harmonic coefficients as given by
             the format YYYYMMDD.DD.
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         Notes
         -----
@@ -1066,6 +1074,9 @@ def from_shape(self, shape, rho, gm, nmax=7, lmax=None, lmax_grid=None,
         """
         mass = gm / _G.value
 
+        if backend is None:
+            backend = preferred_backend()
+
         if type(shape) is not _SHRealCoeffs and type(shape) is not _DHRealGrid:
             raise ValueError('shape must be of type SHRealCoeffs '
                              'or DHRealGrid. Input type is {:s}.'
@@ -1080,10 +1091,12 @@ def from_shape(self, shape, rho, gm, nmax=7, lmax=None, lmax_grid=None,
                              .format(repr(type(rho))))
 
         if type(shape) is _SHRealCoeffs:
-            shape = shape.expand(lmax=lmax_grid, lmax_calc=lmax_calc)
+            shape = shape.expand(lmax=lmax_grid, lmax_calc=lmax_calc,
+                                 backend=backend, nthreads=nthreads)
 
         if type(rho) is _SHRealCoeffs:
-            rho = rho.expand(lmax=lmax_grid, lmax_calc=lmax_calc)
+            rho = rho.expand(lmax=lmax_grid, lmax_calc=lmax_calc,
+                             backend=backend, nthreads=nthreads)
 
         if type(rho) is _DHRealGrid:
             if shape.lmax != rho.lmax:
@@ -2049,7 +2062,7 @@ def admitcorr(self, hlm, errors=True, function='radial', lmax=None):
 
     # ---- Operations that return a new SHGravCoeffs class instance ----
     def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
-               body=False, dj_matrix=None):
+               body=False, dj_matrix=None, backend=None, nthreads=None):
         """
         Rotate either the coordinate system used to express the spherical
         harmonic coefficients or the physical body, and return a new class
@@ -2058,7 +2071,7 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
         Usage
         -----
         x_rotated = x.rotate(alpha, beta, gamma, [degrees, convention,
-                             body, dj_matrix])
+                             body, dj_matrix, backend, nthreads])
 
         Returns
         -------
@@ -2078,7 +2091,14 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
         body : bool, optional, default = False
             If true, rotate the physical body and not the coordinate system.
         dj_matrix : ndarray, optional, default = None
-            The djpi2 rotation matrix computed by a call to djpi2.
+            The djpi2 rotation matrix computed by a call to djpi2 (not used if
+            the backend is 'ducc').
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         Notes
         -----
@@ -2146,15 +2166,12 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
 
         if degrees:
             angles = _np.radians(angles)
-
-        if self.lmax > 1200:
-            _warnings.warn("The rotate() method is accurate only to about"
-                           " spherical harmonic degree 1200. "
-                           "lmax = {:d}".format(self.lmax),
-                           category=RuntimeWarning)
+        if backend is None:
+            backend = preferred_backend()
 
         rot = self._rotate(angles, dj_matrix, gm=self.gm, r0=self.r0,
-                           omega=self.omega)
+                           omega=self.omega, backend=backend,
+                           nthreads=nthreads)
         return rot
 
     def convert(self, normalization=None, csphase=None, lmax=None):
@@ -3754,15 +3771,23 @@ def __repr__(self):
                         repr(self.header), repr(self.header2),
                         repr(self.name), repr(self.epoch)))
 
-    def _rotate(self, angles, dj_matrix, gm=None, r0=None, omega=None):
+    def _rotate(self, angles, dj_matrix, gm=None, r0=None, omega=None,
+                backend=None, nthreads=None):
         """Rotate the coefficients by the Euler angles alpha, beta, gamma."""
-        if dj_matrix is None:
+        if self.lmax > 1200 and backend.lower() == "shtools":
+            _warnings.warn("The rotate() method is accurate only to about" +
+                           " spherical harmonic degree 1200 when using the" +
+                           " shtools backend. " +
+                           "lmax = {:d}".format(self.lmax),
+                           category=RuntimeWarning)
+        if backend == "shtools" and dj_matrix is None:
             dj_matrix = _djpi2(self.lmax + 1)
 
         # The coefficients need to be 4pi normalized with csphase = 1
-        coeffs = _SHRotateRealCoef(
-            self.to_array(normalization='4pi', csphase=1, errors=False),
-            angles, dj_matrix)
+        coeffs = backend_module(
+            backend=backend, nthreads=nthreads).SHRotateRealCoef(
+                self.to_array(normalization='4pi', csphase=1, errors=False),
+                angles, dj_matrix)
 
         # Convert 4pi normalized coefficients to the same normalization
         # as the unrotated coefficients.
diff --git a/pyshtools/shclasses/shgrid.py b/pyshtools/shclasses/shgrid.py
index 74d078ce1..70b0004c9 100644
--- a/pyshtools/shclasses/shgrid.py
+++ b/pyshtools/shclasses/shgrid.py
@@ -8,7 +8,8 @@
 import copy as _copy
 import xarray as _xr
 import tempfile as _tempfile
-
+from ..backends import backend_module
+from ..backends import preferred_backend
 from .. import shtools as _shtools
 
 try:
@@ -951,13 +952,14 @@ def histogram(self, bins=10, range=None):
         """
         return self._histogram(bins=bins, range=range)
 
-    def expand(self, normalization='4pi', csphase=1, **kwargs):
+    def expand(self, normalization='4pi', csphase=1, lmax_calc=None,
+               backend=None, nthreads=None):
         """
         Expand the grid into spherical harmonics.
 
         Usage
         -----
-        clm = x.expand([normalization, csphase, lmax_calc])
+        clm = x.expand([normalization, csphase, lmax_calc, backend, nthreads])
 
         Returns
         -------
@@ -974,6 +976,12 @@ def expand(self, normalization='4pi', csphase=1, **kwargs):
             or -1 to include it.
         lmax_calc : int, optional, default = x.lmax
             Maximum spherical harmonic degree to return.
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         Notes
         -----
@@ -1000,8 +1008,12 @@ def expand(self, normalization='4pi', csphase=1, **kwargs):
                 .format(repr(csphase))
                 )
 
+        if backend is None:
+            backend = preferred_backend()
+
         return self._expand(normalization=normalization, csphase=csphase,
-                            **kwargs)
+                            lmax_calc=lmax_calc, backend=backend,
+                            nthreads=nthreads)
 
     # ---- Plotting routines ----
     def plot3d(self, elevation=20, azimuth=30, cmap='viridis',
@@ -1841,7 +1853,7 @@ def _histogram(self, bins=None, range=None):
                              weights=da[:, :self.nlon-self.extend],
                              density=True, range=range)
 
-    def _expand(self, normalization, csphase, **kwargs):
+    def _expand(self, normalization, csphase, lmax_calc, backend, nthreads):
         """Expand the grid into real spherical harmonics."""
         from .shcoeffs import SHCoeffs
         if normalization.lower() == '4pi':
@@ -1859,10 +1871,11 @@ def _expand(self, normalization, csphase, **kwargs):
                 .format(repr(normalization))
                 )
 
-        cilm = _shtools.SHExpandDH(self.data[:self.nlat-self.extend,
-                                             :self.nlon-self.extend],
-                                   norm=norm, csphase=csphase,
-                                   sampling=self.sampling, **kwargs)
+        cilm = backend_module(
+            backend=backend, nthreads=nthreads).SHExpandDH(
+                self.data[:self.nlat-self.extend, :self.nlon-self.extend],
+                norm=norm, csphase=csphase, sampling=self.sampling,
+                lmax_calc=lmax_calc)
         coeffs = SHCoeffs.from_array(cilm,
                                      normalization=normalization.lower(),
                                      csphase=csphase, units=self.units,
@@ -2390,7 +2403,7 @@ def _histogram(self, **kwargs):
         """Return an area-weighted histogram normalized to unity."""
         raise NotImplementedError('histogram() does not support complex data.')
 
-    def _expand(self, normalization, csphase, **kwargs):
+    def _expand(self, normalization, csphase, lmax_calc, backend, nthreads):
         """Expand the grid into real spherical harmonics."""
         from .shcoeffs import SHCoeffs
         if normalization.lower() == '4pi':
@@ -2408,10 +2421,11 @@ def _expand(self, normalization, csphase, **kwargs):
                 .format(repr(normalization))
                 )
 
-        cilm = _shtools.SHExpandDHC(self.data[:self.nlat-self.extend,
-                                              :self.nlon-self.extend],
-                                    norm=norm, csphase=csphase,
-                                    sampling=self.sampling, **kwargs)
+        cilm = backend_module(
+            backend=backend, nthreads=nthreads).SHExpandDHC(
+                self.data[:self.nlat-self.extend, :self.nlon-self.extend],
+                norm=norm, csphase=csphase, sampling=self.sampling,
+                lmax_calc=lmax_calc)
         coeffs = SHCoeffs.from_array(cilm, normalization=normalization.lower(),
                                      csphase=csphase, units=self.units,
                                      copy=False)
@@ -2629,7 +2643,7 @@ def _histogram(self, bins=None, range=None):
                              weights=da[:, :self.nlon-self.extend],
                              density=True, range=range)
 
-    def _expand(self, normalization, csphase, **kwargs):
+    def _expand(self, normalization, csphase, lmax_calc, backend, nthreads):
         """Expand the grid into real spherical harmonics."""
         from .shcoeffs import SHCoeffs
         if normalization.lower() == '4pi':
@@ -2647,9 +2661,10 @@ def _expand(self, normalization, csphase, **kwargs):
                 .format(repr(normalization))
                 )
 
-        cilm = _shtools.SHExpandGLQ(self.data[:, :self.nlon-self.extend],
-                                    self.weights, self.zeros, norm=norm,
-                                    csphase=csphase, **kwargs)
+        cilm = backend_module(
+            backend=backend, nthreads=nthreads).SHExpandGLQ(
+                self.data[:, :self.nlon-self.extend], self.weights, self.zeros,
+                norm=norm, csphase=csphase, lmax_calc=lmax_calc)
         coeffs = SHCoeffs.from_array(cilm, normalization=normalization.lower(),
                                      csphase=csphase, units=self.units,
                                      copy=False)
@@ -2940,7 +2955,7 @@ def _histogram(self, **kwargs):
         """Return an area-weighted histogram normalized to unity."""
         raise NotImplementedError('histogram() does not support complex data.')
 
-    def _expand(self, normalization, csphase, **kwargs):
+    def _expand(self, normalization, csphase, lmax_calc, backend, nthreads):
         """Expand the grid into real spherical harmonics."""
         from .shcoeffs import SHCoeffs
         if normalization.lower() == '4pi':
@@ -2958,9 +2973,10 @@ def _expand(self, normalization, csphase, **kwargs):
                 .format(repr(normalization))
                 )
 
-        cilm = _shtools.SHExpandGLQC(self.data[:, :self.nlon-self.extend],
-                                     self.weights, self.zeros, norm=norm,
-                                     csphase=csphase, **kwargs)
+        cilm = backend_module(
+            backend=backend, nthreads=nthreads).SHExpandGLQC(
+                self.data[:, :self.nlon-self.extend], self.weights, self.zeros,
+                norm=norm, csphase=csphase, lmax_calc=lmax_calc)
         coeffs = SHCoeffs.from_array(cilm, normalization=normalization.lower(),
                                      csphase=csphase, units=self.units,
                                      copy=False)
diff --git a/pyshtools/shclasses/shmagcoeffs.py b/pyshtools/shclasses/shmagcoeffs.py
index 1dd49ac4c..1571152e5 100644
--- a/pyshtools/shclasses/shmagcoeffs.py
+++ b/pyshtools/shclasses/shmagcoeffs.py
@@ -29,8 +29,9 @@
 from ..shtools import MakeMagGridDH as _MakeMagGridDH
 from ..shtools import MakeMagGradGridDH as _MakeMagGradGridDH
 from ..shtools import djpi2 as _djpi2
-from ..shtools import SHRotateRealCoef as _SHRotateRealCoef
 from ..shtools import MakeMagGridPoint as _MakeMagGridPoint
+from ..backends import backend_module
+from ..backends import preferred_backend
 
 
 class SHMagCoeffs(object):
@@ -1538,7 +1539,7 @@ def correlation(self, hlm, lmax=None):
 
     # ---- Operations that return a new SHMagCoeffs class instance ----
     def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
-               body=False, dj_matrix=None):
+               body=False, dj_matrix=None, backend=None, nthreads=None):
         """
         Rotate either the coordinate system used to express the spherical
         harmonic coefficients or the physical body, and return a new class
@@ -1547,7 +1548,7 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
         Usage
         -----
         x_rotated = x.rotate(alpha, beta, gamma, [degrees, convention,
-                             body, dj_matrix])
+                             body, dj_matrix, backend, nthreads])
 
         Returns
         -------
@@ -1567,7 +1568,14 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
         body : bool, optional, default = False
             If true, rotate the physical body and not the coordinate system.
         dj_matrix : ndarray, optional, default = None
-            The djpi2 rotation matrix computed by a call to djpi2.
+            The djpi2 rotation matrix computed by a call to djpi2 (not used if
+            the backend is 'ducc').
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
 
         Notes
         -----
@@ -1636,6 +1644,8 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
 
         if degrees:
             angles = _np.radians(angles)
+        if backend is None:
+            backend = preferred_backend()
 
         if self.lmax > 1200:
             _warnings.warn("The rotate() method is accurate only to about"
@@ -1643,7 +1653,8 @@ def rotate(self, alpha, beta, gamma, degrees=True, convention='y',
                            "lmax = {:d}.".format(self.lmax),
                            category=RuntimeWarning)
 
-        rot = self._rotate(angles, dj_matrix, r0=self.r0)
+        rot = self._rotate(angles, dj_matrix, r0=self.r0, backend=backend,
+                           nthread=nthreads)
         return rot
 
     def convert(self, normalization=None, csphase=None, lmax=None):
@@ -2848,15 +2859,22 @@ def __repr__(self):
                         repr(self.header2), repr(self.name), repr(self.units),
                         repr(self.year)))
 
-    def _rotate(self, angles, dj_matrix, r0=None):
+    def _rotate(self, angles, dj_matrix, r0=None, backend=None, nthreads=None):
         """Rotate the coefficients by the Euler angles alpha, beta, gamma."""
-        if dj_matrix is None:
+        if self.lmax > 1200 and backend.lower() == "shtools":
+            _warnings.warn("The rotate() method is accurate only to about" +
+                           " spherical harmonic degree 1200 when using the" +
+                           " shtools backend. " +
+                           "lmax = {:d}".format(self.lmax),
+                           category=RuntimeWarning)
+        if backend == "shtools" and dj_matrix is None:
             dj_matrix = _djpi2(self.lmax + 1)
 
         # The coefficients need to be 4pi normalized with csphase = 1
-        coeffs = _SHRotateRealCoef(
-            self.to_array(normalization='4pi', csphase=1, errors=False),
-            angles, dj_matrix)
+        coeffs = backend_module(
+            backend=backend, nthreads=nthreads).SHRotateRealCoef(
+                self.to_array(normalization='4pi', csphase=1, errors=False),
+                angles, dj_matrix)
 
         # Convert 4pi normalized coefficients to the same normalization
         # as the unrotated coefficients.
diff --git a/pyshtools/shclasses/slepiancoeffs.py b/pyshtools/shclasses/slepiancoeffs.py
index 0a4c0fdfb..f42a9bbf3 100644
--- a/pyshtools/shclasses/slepiancoeffs.py
+++ b/pyshtools/shclasses/slepiancoeffs.py
@@ -4,6 +4,8 @@
 import numpy as _np
 import copy as _copy
 
+from ..backends import backend_module
+from ..backends import preferred_backend
 from .. import shtools as _shtools
 
 from .shcoeffs import SHCoeffs
@@ -83,7 +85,8 @@ def __repr__(self):
         str += '\nSlepian functions:\n' + self.galpha.__repr__()
         return str
 
-    def expand(self, nmax=None, grid='DH2', zeros=None, extend=True):
+    def expand(self, nmax=None, grid='DH2', zeros=None, extend=True,
+               backend=None, nthreads=None):
         """
         Expand the function on a grid using the first n Slepian coefficients.
 
@@ -110,6 +113,12 @@ def expand(self, nmax=None, grid='DH2', zeros=None, extend=True):
         extend : bool, optional, default = True
             If True, compute the longitudinal band for 360 E (DH and GLQ grids)
             and the latitudinal band for 90 S (DH grids only).
+        backend : str, optional, default = preferred_backend()
+            Name of the preferred backend, either 'shtools' or 'ducc'.
+        nthreads : int, optional, default = 1
+            Number of threads to use for the 'ducc' backend. Setting this
+            parameter to 0 will use as many threads as there are hardware
+            threads on the system.
         """
         if type(grid) != str:
             raise ValueError('grid must be a string. ' +
@@ -118,6 +127,8 @@ def expand(self, nmax=None, grid='DH2', zeros=None, extend=True):
 
         if nmax is None:
             nmax = self.nmax
+        if backend is None:
+            backend = preferred_backend()
 
         if self.galpha.kind == 'cap':
             shcoeffs = _shtools.SlepianCoeffsToSH(self.falpha,
@@ -127,18 +138,21 @@ def expand(self, nmax=None, grid='DH2', zeros=None, extend=True):
                                                   self.galpha.tapers, nmax)
 
         if grid.upper() in ('DH', 'DH1'):
-            gridout = _shtools.MakeGridDH(shcoeffs, sampling=1,
-                                          norm=1, csphase=1, extend=extend)
+            gridout = backend_module(
+                backend=backend, nthreads=nthreads).MakeGridDH(
+                    shcoeffs, sampling=1, norm=1, csphase=1, extend=extend)
             return SHGrid.from_array(gridout, grid='DH', copy=False)
         elif grid.upper() == 'DH2':
-            gridout = _shtools.MakeGridDH(shcoeffs, sampling=2,
-                                          norm=1, csphase=1, extend=extend)
+            gridout = backend_module(
+                backend=backend, nthreads=nthreads).MakeGridDH(
+                    shcoeffs, sampling=2, norm=1, csphase=1, extend=extend)
             return SHGrid.from_array(gridout, grid='DH', copy=False)
         elif grid.upper() == 'GLQ':
-            if zeros is None:
+            if backend == "shtools" and zeros is None:
                 zeros, weights = _shtools.SHGLQ(self.galpha.lmax)
-            gridout = _shtools.MakeGridGLQ(shcoeffs, zeros, norm=1, csphase=1,
-                                           extend=extend)
+            gridout = backend_module(
+                backend=backend, nthreads=nthreads).MakeGridGLQ(
+                    shcoeffs, zeros=zeros, norm=1, csphase=1, extend=extend)
             return SHGrid.from_array(gridout, grid='GLQ', copy=False)
         else:
             raise ValueError(
diff --git a/pyshtools/wrappers/__init__.py b/pyshtools/wrappers/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/pyshtools/wrappers/ducc0_wrapper.py b/pyshtools/wrappers/ducc0_wrapper.py
new file mode 100644
index 000000000..6cb1fc769
--- /dev/null
+++ b/pyshtools/wrappers/ducc0_wrapper.py
@@ -0,0 +1,443 @@
+"""
+    Distinctly Useful Code Collection (DUCC)
+
+    DUCC wrapper functions for use in pyshtools.
+"""
+import numpy as _np
+
+try:
+    import ducc0
+
+    major, minor, patch = ducc0.__version__.split(".")
+    if int(major) < 1 and int(minor) < 15:
+        raise RuntimeError
+except:
+    ducc0 = None
+
+# setup a few required variables
+if ducc0 is not None:
+    import os as _os
+
+    try:
+        nthreads = int(_os.environ["OMP_NUM_THREADS"])
+    except:
+        nthreads = 0
+
+
+def set_nthreads(ntnew):
+    global nthreads
+    nthreads = ntnew
+
+
+def available():
+    return ducc0 is not None
+
+
+def _nalm(lmax, mmax):
+    return ((mmax + 1) * (mmax + 2)) // 2 + (mmax + 1) * (lmax - mmax)
+
+
+# ducc0's only accepted conventions are
+#   normalization = 'ortho'
+#   csphase = 1
+# so we need to make the required adjustments
+
+
+def _get_norm(lmax, norm):
+    if norm == 1:
+        return _np.full(lmax + 1, _np.sqrt(4 * _np.pi))
+    if norm == 2:
+        return _np.sqrt(4 * _np.pi / (2 * _np.arange(lmax + 1) + 1.0))
+    if norm == 3:
+        return _np.sqrt(2 * _np.pi / (2 * _np.arange(lmax + 1) + 1.0))
+    if norm == 4:
+        return _np.ones(lmax + 1)
+    raise RuntimeError("unsupported normalization")
+
+
+def _rcilm2alm(cilm, lmax):
+    alm = _np.empty((_nalm(lmax, lmax),), dtype=_np.complex128)
+    alm[0:lmax + 1] = cilm[0, :, 0]
+    ofs = lmax + 1
+    for m in range(1, lmax + 1):
+        alm[ofs:ofs + lmax + 1 - m].real = cilm[0, m:, m]
+        alm[ofs:ofs + lmax + 1 - m].imag = cilm[1, m:, m]
+        ofs += lmax + 1 - m
+    return alm
+
+
+def _ralm2cilm(alm, lmax):
+    cilm = _np.zeros((2, lmax + 1, lmax + 1), dtype=_np.float64)
+    cilm[0, :, 0] = alm[0:lmax + 1].real
+    ofs = lmax + 1
+    for m in range(1, lmax + 1):
+        cilm[0, m:, m] = alm[ofs:ofs + lmax + 1 - m].real
+        cilm[1, m:, m] = alm[ofs:ofs + lmax + 1 - m].imag
+        ofs += lmax + 1 - m
+    return cilm
+
+
+def _apply_norm(alm, lmax, norm, csphase, reverse):
+    lnorm = _get_norm(lmax, norm)
+    if reverse:
+        lnorm = 1.0 / lnorm
+    alm[0:lmax + 1] *= lnorm[0:lmax + 1]
+    lnorm *= _np.sqrt(2.0) if reverse else (1.0 / _np.sqrt(2.0))
+    mlnorm = -lnorm
+    ofs = lmax + 1
+    for m in range(1, lmax + 1):
+        if csphase == 1:
+            if m & 1:
+                alm[ofs:ofs + lmax + 1 - m].real *= mlnorm[m:]
+                alm[ofs:ofs + lmax + 1 - m].imag *= lnorm[m:]
+            else:
+                alm[ofs:ofs + lmax + 1 - m].real *= lnorm[m:]
+                alm[ofs:ofs + lmax + 1 - m].imag *= mlnorm[m:]
+        else:
+            alm[ofs:ofs + lmax + 1 - m].real *= lnorm[m:]
+            alm[ofs:ofs + lmax + 1 - m].imag *= mlnorm[m:]
+        ofs += lmax + 1 - m
+    if norm == 3:  # special treatment for unnormalized a_lm
+        r = _np.arange(lmax + 1)
+        fct = _np.ones(lmax + 1)
+        ofs = lmax + 1
+        if reverse:
+            alm[0:lmax + 1] /= _np.sqrt(2)
+            for m in range(1, lmax + 1):
+                fct[m:] *= _np.sqrt((r[m:] + m) * (r[m:] - m + 1))
+                alm[ofs:ofs + lmax + 1 - m] /= fct[m:]
+                ofs += lmax + 1 - m
+        else:
+            alm[0:lmax + 1] *= _np.sqrt(2)
+            for m in range(1, lmax + 1):
+                fct[m:] *= _np.sqrt((r[m:] + m) * (r[m:] - m + 1))
+                alm[ofs:ofs + lmax + 1 - m] *= fct[m:]
+                ofs += lmax + 1 - m
+    return alm
+
+
+def _make_alm(cilm, lmax, norm, csphase):
+    alm = _rcilm2alm(cilm, lmax)
+    return _apply_norm(alm, lmax, norm, csphase, False)
+
+
+def _extract_alm(alm, lmax, norm, csphase):
+    _apply_norm(alm, lmax, norm, csphase, True)
+    return _ralm2cilm(alm, lmax)
+
+
+def _synthesize_DH(alm, lmax, extend, out):
+    ducc0.sht.experimental.synthesis_2d(
+        alm=alm.reshape((1, -1)),
+        map=out[:, : out.shape[1] - extend].reshape(
+            (1, out.shape[0], out.shape[1] - extend)
+        ),
+        spin=0,
+        lmax=lmax,
+        geometry="CC" if extend else "DH",
+        nthreads=nthreads,
+    )
+    if extend:
+        out[:, -1] = out[:, 0]
+    return out
+
+
+def _synthesize_DH_deriv1(alm, lmax, extend, out):
+    ducc0.sht.experimental.synthesis_2d_deriv1(
+        alm=alm.reshape((1, -1)),
+        map=out[:, :, : out.shape[2] - extend],
+        lmax=lmax,
+        geometry="CC" if extend else "DH",
+        nthreads=nthreads,
+    )
+    out[:, 0, :] = 0.0
+    if extend:
+        out[:, -1, :] = 0.0
+        out[:, :, -1] = out[:, :, 0]
+    return out
+
+
+def _synthesize_GLQ(alm, lmax, extend, out):
+    ducc0.sht.experimental.synthesis_2d(
+        alm=alm.reshape((1, -1)),
+        map=out[:, : out.shape[1] - extend].reshape(
+            (1, out.shape[0], out.shape[1] - extend)
+        ),
+        spin=0,
+        lmax=lmax,
+        geometry="GL",
+        nthreads=nthreads,
+    )
+    if extend:
+        out[:, -1] = out[:, 0]
+    return out
+
+
+def _analyze_DH(map, lmax):
+    alm = ducc0.sht.experimental.analysis_2d(
+        map=map.reshape((1, map.shape[0], map.shape[1])),
+        spin=0,
+        lmax=lmax,
+        geometry="DH",
+        nthreads=nthreads,
+    )
+    return alm[0]
+
+
+def _analyze_GLQ(map, lmax):
+    alm = ducc0.sht.experimental.analysis_2d(
+        map=map.reshape((1, map.shape[0], map.shape[1])),
+        spin=0,
+        lmax=lmax,
+        geometry="GL",
+        nthreads=nthreads,
+    )
+    return alm[0]
+
+
+def _ccilm2almr(cilm):
+    lmax = cilm.shape[1] - 1
+    alm = _np.empty((_nalm(lmax, lmax),), dtype=_np.complex128)
+    fct = (-1) ** _np.arange(lmax + 1)
+    alm[0:lmax + 1] = cilm[0, :, 0].real
+    ofs = lmax + 1
+    for m in range(1, lmax + 1):
+        tmp = _np.conj(cilm[1, m:, m])
+        tmp *= fct[m]
+        tmp += cilm[0, m:, m]
+        tmp *= 1.0 / _np.sqrt(2.0)
+        alm[ofs:ofs + lmax + 1 - m] = _np.conj(tmp)
+        ofs += lmax + 1 - m
+    return alm
+
+
+def _ccilm2almi(cilm):
+    lmax = cilm.shape[1] - 1
+    alm = _np.empty((_nalm(lmax, lmax),), dtype=_np.complex128)
+    fct = (-1) ** _np.arange(lmax + 1)
+    alm[0:lmax + 1] = cilm[0, :, 0].imag
+    ofs = lmax + 1
+    for m in range(1, lmax + 1):
+        tmp = _np.conj(cilm[1, m:, m])
+        tmp *= -fct[m]
+        tmp += cilm[0, m:, m]
+        tmp *= 1.0 / _np.sqrt(2.0)
+        alm[ofs:ofs + lmax + 1 - m] = tmp.imag + 1j * tmp.real
+        ofs += lmax + 1 - m
+    return alm
+
+
+def _addRealpart(cilm, alm):
+    lmax = cilm.shape[1] - 1
+    cilm[0, :, 0].real += alm[0:lmax + 1].real
+    ofs = lmax + 1
+    for m in range(1, lmax + 1):
+        tmp = alm[ofs:ofs + lmax + 1 - m] / _np.sqrt(2.0)
+        cilm[0, m:, m].real += tmp.real
+        cilm[0, m:, m].imag -= tmp.imag
+        if m & 1:
+            cilm[1, m:, m] -= tmp
+        else:
+            cilm[1, m:, m] += tmp
+        ofs += lmax + 1 - m
+    return cilm
+
+
+def _addImagpart(cilm, alm):
+    lmax = cilm.shape[1] - 1
+    cilm[0, :, 0].imag += alm[0:lmax + 1].real
+    ofs = lmax + 1
+    for m in range(1, lmax + 1):
+        tmp = alm[ofs:ofs + lmax + 1 - m] / _np.sqrt(2.0)
+        cilm[0, m:, m].real += tmp.imag
+        cilm[0, m:, m].imag += tmp.real
+        if m & 1:
+            cilm[1, m:, m].real += tmp.imag
+            cilm[1, m:, m].imag -= tmp.real
+        else:
+            cilm[1, m:, m].real -= tmp.imag
+            cilm[1, m:, m].imag += tmp.real
+        ofs += lmax + 1 - m
+    return cilm
+
+
+def _prep_lmax(lmax, lmax_calc, cilm):
+    if lmax is None:
+        lmax = cilm.shape[1] - 1
+    if lmax_calc is None:
+        lmax_calc = cilm.shape[1] - 1
+    if lmax_calc > lmax:
+        raise RuntimeError(
+            "lmax_calc ({}) must be less than or equal to lmax ({})"
+                .format(lmax_calc, lmax))
+    # lmax_calc need must not be higher than cilm.shape[1] - 1.
+    lmax_calc = min(cilm.shape[1] - 1, lmax_calc)
+    return lmax, lmax_calc, cilm[:, : lmax_calc + 1, : lmax_calc + 1]
+
+
+# dj_matrix is ignored
+def SHRotateRealCoef(rcoeffs, angles, dj_matrix=None):
+    lmax = rcoeffs.shape[1] - 1
+    alm = _make_alm(rcoeffs, lmax, 1, 1)
+    alm = ducc0.sht.rotate_alm(
+        alm, lmax, -angles[0], -angles[1], -angles[2], nthreads=nthreads
+    )
+    return _extract_alm(alm, lmax, 1, 1)
+
+
+# dj_matrix is ignored
+def SHRotateComplexCoef(ccoeffs, angles, dj_matrix=None):
+    lmax = ccoeffs.shape[1] - 1
+    alm = _ccilm2almr(ccoeffs)
+    alm = _apply_norm(alm, lmax, 1, 1, False)
+    alm = ducc0.sht.rotate_alm(
+        alm, lmax, -angles[0], -angles[1], -angles[2], nthreads=nthreads
+    )
+    alm = _apply_norm(alm, lmax, 1, 1, True)
+    res = _np.zeros((2, lmax + 1, lmax + 1), dtype=_np.complex128)
+    _addRealpart(res, alm)
+    alm = _ccilm2almi(ccoeffs)
+    alm = _apply_norm(alm, lmax, 1, 1, False)
+    alm = ducc0.sht.rotate_alm(
+        alm, lmax, -angles[0], -angles[1], -angles[2], nthreads=nthreads
+    )
+    alm = _apply_norm(alm, lmax, 1, 1, True)
+    _addImagpart(res, alm)
+    return res
+
+
+def MakeGridDH(
+    cilm,
+    lmax=None,
+    norm=1,
+    sampling=1,
+    csphase=1,
+    lmax_calc=None,
+    extend=False,
+):
+    lmax, lmax_calc, cilm = _prep_lmax(lmax, lmax_calc, cilm)
+    alm = _make_alm(cilm, lmax_calc, norm, csphase)
+    out = _np.empty(
+        [2 * lmax + 2 + extend, sampling * (2 * lmax + 2) + extend])
+    return _synthesize_DH(alm, lmax_calc, extend, out)
+
+
+def MakeGridDHC(
+    cilm,
+    lmax=None,
+    norm=1,
+    sampling=1,
+    csphase=1,
+    lmax_calc=None,
+    extend=False,
+):
+    lmax, lmax_calc, cilm = _prep_lmax(lmax, lmax_calc, cilm)
+    alm = _ccilm2almi(cilm)
+    alm = _apply_norm(alm, lmax, norm, csphase, False)
+    res = _np.empty(
+        [2 * lmax + 2 + extend, sampling * (2 * lmax + 2) + extend],
+        dtype=_np.complex128,
+    )
+    _synthesize_DH(alm, lmax_calc, extend, res.imag)
+    alm = _ccilm2almr(cilm)
+    alm = _apply_norm(alm, lmax, norm, csphase, False)
+    _synthesize_DH(alm, lmax_calc, extend, res.real)
+    return res
+
+
+def SHExpandDH(grid, norm=1, sampling=1, csphase=1, lmax_calc=None):
+    if grid.shape[1] != sampling * grid.shape[0]:
+        raise RuntimeError("grid resolution mismatch")
+    if lmax_calc is None:
+        lmax_calc = grid.shape[0] // 2 - 1
+    if lmax_calc > (grid.shape[0] // 2 - 1):
+        raise RuntimeError("lmax_calc too high")
+    alm = _analyze_DH(grid, lmax_calc)
+    return _extract_alm(alm, lmax_calc, norm, csphase)
+
+
+def SHExpandDHC(grid, norm=1, sampling=1, csphase=1, lmax_calc=None):
+    if grid.shape[1] != sampling * grid.shape[0]:
+        raise RuntimeError("grid resolution mismatch")
+    if lmax_calc is None:
+        lmax_calc = grid.shape[0] // 2 - 1
+    if lmax_calc > (grid.shape[0] // 2 - 1):
+        raise RuntimeError("lmax_calc too high")
+    lmax = grid.shape[0] // 2 - 1 if lmax_calc is None else lmax_calc
+    res = _np.zeros((2, lmax + 1, lmax + 1), dtype=_np.complex128)
+    alm = _analyze_DH(grid.real, lmax_calc)
+    alm = _apply_norm(alm, lmax, norm, csphase, True)
+    _addRealpart(res, alm)
+    alm = _analyze_DH(grid.imag, lmax_calc)
+    alm = _apply_norm(alm, lmax, norm, csphase, True)
+    _addImagpart(res, alm)
+    return res
+
+
+# zero is ignored (they are computed internally)
+def MakeGridGLQ(
+    cilm, lmax=None, zero=None, norm=1, csphase=1, lmax_calc=None, extend=False
+):
+    lmax, lmax_calc, cilm = _prep_lmax(lmax, lmax_calc, cilm)
+    alm = _make_alm(cilm, lmax_calc, norm, csphase)
+    out = _np.empty([lmax + 1, (2 * lmax + 1) + extend])
+    return _synthesize_GLQ(alm, lmax_calc, extend, out)
+
+
+# zero is ignored (they are computed internally)
+def MakeGridGLQC(
+    cilm, lmax=None, zero=None, norm=1, csphase=1, lmax_calc=None, extend=False
+):
+    lmax, lmax_calc, cilm = _prep_lmax(lmax, lmax_calc, cilm)
+    alm = _ccilm2almi(cilm)
+    alm = _apply_norm(alm, lmax, norm, csphase, False)
+    res = _np.empty([lmax + 1, 2 * lmax + 1 + extend], dtype=_np.complex128)
+    _synthesize_GLQ(alm, lmax_calc, extend, res.imag)
+    alm = _ccilm2almr(cilm)
+    alm = _apply_norm(alm, lmax, norm, csphase, False)
+    _synthesize_GLQ(alm, lmax_calc, extend, res.real)
+    return res
+
+
+# weights and zeros are ignored (they are computed internally)
+def SHExpandGLQ(
+    grid, weights=None, zeros=None, norm=1, csphase=1, lmax_calc=None
+):
+    if lmax_calc is None:
+        lmax_calc = grid.shape[0] - 1
+    if lmax_calc > (grid.shape[0] - 1):
+        raise RuntimeError("lmax_calc too high")
+    alm = _analyze_GLQ(grid, lmax_calc)
+    return _extract_alm(alm, lmax_calc, norm, csphase)
+
+
+# weights and zeros are ignored (they are computed internally)
+def SHExpandGLQC(
+    grid, weights=None, zeros=None, norm=1, csphase=1, lmax_calc=None
+):
+    if lmax_calc is None:
+        lmax_calc = grid.shape[0] - 1
+    if lmax_calc > (grid.shape[0] - 1):
+        raise RuntimeError("lmax_calc too high")
+    res = _np.zeros((2, lmax_calc + 1, lmax_calc + 1), dtype=_np.complex128)
+    alm = _analyze_GLQ(grid.real, lmax_calc)
+    alm = _apply_norm(alm, lmax_calc, norm, csphase, True)
+    _addRealpart(res, alm)
+    alm = _analyze_GLQ(grid.imag, lmax_calc)
+    alm = _apply_norm(alm, lmax_calc, norm, csphase, True)
+    _addImagpart(res, alm)
+    return res
+
+
+def MakeGradientDH(
+    cilm, lmax=None, sampling=1, lmax_calc=None, extend=False, radius=None
+):
+    lmax, lmax_calc, cilm = _prep_lmax(lmax, lmax_calc, cilm)
+    alm = _make_alm(cilm, lmax_calc, norm=1, csphase=1)
+    res = _np.empty(
+        (2, 2 * lmax + 2 + extend, sampling * (2 * lmax + 2) + extend)
+    )
+    res = _synthesize_DH_deriv1(alm, lmax_calc, extend, res)
+    if radius is not None:
+        res *= 1.0 / radius
+    return res
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 02ea4ac1e..5fdbd1c11 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -8,4 +8,5 @@ six
 shapely
 cartopy>=0.18.0
 pygmt>=0.3
+ducc0>=0.15
 flake8
diff --git a/setup.py b/setup.py
index b656c47a8..0a77d95a6 100644
--- a/setup.py
+++ b/setup.py
@@ -81,7 +81,10 @@
 ]
 
 EXTRAS_REQUIRE = {
-    'extras': ['cartopy>=0.18.0', 'pygmt>=0.2', 'palettable>=3.3']
+    'cartopy': ['cython', 'pyshp', 'six', 'shapely', 'cartopy>=0.18.0'],
+    'pygmt': ['pygmt==0.3'],
+    'palettable': ['palettable>=3.3'],
+    'ducc': ['ducc0>=0.15']
 }
 
 VERSION = versioneer.get_version()
diff --git a/src/pydoc/backend_module.md b/src/pydoc/backend_module.md
new file mode 100644
index 000000000..f577000da
--- /dev/null
+++ b/src/pydoc/backend_module.md
@@ -0,0 +1,21 @@
+# backend_module()
+
+Return the specified backend module used for the spherical harmonic
+transforms in pyshtools.
+
+# Usage
+
+```python
+module = backend_module([backend, nthreads])
+```
+
+# Parameters
+
+**backend : str, optional, default = preferred_backend()**
+:   Name of the preferred backend, either 'shtools' or 'ducc'.
+
+**nthreads : int, optional, default = 1**
+:   Number of threads to use for the 'ducc' backend. Setting this parameter
+        to 0 will use as many threads as there are hardware threads on the
+        system.
+    
\ No newline at end of file
diff --git a/src/pydoc/doc2md.py b/src/pydoc/doc2md.py
index 8bbae8b34..9f7abe404 100644
--- a/src/pydoc/doc2md.py
+++ b/src/pydoc/doc2md.py
@@ -24,7 +24,11 @@
     ['pysh.shio.YilmIndexVector', 'pyyilmindexvector.md'],
     ['pysh.spectralanalysis.cross_spectrum', 'cross_spectrum.md'],
     ['pysh.spectralanalysis.spectrum', 'spectrum.md'],
-    ['pysh.utils.figstyle', 'figstyle.md']
+    ['pysh.utils.figstyle', 'figstyle.md'],
+    ['pysh.backends.select_preferred_backend', 'select_preferred_backend.md'],
+    ['pysh.backends.preferred_backend', 'preferred_backend.md'],
+    ['pysh.backends.preferred_backend_module', 'preferred_backend_module.md'],
+    ['pysh.backends.backend_module', 'backend_module.md']
 ]
 
 for func, file in functions:
diff --git a/src/pydoc/preferred_backend.md b/src/pydoc/preferred_backend.md
new file mode 100644
index 000000000..61e43390c
--- /dev/null
+++ b/src/pydoc/preferred_backend.md
@@ -0,0 +1,11 @@
+# preferred_backend()
+
+Return the name of the current preferred backend used for the spherical
+harmonic transforms in pyshtools as a string.
+
+# Usage
+
+```python
+preferred_backend()
+```
+    
\ No newline at end of file
diff --git a/src/pydoc/preferred_backend_module.md b/src/pydoc/preferred_backend_module.md
new file mode 100644
index 000000000..6c807967f
--- /dev/null
+++ b/src/pydoc/preferred_backend_module.md
@@ -0,0 +1,11 @@
+# preferred_backend_module()
+
+Return a reference to the current preferred backend module used for the
+spherical harmonic transforms in pyshtools.
+
+# Usage
+
+```python
+module = preferred_backend_module()
+```
+    
\ No newline at end of file
diff --git a/src/pydoc/pyshread.md b/src/pydoc/pyshread.md
index 411b9c2d2..afb6d8819 100644
--- a/src/pydoc/pyshread.md
+++ b/src/pydoc/pyshread.md
@@ -6,8 +6,9 @@ Read shtools-formatted spherical harmonic coefficients from a text file.
 
 ```python
 coeffs, [errors], lmaxout, [header], [header2] = shread(
-    filename, [error=True, header=True, header2=True, lmax, skip])
+    filename, [error=True, header=True, header2=True, lmax, skip,
 ```
+        encoding])
 
 # Returns
 
@@ -55,6 +56,9 @@ coeffs, [errors], lmaxout, [header], [header2] = shread(
 **skip : int, optional, default = 0**
 :   The number of lines to skip before parsing the file.
 
+**encoding : str, optional, default = None**
+:   Encoding of the input file. The default is to use the system default.
+
 # Notes
 
 This function will read spherical harmonic coefficients from an
diff --git a/src/pydoc/read_dov.md b/src/pydoc/read_dov.md
index cc0298a00..c200fc20a 100644
--- a/src/pydoc/read_dov.md
+++ b/src/pydoc/read_dov.md
@@ -7,8 +7,9 @@ Read spherical harmonic coefficients from a text file formatted as
 
 ```python
 coeffs, [errors], lmaxout, [header], [header2] = read_dov(
-    filename, [error=True, header=True, header2=True, lmax, skip])
+    filename, [error=True, header=True, header2=True, lmax, skip,
 ```
+        encoding])
 
 # Returns
 
@@ -56,6 +57,9 @@ coeffs, [errors], lmaxout, [header], [header2] = read_dov(
 **skip : int, optional, default = 0**
 :   The number of lines to skip before parsing the file.
 
+**encoding : str, optional, default = None**
+:   Encoding of the input file. The default is to use the system default.
+
 # Notes
 
 This function will read spherical harmonic coefficients from a 'dov'-
diff --git a/src/pydoc/read_icgem_gfc.md b/src/pydoc/read_icgem_gfc.md
index d672c2d5f..562c153ff 100644
--- a/src/pydoc/read_icgem_gfc.md
+++ b/src/pydoc/read_icgem_gfc.md
@@ -7,8 +7,9 @@ file.
 
 ```python
 cilm, gm, r0, [errors] = read_icgem_gfc(filename,
-    [errors, lmax, epoch, encoding)
+    [errors, lmax, epoch, encoding,
 ```
+                                             quiet)
 
 # Returns
 
@@ -47,7 +48,40 @@ cilm, gm, r0, [errors] = read_icgem_gfc(filename,
         If the format of the file is 'icgem2.0' then the epoch must be
         specified.
 
-**encoding : str, optional**
-:   Encoding of the input file. Try to use 'iso-8859-1' if the default
-        (UTF-8) fails.
+**encoding : str, optional, default = None**
+:   Encoding of the input file. The default is to use the system default.
+        If the default encoding doesn't work, try 'iso-8859-1'.
+
+**quiet : bool, default = False**
+:   If True, suppress warnings about undefined keywords when reading the
+        file.
+
+# Notes
+
+This routine reads ICGEM formatted files of gravitational potential models
+and outputs arrays of the gravitational potential coefficients, errors, GM,
+and the reference radius. If epoch is specified, the coefficients will make
+use of the time variable terms in order to compute and return the potential
+coefficients for the specified epoch. Otherwise, the coefficients will be
+returned for the reference epoch of the model.
+
+Valid keys in the header section include:
+        modelname (not used)
+        product_type (only 'gravity_field' is allowed)
+        earth_gravity_constant or gravity_constant
+        radius
+        max_degree
+        errors ('unknown', 'formal', 'calibrated' or 'calibrated_and_formal')
+        tide_system (not used)
+        norm (not used)
+        format (either None or 'icgem2.0')
+
+Valid keys in the data section include:
+        gfc
+        gfct
+        trnd or dot
+        asin
+        acos
+
+Data lines starting with an unknown key are ignored.
     
\ No newline at end of file
diff --git a/src/pydoc/read_igrf.md b/src/pydoc/read_igrf.md
index 36fdeaabe..dc0c8ea57 100644
--- a/src/pydoc/read_igrf.md
+++ b/src/pydoc/read_igrf.md
@@ -25,6 +25,9 @@ read_igrm(filename, [year])
 **year : float, optional, default = 2020.**
 :   The year to compute the coefficients.
 
+**encoding : str, optional, default = None**
+:   Encoding of the input file. The default is to use the system default.
+
 # Notes
 
 The current International Geomagnetic Reference Field (IGRF-13) is a
diff --git a/src/pydoc/select_preferred_backend.md b/src/pydoc/select_preferred_backend.md
new file mode 100644
index 000000000..1d2f924a4
--- /dev/null
+++ b/src/pydoc/select_preferred_backend.md
@@ -0,0 +1,21 @@
+# select_preferred_backend()
+
+Select the preferred backend module used for the spherical harmonic
+transforms in pyshtools.
+
+# Usage
+
+```python
+select_preferred_backend_module([backend, nthreads])
+```
+
+# Parameters
+
+**backend : str, optional, default = 'shtools'**
+:   Name of the preferred backend, either 'shtools' or 'ducc'.
+
+**nthreads : int, optional, default = 1**
+:   Number of threads to use for the 'ducc' backend. Setting this parameter
+        to 0 will use as many threads as there are hardware threads on the
+        system.
+    
\ No newline at end of file
diff --git a/src/pydoc/shwrite.md b/src/pydoc/shwrite.md
index 23894b99f..cf0f5348d 100644
--- a/src/pydoc/shwrite.md
+++ b/src/pydoc/shwrite.md
@@ -5,7 +5,7 @@ Write shtools-formatted spherical harmonic coefficients to a text file.
 # Usage
 
 ```python
-shwrite(filename, coeffs, [errors, header, header2, lmax])
+shwrite(filename, coeffs, [errors, header, header2, lmax, encoding])
 ```
 
 # Parameters
@@ -32,6 +32,9 @@ shwrite(filename, coeffs, [errors, header, header2, lmax])
 **lmax : int, optional, default = None**
 :   The maximum spherical harmonic degree to write to the file.
 
+**encoding : str, optional, default = None**
+:   Encoding of the output file. The default is to use the system default.
+
 # Notes
 
 This function will write spherical harmonic coefficients (and optionally
diff --git a/src/pydoc/write_dov.md b/src/pydoc/write_dov.md
index e86671b78..e2ec9c4bf 100644
--- a/src/pydoc/write_dov.md
+++ b/src/pydoc/write_dov.md
@@ -6,7 +6,7 @@ Write spherical harmonic coefficients to a text file formatted as
 # Usage
 
 ```python
-write_dov(filename, coeffs, [errors, header, header2, lmax])
+write_dov(filename, coeffs, [errors, header, header2, lmax, encoding])
 ```
 
 # Parameters
@@ -33,6 +33,9 @@ write_dov(filename, coeffs, [errors, header, header2, lmax])
 **lmax : int, optional, default = None**
 :   The maximum spherical harmonic degree to write to the file.
 
+**encoding : str, optional, default = None**
+:   Encoding of the output file. The default is to use the system default.
+
 # Notes
 
 This function will write spherical harmonic coefficients (and optionally
diff --git a/src/pydoc/write_icgem_gfc.md b/src/pydoc/write_icgem_gfc.md
index b0c9e7439..029b608ba 100644
--- a/src/pydoc/write_icgem_gfc.md
+++ b/src/pydoc/write_icgem_gfc.md
@@ -7,8 +7,9 @@ file.
 
 ```python
 write_icgem_gfc(filename, coeffs, [errors, header, lmax, modelname, gm, r0,
-    product_type, earth_gm, error_kind, tide_system, normalization, format)
+    product_type, earth_gm, error_kind, tide_system, normalization, format,
 ```
+        encoding)
 
 # Parameters
 
@@ -46,7 +47,7 @@ write_icgem_gfc(filename, coeffs, [errors, header, lmax, modelname, gm, r0,
 :   Reference radius of the model, in meters.
 
 **error_kind : str, optional, default = None**
-:   Which errors to read. Can be either 'unknown', 'calibrated', or
+:   Which errors to write. Can be either 'unknown', 'calibrated', or
         'formal'.
 
 **tide_system : str, optional, default = 'unknown'**
@@ -58,4 +59,7 @@ write_icgem_gfc(filename, coeffs, [errors, header, lmax, modelname, gm, r0,
 
 **format : str, optional, default = None**
 :   The format of the ICGEM spherical harmonic coefficients.
+
+**encoding : str, optional, default = None**
+:   Encoding of the output file. The default is to use the system default.
     
\ No newline at end of file


Test patch:
diff --git a/docs/images/backend-timing-tests.png b/docs/images/backend-timing-tests.png
new file mode 100644
index 000000000..419469b29
Binary files /dev/null and b/docs/images/backend-timing-tests.png differ
diff --git a/examples/python/ducctest.py b/examples/python/ducctest.py
new file mode 100644
index 000000000..6a27a94fe
--- /dev/null
+++ b/examples/python/ducctest.py
@@ -0,0 +1,281 @@
+import numpy as np
+import pyshtools as pysh
+from time import time
+
+
+nthreads = 1
+
+
+def _l2error(a, b):
+    return np.sqrt(np.sum(np.abs(a - b) ** 2) / np.sum(np.abs(a) ** 2))
+
+
+# force SHTOOLS to deallocate temporary buffers
+def flush_buffers(grd):
+    degrees = np.arange(1, dtype=float)
+    degrees[0] = np.inf
+    power = degrees ** (-2)
+
+    pysh.backends.select_preferred_backend("shtools")
+    clm = pysh.SHCoeffs.from_random(power, seed=12345)
+    grid2 = clm.expand(grid=grd)
+    _ = grid2.expand()
+    clm = pysh.SHCoeffs.from_random(power, seed=12345, kind="complex")
+    grid2 = clm.expand(grid=grd)
+    _ = grid2.expand()
+
+
+def test_SHT(lmax, grd, csphase, normalization, extend):
+    degrees = np.arange(lmax + 1, dtype=float)
+    degrees[0] = np.inf
+    power = degrees ** (-2)
+
+    clm = pysh.SHCoeffs.from_random(power, seed=12345)
+
+    clm = clm.convert(normalization=normalization, csphase=csphase, lmax=lmax)
+
+    pysh.backends.select_preferred_backend("ducc", nthreads=nthreads)
+    t0 = time()
+    grid = clm.expand(grid=grd, extend=extend)
+    cilm = grid.expand()
+    tducc = time() - t0
+    pysh.backends.select_preferred_backend("shtools")
+    t0 = time()
+    grid2 = clm.expand(grid=grd, extend=extend)
+    cilm2 = grid2.expand()
+    tshtools = time() - t0
+
+    flush_buffers(grd)
+
+    return (
+        _l2error(grid.to_array(), grid2.to_array())
+        + _l2error(cilm.to_array(), cilm2.to_array()),
+        tshtools / tducc,
+    )
+
+
+def test_SHTducc(lmax, grd, nthreads):
+    degrees = np.arange(lmax + 1, dtype=float)
+    degrees[0] = np.inf
+    power = degrees ** (-2)
+
+    clm = pysh.SHCoeffs.from_random(power, seed=12345)
+
+    clm = clm.convert(normalization="ortho", csphase=1, lmax=lmax)
+
+    pysh.backends.select_preferred_backend("ducc", nthreads=nthreads)
+    t0 = time()
+    grid = clm.expand(grid=grd, extend=False)
+    cilm = grid.expand(normalization="ortho", csphase=1)
+    tducc = time() - t0
+
+    return _l2error(clm.to_array(), cilm.to_array()), tducc
+
+
+def test_SHTC(lmax, grd, csphase, normalization, extend):
+    degrees = np.arange(lmax + 1, dtype=float)
+    degrees[0] = np.inf
+    power = degrees ** (-2)
+
+    clm = pysh.SHCoeffs.from_random(power, seed=12345, kind="complex")
+
+    clm = clm.convert(normalization=normalization, csphase=csphase, lmax=lmax)
+
+    pysh.backends.select_preferred_backend("ducc", nthreads=nthreads)
+    t0 = time()
+    grid = clm.expand(grid=grd, extend=extend)
+    cilm = grid.expand(normalization=normalization, csphase=csphase)
+    tducc = time() - t0
+    pysh.backends.select_preferred_backend("shtools")
+    t0 = time()
+    grid2 = clm.expand(grid=grd, extend=extend)
+    cilm2 = grid2.expand(normalization=normalization, csphase=csphase)
+    tshtools = time() - t0
+
+    flush_buffers(grd)
+
+    return (
+        _l2error(grid.to_array(), grid2.to_array())
+        + _l2error(cilm.to_array(), cilm2.to_array()),
+        tshtools / tducc,
+    )
+
+
+def test_SHT_deriv(lmax, grd, csphase, extend):
+    degrees = np.arange(lmax + 1, dtype=float)
+    degrees[0] = 1.0
+    power = degrees ** (-2)
+
+    clm = pysh.SHCoeffs.from_random(power, seed=12345)
+
+    clm = clm.convert(csphase=csphase, lmax=lmax)
+
+    pysh.backends.select_preferred_backend("ducc", nthreads=nthreads)
+    t0 = time()
+    grad = clm.gradient(extend=extend, radius=3.4)
+    tducc = time() - t0
+    pysh.backends.select_preferred_backend("shtools")
+    t0 = time()
+    grad2 = clm.gradient(extend=extend, radius=1.0)
+    tshtools = time() - t0
+
+    flush_buffers(grd)
+
+    return (
+        _l2error(3.4 * grad.phi.to_array(), grad2.phi.to_array())
+        + _l2error(3.4 * grad.theta.to_array(), grad2.theta.to_array()),
+        tshtools / tducc,
+    )
+
+
+def test_rot(lmax, alpha, beta, gamma):
+    degrees = np.arange(lmax + 1, dtype=float)
+    degrees[0] = np.inf
+    power = degrees ** (-2)
+
+    clm = pysh.SHCoeffs.from_random(power, seed=12345)
+
+    pysh.backends.select_preferred_backend("ducc", nthreads=nthreads)
+    t0 = time()
+    clm_rotated = clm.rotate(alpha, beta, gamma, degrees=True)
+    tducc = time() - t0
+    pysh.backends.select_preferred_backend("shtools")
+    t0 = time()
+    clm_rotated2 = clm.rotate(alpha, beta, gamma, degrees=True)
+    tshtools = time() - t0
+
+    return (
+        _l2error(clm_rotated.to_array(), clm_rotated2.to_array()),
+        tshtools / tducc,
+    )
+
+
+def test_rotc(lmax, alpha, beta, gamma):
+    degrees = np.arange(lmax + 1, dtype=float)
+    degrees[0] = np.inf
+    power = degrees ** (-2)
+
+    clm = pysh.SHCoeffs.from_random(power, seed=12345, kind="complex")
+
+    pysh.backends.select_preferred_backend("ducc", nthreads=nthreads)
+    t0 = time()
+    clm_rotated = clm.rotate(alpha, beta, gamma, degrees=True)
+    tducc = time() - t0
+    pysh.backends.select_preferred_backend("shtools")
+    t0 = time()
+    clm_rotated2 = clm.rotate(alpha, beta, gamma, degrees=True)
+    tshtools = time() - t0
+
+    return (
+        _l2error(clm_rotated.to_array(), clm_rotated2.to_array()),
+        tshtools / tducc,
+    )
+
+
+def test_rot2(lmax, alpha, beta, gamma):
+    degrees = np.arange(lmax + 1, dtype=float)
+    degrees[0] = np.inf
+    power = degrees ** (-2)
+
+    clm = pysh.SHCoeffs.from_random(power, seed=12345)
+
+    pysh.backends.select_preferred_backend("ducc", nthreads=nthreads)
+    t0 = time()
+    clm_rotated = clm.rotate(alpha, beta, gamma, degrees=True)
+    clm_rotated = clm_rotated.rotate(-gamma, -beta, -alpha, degrees=True)
+    tducc = time() - t0
+
+    return _l2error(clm.to_array(), clm_rotated.to_array()), tducc
+
+
+lmax_list = [127, 255, 511, 1023]
+
+print("SHRealCoeff rotation tests:")
+for lmax in lmax_list:
+    for alpha in [47]:
+        res = test_rot(lmax, alpha, 27, 59)
+        print(
+            "lmax={:4}: L2 error={:e}, speedup factor={:f}".format(
+                lmax, res[0], res[1]
+            )
+        )
+
+print("SHComplexCoeff rotation tests:")
+for lmax in lmax_list:
+    for alpha in [47]:
+        res = test_rotc(lmax, alpha, 27, 59)
+        print(
+            "lmax={:4}: L2 error={:e}, speedup factor={:f}".format(
+                lmax, res[0], res[1]
+            )
+        )
+
+lmax_list = [80]
+
+print("SHT tests unnorm:")
+for grid in ["GLQ", "DH", "DH2"]:
+    for csphase in [-1, 1]:
+        for norm in ["unnorm"]:
+            for extend in [True, False]:
+                for lmax in [5, 10, 20, 85]:
+                    res = test_SHT(lmax, grid, csphase, norm, extend)
+                    print(
+                        "{:3}, CS={:2}, norm={:7}, extend={:5}, lmax={:4}:  "
+                        "L2 error={:e}, speedup factor={:f}".format(
+                            grid, csphase, norm, extend, lmax, res[0], res[1]
+                        )
+                    )
+
+lmax_list = [127, 255, 511, 1023, 2047]
+
+print("SHT tests:")
+for grid in ["GLQ", "DH", "DH2"]:
+    for csphase in [-1, 1]:
+        for norm in ["ortho", "4pi", "schmidt"]:
+            for extend in [True, False]:
+                for lmax in lmax_list:
+                    res = test_SHT(lmax, grid, csphase, norm, extend)
+                    print(
+                        "{:3}, CS={:2}, norm={:7}, extend={:5}, lmax={:4}:  "
+                        "L2 error={:e}, speedup factor={:f}".format(
+                            grid, csphase, norm, extend, lmax, res[0], res[1]
+                        )
+                    )
+print("SHTC tests:")
+for grid in ["GLQ", "DH", "DH2"]:
+    for csphase in [-1, 1]:
+        for norm in ["ortho", "4pi", "schmidt"]:
+            for extend in [True, False]:
+                for lmax in lmax_list:
+                    res = test_SHTC(lmax, grid, csphase, norm, extend)
+                    print(
+                        "{:3}, CS={:2}, norm={:7}, extend={:5}, lmax={:4}:  "
+                        "L2 error={:e}, speedup factor={:f}".format(
+                            grid, csphase, norm, extend, lmax, res[0], res[1]
+                        )
+                    )
+print("SHT gradient tests:")
+for grid in ["DH", "DH2"]:
+    for csphase in [-1, 1]:
+        for extend in [True, False]:
+            for lmax in lmax_list:
+                res = test_SHT_deriv(lmax, grid, csphase, extend)
+                print(
+                    "{:3}, CS={:2}, extend={:5}, lmax={:4}:  L2 error={:e}, "
+                    "speedup factor={:f}".format(
+                        grid, csphase, extend, lmax, res[0], res[1]
+                    )
+                )
+
+print("DUCC: forward/backward rotation with high band limits:")
+for lmax in [4095]:
+    for alpha in [47]:
+        res = test_rot2(lmax, alpha, 27, 59)
+        print(
+            "lmax={:4}: L2 error={:e}, time={:f}".format(lmax, res[0], res[1])
+        )
+
+print("DUCC: forward/backward SHT with high band limits:")
+for lmax in [8191]:
+    res = test_SHTducc(lmax, "GLQ", nthreads=nthreads)
+    print("lmax={:4}: L2 error={:e}, time={:f}".format(lmax, res[0], res[1]))