diff --git a/buildscripts/build_shtools.sh b/buildscripts/build_shtools.sh index 01964da7b..7e82b8846 100755 --- a/buildscripts/build_shtools.sh +++ b/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" diff --git a/environment.yml b/environment.yml index 34c6c6c3a..933cb3fe1 100644 --- a/environment.yml +++ b/environment.yml @@ -21,4 +21,3 @@ dependencies: - x264>=1!157.20191217 - ffmpeg>=4.3.2 - cython>=3.0.0 - - pyshtools diff --git a/pyproject.toml b/pyproject.toml index 413681303..5fde1057a 100644 --- a/pyproject.toml +++ b/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" diff --git a/requirements.txt b/requirements.txt index 9519f962d..79d85b58b 100644 --- a/requirements.txt +++ b/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 \ No newline at end of file +cython>=3.0.0 \ No newline at end of file diff --git a/swiftest/sph_harmonics.py b/swiftest/sph_harmonics.py index 2f192e86d..140d35e2b 100644 --- a/swiftest/sph_harmonics.py +++ b/swiftest/sph_harmonics.py @@ -20,7 +20,6 @@ from astropy.coordinates import SkyCoord import datetime import xarray as xr -import pyshtools as pysh from typing import ( Literal, Dict, @@ -28,121 +27,136 @@ 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.") + \ No newline at end of file