Travel Time Calculations

# Import from obspy.taup
class TauPyModel:
    def __init__(self, model: str = "iasp91", verbose: bool = False, **kwargs):
        """
        Initialize TauP model for travel time calculations.
        
        Args:
            model: Earth model name 
            verbose: Enable verbose output
            **kwargs: Additional model options
            
        Available Models:
            '1066a', '1066b', 'ak135', 'ak135f', 'ak135f_no_mud', 'herrin',
            'iasp91', 'jb', 'prem', 'pwdk', 'sp6'
        """
    
    def get_travel_times(self, source_depth_in_km: float, distance_in_degree: float = None,
                        distance_in_km: float = None, phase_list: list[str] = ['p', 's'],
                        receiver_depth_in_km: float = 0.0) -> list:
        """
        Calculate seismic phase travel times.
        
        Args:
            source_depth_in_km: Earthquake depth in kilometers
            distance_in_degree: Epicentral distance in degrees (or use distance_in_km)
            distance_in_km: Epicentral distance in kilometers (or use distance_in_degree)
            phase_list: List of seismic phases to calculate
            receiver_depth_in_km: Receiver depth in kilometers (default: surface)
            
        Returns:
            List of Arrival objects with travel time information
            
        Common Phases:
            P, S - Direct body waves
            PP, SS - Surface reflected waves  
            PcP, ScS - Core reflected waves
            PKP, SKS - Core transmitted waves
            Pdiff, Sdiff - Diffracted waves
            PmP, SmS - Moho reflected waves
            pP, sS - Depth phases
        """
    
    def get_ray_paths(self, source_depth_in_km: float, distance_in_degree: float = None,
                     distance_in_km: float = None, phase_list: list[str] = ['p', 's'],
                     receiver_depth_in_km: float = 0.0) -> list:
        """
        Calculate seismic ray paths through Earth model.
        
        Args:
            source_depth_in_km: Source depth in km
            distance_in_degree: Distance in degrees (or use distance_in_km)
            distance_in_km: Distance in km (or use distance_in_degree)
            phase_list: Seismic phases to calculate
            receiver_depth_in_km: Receiver depth in km
            
        Returns:
            List of RayPath objects with ray geometry information
        """
    
    def get_pierce_points(self, source_depth_in_km: float, distance_in_degree: float = None,
                         distance_in_km: float = None, phase_list: list[str] = ['p', 's'],
                         receiver_depth_in_km: float = 0.0) -> list:
        """
        Calculate ray pierce points at discontinuities.
        
        Args:
            source_depth_in_km: Source depth in km
            distance_in_degree: Distance in degrees (or use distance_in_km) 
            distance_in_km: Distance in km (or use distance_in_degree)
            phase_list: Seismic phases to calculate
            receiver_depth_in_km: Receiver depth in km
            
        Returns:
            List of PiercePoint objects with pierce point information
        """
    
    def get_travel_times_geo(self, source_latitude: float, source_longitude: float,
                           source_depth_in_km: float, receiver_latitude: float,
                           receiver_longitude: float, phase_list: list[str] = ['p', 's'],
                           receiver_depth_in_km: float = 0.0) -> list:
        """
        Calculate travel times using geographic coordinates.
        
        Args:
            source_latitude: Event latitude in degrees
            source_longitude: Event longitude in degrees
            source_depth_in_km: Event depth in km
            receiver_latitude: Station latitude in degrees
            receiver_longitude: Station longitude in degrees
            phase_list: Seismic phases to calculate
            receiver_depth_in_km: Station elevation (negative for depth) in km
            
        Returns:
            List of Arrival objects with travel times and ray parameters
        """
    
    def get_ray_paths_geo(self, source_latitude: float, source_longitude: float,
                         source_depth_in_km: float, receiver_latitude: float,
                         receiver_longitude: float, phase_list: list[str] = ['p', 's'],
                         receiver_depth_in_km: float = 0.0) -> list:
        """
        Calculate ray paths using geographic coordinates.
        
        Args:
            source_latitude: Event latitude
            source_longitude: Event longitude
            source_depth_in_km: Event depth in km
            receiver_latitude: Station latitude
            receiver_longitude: Station longitude  
            phase_list: Seismic phases
            receiver_depth_in_km: Station depth in km
            
        Returns:
            List of RayPath objects with geographic ray paths
        """

def plot_travel_times(source_depth: float, ax=None, fig=None, phase_list: list[str] = None,
                     npoints: int = 100, model: str = 'iasp91', legend: bool = True,
                     **kwargs):
    """
    Plot travel time curves for seismic phases.
    
    Args:
        source_depth: Source depth in km
        ax: Matplotlib axes object (created if None)
        fig: Matplotlib figure object (created if None)  
        phase_list: Phases to plot (common phases if None)
        npoints: Number of distance points
        model: Earth model to use
        legend: Show phase legend
        **kwargs: Additional plotting options
        
    Returns:
        Matplotlib figure and axes objects
    """

def plot_ray_paths(source_depth: float, distance: float, ax=None, fig=None,
                  phase_list: list[str] = None, model: str = 'iasp91',
                  legend: bool = True, **kwargs):
    """
    Plot ray paths through Earth model.
    
    Args:
        source_depth: Source depth in km
        distance: Epicentral distance in degrees
        ax: Matplotlib axes object
        fig: Matplotlib figure object
        phase_list: Phases to plot
        model: Earth model
        legend: Show legend
        **kwargs: Plotting options
        
    Returns:
        Figure and axes objects
    """

class Arrival:
    """Travel time calculation result for single phase."""
    
    @property
    def name(self) -> str:
        """Phase name (e.g., 'P', 'S', 'PKP')."""
    
    @property
    def time(self) -> float:
        """Travel time in seconds."""
    
    @property
    def distance(self) -> float:
        """Epicentral distance in degrees."""
    
    @property
    def ray_param(self) -> float:
        """Ray parameter in s/degree."""
    
    @property
    def ray_param_sec_degree(self) -> float:
        """Ray parameter in s/degree."""
    
    @property
    def takeoff_angle(self) -> float:
        """Takeoff angle from source in degrees."""
    
    @property
    def incident_angle(self) -> float:
        """Incident angle at receiver in degrees."""
    
    @property
    def purist_distance(self) -> float:
        """Distance for which phase is pure."""
    
    @property
    def purist_name(self) -> str:
        """Purist phase name."""

class RayPath:
    """Ray path through Earth model."""
    
    @property
    def name(self) -> str:
        """Phase name."""
    
    @property
    def time(self) -> float:
        """Travel time in seconds."""
    
    @property
    def distance(self) -> float:
        """Distance in degrees."""
    
    @property 
    def ray_param(self) -> float:
        """Ray parameter."""
    
    @property
    def path(self) -> dict:
        """Ray path coordinates.
        
        Returns:
            Dictionary with 'dist' and 'depth' arrays for plotting
        """

class PiercePoint:
    """Ray pierce point at discontinuity."""
    
    @property
    def name(self) -> str:
        """Phase name."""
    
    @property
    def time(self) -> float:
        """Time to pierce point."""
    
    @property
    def distance(self) -> float:
        """Distance to pierce point."""
    
    @property
    def depth(self) -> float:
        """Depth of pierce point."""

from obspy.taup import TauPyModel
from obspy import UTCDateTime

# Initialize model
model = TauPyModel(model="iasp91")

# Calculate P and S arrival times
arrivals = model.get_travel_times(source_depth_in_km=10.0,
                                 distance_in_degree=45.0,
                                 phase_list=['P', 'S'])

print("Phase Arrivals:")
for arrival in arrivals:
    print(f"{arrival.name}: {arrival.time:.2f} s")
    print(f"  Takeoff angle: {arrival.takeoff_angle:.1f}°")
    print(f"  Ray parameter: {arrival.ray_param:.3f} s/°")

# Calculate more phases including core phases
arrivals = model.get_travel_times(source_depth_in_km=150.0,
                                 distance_in_degree=85.0,
                                 phase_list=['P', 'S', 'PP', 'SS', 'PKP', 'SKS'])

print(f"\nTotal arrivals found: {len(arrivals)}")

from obspy.taup import TauPyModel
from obspy.geodetics import gps2dist_azimuth
from obspy import UTCDateTime

# Event and station information
event_lat, event_lon = 35.0, 140.0  # Japan
event_depth = 10.0  # km
origin_time = UTCDateTime("2023-01-15T10:30:00")

stations = [
    ("Tokyo", 35.7, 139.7),
    ("Osaka", 34.7, 135.5),
    ("Sendai", 38.3, 140.9)
]

model = TauPyModel(model="ak135")

print("Predicted P-wave arrivals:")
for station_name, sta_lat, sta_lon in stations:
    # Calculate distance and azimuth
    dist_m, az, baz = gps2dist_azimuth(event_lat, event_lon, sta_lat, sta_lon)
    dist_deg = dist_m / 111319.9  # Convert to degrees
    
    # Get P-wave travel time
    arrivals = model.get_travel_times(source_depth_in_km=event_depth,
                                     distance_in_degree=dist_deg,
                                     phase_list=['P'])
    
    if arrivals:
        p_arrival = arrivals[0]
        predicted_time = origin_time + p_arrival.time
        print(f"{station_name}: {predicted_time} ({p_arrival.time:.2f} s)")

from obspy.taup import TauPyModel, plot_ray_paths
import matplotlib.pyplot as plt

model = TauPyModel(model="prem")

# Get ray paths for different phases
source_depth = 600.0  # Deep earthquake
distance = 75.0       # Teleseismic distance

ray_paths = model.get_ray_paths(source_depth_in_km=source_depth,
                               distance_in_degree=distance,
                               phase_list=['P', 'PP', 'PKP', 'PKIKP'])

print("Ray path information:")
for path in ray_paths:
    print(f"{path.name}: {path.time:.2f} s, {path.ray_param:.3f} s/°")
    
    # Get path coordinates
    coords = path.path
    max_depth = max(coords['depth'])
    print(f"  Maximum depth: {max_depth:.1f} km")

# Plot ray paths
plot_ray_paths(source_depth, distance, 
               phase_list=['P', 'PP', 'PKP', 'PKIKP'],
               model='prem')
plt.title(f"Ray Paths: {source_depth} km depth, {distance}° distance")
plt.show()

from obspy.taup import plot_travel_times
import matplotlib.pyplot as plt

# Plot travel time curves for shallow earthquake
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))

# Shallow earthquake (crustal phases important)
plot_travel_times(source_depth=10.0, ax=ax1, 
                 phase_list=['P', 'S', 'Pg', 'Sg', 'Pn', 'Sn'],
                 model='iasp91')
ax1.set_title("Shallow Earthquake (10 km depth)")
ax1.set_xlim(0, 30)  # Regional distances

# Deep earthquake (core phases visible)
plot_travel_times(source_depth=600.0, ax=ax2,
                 phase_list=['P', 'S', 'PP', 'SS', 'PKP', 'SKS'],
                 model='iasp91')
ax2.set_title("Deep Earthquake (600 km depth)")
ax2.set_xlim(0, 180)  # Global distances

plt.tight_layout()
plt.show()

from obspy.taup import TauPyModel

model = TauPyModel(model="ak135")

# Get pierce points for P wave
pierce_points = model.get_pierce_points(source_depth_in_km=50.0,
                                       distance_in_degree=60.0,
                                       phase_list=['P'])

if pierce_points:
    p_pierce = pierce_points[0]
    print("P-wave pierce points:")
    
    # Pierce points are stored in path attribute
    for i, (dist, depth) in enumerate(zip(p_pierce.pierce['dist'], 
                                         p_pierce.pierce['depth'])):
        print(f"Point {i}: {dist:.2f}° distance, {depth:.1f} km depth")

# Analyze core phase pierce points
pierce_points = model.get_pierce_points(source_depth_in_km=100.0,
                                       distance_in_degree=120.0,
                                       phase_list=['PKP'])

if pierce_points:
    pkp_pierce = pierce_points[0]
    print(f"\nPKP phase ({pkp_pierce.name}):")
    print(f"Time: {pkp_pierce.time:.2f} s")
    print(f"Ray parameter: {pkp_pierce.ray_param:.3f} s/°")

from obspy.taup import TauPyModel
import matplotlib.pyplot as plt
import numpy as np

# Compare travel times between models
models = ['iasp91', 'ak135', 'prem']
source_depth = 100.0
distances = np.linspace(10, 90, 9)

fig, ax = plt.subplots(figsize=(10, 6))

for model_name in models:
    model = TauPyModel(model=model_name)
    p_times = []
    
    for dist in distances:
        arrivals = model.get_travel_times(source_depth_in_km=source_depth,
                                        distance_in_degree=dist,
                                        phase_list=['P'])
        if arrivals:
            p_times.append(arrivals[0].time)
        else:
            p_times.append(np.nan)
    
    ax.plot(distances, p_times, 'o-', label=f'{model_name} P-wave')

ax.set_xlabel('Distance (degrees)')
ax.set_ylabel('Travel Time (seconds)')
ax.set_title(f'P-wave Travel Times - {source_depth} km depth')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()

# Print model differences
print("P-wave time differences at 45° distance:")
ref_model = TauPyModel(model='iasp91')
ref_time = ref_model.get_travel_times(100.0, 45.0, ['P'])[0].time

for model_name in ['ak135', 'prem']:
    test_model = TauPyModel(model=model_name)
    test_time = test_model.get_travel_times(100.0, 45.0, ['P'])[0].time
    diff = test_time - ref_time
    print(f"{model_name} vs iasp91: {diff:.2f} s")

# Earth model information
ModelInfo = {
    'name': str,              # Model name  
    'description': str,       # Model description
    'min_radius': float,      # Minimum radius (km)
    'max_radius': float,      # Maximum radius (km)
    'discontinuities': list[float],  # Major discontinuity depths
    'phases': list[str]       # Supported seismic phases
}

# Travel time result structure
TravelTimeResult = {
    'phase': str,             # Phase name
    'time': float,            # Travel time (s)
    'distance': float,        # Distance (degrees)
    'ray_param': float,       # Ray parameter (s/degree)
    'takeoff_angle': float,   # Takeoff angle (degrees)
    'incident_angle': float,  # Incident angle (degrees)
    'purist_distance': float, # Purist distance (degrees)
    'purist_name': str        # Purist phase name
}