Skip to content

Refactor tabulate #24

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
andreww wants to merge 2 commits into from
Closed

Refactor tabulate #24

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
02eef62
refactor: pull out and remove duplication for model tabulation
auggiemarignier Nov 19, 2025
f41802e
Use new tabulator in notebook and init
andreww Dec 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
174 changes: 30 additions & 144 deletions notebooks/PREM_travel_times_example.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions src/prem4derg/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
"""PREM-like Earth models"""

from .earth_model import OneDModel
from .earth_model import OneDModel, tabulate_model
from .const import R_EARTH
from .PREM import PREM

__all__ = ["OneDModel", "R_EARTH", "PREM"]
__all__ = ["OneDModel", "R_EARTH", "PREM", "tabulate_model"]
271 changes: 88 additions & 183 deletions src/prem4derg/earth_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,152 +163,97 @@ def _require_polynomial(self, attr: str) -> PP:
raise ValueError(f"{attr} polynomial is not defined.")
return poly

def tabulate_model_inwards(self, min_step):
"""
Return a record array representing the model handling discontiuities

This method creates a numpy record array with the model evaulated
at all depths with a minimum spacing of min_step km. All breakpoints
are also included in the output. If the densioty is discontinuoius,
the depth is represented twice, first with the value above the
discontiuity, then with the value below it. This representation can
be used to construct travel time curves (for examople).

The record array contains fields:

depth (in km)
radius (in km)
density (in kg/m^3)
qkappa (dimensionless quality factor)
qshear (dimensionless quality factor)

and is ordered such that element 0 is at the surface and the last
element (element -1) is at the center of the planet.
"""
# Keep the data as we get it
radii = np.array([])
depths = np.array([])
densities = np.array([])
vps = np.array([])
vss = np.array([])
qks = np.array([])
qms = np.array([])

nbps = len(self.density_poly.breakpoints) - 1
for i in range(nbps):
j = nbps - i
k = j - 1
rs = np.arange(
self.density_poly.breakpoints[j],
self.density_poly.breakpoints[k],
-min_step,
)
ds = self.r_earth - rs
dens = self.density(rs, break_down=True) # As we go inwards
vp = self.vp(rs, break_down=True) # As we go inwards
vs = self.vs(rs, break_down=True) # As we go inwards
qk = self.qkappa(rs, break_down=True) # As we go inwards
qm = self.qshear(rs, break_down=True) # As we go inwards
def tabulate_model(
model: OneDModel, min_step: float, outwards: bool = True
) -> np.recarray:
"""
Return a record array representing the model handling discontiuities

This method creates a numpy record array with the model evaulated
at all depths with a minimum spacing of min_step km. All breakpoints
are also included in the output. If the densioty is discontinuoius,
the depth is represented twice, first with the value above the
discontiuity, then with the value below it. This representation can
be used to construct travel time curves (for examople).

The record array contains fields:

depth (in km)
radius (in km)
density (in kg/m^3)
qkappa (dimensionless quality factor)
qshear (dimensionless quality factor)

If outwards=True, element 0 is at the centre of the planet and element -1
is at the surface. If outwards=False, element 0 is at the surface and
element -1 is at the centre of the planet.
"""
# Keep the data as we get it
radii = np.array([])
depths = np.array([])
densities = np.array([])
vps = np.array([])
vss = np.array([])
qks = np.array([])
qms = np.array([])

nbps = len(model.density_poly.breakpoints) - 1
for i in range(nbps):
j = i if outwards else nbps - i
k = j + 1 if outwards else j - 1
rs = np.arange(
model.density_poly.breakpoints[j],
model.density_poly.breakpoints[k],
min_step if outwards else -min_step,
)
ds = model.r_earth - rs
dens = model.density(rs, break_down=not outwards)
vp = model.vp(rs, break_down=not outwards)
vs = model.vs(rs, break_down=not outwards)
qk = model.qkappa(rs, break_down=not outwards)
qm = model.qshear(rs, break_down=not outwards)
radii = np.append(radii, rs)
depths = np.append(depths, ds)
densities = np.append(densities, dens)
vps = np.append(vps, vp)
vss = np.append(vss, vs)
qks = np.append(qks, qk)
qms = np.append(qms, qm)

# Look at the breakpoint. If it is discontinous in
# value put add it here (i.e. so we have above followed
# by below for the next step). Othersie we can skip it
# (and it gets adder in the next iteration). But we need
# to hadle k = 0 carefully (always stick in the origin)
cond = (k == nbps + 1) if outwards else (k == 0)
if cond:
# Add the value boundary
rs = model.density_poly.breakpoints[k]
ds = model.r_earth - rs
dens = model.density(rs)
vp = model.vp(rs)
vs = model.vs(rs)
qk = model.qkappa(rs)
qm = model.qshear(rs)
radii = np.append(radii, rs)
depths = np.append(depths, ds)
densities = np.append(densities, dens)
vps = np.append(vps, vp)
vss = np.append(vss, vs)
qks = np.append(qks, qk)
qms = np.append(qms, qm)

# Look at the breakpoint. If it is discontinous in
# value put add it here (i.e. so we have above followed
# by below for the next step). Othersie we can skip it
# (and it gets adder in the next iteration). But we need
# to hadle k = 0 carefully (always stick in the origin)
if k == 0:
# Add the value at r=0
rs = self.density_poly.breakpoints[k]
ds = self.r_earth - rs
dens = self.density(rs)
vp = self.vp(rs)
vs = self.vs(rs)
qk = self.qkappa(rs)
qm = self.qshear(rs)
radii = np.append(radii, rs)
depths = np.append(depths, ds)
densities = np.append(densities, dens)
vps = np.append(vps, vp)
vss = np.append(vss, vs)
qks = np.append(qks, qk)
qms = np.append(qms, qm)
elif self.density(self.density_poly.breakpoints[k]) != self.density(
self.density_poly.breakpoints[k], break_down=True
):
# Add the value above the inner boundary of this layer
rs = self.density_poly.breakpoints[k]
ds = self.r_earth - rs
dens = self.density(rs)
vp = self.vp(rs)
vs = self.vs(rs)
qk = self.qkappa(rs)
qm = self.qshear(rs)
radii = np.append(radii, rs)
depths = np.append(depths, ds)
densities = np.append(densities, dens)
vps = np.append(vps, vp)
vss = np.append(vss, vs)
qks = np.append(qks, qk)
qms = np.append(qms, qm)

result = np.rec.fromarrays(
[depths, radii, densities, vps, vss, qks, qms],
names="depth, radius, density, vp, vs, qkappa, qshear",
)
return result

def tabulate_model_outwards(self, min_step):
"""
Return a record array representing the model handling discontiuities

This method creates a numpy record array with the model evaulated
at all depths with a minimum spacing of min_step km. All breakpoints
are also included in the output. If the densioty is discontinuoius,
the depth is represented twice, first with the value above the
discontiuity, then with the value below it. This representation can
be used to construct travel time curves (for examople).

The record array contains fields:

depth (in km)
radius (in km)
density (in kg/m^3)
qkappa (dimensionless quality factor)
qshear (dimensionless quality factor)

and is ordered such that element 0 is at the center of the planet
and the last element (element -1) is at the surface.
"""
# Keep the data as we get it
radii = np.array([])
depths = np.array([])
densities = np.array([])
vps = np.array([])
vss = np.array([])
qks = np.array([])
qms = np.array([])

nbps = len(self.density_poly.breakpoints) - 1
for i in range(nbps):
j = i
k = j + 1
rs = np.arange(
self.density_poly.breakpoints[j],
self.density_poly.breakpoints[k],
min_step,
)
ds = self.r_earth - rs
dens = self.density(rs)
vp = self.vp(rs)
vs = self.vs(rs)
qk = self.qkappa(rs)
qm = self.qshear(rs)
elif model.density(model.density_poly.breakpoints[k]) != model.density(
model.density_poly.breakpoints[k], break_down=True
):
# Add the value above the inner boundary of this layer
rs = model.density_poly.breakpoints[k]
ds = model.r_earth - rs
dens = model.density(rs, break_down=outwards)
vp = model.vp(rs, break_down=outwards)
vs = model.vs(rs, break_down=outwards)
qk = model.qkappa(rs, break_down=outwards)
qm = model.qshear(rs, break_down=outwards)
radii = np.append(radii, rs)
depths = np.append(depths, ds)
densities = np.append(densities, dens)
Expand All @@ -317,48 +262,8 @@ def tabulate_model_outwards(self, min_step):
qks = np.append(qks, qk)
qms = np.append(qms, qm)

# Look at the breakpoint. If it is discontinous in
# value put add it here (i.e. so we have above followed
# by below for the next step). Othersie we can skip it
# (and it gets adder in the next iteration). But we need
# to hadle k = 0 carefully (always stick in the origin)
if k == nbps + 1:
# Add the value surface
rs = self.density_poly.breakpoints[k]
ds = self.r_earth - rs
dens = self.density(rs)
vp = self.vp(rs)
vs = self.vs(rs)
qk = self.qkappa(rs)
qm = self.qshear(rs)
radii = np.append(radii, rs)
depths = np.append(depths, ds)
densities = np.append(densities, dens)
vps = np.append(vps, vp)
vss = np.append(vss, vs)
qks = np.append(qks, qk)
qms = np.append(qms, qm)
elif self.density(self.density_poly.breakpoints[k]) != self.density(
self.density_poly.breakpoints[k], break_down=True
):
# Add the value above the inner boundary of this layer
rs = self.density_poly.breakpoints[k]
ds = self.r_earth - rs
dens = self.density(rs, break_down=True)
vp = self.vp(rs, break_down=True)
vs = self.vs(rs, break_down=True)
qk = self.qkappa(rs, break_down=True)
qm = self.qshear(rs, break_down=True)
radii = np.append(radii, rs)
depths = np.append(depths, ds)
densities = np.append(densities, dens)
vps = np.append(vps, vp)
vss = np.append(vss, vs)
qks = np.append(qks, qk)
qms = np.append(qms, qm)

result = np.rec.fromarrays(
[depths, radii, densities, vps, vss, qks, qms],
names="depth, radius, density, vp, vs, qkappa, qshear",
)
return result
result = np.rec.fromarrays(
[depths, radii, densities, vps, vss, qks, qms],
names="depth, radius, density, vp, vs, qkappa, qshear",
)
return result