Coverage for bilby/gw/utils.py: 48%

1import json 

2import os 

3from functools import lru_cache 

4 

5import numpy as np 

6from scipy.interpolate import interp1d 

7from scipy.special import i0e 

8from bilby_cython.geometry import ( 

9 zenith_azimuth_to_theta_phi as _zenith_azimuth_to_theta_phi, 

10) 

11from bilby_cython.time import greenwich_mean_sidereal_time 

12 

13from ..core.utils import (logger, run_commandline, 

14 check_directory_exists_and_if_not_mkdir, 

15 SamplesSummary, theta_phi_to_ra_dec) 

16from ..core.utils.constants import solar_mass 

17 

18 

19def asd_from_freq_series(freq_data, df): 

20 """ 

21 Calculate the ASD from the frequency domain output of gaussian_noise() 

22 

23 Parameters 

24 ========== 

25 freq_data: array_like 

26 Array of complex frequency domain data 

27 df: float 

28 Spacing of freq_data, 1/(segment length) used to generate the gaussian noise 

29 

30 Returns 

31 ======= 

32 array_like: array of real-valued normalized frequency domain ASD data 

33 

34 """ 

35 return np.absolute(freq_data) * 2 * df**0.5 

36 

37 

38def psd_from_freq_series(freq_data, df): 

39 """ 

40 Calculate the PSD from the frequency domain output of gaussian_noise() 

41 Calls asd_from_freq_series() and squares the output 

42 

43 Parameters 

44 ========== 

45 freq_data: array_like 

46 Array of complex frequency domain data 

47 df: float 

48 Spacing of freq_data, 1/(segment length) used to generate the gaussian noise 

49 

50 Returns 

51 ======= 

52 array_like: Real-valued normalized frequency domain PSD data 

53 

54 """ 

55 return np.power(asd_from_freq_series(freq_data, df), 2) 

56 

57 

58def get_vertex_position_geocentric(latitude, longitude, elevation): 

59 """ 

60 Calculate the position of the IFO vertex in geocentric coordinates in meters. 

61 

62 Based on arXiv:gr-qc/0008066 Eqs. B11-B13 except for the typo in the definition of the local radius. 

63 See Section 2.1 of LIGO-T980044-10 for the correct expression 

64 

65 Parameters 

66 ========== 

67 latitude: float 

68 Latitude in radians 

69 longitude: 

70 Longitude in radians 

71 elevation: 

72 Elevation in meters 

73 

74 Returns 

75 ======= 

76 array_like: A 3D representation of the geocentric vertex position 

77 

78 """ 

79 semi_major_axis = 6378137 # for ellipsoid model of Earth, in m 

80 semi_minor_axis = 6356752.314 # in m 

81 radius = semi_major_axis**2 * (semi_major_axis**2 * np.cos(latitude)**2 + 

82 semi_minor_axis**2 * np.sin(latitude)**2)**(-0.5) 

83 x_comp = (radius + elevation) * np.cos(latitude) * np.cos(longitude) 

84 y_comp = (radius + elevation) * np.cos(latitude) * np.sin(longitude) 

85 z_comp = ((semi_minor_axis / semi_major_axis)**2 * radius + elevation) * np.sin(latitude) 

86 return np.array([x_comp, y_comp, z_comp]) 

87 

88 

89def inner_product(aa, bb, frequency, PSD): 

90 """ 

91 Calculate the inner product defined in the matched filter statistic 

92 

93 Parameters 

94 ========== 

95 aa, bb: array_like 

96 Single-sided Fourier transform, created, e.g., by the nfft function above 

97 frequency: array_like 

98 An array of frequencies associated with aa, bb, also returned by nfft 

99 PSD: bilby.gw.detector.PowerSpectralDensity 

100 

101 Returns 

102 ======= 

103 The matched filter inner product for aa and bb 

104 

105 """ 

106 psd_interp = PSD.power_spectral_density_interpolated(frequency) 

107 

108 # calculate the inner product 

109 integrand = np.conj(aa) * bb / psd_interp 

110 

111 df = frequency[1] - frequency[0] 

112 integral = np.sum(integrand) * df 

113 return 4. * np.real(integral) 

114 

115 

116def noise_weighted_inner_product(aa, bb, power_spectral_density, duration): 

117 """ 

118 Calculate the noise weighted inner product between two arrays. 

119 

120 Parameters 

121 ========== 

122 aa: array_like 

123 Array to be complex conjugated 

124 bb: array_like 

125 Array not to be complex conjugated 

126 power_spectral_density: array_like 

127 Power spectral density of the noise 

128 duration: float 

129 duration of the data 

130 

131 Returns 

132 ======= 

133 Noise-weighted inner product. 

134 """ 

135 

136 integrand = np.conj(aa) * bb / power_spectral_density 

137 return 4 / duration * np.sum(integrand) 

138 

139 

140def matched_filter_snr(signal, frequency_domain_strain, power_spectral_density, duration): 

141 """ 

142 Calculate the _complex_ matched filter snr of a signal. 

143 This is <signal|frequency_domain_strain> / optimal_snr 

144 

145 Parameters 

146 ========== 

147 signal: array_like 

148 Array containing the signal 

149 frequency_domain_strain: array_like 

150 

151 power_spectral_density: array_like 

152 

153 duration: float 

154 Time duration of the signal 

155 

156 Returns 

157 ======= 

158 float: The matched filter signal to noise ratio squared 

159 

160 """ 

161 rho_mf = noise_weighted_inner_product( 

162 aa=signal, bb=frequency_domain_strain, 

163 power_spectral_density=power_spectral_density, duration=duration) 

164 rho_mf /= optimal_snr_squared( 

165 signal=signal, power_spectral_density=power_spectral_density, 

166 duration=duration)**0.5 

167 return rho_mf 

168 

169 

170def optimal_snr_squared(signal, power_spectral_density, duration): 

171 """ 

172 Compute the square of the optimal matched filter SNR for the provided 

173 signal. 

174 

175 

176 Parameters 

177 ========== 

178 signal: array_like 

179 Array containing the signal 

180 power_spectral_density: array_like 

181 

182 duration: float 

183 Time duration of the signal 

184 

185 Returns 

186 ======= 

187 float: The matched filter signal to noise ratio squared 

188 

189 """ 

190 return noise_weighted_inner_product(signal, signal, power_spectral_density, duration) 

191 

192 

193def overlap(signal_a, signal_b, power_spectral_density=None, delta_frequency=None, 

194 lower_cut_off=None, upper_cut_off=None, norm_a=None, norm_b=None): 

195 r""" 

196 Compute the overlap between two signals 

197 

198 .. math:: 

199 

200 {\cal O} = \frac{4 \Delta f}{N_{a} N_{b}} \sum_{i} \frac{h^{*}_{a,i} h_{b,i}}{S_{i}} 

201 

202 Parameters 

203 ---------- 

204 signal_a: array-like 

205 signal_b: array-like 

206 power_spectral_density : array-like 

207 delta_frequency: float 

208 Frequency spacing of the signals 

209 lower_cut_off: float 

210 Minimum frequency for the integral 

211 upper_cut_off: float 

212 Maximum frequency for the integral 

213 norm_a: float 

214 Normalizing factor for signal_a 

215 norm_b: float 

216 Normalizing factor for signal_b 

217 

218 Returns 

219 ------- 

220 float 

221 The overlap 

222 """ 

223 low_index = int(lower_cut_off / delta_frequency) 

224 up_index = int(upper_cut_off / delta_frequency) 

225 integrand = np.conj(signal_a) * signal_b 

226 integrand = integrand[low_index:up_index] / power_spectral_density[low_index:up_index] 

227 integral = (4 * delta_frequency * integrand) / norm_a / norm_b 

228 return sum(integral).real 

229 

230 

231def zenith_azimuth_to_theta_phi(zenith, azimuth, ifos): 

232 """ 

233 Convert from the 'detector frame' to the Earth frame. 

234 

235 Parameters 

236 ========== 

237 kappa: float 

238 The zenith angle in the detector frame 

239 eta: float 

240 The azimuthal angle in the detector frame 

241 ifos: list 

242 List of Interferometer objects defining the detector frame 

243 

244 Returns 

245 ======= 

246 theta, phi: float 

247 The zenith and azimuthal angles in the earth frame. 

248 """ 

249 delta_x = ifos[0].geometry.vertex - ifos[1].geometry.vertex 

250 return _zenith_azimuth_to_theta_phi(zenith, azimuth, delta_x) 

251 

252 

253def zenith_azimuth_to_ra_dec(zenith, azimuth, geocent_time, ifos): 

254 """ 

255 Convert from the 'detector frame' to the Earth frame. 

256 

257 Parameters 

258 ========== 

259 kappa: float 

260 The zenith angle in the detector frame 

261 eta: float 

262 The azimuthal angle in the detector frame 

263 geocent_time: float 

264 GPS time at geocenter 

265 ifos: list 

266 List of Interferometer objects defining the detector frame 

267 

268 Returns 

269 ======= 

270 ra, dec: float 

271 The zenith and azimuthal angles in the sky frame. 

272 """ 

273 theta, phi = zenith_azimuth_to_theta_phi(zenith, azimuth, ifos) 

274 gmst = greenwich_mean_sidereal_time(geocent_time) 

275 ra, dec = theta_phi_to_ra_dec(theta, phi, gmst) 

276 ra = ra % (2 * np.pi) 

277 return ra, dec 

278 

279 

280def get_event_time(event): 

281 """ 

282 Get the merger time for known GW events using the gwosc package 

283 

284 Parameters 

285 ---------- 

286 event: str 

287 Event name, e.g. GW150914 

288 

289 Returns 

290 ------- 

291 event_time: float 

292 Merger time 

293 

294 Raises 

295 ------ 

296 ImportError 

297 If the gwosc package is not installed 

298 ValueError 

299 If the event is not in the gwosc dataset 

300 """ 

301 try: 

302 from gwosc import datasets 

303 except ImportError: 

304 raise ImportError("You do not have the gwosc package installed") 

305 

306 return datasets.event_gps(event) 

307 

308 

309def get_open_strain_data( 

310 name, start_time, end_time, outdir, cache=False, buffer_time=0, **kwargs): 

311 """ A function which accesses the open strain data 

312 

313 This uses `gwpy` to download the open data and then saves a cached copy for 

314 later use 

315 

316 Parameters 

317 ========== 

318 name: str 

319 The name of the detector to get data for 

320 start_time, end_time: float 

321 The GPS time of the start and end of the data 

322 outdir: str 

323 The output directory to place data in 

324 cache: bool 

325 If true, cache the data 

326 buffer_time: float 

327 Time to add to the beginning and end of the segment. 

328 **kwargs: 

329 Passed to `gwpy.timeseries.TimeSeries.fetch_open_data` 

330 

331 Returns 

332 ======= 

333 strain: gwpy.timeseries.TimeSeries 

334 The object containing the strain data. If the connection to the open-data server 

335 fails, this function returns `None`. 

336 

337 """ 

338 from gwpy.timeseries import TimeSeries 

339 filename = '{}/{}_{}_{}.txt'.format(outdir, name, start_time, end_time) 

340 

341 if buffer_time < 0: 

342 raise ValueError("buffer_time < 0") 

343 start_time = start_time - buffer_time 

344 end_time = end_time + buffer_time 

345 

346 if os.path.isfile(filename) and cache: 

347 logger.info('Using cached data from {}'.format(filename)) 

348 strain = TimeSeries.read(filename) 

349 else: 

350 logger.info('Fetching open data from {} to {} with buffer time {}' 

351 .format(start_time, end_time, buffer_time)) 

352 try: 

353 strain = TimeSeries.fetch_open_data(name, start_time, end_time, **kwargs) 

354 logger.info('Saving cache of data to {}'.format(filename)) 

355 strain.write(filename) 

356 except Exception as e: 

357 logger.info("Unable to fetch open data, see debug for detailed info") 

358 logger.info("Call to gwpy.timeseries.TimeSeries.fetch_open_data returned {}" 

359 .format(e)) 

360 strain = None 

361 

362 return strain 

363 

364 

365def read_frame_file(file_name, start_time, end_time, resample=None, channel=None, buffer_time=0, **kwargs): 

366 """ A function which accesses the open strain data 

367 

368 This uses `gwpy` to download the open data and then saves a cached copy for 

369 later use 

370 

371 Parameters 

372 ========== 

373 file_name: str 

374 The name of the frame to read 

375 start_time, end_time: float 

376 The GPS time of the start and end of the data 

377 buffer_time: float 

378 Read in data with `t1-buffer_time` and `end_time+buffer_time` 

379 channel: str 

380 The name of the channel being searched for, some standard channel names are attempted 

381 if channel is not specified or if specified channel is not found. 

382 resample: float 

383 The sampling frequency to use for the TimeSeries in Hz 

384 **kwargs: 

385 Passed to `gwpy.timeseries.TimeSeries.fetch_open_data` 

386 

387 Returns 

388 ======= 

389 strain: gwpy.timeseries.TimeSeries 

390 

391 """ 

392 from gwpy.timeseries import TimeSeries 

393 loaded = False 

394 strain = None 

395 

396 if channel is not None: 

397 try: 

398 strain = TimeSeries.read(source=file_name, channel=channel, start=start_time, end=end_time, **kwargs) 

399 loaded = True 

400 logger.info('Successfully loaded {}.'.format(channel)) 

401 except (RuntimeError, ValueError): 

402 logger.warning('Channel {} not found. Trying preset channel names'.format(channel)) 

403 

404 if loaded is False: 

405 ligo_channel_types = ['GDS-CALIB_STRAIN', 'DCS-CALIB_STRAIN_C01', 'DCS-CALIB_STRAIN_C02', 

406 'DCH-CLEAN_STRAIN_C02', 'GWOSC-16KHZ_R1_STRAIN', 

407 'GWOSC-4KHZ_R1_STRAIN'] 

408 virgo_channel_types = ['Hrec_hoft_V1O2Repro2A_16384Hz', 'FAKE_h_16384Hz_4R', 

409 'GWOSC-16KHZ_R1_STRAIN', 'GWOSC-4KHZ_R1_STRAIN'] 

410 channel_types = dict(H1=ligo_channel_types, L1=ligo_channel_types, V1=virgo_channel_types) 

411 for detector in channel_types.keys(): 

412 for channel_type in channel_types[detector]: 

413 if loaded: 

414 break 

415 channel = '{}:{}'.format(detector, channel_type) 

416 try: 

417 strain = TimeSeries.read(source=file_name, channel=channel, start=start_time, end=end_time, 

418 **kwargs) 

419 loaded = True 

420 logger.info('Successfully read strain data for channel {}.'.format(channel)) 

421 except (RuntimeError, ValueError): 

422 pass 

423 

424 if loaded: 

425 if resample and (strain.sample_rate.value != resample): 

426 strain = strain.resample(resample) 

427 return strain 

428 else: 

429 logger.warning('No data loaded.') 

430 return None 

431 

432 

433def get_gracedb(gracedb, outdir, duration, calibration, detectors, query_types=None, server=None): 

434 """ 

435 Download information about a trigger from GraceDb and create cache files 

436 for finding gravitational-wave strain data. 

437 

438 Parameters 

439 ---------- 

440 gracedb: str 

441 The GraceDb ID of the trigger. 

442 outdir: str 

443 Directory to write output. 

444 duration: int 

445 Duration of data to find about the trigger time (units: :code:`s`). 

446 calibration: int 

447 Calibration label of the data, should be one of :code:`1, 2`. 

448 detectors: list 

449 The detectors to look for data for. 

450 query_types: str 

451 The LDR query type 

452 server: str 

453 The LIGO datafind server to look for data on. 

454 

455 Returns 

456 ------- 

457 candidate: dict 

458 Information downloaded from GraceDb about the trigger. 

459 cache_files: list 

460 List of cache filenames, one per interferometer. 

461 """ 

462 candidate = gracedb_to_json(gracedb, outdir=outdir) 

463 trigger_time = candidate['gpstime'] 

464 gps_start_time = trigger_time - duration 

465 cache_files = [] 

466 if query_types is None: 

467 query_types = [None] * len(detectors) 

468 for i, det in enumerate(detectors): 

469 output_cache_file = gw_data_find( 

470 det, gps_start_time, duration, calibration, 

471 outdir=outdir, query_type=query_types[i], server=server) 

472 cache_files.append(output_cache_file) 

473 return candidate, cache_files 

474 

475 

476def gracedb_to_json(gracedb, cred=None, service_url='https://gracedb.ligo.org/api/', outdir=None): 

477 """ Script to download a GraceDB candidate 

478 

479 Parameters 

480 ========== 

481 gracedb: str 

482 The UID of the GraceDB candidate 

483 cred: 

484 Credentials for authentications, see ligo.gracedb.rest.GraceDb 

485 service_url: 

486 The url of the GraceDB candidate 

487 GraceDB 'https://gracedb.ligo.org/api/' (default) 

488 GraceDB-playground 'https://gracedb-playground.ligo.org/api/' 

489 outdir: str, optional 

490 If given, a string identfying the location in which to store the json 

491 """ 

492 logger.info( 

493 'Starting routine to download GraceDb candidate {}'.format(gracedb)) 

494 from ligo.gracedb.rest import GraceDb 

495 

496 logger.info('Initialise client and attempt to download') 

497 logger.info('Fetching from {}'.format(service_url)) 

498 try: 

499 client = GraceDb(cred=cred, service_url=service_url) 

500 except IOError: 

501 raise ValueError( 

502 'Failed to authenticate with gracedb: check your X509 ' 

503 'certificate is accessible and valid') 

504 try: 

505 candidate = client.event(gracedb) 

506 logger.info('Successfully downloaded candidate') 

507 except Exception as e: 

508 raise ValueError("Unable to obtain GraceDB candidate, exception: {}".format(e)) 

509 

510 json_output = candidate.json() 

511 

512 if outdir is not None: 

513 check_directory_exists_and_if_not_mkdir(outdir) 

514 outfilepath = os.path.join(outdir, '{}.json'.format(gracedb)) 

515 logger.info('Writing candidate to {}'.format(outfilepath)) 

516 with open(outfilepath, 'w') as outfile: 

517 json.dump(json_output, outfile, indent=2) 

518 

519 return json_output 

520 

521 

522def gw_data_find(observatory, gps_start_time, duration, calibration, 

523 outdir='.', query_type=None, server=None): 

524 """ Builds a gw_data_find call and process output 

525 

526 Parameters 

527 ========== 

528 observatory: str, {H1, L1, V1} 

529 Observatory description 

530 gps_start_time: float 

531 The start time in gps to look for data 

532 duration: int 

533 The duration (integer) in s 

534 calibration: int {1, 2} 

535 Use C01 or C02 calibration 

536 outdir: string 

537 A path to the directory where output is stored 

538 query_type: string 

539 The LDRDataFind query type 

540 

541 Returns 

542 ======= 

543 output_cache_file: str 

544 Path to the output cache file 

545 

546 """ 

547 logger.info('Building gw_data_find command line') 

548 

549 observatory_lookup = dict(H1='H', L1='L', V1='V') 

550 observatory_code = observatory_lookup[observatory] 

551 

552 if query_type is None: 

553 logger.warning('No query type provided. This may prevent data from being read.') 

554 if observatory_code == 'V': 

555 query_type = 'V1Online' 

556 else: 

557 query_type = '{}_HOFT_C0{}'.format(observatory, calibration) 

558 

559 logger.info('Using LDRDataFind query type {}'.format(query_type)) 

560 

561 cache_file = '{}-{}_CACHE-{}-{}.lcf'.format( 

562 observatory, query_type, gps_start_time, duration) 

563 output_cache_file = os.path.join(outdir, cache_file) 

564 

565 gps_end_time = gps_start_time + duration 

566 if server is None: 

567 server = os.environ.get("LIGO_DATAFIND_SERVER") 

568 if server is None: 

569 logger.warning('No data_find server found, defaulting to CIT server. ' 

570 'To run on other clusters, pass the output of ' 

571 '`echo $LIGO_DATAFIND_SERVER`') 

572 server = 'ldr.ldas.cit:80' 

573 

574 cl_list = ['gw_data_find'] 

575 cl_list.append('--observatory {}'.format(observatory_code)) 

576 cl_list.append('--gps-start-time {}'.format(int(np.floor(gps_start_time)))) 

577 cl_list.append('--gps-end-time {}'.format(int(np.ceil(gps_end_time)))) 

578 cl_list.append('--type {}'.format(query_type)) 

579 cl_list.append('--output {}'.format(output_cache_file)) 

580 cl_list.append('--server {}'.format(server)) 

581 cl_list.append('--url-type file') 

582 cl_list.append('--lal-cache') 

583 cl = ' '.join(cl_list) 

584 run_commandline(cl) 

585 return output_cache_file 

586 

587 

588def convert_args_list_to_float(*args_list): 

589 """ Converts inputs to floats, returns a list in the same order as the input""" 

590 try: 

591 args_list = [float(arg) for arg in args_list] 

592 except ValueError: 

593 raise ValueError("Unable to convert inputs to floats") 

594 return args_list 

595 

596 

597def lalsim_SimInspiralTransformPrecessingNewInitialConditions( 

598 theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2, 

599 reference_frequency, phase): 

600 from lalsimulation import SimInspiralTransformPrecessingNewInitialConditions 

601 

602 args_list = convert_args_list_to_float( 

603 theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2, 

604 reference_frequency, phase) 

605 

606 return SimInspiralTransformPrecessingNewInitialConditions(*args_list) 

607 

608 

609@lru_cache(maxsize=10) 

610def lalsim_GetApproximantFromString(waveform_approximant): 

611 from lalsimulation import GetApproximantFromString 

612 if isinstance(waveform_approximant, str): 

613 return GetApproximantFromString(waveform_approximant) 

614 else: 

615 raise ValueError("waveform_approximant must be of type str") 

616 

617 

618def lalsim_SimInspiralFD( 

619 mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, 

620 spin_2z, luminosity_distance, iota, phase, 

621 longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency, 

622 minimum_frequency, maximum_frequency, reference_frequency, 

623 waveform_dictionary, approximant): 

624 """ 

625 Safely call lalsimulation.SimInspiralFD 

626 

627 Parameters 

628 ========== 

629 phase: float, int 

630 mass_1: float, int 

631 mass_2: float, int 

632 spin_1x: float, int 

633 spin_1y: float, int 

634 spin_1z: float, int 

635 spin_2x: float, int 

636 spin_2y: float, int 

637 spin_2z: float, int 

638 reference_frequency: float, int 

639 luminosity_distance: float, int 

640 iota: float, int 

641 longitude_ascending_nodes: float, int 

642 eccentricity: float, int 

643 mean_per_ano: float, int 

644 delta_frequency: float, int 

645 minimum_frequency: float, int 

646 maximum_frequency: float, int 

647 waveform_dictionary: None, lal.Dict 

648 approximant: int, str 

649 """ 

650 from lalsimulation import SimInspiralFD 

651 

652 args = convert_args_list_to_float( 

653 mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, 

654 luminosity_distance, iota, phase, longitude_ascending_nodes, 

655 eccentricity, mean_per_ano, delta_frequency, minimum_frequency, 

656 maximum_frequency, reference_frequency) 

657 

658 approximant = _get_lalsim_approximant(approximant) 

659 

660 return SimInspiralFD(*args, waveform_dictionary, approximant) 

661 

662 

663def lalsim_SimInspiralChooseFDWaveform( 

664 mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, 

665 spin_2z, luminosity_distance, iota, phase, 

666 longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency, 

667 minimum_frequency, maximum_frequency, reference_frequency, 

668 waveform_dictionary, approximant): 

669 """ 

670 Safely call lalsimulation.SimInspiralChooseFDWaveform 

671 

672 Parameters 

673 ========== 

674 phase: float, int 

675 mass_1: float, int 

676 mass_2: float, int 

677 spin_1x: float, int 

678 spin_1y: float, int 

679 spin_1z: float, int 

680 spin_2x: float, int 

681 spin_2y: float, int 

682 spin_2z: float, int 

683 reference_frequency: float, int 

684 luminosity_distance: float, int 

685 iota: float, int 

686 longitude_ascending_nodes: float, int 

687 eccentricity: float, int 

688 mean_per_ano: float, int 

689 delta_frequency: float, int 

690 minimum_frequency: float, int 

691 maximum_frequency: float, int 

692 waveform_dictionary: None, lal.Dict 

693 approximant: int, str 

694 """ 

695 from lalsimulation import SimInspiralChooseFDWaveform 

696 

697 args = convert_args_list_to_float( 

698 mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, 

699 luminosity_distance, iota, phase, longitude_ascending_nodes, 

700 eccentricity, mean_per_ano, delta_frequency, minimum_frequency, 

701 maximum_frequency, reference_frequency) 

702 

703 approximant = _get_lalsim_approximant(approximant) 

704 

705 return SimInspiralChooseFDWaveform(*args, waveform_dictionary, approximant) 

706 

707 

708@lru_cache(maxsize=10) 

709def _get_lalsim_approximant(approximant): 

710 if isinstance(approximant, int): 

711 pass 

712 elif isinstance(approximant, str): 

713 approximant = lalsim_GetApproximantFromString(approximant) 

714 else: 

715 raise ValueError("approximant not an int") 

716 return approximant 

717 

718 

719def lalsim_SimInspiralChooseFDWaveformSequence( 

720 phase, mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, 

721 spin_2z, reference_frequency, luminosity_distance, iota, 

722 waveform_dictionary, approximant, frequency_array): 

723 """ 

724 Safely call lalsimulation.SimInspiralChooseFDWaveformSequence 

725 

726 Parameters 

727 ========== 

728 phase: float, int 

729 mass_1: float, int 

730 mass_2: float, int 

731 spin_1x: float, int 

732 spin_1y: float, int 

733 spin_1z: float, int 

734 spin_2x: float, int 

735 spin_2y: float, int 

736 spin_2z: float, int 

737 reference_frequency: float, int 

738 luminosity_distance: float, int 

739 iota: float, int 

740 waveform_dictionary: None, lal.Dict 

741 approximant: int, str 

742 frequency_array: np.ndarray, lal.REAL8Vector 

743 """ 

744 from lal import REAL8Vector, CreateREAL8Vector 

745 from lalsimulation import SimInspiralChooseFDWaveformSequence 

746 

747 [mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, 

748 luminosity_distance, iota, phase, reference_frequency] = convert_args_list_to_float( 

749 mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, 

750 luminosity_distance, iota, phase, reference_frequency) 

751 

752 if isinstance(approximant, int): 

753 pass 

754 elif isinstance(approximant, str): 

755 approximant = lalsim_GetApproximantFromString(approximant) 

756 else: 

757 raise ValueError("approximant not an int") 

758 

759 if not isinstance(frequency_array, REAL8Vector): 

760 old_frequency_array = frequency_array.copy() 

761 frequency_array = CreateREAL8Vector(len(old_frequency_array)) 

762 frequency_array.data = old_frequency_array 

763 

764 return SimInspiralChooseFDWaveformSequence( 

765 phase, mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, 

766 spin_2z, reference_frequency, luminosity_distance, iota, 

767 waveform_dictionary, approximant, frequency_array) 

768 

769 

770def lalsim_SimInspiralWaveformParamsInsertTidalLambda1( 

771 waveform_dictionary, lambda_1): 

772 from lalsimulation import SimInspiralWaveformParamsInsertTidalLambda1 

773 try: 

774 lambda_1 = float(lambda_1) 

775 except ValueError: 

776 raise ValueError("Unable to convert lambda_1 to float") 

777 

778 return SimInspiralWaveformParamsInsertTidalLambda1( 

779 waveform_dictionary, lambda_1) 

780 

781 

782def lalsim_SimInspiralWaveformParamsInsertTidalLambda2( 

783 waveform_dictionary, lambda_2): 

784 from lalsimulation import SimInspiralWaveformParamsInsertTidalLambda2 

785 try: 

786 lambda_2 = float(lambda_2) 

787 except ValueError: 

788 raise ValueError("Unable to convert lambda_2 to float") 

789 

790 return SimInspiralWaveformParamsInsertTidalLambda2( 

791 waveform_dictionary, lambda_2) 

792 

793 

794def lalsim_SimNeutronStarEOS4ParamSDGammaCheck(g0, g1, g2, g3): 

795 from lalsimulation import SimNeutronStarEOS4ParamSDGammaCheck 

796 try: 

797 g0 = float(g0) 

798 g1 = float(g1) 

799 g2 = float(g2) 

800 g3 = float(g3) 

801 except ValueError: 

802 raise ValueError("Unable to convert EOS spectral parameters to floats") 

803 except TypeError: 

804 raise TypeError("Unable to convert EOS spectral parameters to floats") 

805 

806 return SimNeutronStarEOS4ParamSDGammaCheck(g0, g1, g2, g3) 

807 

808 

809def lalsim_SimNeutronStarEOS4ParameterSpectralDecomposition(g0, g1, g2, g3): 

810 from lalsimulation import SimNeutronStarEOS4ParameterSpectralDecomposition 

811 try: 

812 g0 = float(g0) 

813 g1 = float(g1) 

814 g2 = float(g2) 

815 g3 = float(g3) 

816 except ValueError: 

817 raise ValueError("Unable to convert EOS spectral parameters to floats") 

818 except TypeError: 

819 raise TypeError("Unable to convert EOS spectral parameters to floats") 

820 

821 return SimNeutronStarEOS4ParameterSpectralDecomposition(g0, g1, g2, g3) 

822 

823 

824def lalsim_SimNeutronStarEOS4ParamSDViableFamilyCheck(g0, g1, g2, g3): 

825 from lalsimulation import SimNeutronStarEOS4ParamSDViableFamilyCheck 

826 try: 

827 g0 = float(g0) 

828 g1 = float(g1) 

829 g2 = float(g2) 

830 g3 = float(g3) 

831 except ValueError: 

832 raise ValueError("Unable to convert EOS spectral parameters to floats") 

833 except TypeError: 

834 raise TypeError("Unable to convert EOS spectral parameters to floats") 

835 

836 return SimNeutronStarEOS4ParamSDViableFamilyCheck(g0, g1, g2, g3) 

837 

838 

839def lalsim_SimNeutronStarEOS3PieceDynamicPolytrope(g0, log10p1_si, g1, log10p2_si, g2): 

840 from lalsimulation import SimNeutronStarEOS3PieceDynamicPolytrope 

841 try: 

842 g0 = float(g0) 

843 g1 = float(g1) 

844 g2 = float(g2) 

845 log10p1_si = float(log10p1_si) 

846 log10p2_si = float(log10p2_si) 

847 except ValueError: 

848 raise ValueError("Unable to convert EOS polytrope parameters to floats") 

849 except TypeError: 

850 raise TypeError("Unable to convert EOS polytrope parameters to floats") 

851 

852 return SimNeutronStarEOS3PieceDynamicPolytrope(g0, log10p1_si, g1, log10p2_si, g2) 

853 

854 

855def lalsim_SimNeutronStarEOS3PieceCausalAnalytic(v1, log10p1_si, v2, log10p2_si, v3): 

856 from lalsimulation import SimNeutronStarEOS3PieceCausalAnalytic 

857 try: 

858 v1 = float(v1) 

859 v2 = float(v2) 

860 v3 = float(v3) 

861 log10p1_si = float(log10p1_si) 

862 log10p2_si = float(log10p2_si) 

863 except ValueError: 

864 raise ValueError("Unable to convert EOS causal parameters to floats") 

865 except TypeError: 

866 raise TypeError("Unable to convert EOS causal parameters to floats") 

867 

868 return SimNeutronStarEOS3PieceCausalAnalytic(v1, log10p1_si, v2, log10p2_si, v3) 

869 

870 

871def lalsim_SimNeutronStarEOS3PDViableFamilyCheck(p0, log10p1_si, p1, log10p2_si, p2, causal): 

872 from lalsimulation import SimNeutronStarEOS3PDViableFamilyCheck 

873 try: 

874 p0 = float(p0) 

875 p1 = float(p1) 

876 p2 = float(p2) 

877 log10p1_si = float(log10p1_si) 

878 log10p2_si = float(log10p2_si) 

879 causal = int(causal) 

880 except ValueError: 

881 raise ValueError("Unable to convert EOS parameters to floats or int") 

882 except TypeError: 

883 raise TypeError("Unable to convert EOS parameters to floats or int") 

884 

885 return SimNeutronStarEOS3PDViableFamilyCheck(p0, log10p1_si, p1, log10p2_si, p2, causal) 

886 

887 

888def lalsim_CreateSimNeutronStarFamily(eos): 

889 from lalsimulation import CreateSimNeutronStarFamily 

890 

891 return CreateSimNeutronStarFamily(eos) 

892 

893 

894def lalsim_SimNeutronStarEOSMaxPseudoEnthalpy(eos): 

895 from lalsimulation import SimNeutronStarEOSMaxPseudoEnthalpy 

896 

897 return SimNeutronStarEOSMaxPseudoEnthalpy(eos) 

898 

899 

900def lalsim_SimNeutronStarEOSSpeedOfSoundGeometerized(max_pseudo_enthalpy, eos): 

901 from lalsimulation import SimNeutronStarEOSSpeedOfSoundGeometerized 

902 try: 

903 max_pseudo_enthalpy = float(max_pseudo_enthalpy) 

904 except ValueError: 

905 raise ValueError("Unable to convert max_pseudo_enthalpy to float.") 

906 except TypeError: 

907 raise TypeError("Unable to convert max_pseudo_enthalpy to float.") 

908 

909 return SimNeutronStarEOSSpeedOfSoundGeometerized(max_pseudo_enthalpy, eos) 

910 

911 

912def lalsim_SimNeutronStarFamMinimumMass(fam): 

913 from lalsimulation import SimNeutronStarFamMinimumMass 

914 

915 return SimNeutronStarFamMinimumMass(fam) 

916 

917 

918def lalsim_SimNeutronStarMaximumMass(fam): 

919 from lalsimulation import SimNeutronStarMaximumMass 

920 

921 return SimNeutronStarMaximumMass(fam) 

922 

923 

924def lalsim_SimNeutronStarRadius(mass_in_SI, fam): 

925 from lalsimulation import SimNeutronStarRadius 

926 try: 

927 mass_in_SI = float(mass_in_SI) 

928 except ValueError: 

929 raise ValueError("Unable to convert mass_in_SI to float.") 

930 except TypeError: 

931 raise TypeError("Unable to convert mass_in_SI to float.") 

932 

933 return SimNeutronStarRadius(mass_in_SI, fam) 

934 

935 

936def lalsim_SimNeutronStarLoveNumberK2(mass_in_SI, fam): 

937 from lalsimulation import SimNeutronStarLoveNumberK2 

938 try: 

939 mass_in_SI = float(mass_in_SI) 

940 except ValueError: 

941 raise ValueError("Unable to convert mass_in_SI to float.") 

942 except TypeError: 

943 raise TypeError("Unable to convert mass_in_SI to float.") 

944 

945 return SimNeutronStarLoveNumberK2(mass_in_SI, fam) 

946 

947 

948def spline_angle_xform(delta_psi):