0
# Astrophysics Calculations
1

2
Sensitivity and detection range calculations for gravitational-wave sources including binary inspirals and burst sources. These functions provide essential tools for understanding detector performance and planning astrophysical searches.
3

4
## Capabilities
5

6
### Range Calculations
7

8
Functions for calculating detection ranges for different types of gravitational-wave sources.
9

10
```python { .api }
11
from gwpy.astro import (inspiral_range, inspiral_range_psd, 
12
                       burst_range, burst_range_spectrum,
13
                       sensemon_range, sensemon_range_psd,
14
                       range_timeseries, range_spectrogram)
15

16
def inspiral_range(psd, mass1=1.4, mass2=1.4, **kwargs):
17
    """
18
    Calculate binary inspiral detection range.
19
    
20
    Parameters:
21
    - psd: FrequencySeries, detector power spectral density
22
    - mass1: float, primary mass in solar masses  
23
    - mass2: float, secondary mass in solar masses
24
    - snr: float, signal-to-noise ratio threshold (default: 8)
25
    - fmin: float, minimum frequency for integration (Hz)
26
    - fmax: float, maximum frequency for integration (Hz)
27
    
28
    Returns:
29
    Quantity, detection range in Mpc
30
    """
31

32
def inspiral_range_psd(psd, **kwargs):
33
    """
34
    Calculate inspiral range from PSD (alias for inspiral_range).
35
    
36
    Parameters:
37
    - psd: FrequencySeries, power spectral density
38
    - **kwargs: additional parameters for inspiral_range
39
    
40
    Returns:
41
    Quantity, detection range in Mpc
42
    """
43

44
def burst_range(spectrum, **kwargs):
45
    """
46
    Calculate burst detection range.
47
    
48
    Parameters:
49
    - spectrum: FrequencySeries, amplitude spectral density
50
    - energy: float, gravitational-wave energy (solar masses × c²)
51
    - frequency: float, characteristic frequency (Hz)
52
    - snr: float, SNR threshold (default: 8)
53
    
54
    Returns:
55
    Quantity, detection range in Mpc
56
    """
57

58
def burst_range_spectrum(spectrum, **kwargs):
59
    """
60
    Calculate burst range from spectrum (alias for burst_range).
61
    
62
    Parameters:
63
    - spectrum: FrequencySeries, amplitude spectral density  
64
    - **kwargs: additional parameters for burst_range
65
    
66
    Returns:
67
    Quantity, detection range in Mpc
68
    """
69

70
def sensemon_range(psd, **kwargs):
71
    """
72
    Calculate SenseMon-style detection range.
73
    
74
    Parameters:
75
    - psd: FrequencySeries, power spectral density
76
    - mass1: float, primary mass (default: 1.4 solar masses)
77
    - mass2: float, secondary mass (default: 1.4 solar masses)
78
    
79
    Returns:
80
    Quantity, SenseMon range in Mpc
81
    """
82

83
def sensemon_range_psd(psd, **kwargs):
84
    """
85
    Calculate SenseMon range from PSD (alias for sensemon_range).
86
    
87
    Parameters:
88
    - psd: FrequencySeries, power spectral density
89
    - **kwargs: additional parameters for sensemon_range
90
    
91
    Returns:  
92
    Quantity, SenseMon range in Mpc
93
    """
94

95
def range_timeseries(timeseries, stride, **kwargs):
96
    """
97
    Calculate detection range as a time series.
98
    
99
    Parameters:
100
    - timeseries: TimeSeries, input strain data
101
    - stride: float, time between range calculations (seconds)
102
    - fftlength: float, FFT length for PSD calculation
103
    - overlap: float, overlap between FFT segments
104
    - **kwargs: additional range calculation parameters
105
    
106
    Returns:
107
    TimeSeries, detection range over time
108
    """
109

110
def range_spectrogram(timeseries, stride, **kwargs):
111
    """
112
    Calculate range as a time-frequency spectrogram.
113
    
114
    Parameters:
115
    - timeseries: TimeSeries, input strain data
116
    - stride: float, time between calculations
117
    - **kwargs: additional parameters
118
    
119
    Returns:
120
    Spectrogram, range evolution over time and frequency
121
    """
122
```
123

124
### Usage Examples
125

126
#### Basic Range Calculations
127

128
```python
129
from gwpy.timeseries import TimeSeries
130
from gwpy.astro import inspiral_range, burst_range
131

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

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

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

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

150
#### Range Time Series Analysis
151

152
```python
153
from gwpy.astro import range_timeseries
154
import matplotlib.pyplot as plt
155

156
# Calculate range evolution over time
157
range_ts = range_timeseries(strain, 
158
                           stride=60,        # 1-minute intervals
159
                           fftlength=8,      # 8-second FFTs
160
                           overlap=4,        # 4-second overlap
161
                           mass1=1.4,        # BNS masses
162
                           mass2=1.4)
163

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

170
# Calculate statistics
171
print(f"Mean range: {range_ts.mean():.1f}")
172
print(f"Range std: {range_ts.std():.1f}")
173
print(f"Best range: {range_ts.max():.1f}")
174
print(f"Worst range: {range_ts.min():.1f}")
175
```
176

177
#### Multi-Detector Range Comparison
178

179
```python
180
from gwpy.timeseries import TimeSeriesDict
181
from gwpy.astro import inspiral_range
182

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

188
# Calculate PSDs and ranges
189
ranges = {}
190
psds = {}
191

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

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

201
# Plot PSD comparison
202
plot = psds['H1'].plot(label='H1', alpha=0.8)
203
plot.add_frequencyseries(psds['L1'], label='L1', alpha=0.8)
204
plot.set_xlabel('Frequency [Hz]')
205
plot.set_ylabel('PSD [strain²/Hz]')
206
plot.set_xlim(10, 2000)
207
plot.set_yscale('log')
208
plot.legend()
209
plot.title('Detector Sensitivity Comparison')
210
plot.show()
211
```
212

213
#### Mass-Dependent Range Analysis
214

215
```python
216
import numpy as np
217

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

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

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

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

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

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

252
#### SenseMon Range Monitoring
253

254
```python
255
from gwpy.astro import sensemon_range
256

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

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

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

273
# Plot for monitoring dashboard
274
plot = sensemon_ts.plot(figsize=(12, 4))
275
plot.set_ylabel('SenseMon Range [Mpc]')
276
plot.set_title('H1 Detector Range Monitor')
277
plot.axhline(sensemon_ts.mean().value, color='red', linestyle='--', 
278
            label=f'Mean: {sensemon_ts.mean():.1f}')
279
plot.legend()
280
plot.show()
281
```
282

283
#### Burst Source Range Analysis
284

285
```python
286
# Calculate ranges for different burst sources
287
burst_sources = {
288
    'Core Collapse Supernova': {'energy': 1e-8, 'frequency': 235},
289
    'Stellar Collapse': {'energy': 1e-7, 'frequency': 150},
290
    'Cosmic String Cusp': {'energy': 1e-9, 'frequency': 300}
291
}
292

293
asd = psd.sqrt()
294
burst_ranges = {}
295

296
print("Burst Source Detection Ranges:")
297
for source, params in burst_sources.items():
298
    range_val = burst_range(asd, 
299
                          energy=params['energy'],
300
                          frequency=params['frequency'],
301
                          snr=8)
302
    burst_ranges[source] = range_val
303
    print(f"  {source}: {range_val:.1f}")
304

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

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

315
plt.xlabel('Characteristic Frequency [Hz]')
316
plt.ylabel('Detection Range [Mpc]')
317
plt.title('Burst Source Detection Ranges')
318
plt.grid(True, alpha=0.3)
319
plt.show()
320
```
321

322
#### Range Spectrogram Analysis
323

324
```python
325
from gwpy.astro import range_spectrogram
326

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

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

337
# Plot range evolution
338
plot = range_spec.plot(figsize=(12, 8))
339
plot.set_xlabel('Time')
340
plot.set_ylabel('Frequency [Hz]')
341
plot.set_title('BNS Range Evolution')
342
plot.colorbar(label='Range [Mpc]')
343
plot.show()
344
```
345

346
#### Horizon Distance Calculations
347

348
```python
349
# Calculate horizon distances (maximum possible range)
350
def horizon_distance(psd, mass1=1.4, mass2=1.4, **kwargs):
351
    """Calculate horizon distance (SNR=1 threshold)."""
352
    return inspiral_range(psd, mass1=mass1, mass2=mass2, snr=1, **kwargs)
353

354
# Calculate horizons for different mass combinations
355
mass_combinations = [
356
    (1.4, 1.4),   # BNS
357
    (10, 10),     # Light BBH  
358
    (30, 30),     # Heavy BBH
359
    (1.4, 10),    # NSBH
360
]
361

362
print("Horizon Distances (SNR=1):")
363
for m1, m2 in mass_combinations:
364
    horizon = horizon_distance(psd, mass1=m1, mass2=m2)
365
    detection_range = inspiral_range(psd, mass1=m1, mass2=m2, snr=8)
366
    
367
    print(f"  {m1:.1f}-{m2:.1f} M☉:")
368
    print(f"    Horizon: {horizon:.1f}")
369
    print(f"    Range (SNR=8): {detection_range:.1f}")
370
    print(f"    Ratio: {detection_range/horizon:.2f}")
371
```
372

373
#### Network Range Analysis
374

375
```python
376
# Calculate network detection capabilities
377
def network_range(ranges):
378
    """Estimate network detection range from individual detector ranges."""
379
    # Simple approximation: geometric mean for coincident detection
380
    if len(ranges) == 2:
381
        return np.sqrt(ranges[0] * ranges[1])
382
    elif len(ranges) == 3:
383
        return (ranges[0] * ranges[1] * ranges[2])**(1/3)
384
    else:
385
        return np.mean(ranges)
386

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

391
network_bns_range = network_range(individual_ranges)
392

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

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

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