Examples#
Note
In order to access the example files, it is necessary to download or clone the ctplanet repo from GitHub. The files are located in the examples directory.
Moon#
Moon-Crust.py
#!/usr/bin/env python3
"""
Moon-Crust
Create a crustal thickness map of the Moon from gravity and topography.
This script generates two different crustal thickness maps. The first assumes
that the density of both the crust and mantle are constant, whereas the second
includes the effect of lateral variations in crustal density. This script can
be used to reproduce the results presented in Wieczorek et al. (2013).
Wieczorek, M. A., G. A. Neumann, F. Nimmo, W. S. Kiefer, G. J. Taylor,
H. J. Melosh, R. J. Phillips, S. C. Solomon, J. C. Andrews-Hanna,
S. W. Asmar, A. S. Konopliv, F. G. Lemoine, D. E. Smith, M. M. Watkins,
J. G. Williams, M. T. Zuber (2013), The crust of the Moon as seen by GRAIL,
Science, 339, 671-675, doi:10.1126/science.1231530, 2013.
"""
import os
import pyshtools as pysh
from ctplanet import pyMoho
from ctplanet import pyMohoRho
# ==== MAIN FUNCTION ====
def main():
lmax = 900 # determines spatial resolution of grids
lmax_calc = 600 # maximum degree to use in calculations
densityfile = 'Data/density_no_mare_n3000_f3050_719.sh'
pot = pysh.datasets.Moon.GRGM900C()
try:
os.mkdir('figs')
except Exception:
pass
print('Gravity file = {:s}'.format('GRGM900C'))
print('Lmax of potential coefficients = {:d}'.format(pot.lmax))
print('Reference radius (m) = {:e}'.format(pot.r0))
print('GM = {:e}\n'.format(pot.gm))
topo = pysh.datasets.Moon.LDEM_shape_pa(lmax=lmax)
topo.r0 = topo.coeffs[0, 0, 0]
print('Topography file = {:s}'.format('LDEM_shape_pa'))
print('Lmax of topography coefficients = {:d}'.format(topo.lmax))
print('Reference radius (m) = {:e}\n'.format(topo.r0))
density = pysh.SHCoeffs.from_file(densityfile, lmax=719)
rho_c0 = density.coeffs[0, 0, 0] # average grain density
print('Average grain density of crust (kg/m3) = {:e}'.format(rho_c0))
print('Lmax of density coefficients = {:d}\n'.format(density.lmax))
a_12_14_lat = -3.3450
a_12_14_long = -20.450
# Parameters corresponding to model 1 of Wieczorek et al. (2013),
# with the exception that the average thickness has been increased by 1 km.
thickave = 35.e3
porosity = 0.12
rho_m = 3220.0
filter = 1
half = 80
print('Average thickness of the crust (km) = {:e}'.format(thickave / 1.e3))
print('Porosity (%)= {:e}'.format(porosity * 100))
print('Mantle density (kg/m3)= {:e}'.format(rho_m))
# Constant density model
print('\n=== Constant density crust ===')
rho_c = rho_c0 * (1. - porosity) # assumed constant density
print('Bulk density of the crust(kg/m3)= {:e}'.format(rho_c * 100))
moho = pyMoho(pot, topo, lmax, rho_c, rho_m, thickave,
filter_type=filter, half=half, lmax_calc=lmax_calc,
quiet=False, delta_max=50.)
thick_grid = (topo.pad(lmax) - moho.pad(lmax)).expand(grid='DH2') / 1.e3
thick_grid.plot(show=False, colorbar='bottom',
cb_label='Crustal thickness, km',
fname='figs/Thick-Moon-1.png')
moho.plot_spectrum(show=False, fname='figs/Moho-spectrum-Moon-1.png')
print('Crustal thickness at Apollo 12/14 landing sites (km) = {:e}'
.format((topo.pad(lmax) - moho.pad(lmax))
.expand(lat=a_12_14_lat, lon=a_12_14_long) / 1.e3))
# Model with variable density
print('\n=== Variable density crust ===')
moho = pyMohoRho(pot, topo, density, porosity, lmax, rho_m,
thickave, filter_type=filter, half=half,
lmax_calc=lmax_calc, quiet=False, delta_max=50.)
thick_grid = (topo-moho.pad(topo.lmax)).expand(grid='DH2') / 1.e3
thick_grid.plot(show=False, colorbar='bottom',
cb_label='Crustal thickness, km',
fname='figs/Thick-Moon-2.png')
moho.plot_spectrum(show=False, fname='figs/Moho-spectrum-Moon-2.png')
print('Crustal thickness at Apollo 12/14 landing sites (km) = {:e}'.format(
(topo-moho.pad(topo.lmax)).expand(lat=a_12_14_lat,
lon=a_12_14_long) / 1.e3))
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
A script that demonstrates how to calculate the thickness of the lunar crust using either a constant or variable density crust. The latter can be used to reproduce the results presented in Wieczorek et al. (2013).
Moon-Core.py
#!/usr/bin/env python3
"""
Moon-Core
Calculate the hydrostatic shape of the core of the Moon, as presented in
Wieczorek et al. (2019).
"""
import os
import numpy as np
import pyshtools as pysh
from ctplanet import HydrostaticShapeLith
from ctplanet import HydrostaticShape
from ctplanet import InertiaTensor_from_shape
# ==== MAIN FUNCTION ====
def main():
lmax = 20
lmax_grid = 719
omega = pysh.constants.Moon.angular_velocity.value
rem = pysh.constants.Moon.orbit_semimajor_axis.value
mass_earth = pysh.constants.Earth.egm2008.mass.value
cthick = 34.e3 # 43.e3 or 34.0e3
rho_crust = 2550.
try:
os.mkdir('figs')
except Exception:
pass
out_rc_fc = "figs/rc_fc_34_2550.dat"
out_rc_rhoc = "figs/rc_rhoc_34_2550.dat"
out_rc_beta = "figs/rc_beta_34_2550.dat"
sh_core = "figs/core_34.sh"
sh_core_fluid = "figs/core_fluid_34.sh"
core_shape_wo_d1 = "figs/core_shape_wo_d1_330_34.dat"
core_shape = "figs/core_shape_330_34.dat"
rcore_int = 1.e3
rcore_start = 250.e3
rcore_end = 450.e3
rhocore_start = 5000.
rhocore_end = 8000.
rhocore_int = 1.
potential = pysh.datasets.Moon.GRGM900C()
topo = pysh.datasets.Moon.LDEM_shape_pa(lmax=900)
r0 = topo.coeffs[0, 0, 0]
r_sigma = r0 - cthick
ismr2 = 0.3927280 # Williams et al. (2014)
ismr2 = ismr2 * (1738.e3 / r0)**2
print("Mean planetary radius (km) = {:e}".format(r0 / 1.e3))
print("Is/MR2 (solid Moon using mean radius) = {:e}".format(ismr2))
print("Lmax of Gravitational potential = {:d}".format(potential.lmax))
print("Reference radius of potential model (km) = {:e}"
.format(potential.r0/1.e3))
print("GM = {:e}".format(potential.gm))
mass = potential.gm / pysh.constants.G.value
print("Mass (kg) = {:e}".format(mass))
print("Omega = {:e}".format(omega))
print("Period (days) = {:e}".format(2. * np.pi / omega / 60. / 60. / 24.))
print("Average crustal thickness (km) = {:e}".format(cthick / 1.e3))
print("Crustal density (kg/m3) = {:e}".format(rho_crust))
radius = np.zeros(4)
radius[0] = 0.
radius[2] = r0 - cthick
radius[3] = r0
rho = np.zeros(4)
rho[2] = rho_crust
mass_crust = 4. * np.pi / 3. * rho_crust * (radius[3]**3 - radius[2]**3)
n = 3
i_lith = 2
i_core = 1
# For each core radius, find rho_mantle and rho_core that fit total mass
# and moment of inertia
f_rc_fc = open(out_rc_fc, 'w')
f_rc_rhoc = open(out_rc_rhoc, 'w')
f_rc_beta = open(out_rc_beta, 'w')
for r_core in np.arange(rcore_start, rcore_end + rcore_int, rcore_int,
dtype=float):
radius[1] = r_core
first = True
for rho_core in np.arange(rhocore_start, rhocore_end + rhocore_int,
rhocore_int, dtype=float):
mass_core = 4. * np.pi / 3. * rho_core * r_core**3
mass_mantle = mass - mass_crust - mass_core
rho_mantle = mass_mantle * 3. / np.pi / 4. / (radius[2]**3 -
r_core**3)
rho[0] = rho_core
rho[1] = rho_mantle
if rho_mantle >= rho_core:
continue
ismr2_model = moi_solid(radius, rho, n)
if first is True:
diff_old = ismr2 - ismr2_model
first = False
else:
diff_new = ismr2 - ismr2_model
if diff_new * diff_old <= 0.:
# interpolate to get the best fitting core density
rho_core_final = (rho_core - rhocore_int) - diff_old * \
(rhocore_int) / (diff_new - diff_old)
rho[0] = rho_core_final
mass_core = 4. * np.pi / 3. * rho[0] * r_core**3
mass_mantle = mass - mass_crust - mass_core
rho_mantle = mass_mantle * 3. / np.pi / 4. / \
(radius[2]**3 - r_core**3)
rho[1] = rho_mantle
hlm, clm_hydro, mass_model = \
HydrostaticShapeLith(radius, rho, i_lith,
potential, topo, rho_crust,
r_sigma, omega, lmax,
rp=rem, mp=mass_earth)
# set degree-1 terms to zero
hlm[1].coeffs[:, 1, :] = 0.
a = hlm[1].expand(lat=0., lon=0., lmax_calc=lmax)
b = hlm[1].expand(lat=0., lon=90., lmax_calc=lmax)
c = hlm[1].expand(lat=90., lon=0., lmax_calc=lmax)
f_core = ((a+b)/2. - c) / ((a + b) / 2.)
beta_core = (a**2 - b**2) / (a**2 + b**2)
print(r_core/1.e3, rho[0], rho[1], f_core, beta_core)
f_rc_fc.write('{:e}, {:e}\n'.format(r_core/1.e3, f_core))
f_rc_rhoc.write('{:e}, {:e}\n'.format(r_core/1.e3, rho[0]))
f_rc_beta.write('{:e}, {:e}\n'
.format(r_core/1.e3, beta_core))
if r_core == 330.e3:
print("Rcore (km) = {:e}".format(r_core/1.e3))
print("A (km) = {:e}".format(a/1.e3))
print("B (km) = {:e}".format(b/1.e3))
print("C (km) = {:e}".format(c/1.e3))
print("rho_core (kg/m3) = {:e}".format(rho[0]))
II, AA, BB, CC, mass_model, RR, vec = \
InertiaTensor_from_shape(hlm, rho, 1, quiet=False)
grid = hlm[i_core].expand(lmax=lmax_grid, grid='DH2')
grid.to_file(core_shape_wo_d1)
print("Size of output grid = {:d}, {:d}"
.format(grid.nlat, grid.nlon))
print("Maximum = {:e}\nMinimum = {:e}"
.format(grid.max(), grid.min()))
hlm, clm_hydro, mass_model = \
HydrostaticShapeLith(radius, rho, i_lith,
potential, topo, rho_crust,
r_sigma, omega, lmax,
rp=rem, mp=mass_earth)
hlm_fluid, clm_fluid, mass_model = \
HydrostaticShape(radius, rho, omega, potential.gm,
potential.r0,
rp=rem, mp=mass_earth)
grid = hlm[i_core].expand(lmax=lmax_grid, grid='DH2')
grid_fluid = hlm_fluid[i_core].expand(lmax=lmax_grid,
grid='DH2')
a_fluid = hlm_fluid[1].expand(lat=0., lon=0.)
b_fluid = hlm_fluid[1].expand(lat=0., lon=90.)
c_fluid = hlm_fluid[1].expand(lat=90., lon=0.)
grid_fluid_surface = hlm_fluid[3].expand(
lmax=lmax_grid, grid='DH2')
print("Surface relief: Maximum = {:}, Minimum = {:}"
.format(grid_fluid_surface.max(),
grid_fluid_surface.min()))
f_core_fluid = (((a_fluid+b_fluid)/2. - c_fluid)
/ ((a_fluid + b_fluid) / 2.))
beta_core_fluid = ((a_fluid**2 - b_fluid**2) /
(a_fluid**2 + b_fluid**2))
print('f_core for a fluid planet = ', f_core_fluid)
print('beta_core for a fluid planet = ',
beta_core_fluid)
diff = grid - grid_fluid
print('Maximuim and miniumum core relief for '
'a fluid planet (m) = ',
grid_fluid.max(), grid_fluid.min())
print('Maximuim and miniumum difference with respect '
'to a fluid planet (m) = ',
diff.max(), diff.min())
hlm[i_core].to_file(sh_core)
hlm_fluid[i_core].to_file(sh_core_fluid)
grid.to_file(core_shape)
print("Maximum = {:e}\nMinimum = {:e}"
.format(grid.max(), grid.min()))
for l in range(0, 4):
for m in range(0, l+1):
print(l, m, hlm[i_core].coeffs[0, l, m],
hlm[i_core].coeffs[1, l, m])
diff_old = diff_new
def moi_solid(radius, rho, n):
"""
Calculate the mean, normalized, moment of inertia of the solid portion
of the planet.
The radius and density are discretized into shells as in the hydrostatic
flattening routines:
radius[0] = 0
radius[1] = radius of core
radius[n] = surface
rho[i] = density from radius[i] to radius[i+1]
"""
moi_solid = 0.
mass = 4. * np.pi / 3. * rho[0] * radius[1]**3
for i in range(2, n+1):
mass += 4. * np.pi / 3. * rho[i-1] * (radius[i]**3 - radius[i-1]**3)
moi_solid += 8. * np.pi / 15. * rho[i-1] * (radius[i]**5 -
radius[i-1]**5)
return moi_solid / mass / radius[n]**2
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Calculate the hydrostatic relief of the lunar core accounting for the non-hydrostatic potential that comes from the lithosphere.
Mars#
Mars-Crust.py
#!/usr/bin/env python3
"""
Mars-Crust
Create a crustal thickness map of Mars from gravity and topography.
This script generates two different crustal thickness maps. The first assumes
that the density of both the crust and mantle are constant, whereas the second
includes the effect of different densities on either side of the dichotomy
boundary. The average crustal thickness is iterated in order to obtain
a specified minimum crustal thickness. The routine does not consider the low
densities of the polar deposits.
"""
import os
import pyshtools as pysh
from ctplanet import pyMoho
from ctplanet import pyMohoRho
from ctplanet import HydrostaticShapeLith
from ctplanet import HydrostaticShape
from ctplanet import ReadRefModel
# ==== MAIN FUNCTION ====
def main():
densityfile = 'Data/dichotomy_359.sh'
model_name = ['DWThot', 'DWThotCrust1', 'DWThotCrust1r', 'EH45Tcold',
'EH45TcoldCrust1', 'EH45TcoldCrust1r', 'EH45ThotCrust2',
'EH45ThotCrust2r', 'LFAK', 'SANAK', 'TAYAK', 'DWAK',
'ZG_DW']
spec = 'Data/Mars-reference-interior-models/Smrekar/'
interior_file = [spec + name + '.deck' for name in model_name]
lmax_calc = 90
lmax = lmax_calc * 4
potential = pysh.datasets.Mars.GMM3()
try:
os.mkdir('figs')
except Exception:
pass
print('Gravity file = {:s}'.format('GMM3'))
print('Lmax of potential coefficients = {:d}'.format(potential.lmax))
print('Reference radius (km) = {:f}'.format(potential.r0 / 1.e3))
print('GM = {:e}\n'.format(potential.gm))
topo = pysh.datasets.Mars.MOLA_shape(lmax=lmax)
topo.r0 = topo.coeffs[0, 0, 0]
print('Topography file = {:s}'.format('MOLA_shape'))
print('Lmax of topography coefficients = {:d}'.format(topo.lmax))
print('Reference radius (km) = {:f}\n'.format(topo.r0 / 1.e3))
density = pysh.SHCoeffs.from_file(densityfile).pad(lmax=lmax)
print('Lmax of density coefficients = {:d}\n'.format(density.lmax))
lat_insight = 4.502384
lon_insight = 135.623447
filter = 1
half = 50
nmax = 7
lmax_hydro = 15
t0_sigma = 5. # maximum difference between minimum crustal thickness
omega = pysh.constants.Mars.angular_velocity.value
d_lith = 150.e3
t0 = 1.e3 # minimum crustal thickness
model = 10 # identifier for the interior reference model
# --- read 1D reference interior model ---
radius, rho, i_crust, i_core, i_lith = ReadRefModel(
interior_file[model], depth=d_lith, quiet=False)
rho_mantle = rho[i_crust-1]
# rho_core = rho[i_core-1]
n = len(radius) - 1
# r0_model = radius[n]
# --- Compute gravity contribution from hydrostatic density interfaces ---
thickave = 44.e3 # initial guess of average crustal thickness
r_sigma = topo.r0 - thickave
rho_c = 2900.
if True:
# compute values for a planet that is completely fluid
hlm_fluid, clm_fluid, mass_model = \
HydrostaticShape(radius, rho, omega, potential.gm, potential.r0)
print('--- Hydrostatic potential coefficients for a fluid planet ---')
print('c20 = {:e}\nc40 = {:e}'.format(clm_fluid.coeffs[0, 2, 0],
clm_fluid.coeffs[0, 4, 0]))
print('--- Hydrostatic relief of surface for a fluid planet ---')
print('h20 = {:e}\nh40 = {:e}'.format(hlm_fluid[n].coeffs[0, 2, 0],
hlm_fluid[n].coeffs[0, 4, 0]))
hlm, clm_hydro, mass_model = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo, rho_c,
r_sigma, omega, lmax_hydro)
print('Total mass of model (kg) = {:e}'.format(mass_model))
print('% of J2 arising from beneath lithosphere = {:f}'
.format(clm_hydro.coeffs[0, 2, 0]/potential.coeffs[0, 2, 0] * 100.))
potential.coeffs[:, :lmax_hydro+1, :lmax_hydro+1] -= \
clm_hydro.coeffs[:, :lmax_hydro+1, :lmax_hydro+1]
# --- Constant density model ---
rho_c = 2900.
print('-- Constant density model --\nrho_c = {:f}'.format(rho_c))
tmin = 1.e9
thickave = 44.e3 # initial guess of average crustal thickness
while abs(tmin - t0) > t0_sigma:
# iterate to fit assumed minimum crustal thickness
moho = pyMoho(potential, topo, lmax, rho_c, rho_mantle,
thickave, filter_type=filter, half=half,
lmax_calc=lmax_calc, nmax=nmax, quiet=True)
thick_grid = (topo.pad(lmax) - moho.pad(lmax)).expand(grid='DH2')
print('Average crustal thickness (km) = {:f}'.format(thickave / 1.e3))
print('Crustal thickness at InSight landing sites (km) = {:f}'
.format((topo.pad(lmax) - moho.pad(lmax))
.expand(lat=lat_insight, lon=lon_insight) / 1.e3))
tmin = thick_grid.min()
tmax = thick_grid.max()
print('Minimum thickness (km) = {:e}'.format(tmin / 1.e3))
print('Maximum thickness (km) = {:e}'.format(tmax / 1.e3))
thickave += t0 - tmin
thick_grid.plot(colorbar='bottom', show=False,
cb_label='Crustal thickness, km',
fname='figs/Thick-Mars-1.png')
moho.plot_spectrum(show=False, fname='figs/Moho-spectrum-Mars-1.png')
# --- Model with variable density ---
rho_south = 2900.
rho_north = 2900.
porosity = 0.0
print('-- Variable density model ---\n' +
'rho_south = {:f}\n'.format(rho_south) +
'rho_north = {:f}'.format(rho_north))
density = density * (rho_north - rho_south)
density.coeffs[0, 0, 0] += rho_south
tmin = 1.e9
thickave = 44.e3 # initial guess of average crustal thickness
while abs(tmin - t0) > t0_sigma:
# iterate to fit assumed minimum crustal thickness
moho = pyMohoRho(potential, topo, density, porosity, lmax,
rho_mantle, thickave, filter_type=filter,
half=half, lmax_calc=lmax_calc, quiet=True,
nmax=nmax)
thick_grid = (topo.pad(lmax) - moho.pad(lmax)).expand(grid='DH2')
print('Average crustal thickness (km) = {:e}'.format(thickave / 1.e3))
print('Crustal thickness at InSight landing sites (km) = {:e}'
.format((topo.pad(lmax) - moho.pad(lmax))
.expand(lat=lat_insight, lon=lon_insight) / 1.e3))
tmin = thick_grid.data.min()
tmax = thick_grid.data.max()
print('Minimum thickness (km) = {:e}'.format(tmin / 1.e3))
print('Maximum thickness (km) = {:e}'.format(tmax / 1.e3))
thickave += t0 - tmin
thick_grid.plot(colorbar='bottom', show=False,
cb_label='Crustal thickness, km',
fname='figs/Thick-Mars-2.png')
moho.plot_spectrum(show=False, fname='figs/Moho-spectrum-Mars-2.png')
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
A script that demonstrates how to calculate the thickness of the Martian crust using either a constant or variable density crust. For the variable density crust, the density is assumed to change discontinuously across the dichotomy boundary.
Mars-Crust-hydrostatic-tests.py
#!/usr/bin/env python3
"""
Mars-Crust-hydrostatic-tests
Create a crustal thickness map of Mars from gravity and topography and compare
how results change if hydrostatic interfaces are not taken into account.
"""
import os
import numpy as np
import pyshtools as pysh
from ctplanet import pyMoho
from ctplanet import HydrostaticShapeLith
from ctplanet import HydrostaticShape
# ==== MAIN FUNCTION ====
def main():
model_name = ['DWThot', 'DWThotCrust1', 'DWThotCrust1r', 'EH45Tcold',
'EH45TcoldCrust1', 'EH45TcoldCrust1r', 'EH45ThotCrust2',
'EH45ThotCrust2r', 'LFAK', 'SANAK', 'TAYAK', 'DWAK',
'ZG_DW']
spec = 'Data/Mars-reference-interior-models/Smrekar/'
interior_file = [spec + name + '.deck' for name in model_name]
model = 10
d_lith = 150.e3
lmax_calc = 90
lmax = lmax_calc * 4
lmax_grid = 719
potential = pysh.datasets.Mars.GMM3()
topo = pysh.datasets.Mars.MOLA_shape(lmax=lmax_grid)
topo.r0 = topo.coeffs[0, 0, 0]
try:
os.mkdir('figs')
except Exception:
pass
print('Gravity file = {:s}'.format('GMM3'))
print('Lmax of potential coefficients = {:d}'.format(potential.lmax))
print('Reference radius (km) = {:f}'.format(potential.r0 / 1.e3))
print('GM = {:e}'.format(potential.gm))
print('Mass = {:e}'.format(potential.mass))
print('Topography file = {:s}'.format('MOLA_shape'))
print('Lmax of topography coefficients = {:d}'.format(topo.lmax))
print('Reference radius (km) = {:f}\n'.format(topo.r0 / 1.e3))
lat_insight = 4.43
lon_insight = 135.84
filter = 1
half = 50
nmax = 7
lmax_hydro = 50
t0_sigma = 5. # maximum difference between minimum crustal thickness
rho_c = 2900.
t0 = 1.e3 # minimum crustal thickness
omega = pysh.constants.Mars.angular_velocity.value
print('Omega = {:e}'.format(omega))
# --- read 1D reference interior model ---
with open(interior_file[model], 'r') as f:
lines = f.readlines()
print(lines[0].strip())
data = lines[1].split()
if float(data[2]) != 1:
raise RuntimeError('Program not capable of reading polynomial ' +
'files')
num_file = int(lines[2].split()[0])
crust_index_file = int(lines[2].split()[3])
core_index_file = int(lines[2].split()[2])
i_crust_file = crust_index_file - 1
i_core_file = core_index_file - 1
radius = np.zeros(num_file)
rho = np.zeros(num_file)
num = 0
for i in range(0, num_file-1):
data = lines[i+3].split()
rb = float(data[0])
rhob = float(data[1])
data = lines[i+4].split()
rt = float(data[0])
rhot = float(data[1])
if rb == rt:
if i == i_core_file:
i_core = num
if i == i_crust_file:
i_crust = num
else:
radius[num] = rb
rho[num] = (rhot + rhob) / 2.
num += 1
radius[num] = rt
rho[num] = 0. # the density above the surface is zero
num += 1
n = num - 1
radius = radius[:n+1]
rho = rho[:n+1]
r0_model = radius[n]
print('Surface radius of model (km) = {:f}'.format(r0_model / 1.e3))
for i in range(0, n+1):
if radius[i] <= (r0_model - d_lith) and \
radius[i+1] > (r0_model - d_lith):
if radius[i] == (r0_model - d_lith):
i_lith = i
elif (r0_model - d_lith) - radius[i] <= radius[i+1] -\
(r0_model - d_lith):
i_lith = i
else:
i_lith = i + 1
break
rho_mantle = rho[i_crust-1]
rho_core = rho[i_core-1]
print('Mantle density (kg/m3) = {:f}'.format(rho_mantle))
print('Mantle radius (km) = {:f}'.format(radius[i_crust]/1.e3))
print('Core density (kg/m3) = {:f}'.format(rho_core))
print('Core radius (km) = {:f}'.format(radius[i_core]/1.e3))
print('Assumed depth of lithosphere (km) = {:f}'.format(d_lith / 1.e3))
print('Actual depth of lithosphere in discretized model (km) = {:f}'
.format((r0_model - radius[i_lith]) / 1.e3))
# --- Compute gravity contribution from hydrostatic density interfaces ---
thickave = 44.e3 # initial guess of average crustal thickness
r_sigma = topo.r0 - thickave
if True:
# compute values for a planet that is completely fluid
hlm_fluid, clm_fluid, mass_model = \
HydrostaticShape(radius, rho, omega, potential.gm, potential.r0,
i_clm_hydro=i_lith)
print('--- Hydrostatic potential coefficients for a fluid planet ---')
print('c20 = {:e}\nc40 = {:e}'.format(clm_fluid.coeffs[0, 2, 0],
clm_fluid.coeffs[0, 4, 0]))
print('--- Hydrostatic relief of surface for a fluid planet ---')
print('h20 = {:e}\nh40 = {:e}'.format(hlm_fluid[n].coeffs[0, 2, 0],
hlm_fluid[n].coeffs[0, 4, 0]))
print('% of J2 arising from beneath lithosphere = {:f}'
.format(clm_fluid.coeffs[0, 2, 0]
/ potential.coeffs[0, 2, 0] * 100.))
hlm, clm_hydro, mass_model = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo, rho_c,
r_sigma, omega, lmax_hydro)
print('Total mass of model (kg) = {:e}'.format(mass_model))
print('% of J2 arising from beneath lithosphere = {:f}'
.format(clm_hydro.coeffs[0, 2, 0]/potential.coeffs[0, 2, 0] * 100.))
# --- Constant density model without hydrostatic interfaces ---
print('-- Constant density model --\nrho_c = {:f}'.format(rho_c))
tmin = 1.e9
thickave = 44.e3 # initial guess of average crustal thickness
while abs(tmin - t0) > t0_sigma:
# iterate to fit assumed minimum crustal thickness
moho = pyMoho(potential, topo, lmax, rho_c, rho_mantle,
thickave, filter_type=filter, half=half,
lmax_calc=lmax_calc, nmax=nmax, quiet=True)
thick_grid = (topo.pad(lmax_grid)
- moho.pad(lmax_grid)).expand(grid='DH2')
print('Average crustal thickness (km) = {:f}'.format(thickave / 1.e3))
print('Crustal thickness at InSight landing sites (km) = {:f}'
.format((topo.pad(lmax) - moho.pad(lmax))
.expand(lat=lat_insight, lon=lon_insight) / 1.e3))
tmin = thick_grid.min()
tmax = thick_grid.max()
print('Minimum thickness (km) = {:e}'.format(tmin / 1.e3))
print('Maximum thickness (km) = {:e}'.format(tmax / 1.e3))
thickave += t0 - tmin
(thick_grid/1.e3).plot(show=False, colorbar='bottom',
cb_label='Crustal thickness, km',
fname='figs/Thick-Mars-without-hydro.png')
print('Thickness at north pole (km) = ', thick_grid.data[0, 0] / 1.e3)
# --- Constant density model with hydrostatic interfaces beneath
# --- lithosphere
potential.coeffs[:, :lmax_hydro+1, :lmax_hydro+1] -= \
clm_hydro.coeffs[:, :lmax_hydro+1, :lmax_hydro+1]
print('-- Constant density model with hydrostatic interfaces beneath '
'lithosphere--\nrho_c = {:f}'.format(rho_c))
tmin = 1.e9
thickave = 44.e3 # initial guess of average crustal thickness
while abs(tmin - t0) > t0_sigma:
# iterate to fit assumed minimum crustal thickness
moho = pyMoho(potential, topo, lmax, rho_c, rho_mantle,
thickave, filter_type=filter, half=half,
lmax_calc=lmax_calc, nmax=nmax, quiet=True)
thick2_grid = (topo.pad(lmax_grid)
- moho.pad(lmax_grid)).expand(grid='DH2')
print('Average crustal thickness (km) = {:f}'.format(thickave / 1.e3))
print('Crustal thickness at InSight landing sites (km) = {:f}'
.format((topo.pad(lmax) - moho.pad(lmax))
.expand(lat=lat_insight, lon=lon_insight) / 1.e3))
tmin = thick2_grid.min()
tmax = thick2_grid.max()
print('Minimum thickness (km) = {:e}'.format(tmin / 1.e3))
print('Maximum thickness (km) = {:e}'.format(tmax / 1.e3))
thickave += t0 - tmin
(thick2_grid/1.e3).plot(show=False, colorbar='bottom',
cb_label='Crustal thickness, km',
fname='figs/Thick-Mars-with-hydro-lith.png')
(thick2_grid/1.e3 - thick_grid/1.e3).plot(
show=False, colorbar='bottom', cb_label='Crustal thickness, km',
fname='figs/Thick-Mars-diff-hydro-lith.png')
print('Thickness at north pole (km) = ', thick2_grid.data[0, 0] / 1.e3)
min = (thick2_grid/1.e3 - thick_grid/1.e3).min()
max = (thick2_grid/1.e3 - thick_grid/1.e3).max()
print('Minimum and maximum difference (km) = ', min, max)
(thick2_grid/1.e3).to_file('figs/Thick-Mars-hydro-lith.dat')
(thick2_grid/1.e3 - thick_grid/1.e3).to_file(
'figs/Thick-Mars-diff-hydro-lith.dat')
# --- Constant density model with hydrostatic interfaces of fluid planet
potential.coeffs[:, :lmax_hydro+1, :lmax_hydro+1] += \
clm_hydro.coeffs[:, :lmax_hydro+1, :lmax_hydro+1]
potential.coeffs[:, :5, :5] -= clm_fluid.coeffs[:, :5, :5]
print('-- Constant density model with hydrostatic interfaces of a fluid '
'planet--\nrho_c = {:f}'.format(rho_c))
tmin = 1.e9
thickave = 44.e3 # initial guess of average crustal thickness
while abs(tmin - t0) > t0_sigma:
# iterate to fit assumed minimum crustal thickness
moho = pyMoho(potential, topo, lmax, rho_c, rho_mantle,
thickave, filter_type=filter, half=half,
lmax_calc=lmax_calc, nmax=nmax, quiet=True)
thick3_grid = (topo.pad(lmax_grid)
- moho.pad(lmax_grid)).expand(grid='DH2')
print('Average crustal thickness (km) = {:f}'.format(thickave / 1.e3))
print('Crustal thickness at InSight landing sites (km) = {:f}'
.format((topo.pad(lmax) - moho.pad(lmax))
.expand(lat=lat_insight, lon=lon_insight) / 1.e3))
tmin = thick3_grid.min()
tmax = thick3_grid.max()
print('Minimum thickness (km) = {:e}'.format(tmin / 1.e3))
print('Maximum thickness (km) = {:e}'.format(tmax / 1.e3))
thickave += t0 - tmin
(thick3_grid/1.e3).plot(show=False, colorbar='bottom',
cb_label='Crustal thickness, km',
fname='figs/Thick-Mars-with-fluid.png')
print('Thickness at north pole (km) = ', thick3_grid.data[0, 0] / 1.e3)
(thick2_grid/1.e3 - thick3_grid/1.e3).plot(
show=False, colorbar='bottom', cb_label='Crustal thickness, km',
fname='figs/Thick-Mars-diff-lith-fluid.png')
min = (thick2_grid/1.e3 - thick3_grid/1.e3).min()
max = (thick2_grid/1.e3 - thick3_grid/1.e3).max()
print('Minimum and maximum difference (km) = ', min, max)
(thick2_grid/1.e3 - thick3_grid/1.e3).to_file(
'figs/Thick-Mars-diff-lith-fluid.dat')
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Create a crustal thickness map of Mars from gravity and topography and compare how results change if hydrostatic interfaces are not taken into account.
Mars-Crust-InSight.py
#!/usr/bin/env python3
"""
Mars-Crust-InSight
Create a crustal thickness map of Mars from gravity and topography, using
the InSight crustal thickness constraint.
The script assumes that the density of both the crust and mantle are constant,
and that porosity in present in a constant thickness surficial layer (excluding
the polar caps). The gravitational contribution of the polar caps are
explicitly accounted for, and the average crustal thickness is iterated in
order to fit the specified thickness at the InSight landing site.
"""
import os
import matplotlib.pyplot as plt
import pyshtools as pysh
from ctplanet import pyMoho
from ctplanet import HydrostaticShapeLith
from ctplanet import ReadRefModel
# ==== MAIN FUNCTION ====
def main():
northpolarcap = 'Data/Mars_NorthPolarCapThickness719.sh.gz'
southpolarcap = 'Data/Mars_SouthPolarCapThickness719.sh.gz'
model_name1 = ['DWThot', 'DWThotCrust1', 'DWThotCrust1r', 'EH45Tcold',
'EH45TcoldCrust1', 'EH45TcoldCrust1r', 'EH45ThotCrust2',
'EH45ThotCrust2r', 'LFAK', 'SANAK', 'TAYAK', 'DWAK',
'ZG_DW', 'YOTHotRc1760kmDc40km', 'YOTHotRc1810kmDc40km']
spec1 = 'Data/Mars-reference-interior-models/Smrekar/'
model_name2 = ['Khan2022']
spec2 = 'Data/Mars-reference-interior-models/'
interior_file = [spec1 + name + '.deck' for name in model_name1]
interior_file += [spec2 + name + '.deck' for name in model_name2]
model_name = model_name1 + model_name2
lmax_calc = 90
lmax = lmax_calc * 4
potential = pysh.datasets.Mars.GMM3()
print('Gravity file = {:s}'.format('GMM3'))
print('Lmax of potential coefficients = {:d}'.format(potential.lmax))
print('Reference radius (km) = {:f}'.format(potential.r0 / 1.e3))
print('GM = {:e}\n'.format(potential.gm))
# read topo up to degree lmax
topo = pysh.datasets.Mars.MOLA_shape(lmax=lmax)
topo_grid = topo.pad(lmax).expand(grid='DH2')
print('Topography file = {:s}'.format('MOLA_shape'))
print('Lmax of topography coefficients = {:d}'.format(topo.lmax))
print('Reference radius (km) = {:f}\n'.format(topo.coeffs[0, 0, 0] / 1.e3))
# read polar cap thicknesses up to lmax
npc = pysh.SHCoeffs.from_file(northpolarcap, lmax=lmax)
spc = pysh.SHCoeffs.from_file(southpolarcap, lmax=lmax)
rho_npc = 1250.
rho_spc = 1300.
print('North polar cap thickness file = {:s}'.format(northpolarcap))
print('Lmax of North polar cap coefficients = {:d}'.format(npc.lmax))
print('Assumed density of the North polar cap (kg m-3) = {:e}'
.format(rho_npc))
print('South polar cap thickness file = {:s}'.format(southpolarcap))
print('Lmax of South polar cap coefficients = {:d}'.format(spc.lmax))
print('Assumed density of the South polar cap (kg m-3) = {:e}'
.format(rho_npc))
# Topography excluding the polar caps
topo_no_pc = topo - npc - spc
topo_no_pc_grid = topo_no_pc.expand(grid='DH2')
# compute gravitational contribution of the polar caps
npc_correction = pysh.SHGravCoeffs.from_shape(shape=topo_no_pc+npc,
rho=1.0,
gm=potential.gm,
nmax=8,
lmax=lmax_calc,
lmax_grid=lmax,
lmax_calc=lmax)
spc_correction = pysh.SHGravCoeffs.from_shape(shape=topo_no_pc + spc,
rho=1.0,
gm=potential.gm,
nmax=8,
lmax=lmax_calc,
lmax_grid=lmax,
lmax_calc=lmax)
base_correction = pysh.SHGravCoeffs.from_shape(shape=topo_no_pc,
rho=1.0,
gm=potential.gm,
nmax=8,
lmax=lmax_calc,
lmax_grid=lmax,
lmax_calc=lmax)
pc_correction = rho_npc * npc_correction.change_ref(r0=potential.r0)
pc_correction += rho_spc * spc_correction.change_ref(r0=potential.r0)
pc_correction -= (rho_spc + rho_npc) * \
base_correction.change_ref(r0=potential.r0)
filter = 1
half = 50
nmax = 7
lmax_hydro = 15
t_sigma = 5. # maximum difference between crustal thickness iterations
omega = pysh.constants.Mars.angular_velocity.value
d_lith = 150.e3
# Get input model and InSight crustal thickness
lat_insight = 4.502384
lon_insight = 135.623447
print('Input crustal thickness at InSight landing site (km) >')
t_insight_ref = float(input())
t_insight_ref *= 1.e3
ident = str(int(t_insight_ref / 1.e3))
print('Input thickness of porous layer (km) >')
t_porosity = float(input())
t_porosity *= 1.e3
porosity = 0.
if t_porosity != 0.:
directory = 'InSight/constant-porosity/'
print('Input porosity of porous layer in percent >')
porosity = float(input())
porosity /= 100.
else:
directory = 'InSight/constant/'
if porosity != 0.:
ident += '-' + str(int(t_porosity / 1.e3)) + '-' \
+ str(int(porosity * 100))
try:
os.mkdir(directory)
except Exception:
pass
rho_c_min = 2550.
rho_c_max = 3300.
rho_c_int = 50.
f_summary = open(directory + 'summary-' + ident + '.txt', 'w')
f_summary.write('Model rho_c rho_mantle t_insight t_ave '
+ ' t_min t_max\n')
# Compute contribution to the gravitational potential resulting
# from porosity. The density contrast at the surface of this layer
# is: - porosity * rho_crust and at the base it is porosity * rho_crust.
# Here we assume that rho_c is 1 and then multiply this later by the
# correct value. Note that no porosity is added to the polar caps.
if porosity != 0.:
pot_porosity = pysh.SHGravCoeffs.from_shape(topo_no_pc, -porosity,
potential.gm,
nmax=nmax,
lmax=lmax_calc,
lmax_grid=lmax,
lmax_calc=lmax
)
pot_porosity = pot_porosity.change_ref(r0=potential.r0)
pot_lower = pysh.SHGravCoeffs.from_shape(topo_no_pc - t_porosity,
porosity,
potential.gm,
nmax=nmax,
lmax=lmax_calc,
lmax_grid=lmax,
lmax_calc=lmax
)
pot_porosity += pot_lower.change_ref(r0=potential.r0)
for model in range(len(model_name)):
print('Working on model : {:s}'.format(model_name[model]))
# --- read 1D reference interior model ---
radius, rho, i_crust, i_core, i_lith = ReadRefModel(
interior_file[model], depth=d_lith, quiet=False)
rho_mantle = rho[i_crust-1]
# rho_core = rho[i_core-1]
# n = len(radius) - 1
# r0_model = radius[n]
for i in range(int((rho_c_max - rho_c_min) / rho_c_int + 1)):
rho_c = rho_c_min + i * rho_c_int
print('rho_crust (kg/m3) = {:f}'.format(rho_c))
ident = model_name[model] + '-' + str(int(t_insight_ref / 1.e3)) \
+ '-' + str(int(rho_c))
if porosity != 0.:
ident += '-' + str(int(t_porosity / 1.e3)) + '-' \
+ str(int(porosity * 100))
# Compute gravity contribution from hydrostatic density interfaces
thickave = 44.e3 # initial guess of average crustal thickness
r_sigma = topo.coeffs[0, 0, 0] - thickave
hlm, clm_hydro, mass_model = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo,
rho_c, r_sigma, omega, lmax_hydro)
print('Total mass of model (kg) = {:e}'.format(mass_model))
print('% of J2 arising from beneath lithosphere = {:f}'
.format(clm_hydro.coeffs[0, 2, 0] /
potential.coeffs[0, 2, 0] * 100.))
pot_lith = potential.copy()
pot_lith.coeffs[:, :lmax_hydro+1, :lmax_hydro+1] -= \
clm_hydro.coeffs[:, :lmax_hydro+1, :lmax_hydro+1]
if porosity != 0.:
pot_lith.coeffs[:, :lmax_calc+1, :lmax_calc+1] -= \
rho_c * pot_porosity.coeffs[:, :lmax_calc+1, :lmax_calc+1]
t_insight = 1.e9
pot_temp = pot_lith.copy()
pot_porosity_correction = pysh.SHGravCoeffs.from_zeros(
lmax_calc, potential.gm, potential.r0)
# iterate to fit assumed InSight crustal thickness
while abs(t_insight_ref - t_insight) > t_sigma:
if porosity != 0.:
pot_temp = pot_lith.pad(lmax=lmax_calc) \
- pot_porosity_correction.pad(lmax=lmax_calc)
moho = pyMoho(pot_temp, topo_no_pc, lmax, rho_c, rho_mantle,
thickave, filter_type=filter, half=half,
lmax_calc=lmax_calc, nmax=nmax,
correction=pc_correction, quiet=True)
moho_grid = moho.pad(lmax).expand(grid='DH2')
thick_grid = topo_grid - moho_grid
t_insight = (topo.pad(lmax) -
moho.pad(lmax)).expand(lat=lat_insight,
lon=lon_insight)
tmin = thick_grid.min()
tmax = thick_grid.max()
thickave += t_insight_ref - t_insight
# Compute correction resulting from porosity in the mantle
if porosity != 0.:
lower_grid = topo_no_pc_grid - t_porosity
upper_grid = lower_grid.copy()
upper_grid.data[moho_grid.data > lower_grid.data] = \
moho_grid.data[moho_grid.data > lower_grid.data]
pot_porosity_correction = \
pysh.SHGravCoeffs.from_shape(
upper_grid, -porosity * (rho_mantle - rho_c),
potential.gm, nmax=nmax, lmax_grid=lmax,
lmax=lmax_calc)
pot_porosity_correction = \
pot_porosity_correction.change_ref(r0=potential.r0)
pot_lower = pysh.SHGravCoeffs.from_shape(
lower_grid, porosity * (rho_mantle - rho_c),
potential.gm, nmax=nmax, lmax_grid=lmax,
lmax=lmax_calc)
pot_porosity_correction += \
pot_lower.change_ref(r0=potential.r0)
print('Average crustal thickness (km) = {:f}'
.format(thickave / 1.e3))
print('Crustal thickness at InSight landing sites (km) = {:f}'
.format(t_insight / 1.e3))
print('Minimum thickness (km) = {:e}'.format(tmin / 1.e3))
print('Maximum thickness (km) = {:e}'.format(tmax / 1.e3))
if tmin < 0:
break
moho.pad(lmax_calc).to_file(directory + 'Moho-Mars-'
+ ident + '.sh')
fig, ax = (thick_grid/1.e3).plot(
show=False, colorbar='bottom',
cb_label='Crustal thickness (km)',
fname=directory + 'Thick-Mars-' + ident + '.png')
# fig2, ax2 = moho.plot_spectrum(
# show=False,
# fname=directory + 'Moho-spectrum-Mars-' + ident + '.png',
# legend=ident)
f_summary.write('{:s} {:f} {:f} {:f} {:f} '
'{:f} {:f}\n'
.format(model_name[model], rho_c, rho_mantle,
t_insight/1.e3, thickave / 1.e3,
tmin / 1.e3, tmax / 1.e3))
plt.close(fig)
# plt.close(fig2)
f_summary.close()
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Create a crustal thickness map of Mars from gravity and topography, using the InSight crustal thickness constraint.
Mars-Crust-InSight-dichotomy.py
#!/usr/bin/env python3
"""
Mars-Crust-InSight-dichotomy
Create a crustal thickness map of Mars from gravity and topography, using
the InSight crustal thickness constraint.
The script assumes that the mantle density is constant and that the density of
the crust differs on either side of the dichotomy boundary. The gravitational
contribution of the polar caps are explicitly accounted for, and the average
crustal thickness is iterated in order to fit the specified thickness at the
InSight landing site.
"""
import os
import matplotlib.pyplot as plt
import pyshtools as pysh
from ctplanet import pyMohoRho
from ctplanet import HydrostaticShapeLith
from ctplanet import ReadRefModel
# ==== MAIN FUNCTION ====
def main():
densityfile = 'Data/dichotomy_359.sh'
northpolarcap = 'Data/Mars_NorthPolarCapThickness719.sh.gz'
southpolarcap = 'Data/Mars_SouthPolarCapThickness719.sh.gz'
model_name1 = ['DWThot', 'DWThotCrust1', 'DWThotCrust1r', 'EH45Tcold',
'EH45TcoldCrust1', 'EH45TcoldCrust1r', 'EH45ThotCrust2',
'EH45ThotCrust2r', 'LFAK', 'SANAK', 'TAYAK', 'DWAK',
'ZG_DW', 'YOTHotRc1760kmDc40km', 'YOTHotRc1810kmDc40km']
spec1 = 'Data/Mars-reference-interior-models/Smrekar/'
model_name2 = ['Khan2022']
spec2 = 'Data/Mars-reference-interior-models/'
interior_file = [spec1 + name + '.deck' for name in model_name1]
interior_file += [spec2 + name + '.deck' for name in model_name2]
model_name = model_name1 + model_name2
lmax_calc = 90
lmax = lmax_calc * 4
potential = pysh.datasets.Mars.GMM3()
directory = 'InSight/dichotomy/'
try:
os.mkdir(directory)
except Exception:
pass
print('Gravity file = {:s}'.format('GMM3'))
print('Lmax of potential coefficients = {:d}'.format(potential.lmax))
print('Reference radius (km) = {:f}'.format(potential.r0 / 1.e3))
print('GM = {:e}\n'.format(potential.gm))
topo = pysh.datasets.Mars.MOLA_shape(lmax=lmax)
topo_grid = topo.pad(lmax).expand(grid='DH2')
print('Topography file = {:s}'.format('MOLA_shape'))
print('Lmax of topography coefficients = {:d}'.format(topo.lmax))
print('Reference radius (km) = {:f}\n'.format(topo.coeffs[0, 0, 0] / 1.e3))
print('Crustal density file = {:s}'.format(densityfile))
dichotomy = pysh.SHCoeffs.from_file(densityfile).pad(lmax=lmax)
print('Lmax of density coefficients = {:d}\n'.format(dichotomy.lmax))
# read polar cap thicknesses up to lmax
npc = pysh.SHCoeffs.from_file(northpolarcap, lmax=lmax)
spc = pysh.SHCoeffs.from_file(southpolarcap, lmax=lmax)
rho_npc = 1250.
rho_spc = 1300.
print('North polar cap thickness file = {:s}'.format(northpolarcap))
print('Lmax of North polar cap coefficients = {:d}'.format(npc.lmax))
print('Assumed density of the North polar cap (kg m-3) = {:e}'
.format(rho_npc))
print('South polar cap thickness file = {:s}'.format(southpolarcap))
print('Lmax of South polar cap coefficients = {:d}'.format(spc.lmax))
print('Assumed density of the South polar cap (kg m-3) = {:e}'
.format(rho_npc))
# Topography excluding the polar caps
topo_no_pc = topo - npc - spc
# compute gravitational contribution of the polar caps
npc_correction = pysh.SHGravCoeffs.from_shape(shape=topo_no_pc+npc,
rho=1.0,
gm=potential.gm,
nmax=8,
lmax=lmax_calc,
lmax_grid=lmax,
lmax_calc=lmax)
spc_correction = pysh.SHGravCoeffs.from_shape(shape=topo_no_pc + spc,
rho=1.0,
gm=potential.gm,
nmax=8,
lmax=lmax_calc,
lmax_grid=lmax,
lmax_calc=lmax)
base_correction = pysh.SHGravCoeffs.from_shape(shape=topo_no_pc,
rho=1.0,
gm=potential.gm,
nmax=8,
lmax=lmax_calc,
lmax_grid=lmax,
lmax_calc=lmax)
pc_correction = rho_npc * npc_correction.change_ref(r0=potential.r0)
pc_correction += rho_spc * spc_correction.change_ref(r0=potential.r0)
pc_correction -= (rho_spc + rho_npc) * \
base_correction.change_ref(r0=potential.r0)
filter = 1
half = 50
nmax = 7
lmax_hydro = 15
t_sigma = 5. # maximum difference between crustal thickness iterations
omega = pysh.constants.Mars.angular_velocity.value
d_lith = 150.e3
# Get input model and InSight crustal thickness
lat_insight = 4.502384
lon_insight = 135.623447
print('Input crustal thickness at InSight landing site (km) >')
t_insight_ref = float(input())
t_insight_ref *= 1.e3
rho_c_min = 2550.
rho_c_max = 3200.
rho_c_int = 50.
ident = str(int(t_insight_ref / 1.e3))
f_summary = open(directory + 'summary-' + ident + '.txt', 'w')
f_summary.write('Model rho_south rho_north rho_mantle '
't_insight t_ave t_min t_max\n')
for model in range(len(model_name)):
print('Working on model : {:s}'.format(model_name[model]))
# --- read 1D reference interior model ---
radius, rho, i_crust, i_core, i_lith = ReadRefModel(
interior_file[model], depth=d_lith, quiet=False)
rho_mantle = rho[i_crust-1]
# rho_core = rho[i_core-1]
# n = len(radius) - 1
# r0_model = radius[n]
for i in range(int((rho_c_max - rho_c_min) / rho_c_int + 1)):
rho_south = rho_c_min + i * rho_c_int
for j in range(int((rho_c_max - rho_c_min) / rho_c_int + 1)):
rho_north = rho_c_min + j * rho_c_int
if rho_north < rho_south:
continue
print('rho_south = {:f}\n'.format(rho_south) +
'rho_north = {:f}'.format(rho_north))
density = dichotomy.copy()
density = density * (rho_north - rho_south)
density.coeffs[0, 0, 0] += rho_south
crustal_porosity = 0.0
ident = model_name[model] + '-' \
+ str(int(t_insight_ref / 1.e3)) \
+ '-' + str(int(rho_south)) + '-' + str(int(rho_north))
# Compute gravity contribution from hydrostatic density
# interfaces
thickave = 44.e3 # initial guess of average crustal thickness
r_sigma = topo.coeffs[0, 0, 0] - thickave
hlm, clm_hydro, mass_model = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo,
(rho_south + rho_north) / 2, r_sigma,
omega, lmax_hydro)
print('Total mass of model (kg) = {:e}'.format(mass_model))
print('% of J2 arising from beneath lithosphere = {:f}'
.format(clm_hydro.coeffs[0, 2, 0] /
potential.coeffs[0, 2, 0] * 100.))
pot_lith = potential.copy()
pot_lith.coeffs[:, :lmax_hydro+1, :lmax_hydro+1] -= \
clm_hydro.coeffs[:, :lmax_hydro+1, :lmax_hydro+1]
t_insight = 1.e9
thickave = 44.e3 # initial guess of average thickness
while abs(t_insight_ref - t_insight) > t_sigma:
# iterate to fit assumed minimum crustal thickness
moho = pyMohoRho(pot_lith, topo_no_pc, density,
crustal_porosity, lmax,
rho_mantle, thickave,
filter_type=filter,
half=half, lmax_calc=lmax_calc,
quiet=True,
correction=pc_correction,
nmax=nmax)
moho_grid = moho.pad(lmax).expand(grid='DH2')
thick_grid = topo_grid - moho_grid
t_insight = (topo.pad(lmax) -
moho.pad(lmax)).expand(lat=lat_insight,
lon=lon_insight)
tmin = thick_grid.min()
tmax = thick_grid.max()
thickave += t_insight_ref - t_insight
print('Average crustal thickness (km) = {:f}'
.format(thickave / 1.e3))
print('Crustal thickness at InSight landing sites (km) = {:f}'
.format(t_insight / 1.e3))
print('Minimum thickness (km) = {:e}'.format(tmin / 1.e3))
print('Maximum thickness (km) = {:e}'.format(tmax / 1.e3))
if tmin < 0:
break
moho.pad(lmax_calc).to_file(
directory + 'Moho-Mars-' + ident + '.sh')
fig, ax = (thick_grid/1.e3).plot(
show=False,
colorbar='bottom',
cb_label='Crustal thickness (km)',
fname=directory + 'Thick-Mars-' + ident + '.png')
# fig2, ax2 = moho.plot_spectrum(
# show=False,
# fname=directory + 'Moho-spectrum-Mars-'
# + ident + '.png',
# legend=ident)
f_summary.write('{:s} {:f} {:f} {:f} {:f} '
'{:f} {:f} {:f}\n'
.format(model_name[model], rho_south,
rho_north, rho_mantle,
t_insight / 1.e3,
thickave / 1.e3, tmin / 1.e3,
tmax / 1.e3))
plt.close(fig)
# plt.close(fig2)
f_summary.close()
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Create a crustal thickness map of Mars from gravity and topography, using the InSight crustal thickness constraint and different densities across the dichotomy boundary.
Mars-fcn.py
#!/usr/bin/env python3
"""
Mars-fcn
Compute the free core nutation period of Mars.
"""
import numpy as np
import pyshtools as pysh
from ctplanet import InertiaTensor_from_C
# ==== MAIN FUNCTION ====
def main():
lmax = 2
potential = pysh.datasets.Mars.GMM3(lmax=lmax)
mass_mars = potential.mass
r0_pot = potential.r0
omega = pysh.constants.Mars.angular_velocity.value
print('Gravity file = {:s}'.format('GMM3'))
print('Reference radius (km) = {:f}'.format(potential.r0 / 1.e3))
print('GM = {:e}'.format(potential.gm))
print('G, m3 / (kg s2) = {:e}'.format(pysh.constants.G.value))
print('Mass = {:e}'.format(potential.mass))
print('Omega = {:e}'.format(omega))
topo = pysh.datasets.Mars.MOLA_shape(lmax=lmax)
r_mars = topo.coeffs[0, 0, 0]
print('Topography file = {:s}'.format('MOLA_shape'))
print('Mean planetary radius (km) = {:f}\n'.format(r_mars / 1.e3))
# Compute polar moment using precession rate and other constants
# defined in Konopliv et al. 2016 (and 2011).
J2 = - potential.to_array(normalization='unnorm', errors=False)[0, 2, 0]
print('J2 = {:e}'.format(J2))
phidot = -7608.3 # +- 2.1 mas / yr (Konopliv et al. 2016)
phidotp = -0.2 # planetary torque correction, from Konopliv et al. 2011
phidotg = 6.7 # relativistic correction, from Konopliv et al. 2011
phi0dot = phidot - phidotp - phidotg
print('Phi0-dot, mas/yr = {:f}'.format(phi0dot))
phi0dot *= np.pi / 180 / 1000 / 60 / 60 / 365.25 / 24 / 60 / 60
e_mars = 0.09341 # Konopliv et al. 2011
print('e_mars = {:f}'.format(e_mars))
obliquity = 25.1893823 # degrees, Konopliv et al. 2016
print('obliquity, degrees = {:f}'.format(obliquity))
obliquity *= np.pi / 180
n0 = 191.408 # degrees per year, Konopliv et al. 2011
print('n0, degrees/yr = {:f}'.format(n0))
n0 *= np.pi / 180 / 24 / 60 / 60 / 365.25
C = (-1 / phi0dot) * (3/2) * (n0**2 / omega) * (1 - e_mars**2)**(-3/2) \
* J2 * np.cos(obliquity)
C *= mass_mars * r0_pot**2
print('C (Konopliv et al. 2016) = {:e}'.format(C))
print('C / (M R^2) (Konopliv et al. 2016) = {:e}'.format(C / mass_mars /
(r0_pot)**2))
print('C / (M R_mpr^2) (Konopliv et al. 2016) = {:e}\n'
.format(C / mass_mars / r_mars**2))
II, AA, BB, CC, angles = InertiaTensor_from_C(C, potential, quiet=False,
normalize=True,
r_norm=r_mars)
beta = 0.00032
Acf = 0.0254088
Bcf = 0.0254088
Ccf = 0.0255172
Ac = 0.0254036
Bc = 0.0254113
Cc = 0.0255199
print('Free core nutation period (days), fluid planet = {:f}'
.format(2 * np.pi / sigma_fcn(omega, AA, Acf, Bcf, Ccf, beta)
/ 60 / 60 / 24))
print('Free core nutation period (days), planet with lithosphere = {:f}'
.format(2 * np.pi / sigma_fcn(omega, AA, Ac, Bc, Cc, beta)
/ 60 / 60 / 24))
def sigma_fcn(omega, A, Ac, Bc, Cc, beta):
'''
Compute the free core nutation frequency.
'''
alpha = (2 * Cc - (Ac + Bc)) / (Ac + Bc)
sigma = - omega * (A / (A - Ac)) * (alpha - beta)
return sigma
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Compute the free core nutation period of Mars.
Mars-shape.py
#!/usr/bin/env python3
"""
Mars-shape
Calculations and images of the shape of density interfaces of Mars as presented
in Wieczorek et al. (2019).
"""
import os
import numpy as np
import pyshtools as pysh
from ctplanet import HydrostaticShapeLith
from ctplanet import HydrostaticShape
from ctplanet import InertiaTensor_from_shape
from ctplanet import ReadRefModel
# ==== MAIN FUNCTION ====
def main():
d_lith = 150.e3
rho_crust = 2900.
d_sigma = 45.e3
lmax_hydro = 50
model_name = ['DWThot', 'DWThotCrust1', 'DWThotCrust1r', 'EH45Tcold',
'EH45TcoldCrust1', 'EH45TcoldCrust1r', 'EH45ThotCrust2',
'EH45ThotCrust2r', 'LFAK', 'SANAK', 'TAYAK', 'DWAK',
'ZG_DW']
spec = 'Data/Mars-reference-interior-models/Smrekar/'
interior_file = [spec + name + '.deck' for name in model_name]
potential = pysh.datasets.Mars.GMM3()
omega = pysh.constants.Mars.angular_velocity.value
potential.omega = omega
mass_mars = np.float64(potential.mass)
try:
os.mkdir('figs')
except Exception:
pass
print('Gravity file = {:s}'.format('GMM3'))
print('Lmax of potential coefficients = {:d}'.format(potential.lmax))
print('Reference radius (km) = {:f}'.format(potential.r0 / 1.e3))
print('GM = {:e}'.format(potential.gm))
print('Mass = {:e}'.format(potential.mass))
print('Omega = {:e}'.format(potential.omega))
lmax_calc = 90
lmax = 359
model = 10
topo = pysh.datasets.Mars.MOLA_shape(lmax=lmax)
topo.r0 = topo.coeffs[0, 0, 0]
r_mars = topo.coeffs[0, 0, 0]
print('Topography file = {:s}'.format('MOLA_shape'))
print('Lmax of topography coefficients = {:d}'.format(topo.lmax))
print('Reference radius (km) = {:f}\n'.format(topo.r0 / 1.e3))
# --- Make geoid map, expand in spherical harmonics
# --- and then remove degree-0 term, for plotting purposes only.
u0 = potential.gm/potential.r0
geoid = potential.geoid(u0, r=topo.r0, order=2, lmax_calc=lmax_calc,
lmax=lmax)
geoidsh = geoid.geoid.expand()
geoidsh.coeffs[0, 0, 0] = 0.
geoid = geoidsh.expand(grid='DH2')
# --- read 1D reference interior model ---
radius, rho, i_crust, i_core, i_lith = ReadRefModel(interior_file[model],
depth=d_lith)
# rho_mantle = rho[i_crust-1]
# rho_core = rho[i_core-1]
n = len(radius)-1
# r0_model = radius[n]
# --- Compute purely hydrostatic relief of all interfaces ---
print('\n=== Fluid planet ===')
if False:
for i in range(1, n+1):
hlm_fluid, clm_fluid, mass_model = HydrostaticShape(
radius, rho, omega, potential.gm, potential.r0, i_clm_hydro=i)
print('i = {:d}, r = {:f}, rho = {:f}, %C20 = {:f}'
.format(i, radius[i], rho[i], clm_fluid.coeffs[0, 2, 0] /
potential.coeffs[0, 2, 0] * 100))
hlm_fluid, clm_fluid, mass_model = \
HydrostaticShape(radius, rho, omega, potential.gm, potential.r0)
print('--- Hydrostatic relief of surface for a fluid planet ---')
print('h20 = {:e}\n'.format(hlm_fluid[n].coeffs[0, 2, 0]) +
'h40 = {:e}'.format(hlm_fluid[n].coeffs[0, 4, 0]))
hydro_surface = hlm_fluid[n].expand()
print('Elevation difference between pole and equator (km)'
' {:e}'.format(hydro_surface.max()/1.e3 - hydro_surface.min()/1.e3))
print('--- Hydrostatic relief of core-mantle boundary for a fluid planet '
'---')
print('h20 = {:e}\n'.format(hlm_fluid[i_core].coeffs[0, 2, 0]) +
'h40 = {:e}'.format(hlm_fluid[i_core].coeffs[0, 4, 0]))
print('Moments of hydrostatic core')
II, AA, BB, CC, mass, RR, vec = InertiaTensor_from_shape(hlm_fluid, rho,
i_core,
quiet=True)
print('I = ', II)
print('A, B, C = ', AA, BB, CC)
print('A, B, C / (mass_mars r0^2) = ', (AA, BB, CC) / mass_mars /
r_mars**2)
print('mass of core (kg) = ', mass)
print('R core (m) = ', RR)
print('Moments of hydrostatic planet')
II, AA, BB, CC, mass, RR, vec = InertiaTensor_from_shape(hlm_fluid, rho, n,
quiet=True)
print('I = ', II)
print('A, B, C = ', AA, BB, CC)
print('A, B, C / (mass_mars r0^2) = ', (AA, BB, CC) / mass_mars /
r_mars**2)
print('mass of planet (kg) = ', mass)
print('R surface (m) = ', RR)
# --- Compute relief along hydrostatic interfaces with a lithosphere ---
print('\n=== Planet with a lithosphere ===')
r_sigma = topo.r0 - d_sigma
hlm, clm_hydro, mass_model = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo, rho_crust,
r_sigma, omega, lmax_hydro)
print('--- Core shape, with lithosphere ---')
for l in range(0, 5):
for m in range(0, l+1):
print(l, m, hlm[i_core].coeffs[0, l, m],
hlm[i_core].coeffs[1, l, m], )
print('Max and min of degree-1 core shape =',
hlm[i_core].expand(lmax=1).max(),
hlm[i_core].expand(lmax=1).min())
II, AA, BB, CC, mass, RR, vec = InertiaTensor_from_shape(hlm, rho, i_core,
quiet=True)
print('I = ', II)
print('A, B, C = ', AA, BB, CC)
print('A, B, C / (mass_mars r0^2) = ', (AA, BB, CC) / mass_mars /
r_mars**2)
print('mass of core (kg) = ', mass)
print('R core (m) = ', RR)
# --- Calculate relief, with respect to hydrostatic solution ---
# --- at i_lith and i_core
print('\n=== Difference between fluid planet and planet with a '
'lithosphere ===')
diff_ilith = hlm[i_lith] - hlm_fluid[i_lith].pad(lmax=lmax_hydro)
grid_ilith = diff_ilith.expand()
print('Maximum and minimum difference at i_lith (km) =- ',
grid_ilith.max()/1.e3, grid_ilith.min()/1.e3)
diff_icore = hlm[i_core] - hlm_fluid[i_core].pad(lmax=lmax_hydro)
grid_icore = diff_icore.expand()
print('Maximum and minimum difference at i_core (km) =- ',
grid_icore.max()/1.e3, grid_icore.min()/1.e3)
# ---- Write data to files ---
print('\n=== Output gridded data for use with GMT ===')
print('Output grid sizes = ', geoid.nlat, geoid.nlon)
(geoid/1000).to_file('figs/Mars_geoid.dat')
diff = (geoid - hlm_fluid[n].pad(lmax=lmax).expand(grid='DH2')
+ radius[n])
(diff/1000).to_file('figs/Mars_geoid_diff.dat')
diff = hlm[i_lith].pad(lmax=lmax).expand(grid='DH2') - radius[i_lith]
(diff/1000).to_file('figs/hydro_ilith.dat')
diff = hlm[i_core].pad(lmax=lmax).expand(grid='DH2') - radius[i_core]
(diff/1000).to_file('figs/hydro_icore.dat')
diff = hlm[i_lith] - hlm_fluid[i_lith].pad(lmax=lmax_hydro)
(diff.expand(grid='DH2', lmax=lmax)/1000).to_file(
'figs/hydro_ilith_diff.dat')
diff = hlm[i_core] - hlm_fluid[i_core].pad(lmax=lmax_hydro)
(diff.expand(grid='DH2', lmax=lmax)/1000).to_file(
'figs/hydro_icore_diff.dat')
# Sensitivity tests
print('\n=== Sensitivity tests ===')
r_sigma = topo.r0
hlm2, clm_hydro2, mass_model2 = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo, rho_crust,
r_sigma, omega, lmax_hydro)
print('d_sigma = 0')
print('Minimum and maximum differences at i_lith (m) = ',
(hlm2[i_lith] - hlm[i_lith]).expand().min(),
(hlm2[i_lith] - hlm[i_lith]).expand().max())
print('Minimum and maximum differences at i_core (m) = ',
(hlm2[i_core] - hlm[i_core]).expand().min(),
(hlm2[i_core] - hlm[i_core]).expand().max())
r_sigma = topo.r0 - 100.e3
hlm2, clm_hydro2, mass_model2 = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo, rho_crust,
r_sigma, omega, lmax_hydro)
print('d_sigma = 100 km')
print('Minimum and maximum differences at i_lith (m) = ',
(hlm2[i_lith] - hlm[i_lith]).expand().min(),
(hlm2[i_lith] - hlm[i_lith]).expand().max())
print('Minimum and maximum differences at i_core (m) = ',
(hlm2[i_core] - hlm[i_core]).expand().min(),
(hlm2[i_core] - hlm[i_core]).expand().max())
r_sigma = topo.r0 - d_sigma
hlm2, clm_hydro2, mass_model2 = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo, 2500.,
r_sigma, omega, lmax_hydro)
print('rho = 2500 kg / m3')
print('Minimum and maximum differences at i_lith (m) = ',
(hlm2[i_lith] - hlm[i_lith]).expand().min(),
(hlm2[i_lith] - hlm[i_lith]).expand().max())
print('Minimum and maximum differences at i_core (m) = ',
(hlm2[i_core] - hlm[i_core]).expand().min(),
(hlm2[i_core] - hlm[i_core]).expand().max())
r_sigma = topo.r0 - d_sigma
hlm2, clm_hydro2, mass_model2 = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo, 3300.,
r_sigma, omega, lmax_hydro)
print('rho = 3300 kg / m3')
print('Minimum and maximum differences at i_lith (m) = ',
(hlm2[i_lith] - hlm[i_lith]).expand().min(),
(hlm2[i_lith] - hlm[i_lith]).expand().max())
print('Minimum and maximum differences at i_core (m) = ',
(hlm2[i_core] - hlm[i_core]).expand().min(),
(hlm2[i_core] - hlm[i_core]).expand().max())
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Create images related to Mars in Wieczorek et al. (2019).
Mars-j2.py
#!/usr/bin/env python3
"""
Mars-j2
Compute the contribution to the gravitational J2 of Mars from hydrostatic
interfaces beneath the lithosphere.
"""
import numpy as np
import pyshtools as pysh
from ctplanet import HydrostaticShapeLith
# ==== MAIN FUNCTION ====
def main():
d_lith = 150.e3
rho_crust = 2900.
d_sigma = 45.e3
lmax_hydro = 2
potential = pysh.datasets.Mars.GMM3()
omega = pysh.constants.Mars.angular_velocity.value
potential.omega = omega
print('Gravity file = {:s}'.format('GMM3'))
print('Lmax of potential coefficients = {:d}'.format(potential.lmax))
print('Reference radius (km) = {:f}'.format(potential.r0 / 1.e3))
print('GM = {:e}'.format(potential.gm))
print('Mass = {:e}'.format(potential.mass))
print('Omega = {:e}'.format(potential.omega))
lmax = 359
topo = pysh.datasets.Mars.MOLA_shape(lmax=lmax)
topo.r0 = topo.coeffs[0, 0, 0]
print('Topography file = {:s}'.format('MOLA_shape'))
print('Lmax of topography coefficients = {:d}'.format(topo.lmax))
print('Reference radius (km) = {:f}\n'.format(topo.r0 / 1.e3))
# --- read 1D reference interior model ---
model_name = ['DWThot', 'DWThotCrust1', 'DWThotCrust1r', 'EH45Tcold',
'EH45TcoldCrust1', 'EH45TcoldCrust1r', 'EH45ThotCrust2',
'EH45ThotCrust2r', 'LFAK', 'SANAK', 'TAYAK', 'DWAK',
'ZG_DW']
spec = 'Data/Mars-reference-interior-models/Smrekar/'
interior_file = [spec + name + '.deck' for name in model_name]
for file in interior_file:
print('=== Reading model {:s} ==='.format(file))
with open(file, 'r') as f:
lines = f.readlines()
print(lines[0].strip())
data = lines[1].split()
if float(data[2]) != 1:
raise RuntimeError('Program not capable of reading '
'polynomial files')
num_file = int(lines[2].split()[0])
crust_index_file = int(lines[2].split()[3])
core_index_file = int(lines[2].split()[2])
i_crust_file = crust_index_file - 1
i_core_file = core_index_file - 1
radius = np.zeros(num_file)
rho = np.zeros(num_file)
num = 0
for i in range(0, num_file-1):
data = lines[i+3].split()
rb = float(data[0])
rhob = float(data[1])
data = lines[i+4].split()
rt = float(data[0])
rhot = float(data[1])
if rb == rt:
if i == i_core_file:
i_core = num
if i == i_crust_file:
i_crust = num
else:
radius[num] = rb
rho[num] = (rhot + rhob) / 2.
num += 1
radius[num] = rt
rho[num] = 0. # the density above the surface is zero
num += 1
n = num - 1
radius = radius[:n+1]
rho = rho[:n+1]
r0_model = radius[n]
print('Surface radius of model (km) = {:f}'
.format(r0_model / 1.e3))
for i in range(0, n+1):
if radius[i] <= (r0_model - d_lith) and \
radius[i+1] > (r0_model - d_lith):
if radius[i] == (r0_model - d_lith):
i_lith = i
elif (r0_model - d_lith) - radius[i] <= radius[i+1] -\
(r0_model - d_lith):
i_lith = i
else:
i_lith = i + 1
break
rho_mantle = rho[i_crust-1]
rho_core = rho[i_core-1]
print('Mantle density (kg/m3) = {:f}'.format(rho_mantle))
print('Mantle radius (km) = {:f}'.format(radius[i_crust]/1.e3))
print('Core density (kg/m3) = {:f}'.format(rho_core))
print('Core radius (km) = {:f}'.format(radius[i_core]/1.e3))
print('Assumed depth of lithosphere (km) = {:f}'
.format(d_lith / 1.e3))
print('Actual depth of lithosphere in discretized '
'model (km) = {:f}'
.format((r0_model - radius[i_lith]) / 1.e3))
r_sigma = topo.r0 - d_sigma
hlm, clm_hydro, mass_model = \
HydrostaticShapeLith(radius, rho, i_lith, potential, topo,
rho_crust, r_sigma, omega, lmax_hydro)
print('Percentage of h20 derived from hydrostatic mantle = '
'{:f}'.format(clm_hydro.coeffs[0, 2, 0] /
potential.coeffs[0, 2, 0]*100))
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Compute the contribution to the gravitational J2 of Mars from hydrostatic interfaces beneath the lithosphere.
Earth#
Earth-shape.py
#!/usr/bin/env python3
"""
Earth-shape
Compute hydrostatic shape of Earth using PREM.
"""
import numpy as np
import pyshtools as pysh
from ctplanet import HydrostaticShape
# ==== MAIN FUNCTION ====
def main():
model_name = ['PREM.dat']
spec = 'Data/'
interior_file = [spec + name for name in model_name]
r_ref = pysh.constants.Earth.wgs84.a.value
gm = pysh.constants.Earth.wgs84.gm.value
mass = gm / pysh.constants.G.value
omega = pysh.constants.Earth.wgs84.omega.value
print('Reference radius (km) = {:f}'.format(r_ref / 1.e3))
print('GM = {:e}'.format(gm))
print('Mass = {:e}'.format(mass))
print('Omega = {:e}'.format(omega))
model = 0
# --- read 1D reference interior model ---
print('=== Reading model {:s} ==='.format(model_name[model]))
with open(interior_file[model], 'r') as f:
lines = f.readlines()
print(lines[0].strip())
data = lines[1].split()
if float(data[2]) != 1:
raise RuntimeError('Program not capable of reading polynomial ' +
'files')
num = int(lines[2].split()[0])
crust_index = int(lines[2].split()[3])
# mantle_index = int(lines[2].split()[3])
radius = np.zeros(num)
rho = np.zeros(num)
for i in range(0, num):
data = lines[i+3].split()
radius[i] = float(data[0])
rho[i] = float(data[1])
r0_model = radius[num-1]
print('Surface radius of model (km) = {:f}'.format(r0_model / 1.e3))
n = num - 1
rho[n] = 0. # the density above the surface is zero
rho_mantle = rho[crust_index-1]
print('Mantle density (kg/m3) = {:f}'.format(rho_mantle))
# --- Compute gravity contribution from hydrostatic density interfaces ---
hlm_fluid, clm_fluid, mass_model = HydrostaticShape(radius, rho, omega,
gm, r_ref)
hlm_surf_unnorm = hlm_fluid[n].convert(normalization='unnorm')
print('--- Hydrostatic relief of surface ---')
print('h20 = {:e}\n'.format(hlm_fluid[n].coeffs[0, 2, 0]) +
'h40 = {:e}'.format(hlm_fluid[n].coeffs[0, 4, 0]))
print('h20 (unnorm)= {:e}\n'
.format(hlm_surf_unnorm.coeffs[0, 2, 0]/r0_model) +
'h40 (unnorm)= {:e}'
.format(hlm_surf_unnorm.coeffs[0, 4, 0]/r0_model))
print(mass, mass_model)
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Compute hydrostatic relief of Earth using PREM.
Ceres#
Ceres-shape.py
#!/usr/bin/env python3
"""
Ceres-shape
Calculate the hydrostatic shape of Ceres, as in Wieczorek et al. (2019).
"""
import numpy as np
import pyshtools
from ctplanet import HydrostaticShape
# ==== MAIN FUNCTION ====
def main():
# Thomas et al. 2005
mass = 9.395e20
gm = mass * pyshtools.constants.G.value
r0 = 476.2e3
omega = 2 * np.pi / (9.076 * 60 * 60)
volume = 4 * np.pi / 3 * r0**3
lmax = 4
# Model C1
radius = np.ndarray(2)
radius[0] = 0.
radius[1] = r0
rho = np.ndarray(2)
rho[0] = mass / volume
rho[1] = 0
print(rho[0])
hlm, clm, mass_model = HydrostaticShape(radius, rho, omega, gm, r0)
a = hlm[1].expand(lat=0., lon=0., lmax_calc=lmax)
b = hlm[1].expand(lat=0., lon=90., lmax_calc=lmax)
c = hlm[1].expand(lat=90., lon=0., lmax_calc=lmax)
print(a, b, c, a-c)
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()
Calculate the hydrostatic shape of Ceres.