0.3.1

Displacement_strain_planet/Displacement_strain.py

@@ -129,26 +129,27 @@ def SH_deriv_store(
 ):
     """
     Compute and store or load spherical harmonic derivatives
-    over the entire sphere (first and second order).
+    over the entire sphere (first and second order). The spherical
+    harmonic degree and order correspond to the index l*(l+1)/2+m
 
     Returns
     -------
-    Y_lm_d1_theta_a : array, size(2,lmax+1,lmax+1)
+    Y_lm_d1_theta_a : array, size(2,(lmax+1)*(lmax+2)/2)
         Array with the first derivative
         of Legendre polynomials with respect to colatitude.
-    Y_lm_d1_phi_a : array, size(2,lmax+1,lmax+1)
+    Y_lm_d1_phi_a : array, size(2,(lmax+1)*(lmax+2)/2)
         Array with the first derivative
         of Legendre polynomials with respect to longitude.
-    Y_lm_d2_theta_a : array, size(2,lmax+1,lmax+1)
+    Y_lm_d2_theta_a : array, size(2,(lmax+1)*(lmax+2)/2)
         Array with the second derivative
         of Legendre polynomials with respect to colatitude.
-    Y_lm_d2_phi_a : array, size(2,lmax+1,lmax+1)
+    Y_lm_d2_phi_a : array, size(2,(lmax+1)*(lmax+2)/2)
         Array with the second derivative
         of Legendre polynomials with respect to longitude.
-    Y_lm_d2_thetaphi_a : array, size(2,lmax+1,lmax+1)
+    Y_lm_d2_thetaphi_a : array, size(2,(lmax+1)*(lmax+2)/2)
         Array with the first derivative
         of Legendre polynomials with respect to colatitude and longitude.
-    y_lm_save : array, size(2,lmax+1,lmax+1)
+    y_lm_save : array, size(2,(lmax+1)*(lmax+2)/2)
         Array of spherical harmonic functions.
 
     Parameters
@@ -216,15 +217,15 @@ def SH_deriv_store(
             )
             print("dtype is %s." % (dtype))
         index_size = int((lmax + 1) * (lmax + 2) / 2)
-        shape_save = (nlat, nlon, 2, lmax + 1, lmax + 1)
+        shape_save = (nlat, nlon, 2, index_size)
         Y_lm_d1_theta_a = np.zeros(shape_save, dtype=dtype)
         Y_lm_d1_phi_a = np.zeros(shape_save, dtype=dtype)
         Y_lm_d2_phi_a = np.zeros(shape_save, dtype=dtype)
         Y_lm_d2_thetaphi_a = np.zeros(shape_save, dtype=dtype)
         Y_lm_d2_theta_a = np.zeros(shape_save, dtype=dtype)
         y_lm_save = np.zeros(shape_save, dtype=dtype)
         phi_ar = np.linspace(0, 360, nlon, endpoint=False) * pi / 180.0
-        y_lm = np.zeros((len(phi_ar), 2, lmax + 1, lmax + 1))
+        y_lm = np.zeros((len(phi_ar), 2, index_size))
 
         t_i = -1
         for theta in np.linspace(0, 180, nlat, endpoint=False) * pi / 180.0:
@@ -257,36 +258,37 @@ def SH_deriv_store(
                         m2cosphi = -(m ** 2) * cosmphi  # Second
                         # cos(m*phi)
                         # derivative
-                        Y_lm_d1_theta_a[t_i, :, 0, l, m] = dp_theta[index] * cosmphi
-                        Y_lm_d1_phi_a[t_i, :, 0, l, m] = p_theta[index] * msinmphi
-                        Y_lm_d2_phi_a[t_i, :, 0, l, m] = p_theta[index] * m2cosphi
-                        Y_lm_d2_thetaphi_a[t_i, :, 0, l, m] = dp_theta[index] * msinmphi
-                        y_lm[:, 0, l, m] = p_theta[index] * cosmphi
+                        Y_lm_d1_theta_a[t_i, :, 0, index] = dp_theta[index] * cosmphi
+                        Y_lm_d1_phi_a[t_i, :, 0, index] = p_theta[index] * msinmphi
+                        Y_lm_d2_phi_a[t_i, :, 0, index] = p_theta[index] * m2cosphi
+                        Y_lm_d2_thetaphi_a[t_i, :, 0, index] = (
+                            dp_theta[index] * msinmphi
+                        )
+                        y_lm[:, 0, index] = p_theta[index] * cosmphi
                     else:
                         mcosmphi = m_abs * cosmphi
                         m2sinphi = -(m_abs ** 2) * sinmphi
-                        Y_lm_d1_theta_a[t_i, :, 1, l, m_abs] = dp_theta[index] * sinmphi
-                        Y_lm_d1_phi_a[t_i, :, 1, l, m_abs] = p_theta[index] * mcosmphi
-                        Y_lm_d2_phi_a[t_i, :, 1, l, m_abs] = p_theta[index] * m2sinphi
-                        Y_lm_d2_thetaphi_a[t_i, :, 1, l, m_abs] = (
+                        Y_lm_d1_theta_a[t_i, :, 1, index] = dp_theta[index] * sinmphi
+                        Y_lm_d1_phi_a[t_i, :, 1, index] = p_theta[index] * mcosmphi
+                        Y_lm_d2_phi_a[t_i, :, 1, index] = p_theta[index] * m2sinphi
+                        Y_lm_d2_thetaphi_a[t_i, :, 1, index] = (
                             dp_theta[index] * mcosmphi
                         )
-                        y_lm[:, 1, l, m_abs] = p_theta[index] * sinmphi
-
-                y_lm_save[t_i, :, :, l, : l + 1] = y_lm[:, :, l, : l + 1]
-
-                if theta == 0:
-                    Y_lm_d2_theta_a[t_i, :, :, l, : l + 1] = 0.0
-                    # Not defined.
-                else:
-                    # Make use of the Laplacian identity to
-                    # estimate last derivative.
-                    Y_lm_d2_theta_a[t_i, :, :, l, : l + 1] = (
-                        lapla * y_lm[:, :, l, : l + 1]
-                        - Y_lm_d1_theta_a[t_i, :, :, l, : l + 1] * costsint
-                        - sintt * Y_lm_d2_phi_a[t_i, :, :, l, : l + 1]
-                    )
+                        y_lm[:, 1, index] = p_theta[index] * sinmphi
 
+                    y_lm_save[t_i, :, :, index] = y_lm[:, :, index]
+
+                    if theta == 0:
+                        Y_lm_d2_theta_a[t_i, :, :, index] = 0.0
+                        # Not defined.
+                    else:
+                        # Make use of the Laplacian identity to
+                        # estimate last derivative.
+                        Y_lm_d2_theta_a[t_i, :, :, index] = (
+                            lapla * y_lm[:, :, index]
+                            - Y_lm_d1_theta_a[t_i, :, :, index] * costsint
+                            - sintt * Y_lm_d2_phi_a[t_i, :, :, index]
+                        )
         if save:
             if quiet is False:
                 print("Saving SH derivatives at: %s" % (path))
@@ -371,7 +373,6 @@ def Displacement_strains(
     lon_max=360,
     lmaxgrid=None,
     grid=None,
-    only_deflec=False,
     Y_lm_d1_t=None,
     Y_lm_d1_p=None,
     Y_lm_d2_t=None,
@@ -450,8 +451,6 @@ def Displacement_strains(
         Either 'DH' or 'GLQ' for Driscoll and Healy grids or Gauss-Legendre
         Quadrature grids following the convention of SHTOOLs.
         If None, the grid is set to 'GLQ'.
-    only_deflec : bool, optional, default = False
-        Output only the displacement grid for all latitude and longitudes.
     Y_lm_d1_t : array, float, size(2,lmax+1,lmax+1), optional, default = None
         Array with the first derivative
         of Legendre polynomials with respect to colatitude.
@@ -568,17 +567,22 @@ def Displacement_strains(
     )
     cotcsc = csc * cot
 
+    # Convert 3-D of SH to 2-D indexed array
+    w_lm = pysh.shio.SHCilmToCindex(w_lm, lmax)
+    A_lm = pysh.shio.SHCilmToCindex(A_lm, lmax)
+
     y_lm = y_lm[mask]
     Y_lm_d2_t = Y_lm_d2_t[mask]
     Y_lm_d2_p = Y_lm_d2_p[mask]
     Y_lm_d1_t = Y_lm_d1_t[mask]
     Y_lm_d1_p = Y_lm_d1_p[mask]
     Y_lm_d2_tp = Y_lm_d2_tp[mask]
 
-    ein_sum = "mijk,ijk->m"
-    ein_sum_mul = "mijk,ijk,m->m"
+    ein_sum = "mij,ij->m"
+    ein_sum_mul = "mik,ik,m->m"
     path_sum = ["einsum_path", (0, 1)]  # Generated from np.einsum_path
     path_mul = ["einsum_path", (0, 1), (0, 1)]  # Generated from np.einsum_path
+
     w_deflec_ylm = R_m1 * np.einsum(ein_sum, y_lm, w_lm, optimize=path_sum)
     eps_theta[mask] = (
         R_m1 * np.einsum(ein_sum, Y_lm_d2_t, A_lm, optimize=path_sum) + w_deflec_ylm
@@ -768,6 +772,7 @@ def Plt_tecto_Mars(
                 ind_fault_check = range(faults[indx - 1] + 1, faults[indx])
                 fault_dat_lon = dat[ind_fault_check][::2]
                 fault_dat_lat = dat[ind_fault_check][1::2]
+
             split = (
                 np.argwhere((fault_dat_lon[:-1] * fault_dat_lon[1:] < 0)).ravel() + 1
             )
@@ -776,9 +781,9 @@ def Plt_tecto_Mars(
                 fault_lon_split = np.split(fault_dat_lon, split)
                 fault_dat_lat = np.split(fault_dat_lat, split)
                 for fault_lon, fault_lat in zip(fault_lon_split, fault_dat_lat):
-                    ax.plot((fault_lon + 360) % 360, fault_lat, color=col, lw=lw)
+                    ax.plot(fault_lon % 360, fault_lat, color=col, lw=lw)
             else:
-                ax.plot((fault_dat_lon + 360) % 360, fault_dat_lat, color=col, lw=lw)
+                ax.plot(fault_dat_lon % 360, fault_dat_lat, color=col, lw=lw)
 
     if legend_show:
         ax.legend(loc=legend_loc)

README.md

@@ -9,7 +9,7 @@ Planetary gravity, crustal thickness, displacement, stress, and strain calculati
 
 **Displacement_strain_planet** provides several functions and example scripts for generating gravity, crustal thickness, displacement, lateral density variations, stress, and strain maps on a planet given a set of input constraints such as from observed gravity and topography data.
 
-These functions solve the [Banerdt (1986)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB091iB01p00403) system of equations under different assumptions. Various improvements have been made to the model including the possibility to account for finite-amplitude correction and filtering [(Wieczorek & Phillips, 1998)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/97JE03136), lateral density variations at any arbitrary depth and within the surface or moho-relief [(Wieczorek et al., 2013)](https://science.sciencemag.org/content/early/2012/12/04/science.1231530?versioned=true), and density difference between the surface topography and crust [(Broquet & Wieczorek, 2019)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JE005959). 
+These functions solve the [Banerdt (1986)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB091iB01p00403) system of equations under different assumptions. Various improvements have been made to the model including the possibility to account for finite-amplitude correction and filtering [(Wieczorek & Phillips, 1998)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/97JE03136), lateral density variations at any arbitrary depth and within the surface, flexure, or moho reliefs [(Wieczorek et al., 2013)](https://science.sciencemag.org/content/early/2012/12/04/science.1231530?versioned=true), and density difference between the surface topography and crust [(Broquet & Wieczorek, 2019)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JE005959). 
 
 ### Comments
 #### Dependencies
@@ -77,4 +77,4 @@ Alternatively, you can install Displacement-strain-planet via pip
 
 ## Cite
 You can cite the latest release of the package as:
-Adrien Broquet. (2021). AB-Ares/Displacement_strain_planet: 0.3.0 (Version 0.3.0). Zenodo. http://doi.org/10.5281/zenodo.4916799
+Adrien Broquet. (2021). AB-Ares/Displacement_strain_planet: 0.3.1 (Version 0.3.1). Zenodo. http://doi.org/10.5281/zenodo.4916799

docs/_build/html/.buildinfo

@@ -1,4 +1,4 @@
 # Sphinx build info version 1
 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: 57838401250f8e1625587620d8d35764
+config: 4152686da091583470c83bc91a1d8896
 tags: 645f666f9bcd5a90fca523b33c5a78b7