Astrophysics Calculations

from gwpy.astro import (inspiral_range, inspiral_range_psd, 
                       burst_range, burst_range_spectrum,
                       sensemon_range, sensemon_range_psd,
                       range_timeseries, range_spectrogram)

def inspiral_range(psd, mass1=1.4, mass2=1.4, **kwargs):
    """
    Calculate binary inspiral detection range.
    
    Parameters:
    - psd: FrequencySeries, detector power spectral density
    - mass1: float, primary mass in solar masses  
    - mass2: float, secondary mass in solar masses
    - snr: float, signal-to-noise ratio threshold (default: 8)
    - fmin: float, minimum frequency for integration (Hz)
    - fmax: float, maximum frequency for integration (Hz)
    
    Returns:
    Quantity, detection range in Mpc
    """

def inspiral_range_psd(psd, **kwargs):
    """
    Calculate inspiral range from PSD (alias for inspiral_range).
    
    Parameters:
    - psd: FrequencySeries, power spectral density
    - **kwargs: additional parameters for inspiral_range
    
    Returns:
    Quantity, detection range in Mpc
    """

def burst_range(spectrum, **kwargs):
    """
    Calculate burst detection range.
    
    Parameters:
    - spectrum: FrequencySeries, amplitude spectral density
    - energy: float, gravitational-wave energy (solar masses × c²)
    - frequency: float, characteristic frequency (Hz)
    - snr: float, SNR threshold (default: 8)
    
    Returns:
    Quantity, detection range in Mpc
    """

def burst_range_spectrum(spectrum, **kwargs):
    """
    Calculate burst range from spectrum (alias for burst_range).
    
    Parameters:
    - spectrum: FrequencySeries, amplitude spectral density  
    - **kwargs: additional parameters for burst_range
    
    Returns:
    Quantity, detection range in Mpc
    """

def sensemon_range(psd, **kwargs):
    """
    Calculate SenseMon-style detection range.
    
    Parameters:
    - psd: FrequencySeries, power spectral density
    - mass1: float, primary mass (default: 1.4 solar masses)
    - mass2: float, secondary mass (default: 1.4 solar masses)
    
    Returns:
    Quantity, SenseMon range in Mpc
    """

def sensemon_range_psd(psd, **kwargs):
    """
    Calculate SenseMon range from PSD (alias for sensemon_range).
    
    Parameters:
    - psd: FrequencySeries, power spectral density
    - **kwargs: additional parameters for sensemon_range
    
    Returns:  
    Quantity, SenseMon range in Mpc
    """

def range_timeseries(timeseries, stride, **kwargs):
    """
    Calculate detection range as a time series.
    
    Parameters:
    - timeseries: TimeSeries, input strain data
    - stride: float, time between range calculations (seconds)
    - fftlength: float, FFT length for PSD calculation
    - overlap: float, overlap between FFT segments
    - **kwargs: additional range calculation parameters
    
    Returns:
    TimeSeries, detection range over time
    """

def range_spectrogram(timeseries, stride, **kwargs):
    """
    Calculate range as a time-frequency spectrogram.
    
    Parameters:
    - timeseries: TimeSeries, input strain data
    - stride: float, time between calculations
    - **kwargs: additional parameters
    
    Returns:
    Spectrogram, range evolution over time and frequency
    """

from gwpy.timeseries import TimeSeries
from gwpy.astro import inspiral_range, burst_range

# Read strain data and calculate PSD
strain = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)
psd = strain.psd(fftlength=4, overlap=2)

# Calculate binary neutron star inspiral range
bns_range = inspiral_range(psd, mass1=1.4, mass2=1.4, snr=8)
print(f"BNS inspiral range: {bns_range:.1f}")

# Calculate black hole binary range  
bbh_range = inspiral_range(psd, mass1=30, mass2=30, snr=8)
print(f"BBH inspiral range: {bbh_range:.1f}")

# Calculate burst range for supernova-like sources
asd = psd.sqrt()
sn_range = burst_range(asd, energy=1e-8, frequency=235, snr=8)  
print(f"Supernova burst range: {sn_range:.1f}")

from gwpy.astro import range_timeseries
import matplotlib.pyplot as plt

# Calculate range evolution over time
range_ts = range_timeseries(strain, 
                           stride=60,        # 1-minute intervals
                           fftlength=8,      # 8-second FFTs
                           overlap=4,        # 4-second overlap
                           mass1=1.4,        # BNS masses
                           mass2=1.4)

# Plot range time series
plot = range_ts.plot(figsize=(12, 6))
plot.set_ylabel('Inspiral Range [Mpc]')
plot.set_title('H1 BNS Inspiral Range Evolution')
plot.show()

# Calculate statistics
print(f"Mean range: {range_ts.mean():.1f}")
print(f"Range std: {range_ts.std():.1f}")
print(f"Best range: {range_ts.max():.1f}")
print(f"Worst range: {range_ts.min():.1f}")

from gwpy.timeseries import TimeSeriesDict
from gwpy.astro import inspiral_range

# Get data from multiple detectors
data = TimeSeriesDict.fetch_open_data(['H1', 'L1'], 
                                     start=1126259446, 
                                     end=1126259478)

# Calculate PSDs and ranges
ranges = {}
psds = {}

for ifo, timeseries in data.items():
    psds[ifo] = timeseries.psd(fftlength=4, overlap=2)
    ranges[ifo] = inspiral_range(psds[ifo], mass1=1.4, mass2=1.4)

# Compare detector performance
print("BNS Inspiral Ranges:")
for ifo, r in ranges.items():
    print(f"  {ifo}: {r:.1f}")

# Plot PSD comparison
plot = psds['H1'].plot(label='H1', alpha=0.8)
plot.add_frequencyseries(psds['L1'], label='L1', alpha=0.8)
plot.set_xlabel('Frequency [Hz]')
plot.set_ylabel('PSD [strain²/Hz]')
plot.set_xlim(10, 2000)
plot.set_yscale('log')
plot.legend()
plot.title('Detector Sensitivity Comparison')
plot.show()

import numpy as np

# Calculate range for different binary masses
masses = np.logspace(0, 2, 20)  # 1 to 100 solar masses
bns_ranges = []
bbh_ranges = []

for mass in masses:
    # Equal mass binaries
    bns_range = inspiral_range(psd, mass1=mass, mass2=mass, snr=8)
    bns_ranges.append(bns_range.value)

# Convert to arrays for plotting
ranges_array = np.array(bns_ranges)

# Plot mass-range relationship
plt.figure(figsize=(10, 6))
plt.loglog(masses, ranges_array, 'b-', linewidth=2, 
          label='Equal mass binaries')
plt.xlabel('Component Mass [M☉]')
plt.ylabel('Inspiral Range [Mpc]')
plt.title('Detection Range vs Binary Mass')
plt.grid(True, alpha=0.3)
plt.legend()

# Mark interesting mass ranges
plt.axvspan(1, 3, alpha=0.2, color='green', label='Neutron Stars')
plt.axvspan(5, 50, alpha=0.2, color='red', label='Stellar Black Holes')
plt.legend()
plt.show()

# Find optimal masses
max_idx = np.argmax(ranges_array)
print(f"Best range: {ranges_array[max_idx]:.1f} Mpc at {masses[max_idx]:.1f} M☉")

from gwpy.astro import sensemon_range

# Calculate SenseMon range (standard monitoring metric)
sensemon_bns = sensemon_range(psd, mass1=1.4, mass2=1.4)
print(f"SenseMon BNS range: {sensemon_bns:.1f}")

# Compare with standard inspiral range
standard_bns = inspiral_range(psd, mass1=1.4, mass2=1.4)
print(f"Standard BNS range: {standard_bns:.1f}")
print(f"Ratio: {sensemon_bns/standard_bns:.2f}")

# Calculate range time series for monitoring
sensemon_ts = range_timeseries(strain, 
                              stride=300,  # 5-minute intervals
                              fftlength=32, # Long FFTs for stability
                              range_func=sensemon_range,
                              mass1=1.4, mass2=1.4)

# Plot for monitoring dashboard
plot = sensemon_ts.plot(figsize=(12, 4))
plot.set_ylabel('SenseMon Range [Mpc]')
plot.set_title('H1 Detector Range Monitor')
plot.axhline(sensemon_ts.mean().value, color='red', linestyle='--', 
            label=f'Mean: {sensemon_ts.mean():.1f}')
plot.legend()
plot.show()

# Calculate ranges for different burst sources
burst_sources = {
    'Core Collapse Supernova': {'energy': 1e-8, 'frequency': 235},
    'Stellar Collapse': {'energy': 1e-7, 'frequency': 150},
    'Cosmic String Cusp': {'energy': 1e-9, 'frequency': 300}
}

asd = psd.sqrt()
burst_ranges = {}

print("Burst Source Detection Ranges:")
for source, params in burst_sources.items():
    range_val = burst_range(asd, 
                          energy=params['energy'],
                          frequency=params['frequency'],
                          snr=8)
    burst_ranges[source] = range_val
    print(f"  {source}: {range_val:.1f}")

# Plot burst ranges vs frequency
frequencies = [params['frequency'] for params in burst_sources.values()]
ranges = [burst_ranges[source].value for source in burst_sources.keys()]

plt.figure(figsize=(10, 6))
plt.scatter(frequencies, ranges, s=100, alpha=0.7)
for i, source in enumerate(burst_sources.keys()):
    plt.annotate(source, (frequencies[i], ranges[i]), 
                xytext=(10, 10), textcoords='offset points')

plt.xlabel('Characteristic Frequency [Hz]')
plt.ylabel('Detection Range [Mpc]')
plt.title('Burst Source Detection Ranges')
plt.grid(True, alpha=0.3)
plt.show()

from gwpy.astro import range_spectrogram

# Calculate range evolution in time-frequency
# (Note: This is a conceptual example - actual implementation may vary)
long_strain = TimeSeries.fetch_open_data('H1', 1126259446, 1126259446+3600)

# Create range spectrogram
range_spec = range_spectrogram(long_strain,
                              stride=60,     # 1-minute time bins
                              fftlength=8,   # 8-second FFTs
                              mass1=1.4, mass2=1.4)

# Plot range evolution
plot = range_spec.plot(figsize=(12, 8))
plot.set_xlabel('Time')
plot.set_ylabel('Frequency [Hz]')
plot.set_title('BNS Range Evolution')
plot.colorbar(label='Range [Mpc]')
plot.show()

# Calculate horizon distances (maximum possible range)
def horizon_distance(psd, mass1=1.4, mass2=1.4, **kwargs):
    """Calculate horizon distance (SNR=1 threshold)."""
    return inspiral_range(psd, mass1=mass1, mass2=mass2, snr=1, **kwargs)

# Calculate horizons for different mass combinations
mass_combinations = [
    (1.4, 1.4),   # BNS
    (10, 10),     # Light BBH  
    (30, 30),     # Heavy BBH
    (1.4, 10),    # NSBH
]

print("Horizon Distances (SNR=1):")
for m1, m2 in mass_combinations:
    horizon = horizon_distance(psd, mass1=m1, mass2=m2)
    detection_range = inspiral_range(psd, mass1=m1, mass2=m2, snr=8)
    
    print(f"  {m1:.1f}-{m2:.1f} M☉:")
    print(f"    Horizon: {horizon:.1f}")
    print(f"    Range (SNR=8): {detection_range:.1f}")
    print(f"    Ratio: {detection_range/horizon:.2f}")

# Calculate network detection capabilities
def network_range(ranges):
    """Estimate network detection range from individual detector ranges."""
    # Simple approximation: geometric mean for coincident detection
    if len(ranges) == 2:
        return np.sqrt(ranges[0] * ranges[1])
    elif len(ranges) == 3:
        return (ranges[0] * ranges[1] * ranges[2])**(1/3)
    else:
        return np.mean(ranges)

# Multi-detector network analysis
detectors = ['H1', 'L1', 'V1']
individual_ranges = [150, 140, 90]  # Example ranges in Mpc

network_bns_range = network_range(individual_ranges)

print("Network Analysis:")
for ifo, r in zip(detectors, individual_ranges):
    print(f"  {ifo} range: {r} Mpc")
print(f"  Network range: {network_bns_range:.1f} Mpc")

# Calculate detection volume (proportional to range³)
volumes = [r**3 for r in individual_ranges]
network_volume = network_bns_range**3

print("\nDetection Volumes (relative):")
for ifo, v in zip(detectors, volumes):
    print(f"  {ifo}: {v/1000:.1f} (×10³)")
print(f"  Network: {network_volume/1000:.1f} (×10³)")