Improved the robustness of the build by making pyshtools an optional …

…dependency

buildscripts/build_shtools.sh

@@ -15,9 +15,7 @@ ARGS=$@
 . ${SCRIPT_DIR}/_build_getopts.sh ${ARGS}
 . ${SCRIPT_DIR}/set_compilers.sh
 
-
-SHTOOLS_VER="4.11.10"
-
+SHTOOLS_VER="4.9.1"
 
 printf "*********************************************************\n"
 printf "*             FETCHING SHTOOLS SOURCE                      *\n"

environment.yml

@@ -21,4 +21,3 @@ dependencies:
   - x264>=1!157.20191217
   - ffmpeg>=4.3.2
   - cython>=3.0.0
-  - pyshtools

pyproject.toml

@@ -24,8 +24,8 @@ classifiers=[
 ]
 keywords=['astronomy','astrophysics', 'planetary', 'n-body',  'integrator', 'symplectic', 'wisdom-holman', 'symba']
 dependencies = [
-    'numpy>=1.26.3',
-    'scipy>=1.10.1',
+    'numpy>=1.26',
+    'scipy>=1.10',
     'xarray>=2023.1',
     'dask>=2023.5',
     'distributed>=2023.5',
@@ -38,33 +38,32 @@ dependencies = [
     'astroquery>=0.4.6',
     'tqdm>=4.64',
     'cython>=3.0',
-    'meson>=1.3',
-    'meson-python>=0.15',
-    'setuptools_scm',
-    'pyshtools>=4.11.10'
 ]
 
+[project.optional-dependencies]
+pyshtools = ["pyshtools"] 
+
 [project.urls]
 Repository = 'https://github.itap.purdue.edu/MintonGroup/swiftest'
 
 [build-system]
 requires = [
     "scikit-build-core",
-    "cython>=3.0.0",
+    "cython>=3.0",
     "cmake",
     "pyproject_metadata",
     "pytest",
     "pathspec",
     "sphinx",
     "sphinx-autosummary-accessors",
-    "sphinx-book-theme >= 0.3.0",
+    "sphinx-book-theme >= 0.3",
     "sphinx-copybutton",
     "sphinx-design",
     "sphinx-inline-tabs",
     "sphinxext-rediraffe",
     "sphinxext-opengraph",
     "nbsphinx",
-    "ford"
+    "ford",
 ]
 build-backend = "scikit_build_core.build"
 

requirements.txt

@@ -11,5 +11,4 @@ matplotlib>=3.7.1
 astropy>=5.1
 astroquery>=0.4.6
 tqdm>=4.65.0
-cython>=3.0.0
-pyshtools
+cython>=3.0.0

swiftest/sph_harmonics.py

@@ -20,129 +20,143 @@
 from astropy.coordinates import SkyCoord
 import datetime
 import xarray as xr
-import pyshtools as pysh
 from typing import (
     Literal,
     Dict,
     List,
     Any
 )
 
-class Sph_Harmonics(object):
-    def clm_from_ellipsoid(mass, density, a, b = None, c = None, lmax = 6, lref_radius = False, ref_radius = None):
-        """
-        Creates and returns the gravity coefficients for an ellipsoid with principal axes a, b, c upto a certain maximum degree lmax. 
-        Uses pyshtools. No units necessary for a, b, & c. However, they need to be in the same units (DU).
-
-        Parameters
-        ----------
-        mass : float
-            mass of the central body
-        density : float
-            density of the central body
-        a : float 
-            length of the pricipal axis aligned with the x axis
-        b : float, optional, default = a
-            length of the pricipal axis aligned with the y axis
-        c : float, optional, default = b
-            length of the pricipal axis aligned with the z axis
-        lmax : int, optional, default = 6
-            The maximum spherical harmonic degree resolvable by the grid.
-        lref_radius : boolean, optional, default = False
-            Boolean value to return the reference radius calculated by SHTOOLS
-        ref_radius : float, optional, default = None
-            Reference radius to scale the gravitational coefficients to
-
-        Returns
-        -------
-        clm : ndarry, shape (2, lmax+1, lmax+1)
-            numpy ndarray of the gravitational potential spherical harmonic coefficients. 
-            This is a three-dimensional array formatted as coeffs[i, degree, order], 
-            where i=0 corresponds to positive orders and i=1 to negative orders.
-
-        """
-        Gmass = swiftest.constants.GC * mass # SHTOOLS uses an SI G value, and divides it before using the mass; NO NEED TO CHANGE UNITS
-
-        # cap lmax to ensure fast performance without giving up accuracy
-        lmax_limit = 6              # lmax_limit = 6 derived from Jean's Law; characteristic wavelength = the radius of the CB
-        if(lmax > lmax_limit):                           
-            lmax = lmax_limit
-            print(f'Setting maximum spherical harmonic degree to {lmax_limit}')
-
-        # create shape grid 
-        shape_SH = pysh.SHGrid.from_ellipsoid(lmax = lmax, a = a, b = b, c = c)
-
-        # get gravity coefficients
-        clm_class = pysh.SHGravCoeffs.from_shape(shape_SH, rho = density, gm = Gmass) # 4pi normalization
-        clm = clm_class.to_array(normalization = '4pi') # export as array with 4pi normalization and not scaling by 4*pi to match normalisation
-
-        # Return reference radius EQUALS the radius of the Central Body
-        print(f'Ensure that the Central Body radius equals the reference radius.')
-
-        if(lref_radius == True and ref_radius is None):
-            ref_radius = shape_SH.expand(normalization = '4pi').coeffs[0, 0, 0]
-            return clm, ref_radius
-        elif(lref_radius == True and ref_radius is not None):
-            clm_class = clm_class.change_ref(r0 = ref_radius)
-            clm = clm_class.to_array(normalization = '4pi')
-            return clm, ref_radius
-        else:
-            return clm
-
-    def clm_from_relief(mass, density, grid, lmax = 6, lref_radius = False, ref_radius = None):
-        """
-        Creates and returns the gravity coefficients for a body with a given DH grid upto a certain maximum degree lmax. 
-        Uses pyshtools.
-
-        Parameters
-        ----------
-        mass : float
-            mass of the central body
-        density : float
-            density of the central body
-        grid : array, shape []
-            DH grid of the surface relief of the body
-        lmax : int, optional, default = 6
-            The maximum spherical harmonic degree resolvable by the grid.
-        lref_radius : boolean, optional, default = False
-            Boolean value to return the reference radius calculated by SHTOOLS
-        ref_radius : float, optional, default = None
-            Reference radius to scale the gravitational coefficients to
-
-        Returns
-        -------
-        clm : ndarry, shape (2, lmax+1, lmax+1)
-            numpy ndarray of the gravitational potential spherical harmonic coefficients. 
-            This is a three-dimensional array formatted as coeffs[i, degree, order], 
-            where i=0 corresponds to positive orders and i=1 to negative orders.
-
-        """
-
-        Gmass = swiftest.constants.GC * mass # SHTOOLS uses an SI G value, and divides it before using the mass; NO NEED TO CHANGE UNITS
-
-        # cap lmax to 20 to ensure fast performance
-        lmax_limit = 6
-        if(lmax > lmax_limit): # FIND A BETTER WAY to judge this cut off point, i.e., relative change between coefficients
-            lmax = lmax_limit
-            print(f'Setting maximum spherical harmonic degree to {lmax_limit}')
-
-        # convert to spherical harmonics
-        shape_SH = pysh.SHGrid.from_array(grid)
-
-        # get coefficients
-        clm_class = pysh.SHGravcoeffs.from_shape(shape_SH, rho = density, gm = Gmass) # 4pi normalization
-        clm = clm_class.to_array(normalization = '4pi') # export as array with 4pi normalization
-
-        # Return reference radius EQUALS the radius of the Central Body
-
-        print(f'Ensure that the Central Body radius equals the reference radius.')
-
-        if(lref_radius == True and ref_radius is None):
-            ref_radius = shape_SH.expand(normalization = '4pi').coeffs[0, 0, 0]
-            return clm, ref_radius
-        elif(lref_radius == True and ref_radius is not None):
-            clm_class = clm_class.change_ref(r0 = ref_radius)
-            clm = clm_class.to_array(normalization = '4pi')
-            return clm, ref_radius
-        else:
-            return clm
+try:
+    import pyshtools as pysh
+    PYSHTOOLS_AVAILABLE = True
+except ModuleNotFoundError:
+    PYSHTOOLS_AVAILABLE = False
+    print("pyshtools is not installed. Some features will be unavailable.")
+
+if PYSHTOOLS_AVAILABLE:
+
+    class Sph_Harmonics(object):
+        def clm_from_ellipsoid(mass, density, a, b = None, c = None, lmax = 6, lref_radius = False, ref_radius = None):
+            """
+            Creates and returns the gravity coefficients for an ellipsoid with principal axes a, b, c upto a certain maximum degree lmax. 
+            Uses pyshtools. No units necessary for a, b, & c. However, they need to be in the same units (DU).
+
+            Parameters
+            ----------
+            mass : float
+                mass of the central body
+            density : float
+                density of the central body
+            a : float 
+                length of the pricipal axis aligned with the x axis
+            b : float, optional, default = a
+                length of the pricipal axis aligned with the y axis
+            c : float, optional, default = b
+                length of the pricipal axis aligned with the z axis
+            lmax : int, optional, default = 6
+                The maximum spherical harmonic degree resolvable by the grid.
+            lref_radius : boolean, optional, default = False
+                Boolean value to return the reference radius calculated by SHTOOLS
+            ref_radius : float, optional, default = None
+                Reference radius to scale the gravitational coefficients to
+
+            Returns
+            -------
+            clm : ndarry, shape (2, lmax+1, lmax+1)
+                numpy ndarray of the gravitational potential spherical harmonic coefficients. 
+                This is a three-dimensional array formatted as coeffs[i, degree, order], 
+                where i=0 corresponds to positive orders and i=1 to negative orders.
+
+            """
+            Gmass = swiftest.constants.GC * mass # SHTOOLS uses an SI G value, and divides it before using the mass; NO NEED TO CHANGE UNITS
+
+            # cap lmax to ensure fast performance without giving up accuracy
+            lmax_limit = 6              # lmax_limit = 6 derived from Jean's Law; characteristic wavelength = the radius of the CB
+            if(lmax > lmax_limit):                           
+                lmax = lmax_limit
+                print(f'Setting maximum spherical harmonic degree to {lmax_limit}')
+
+            # create shape grid 
+            shape_SH = pysh.SHGrid.from_ellipsoid(lmax = lmax, a = a, b = b, c = c)
+
+            # get gravity coefficients
+            clm_class = pysh.SHGravCoeffs.from_shape(shape_SH, rho = density, gm = Gmass) # 4pi normalization
+            clm = clm_class.to_array(normalization = '4pi') # export as array with 4pi normalization and not scaling by 4*pi to match normalisation
+
+            # Return reference radius EQUALS the radius of the Central Body
+            print(f'Ensure that the Central Body radius equals the reference radius.')
+
+            if(lref_radius == True and ref_radius is None):
+                ref_radius = shape_SH.expand(normalization = '4pi').coeffs[0, 0, 0]
+                return clm, ref_radius
+            elif(lref_radius == True and ref_radius is not None):
+                clm_class = clm_class.change_ref(r0 = ref_radius)
+                clm = clm_class.to_array(normalization = '4pi')
+                return clm, ref_radius
+            else:
+                return clm
+
+        def clm_from_relief(mass, density, grid, lmax = 6, lref_radius = False, ref_radius = None):
+            """
+            Creates and returns the gravity coefficients for a body with a given DH grid upto a certain maximum degree lmax. 
+            Uses pyshtools.
+
+            Parameters
+            ----------
+            mass : float
+                mass of the central body
+            density : float
+                density of the central body
+            grid : array, shape []
+                DH grid of the surface relief of the body
+            lmax : int, optional, default = 6
+                The maximum spherical harmonic degree resolvable by the grid.
+            lref_radius : boolean, optional, default = False
+                Boolean value to return the reference radius calculated by SHTOOLS
+            ref_radius : float, optional, default = None
+                Reference radius to scale the gravitational coefficients to
+
+            Returns
+            -------
+            clm : ndarry, shape (2, lmax+1, lmax+1)
+                numpy ndarray of the gravitational potential spherical harmonic coefficients. 
+                This is a three-dimensional array formatted as coeffs[i, degree, order], 
+                where i=0 corresponds to positive orders and i=1 to negative orders.
+
+            """
+
+            Gmass = swiftest.constants.GC * mass # SHTOOLS uses an SI G value, and divides it before using the mass; NO NEED TO CHANGE UNITS
+
+            # cap lmax to 20 to ensure fast performance
+            lmax_limit = 6
+            if(lmax > lmax_limit): # FIND A BETTER WAY to judge this cut off point, i.e., relative change between coefficients
+                lmax = lmax_limit
+                print(f'Setting maximum spherical harmonic degree to {lmax_limit}')
+
+            # convert to spherical harmonics
+            shape_SH = pysh.SHGrid.from_array(grid)
+
+            # get coefficients
+            clm_class = pysh.SHGravcoeffs.from_shape(shape_SH, rho = density, gm = Gmass) # 4pi normalization
+            clm = clm_class.to_array(normalization = '4pi') # export as array with 4pi normalization
+
+            # Return reference radius EQUALS the radius of the Central Body
+
+            print(f'Ensure that the Central Body radius equals the reference radius.')
+
+            if(lref_radius == True and ref_radius is None):
+                ref_radius = shape_SH.expand(normalization = '4pi').coeffs[0, 0, 0]
+                return clm, ref_radius
+            elif(lref_radius == True and ref_radius is not None):
+                clm_class = clm_class.change_ref(r0 = ref_radius)
+                clm = clm_class.to_array(normalization = '4pi')
+                return clm, ref_radius
+            else:
+                return clm
+
+else:
+    class Sph_Harmonics(object):
+        def clm_from_ellipsoid(*args, **kwargs):
+            raise NotImplementedError("Sph_Harmonics is not available because pyshtools is not installed.")
+