Coverage for /home/runner/work/nr-catalog-tools/nr-catalog-tools/nrcatalogtools/lvc.py: 47%
284 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-06-01 05:18 +0000
1import h5py
2import lal
3import lalsimulation as lalsim
4import numpy as np
5from pycbc.pnutils import mtotal_eta_to_mass1_mass2
6from pycbc.types import TimeSeries
7from scipy.interpolate import interp1d
10def get_lal_mode_dictionary(mode_array):
11 """
12 Get LALDict with all specified modes.
15 Parameters
16 ----------
17 mode_array: list of modes, eg [[2,2], [3,3]]
19 Returns
20 -------
21 waveform_dictionary: LALDict with all modes included
23 """
24 waveform_dictionary = lal.CreateDict()
25 mode_array_lal = lalsim.SimInspiralCreateModeArray()
26 for mode in mode_array:
27 lalsim.SimInspiralModeArrayActivateMode(mode_array_lal, mode[0], mode[1])
28 lalsim.SimInspiralWaveformParamsInsertModeArray(waveform_dictionary, mode_array_lal)
30 return waveform_dictionary
33def get_lal_mode_dictionary_from_lmax(lmax):
34 r"""
35 Get LALDict with modes derived from `lmax`.
38 Parameters
39 ----------
40 lmax: max value of :math:`\ell` to use
42 Returns
43 -------
44 waveform_dictionary: LALDict with all modes upto `lmax` included
46 """
48 mode_array = [
49 [ell, emm] for ell in range(2, lmax + 1) for emm in range(-ell, ell + 1)
50 ]
51 return get_lal_mode_dictionary(mode_array)
54def get_modes_from_lvcnr_file(path_to_file, Mtot, distance, srate, lmax=4, f_low=None):
55 r"""
56 Get individual modes from LVCNR format file.
59 Parameters
60 ==========
61 path_to_file: string
62 Path to LVCNR file
63 Mtot: float
64 Total mass (in units of MSUN) to scale the waveform to
65 distance: float
66 Luminosity Distance (in units of MPc) to scale the waveform to
67 srate: int
68 Sampling rate for the waveform
69 lmax: int
70 Max value of :math:`\ell` to use
71 (Default: 4)
72 f_low: float
73 Value of the low frequency to start waveform generation
74 Uses value given from the LVCNR file if `None` is provided
75 (Default: None)
77 Returns
78 =======
79 mass_ratio: float
80 Mass ratio derived from the LVCNR file
81 spins_args: list
82 List of spins derived from the LVCNR file
83 eccentricity: float
84 Eccentricty derived from the LVCNR file.
85 Returns `None` is eccentricity is not a number.
86 f_low: float
87 Low Frequency derived either from the file, or provided
88 in the call
89 f_ref: float
90 Reference Frequency derived from the file
91 modes: dict of pycbc TimeSeries objects
92 dict containing all the read in modes
93 """
95 with h5py.File(path_to_file) as h5file:
96 waveform_dict = get_lal_mode_dictionary_from_lmax(lmax)
98 f_low_in_file = h5file.attrs["f_lower_at_1MSUN"] / Mtot
99 f_ref = f_low_in_file
101 if f_low is None:
102 f_low = f_low_in_file
104 if h5file.attrs["Format"] < 3:
105 spin_args = lalsim.SimInspiralNRWaveformGetSpinsFromHDF5File(
106 -1, Mtot, path_to_file
107 )
108 else:
109 spin_args = lalsim.SimInspiralNRWaveformGetSpinsFromHDF5File(
110 f_low, Mtot, path_to_file
111 )
113 mass_args = list(
114 mtotal_eta_to_mass1_mass2(Mtot * lal.MSUN_SI, h5file.attrs["eta"])
115 )
117 try:
118 eccentricity = float(h5file.attrs["eccentricity"])
119 except ValueError:
120 eccentricity = None
122 values_mode_array = lalsim.SimInspiralWaveformParamsLookupModeArray(waveform_dict)
123 _, modes = lalsim.SimInspiralNRWaveformGetHlms(
124 1 / srate,
125 *mass_args,
126 distance * 1e6 * lal.PC_SI,
127 f_low,
128 f_ref,
129 *spin_args,
130 path_to_file,
131 values_mode_array,
132 )
133 mode = modes
134 return_dict = dict()
135 while 1 > 0:
136 try:
137 l, m = mode.l, mode.m
138 read_mode = mode.mode.data.data
139 return_dict[(l, m)] = TimeSeries(read_mode, 1 / srate)
140 mode = mode.next
141 except AttributeError:
142 break
144 return (
145 mass_args[1] / mass_args[0],
146 spin_args,
147 eccentricity,
148 f_low,
149 f_ref,
150 return_dict,
151 )
154def get_strain_from_lvcnr_file(
155 path_to_file, Mtot, distance, inclination, phi_ref, srate, mode_array=None
156):
157 """
158 Get full strain from LVCNR format file.
161 Parameters
162 ==========
163 path_to_file: string
164 Path to LVCNR file
165 Mtot: float
166 Total mass (in units of MSUN) to scale the waveform to
167 distance: float
168 Luminosity Distance (in units of MPc) to scale the waveform to
169 srate: int
170 Sampling rate for the waveform
171 mode_array: list
172 list of modes to be included. `None` means all modes are included.
173 (Default:None)
175 Returns
176 =======
177 UNDER CONSTRUCTION
178 """
180 longAscNodes = 0
181 eccentricity = 0
182 meanPerAno = 0
184 fixed_args = [
185 distance * lal.PC_SI * 1e6,
186 inclination,
187 phi_ref,
188 longAscNodes,
189 eccentricity,
190 meanPerAno,
191 1 / srate,
192 ]
194 with h5py.File(path_to_file) as h5file:
195 if mode_array is not None:
196 waveform_dict = get_lal_mode_dictionary(mode_array)
197 else:
198 waveform_dict = lal.CreateDict()
200 lalsim.SimInspiralWaveformParamsInsertNumRelData(waveform_dict, path_to_file)
201 f_low = h5file.attrs["f_lower_at_1MSUN"] / Mtot
203 if h5file.attrs["Format"] < 3:
204 spin_args = lalsim.SimInspiralNRWaveformGetSpinsFromHDF5File(
205 -1, Mtot, path_to_file
206 )
207 else:
208 spin_args = lalsim.SimInspiralNRWaveformGetSpinsFromHDF5File(
209 f_low, Mtot, path_to_file
210 )
212 mass_args = list(
213 mtotal_eta_to_mass1_mass2(Mtot * lal.MSUN_SI, h5file.attrs["eta"])
214 )
216 hp, hc = lalsim.SimInspiralChooseTDWaveform(
217 *mass_args,
218 *spin_args,
219 *fixed_args,
220 f_low,
221 f_low,
222 waveform_dict,
223 lalsim.GetApproximantFromString("NR_hdf5"),
224 )
226 hp_timeseries = TimeSeries(hp.data.data, fixed_args[-1])
227 hc_timeseries = TimeSeries(hc.data.data, fixed_args[-1])
229 return dict(
230 hp=hp_timeseries,
231 hc=hc_timeseries,
232 f_low=f_low,
233 mass_args=mass_args,
234 spin_args=spin_args,
235 fixed_args=fixed_args,
236 )
239def check_interp_req(h5_file=None, metadata=None, ref_time=None, avail_ref_time=None):
240 """Check if the required reference time is different from
241 the available reference time in the NR HDF5 file or the
242 simulation metadata.
244 Parameters
245 ----------
246 h5_file : file object
247 The waveform h5 file handle.
248 metadata :
249 ref_time, avail_ref_time : float
250 The use and available nr reference time.
252 Returns
253 -------
254 interp : bool
255 Whether interpolation across time is required.
256 avail_ref_time: float
257 The ref_time available in the NR HDF5 file.
258 """
260 if avail_ref_time is None:
261 # CCheck for ref time in h5 file
262 if h5_file is not None:
263 keys = list(h5_file.attrs.keys())
265 if "reference_time" in keys:
266 avail_ref_time = h5_file.attrs["reference_time"]
267 elif "ref_time" in keys:
268 avail_ref_time = h5_file.attrs["ref_time"]
269 elif "relaxed_time" in keys:
270 avail_ref_time = h5_file.attrs["relaxed_time"]
271 else:
272 print("Reference time not found in waveform h5 file.")
274 if not avail_ref_time:
275 # If not found, continue search in metadata.
276 if metadata is not None:
277 keys = list(metadata.keys())
279 if "reference_time" in keys:
280 avail_ref_time = metadata["reference_time"]
281 elif "relaxed_time" in keys:
282 avail_ref_time = metadata["relaxed_time"]
283 else:
284 print("Reference time not found in simulation metadata file.")
286 if not avail_ref_time:
287 # Then this is GT simulation!
288 print(
289 "Reference time should be computed from"
290 "the reference orbital frequency!"
291 )
293 interp = True
294 if isinstance(ref_time, float):
295 if abs(avail_ref_time - ref_time) < 1e-5:
296 interp = False
297 else:
298 return interp, avail_ref_time
301def get_ref_freq_from_ref_time(h5_file, ref_time):
302 """Get the reference frequency from reference time
304 Parameters
305 ----------
306 h5_file : file object
307 The waveform h5 file handle.
308 ref_time : float
309 Reference time.
310 Returns
311 -------
312 f_ref : float
313 Reference frequency.
314 """
316 time, freq = h5_file.attrs["Omega-vs-time"]
318 ref_freq = interp1d(time, freq, kind="cubic")[ref_time]
320 return ref_freq
323def get_ref_time_from_ref_freq(h5_file, ref_freq):
324 """Get the reference time from reference frequency
326 Parameters
327 ----------
328 h5_file : file object
329 The waveform h5 file handle.
330 ref_freq : float
331 The reference frequency.
332 Returns
333 -------
334 fTime : float
335 Reference time.
336 """
338 time, freq = h5_file.attrs["Omega-vs-time"]
340 ref_time = interp1d(freq, time, kind="cubic")[ref_freq]
342 return ref_time
345def check_nr_attrs(
346 sim_metadata_object,
347 req_attrs=["LNhatx", "LNhaty", "LNhatz", "nhatx", "nhaty", "nhatz"],
348):
349 """Check if the NR h5 file or a simulation metadata dictionary
350 contains all the attributes required.
352 Parameters
353 ----------
354 sim_metadata_object : h5 file object, dict
355 The NR h5py file handle or simulation metadata.
356 req_attrs : list
357 A list of attribute keys.
358 Returns
359 -------
360 present : bool
361 Whether or not all specified attributes are present.
362 absent_attrs : list
363 The attributes that are absent.
364 """
365 if isinstance(sim_metadata_object, h5py.File):
366 all_attrs = list(sim_metadata_object.attrs.keys())
368 elif isinstance(sim_metadata_object, dict):
369 all_attrs = list(sim_metadata_object.keys())
370 else:
371 raise TypeError("Please supply an open h5py file handle or a dictionary")
373 absent_attrs = []
374 present = True
376 for item in req_attrs:
377 if item not in all_attrs:
378 present = False
379 absent_attrs.append(item)
381 return present, absent_attrs
384def get_interp_ref_values_from_h5_file(h5_file, req_ts_attrs, ref_time):
385 """Get the interpolated reference values at a given reference time
386 from the NR HDF5 File
388 Parameters
389 ----------
390 h5_file : file object
391 The waveform h5 file handle.
392 req_ts_attrs : list
393 A list of attribute keys.
394 ref_time : float
395 Reference time.
396 Returns
397 -------
398 params : dict
399 The parameter values at the reference time.
400 """
402 params = {}
404 for key in req_ts_attrs:
405 time, var = h5_file.attrs[key]
406 RefVal = interp1d(time, var, kind="cubic")[ref_time]
407 params.update({key: RefVal})
408 return params
411def get_ref_vals(
412 sim_metadata_object,
413 req_attrs=["LNhatx", "LNhaty", "LNhatz", "nhatx", "nhaty", "nhatz"],
414):
415 """Get the reference values from a NR HDF5 file
416 or a simulation metadata dictionary.
418 Parameters
419 ----------
420 sim_metadata_object : h5 file object, dict
421 The NR h5py file handle or
422 the simulation metadata.
423 req_attrs : list
424 A list of attribute keys.
425 Returns
426 -------
427 params : dict
428 The parameter values at the reference time.
429 """
430 if isinstance(sim_metadata_object, h5py.File):
431 source = sim_metadata_object.attrs
433 elif isinstance(sim_metadata_object, dict):
434 source = sim_metadata_object
435 else:
436 raise TypeError("Please supply an open h5py file handle or a dictionary")
438 params = {}
440 for key in req_attrs:
441 RefVal = source[key]
442 params.update({key: RefVal})
443 return params
446def compute_lal_source_frame_from_sxs_metadata(sim_metadata):
447 """Compute the LAL source frame vectors at the
448 available reference time from the SXS simulation
449 metadata.
451 Parameters
452 ----------
453 sim_metadata : dict
454 The NR sim_metadata.
455 Returns
456 -------
457 params : dict
458 A dictionary containing the LAL source frame
459 vectors.
460 """
461 # Masses
462 M1 = sim_metadata["reference_mass1"]
463 M2 = sim_metadata["reference_mass2"]
464 M = M1 + M2
466 # Orbital anglular frequency
467 Omega = np.array(sim_metadata["reference_orbital_frequency"])
469 # Positions
470 Pos1 = np.array(sim_metadata["reference_position1"])
471 Pos2 = np.array(sim_metadata["reference_position2"])
473 # CoM position
474 CoMPos = (M1 * Pos1 + M2 * Pos2) / (M)
476 Pos1 = Pos1 - CoMPos
477 Pos2 = Pos2 - CoMPos
479 # This should point from
480 # smaller to larger BH
481 # Assumed here is 1 for smaller
482 # BH
483 DPos = Pos2 - Pos1
485 # Oribital angular momentum, Eucledian
486 P1 = M1 * np.cross(Omega, Pos1)
487 P2 = M2 * np.cross(Omega, Pos2)
489 Lbar = np.cross(Pos1, P1) + np.cross(Pos2, P2)
491 # Orbital normal LNhat
492 LNhat = Lbar / np.linalg.norm(Lbar)
493 LNhatx, LNhaty, LNhatz = LNhat
495 # Position vector nhat
496 nhat = (DPos) / np.linalg.norm(DPos)
497 nhatx, nhaty, nhatz = nhat
499 params = {
500 "LNhatx": LNhatx,
501 "LNhaty": LNhaty,
502 "LNhatz": LNhatz,
503 "nhatx": nhatx,
504 "nhaty": nhaty,
505 "nhatz": nhatz,
506 }
508 return params
511def compute_lal_source_frame_by_interp(h5_file, req_ts_attrs, t_ref):
512 """
513 Compute the LAL source frame vectors at a given reference time
514 by interpolation of time series data.
516 Parameters
517 ----------
518 h5_file : h5 file object
519 The NR h5py file handle that contains the simulation metadata.
520 t_ref : float
521 The reference time.
522 Returns
523 -------
524 params : dict
525 The LAL source frame vectors at the reference time.
527 """
528 # For reference: attributes required for interpolation.
529 # req_ts_attrs = ['LNhatx-vs-time', 'LNhaty-vs-time', 'LNhatz-vs-time', \
530 # 'position1x-vs-time', 'position1y-vs-time', 'position1z-vs-time', \
531 # 'position2x-vs-time', 'position2y-vs-time', 'position2z-vs-time']
533 ref_params = get_interp_ref_values_from_h5_file(h5_file, req_ts_attrs, t_ref)
535 # r21 vec
536 r21_x = ref_params["position1x-vs-time"] - ref_params["position2x-vs-time"]
537 r21_y = ref_params["position1y-vs-time"] - ref_params["position2y-vs-time"]
538 r21_z = ref_params["position1z-vs-time"] - ref_params["position2z-vs-time"]
540 # Position unit vector nhat
541 dPos = np.array([r21_x, r21_y, r21_z])
542 nhat = dPos / np.linalg.norm(dPos)
543 nhatx, nhaty, nhatz = nhat
545 # Orbital unit normal LNhat
546 LNhatx = ref_params["LNhatx-vs-time"]
547 LNhaty = ref_params["LNhaty-vs-time"]
548 LNhatz = ref_params["LNhatz-vs-time"]
550 LNhat = np.array([LNhatx, LNhaty, LNhatz])
551 LNhat = LNhat / np.linalg.norm(LNhat)
552 # LNhatx, LNhaty, LNhatz = LNhat
554 # params = {"LNhat": LNhat, "nhat": nhat}
555 params = {
556 "LNhatx": LNhatx,
557 "LNhaty": LNhaty,
558 "LNhatz": LNhatz,
559 "nhatx": nhatx,
560 "nhaty": nhaty,
561 "nhatz": nhatz,
562 }
564 return params
567def normalize_metadata(sim_metadata):
568 """Ensure that the keys of the metadata are
569 as required.
571 Parameters
572 ----------
573 sim_metadata : dict
574 The NR sim_metadata.
576 Returns
577 -------
578 norm_sim_metadata : dict
579 The normalized simulation metadata
580 """
582 norm_sim_metadata = {}
584 for key, val in sim_metadata.items():
585 norm_sim_metadata.update({key.replace("relaxed-", ""): val})
587 return norm_sim_metadata
590def get_ref_time_from_metadata(sim_metadata):
591 """Get the reference time of definition of the LAL
592 frame from the simulation metadata, if available.
594 Parameters
595 ----------
596 sim_metadata : dict
597 The simulation metadata.
599 Returns
600 -------
601 t_ref : float
602 The reference time
603 """
605 MdataKeys = list(sim_metadata.keys())
607 if "relaxed-time" in MdataKeys:
608 # RIT style
609 t_ref = sim_metadata["relaxed-time"]
610 elif "reference_time" in MdataKeys:
611 # SXS Style
612 t_ref = sim_metadata["reference_time"]
613 else:
614 t_ref = -1
616 return t_ref
619def transform_spins_nr_to_lal(nrSpin1, nrSpin2, n_hat, ln_hat):
620 """Trnasform the spins of the NR simulation from the
621 NR frame to the frame.
622 Parameters
623 ---------
624 nrSpin1, nrSpin2 : list
625 A list of the components of the spins of the objects.
626 nhat, ln_hat : list
627 A list of the components of the unit vectors of the objects,
628 against which the components of the spins are specified.
629 Returns
630 -------
631 S1, S2 : list
632 The transformed spins in LAL frame.
633 """
634 nrSpin1x, nrSpin1y, nrSpin1z = nrSpin1
635 nrSpin2x, nrSpin2y, nrSpin2z = nrSpin2
637 n_hat_x, n_hat_y, n_hat_z = n_hat
638 ln_hat_x, ln_hat_y, ln_hat_z = ln_hat
640 S1x = nrSpin1x * n_hat_x + nrSpin1y * n_hat_y + nrSpin1z * n_hat_z
642 S1y = (
643 nrSpin1x * (-ln_hat_z * n_hat_y + ln_hat_y * n_hat_z)
644 + nrSpin1y * (ln_hat_z * n_hat_x - ln_hat_x * n_hat_z)
645 + nrSpin1z * (-ln_hat_y * n_hat_x + ln_hat_x * n_hat_y)
646 )
648 S1z = nrSpin1x * ln_hat_x + nrSpin1y * ln_hat_y + nrSpin1z * ln_hat_z
650 S2x = nrSpin2x * n_hat_x + nrSpin2y * n_hat_y + nrSpin2z * n_hat_z
651 S2y = (
652 nrSpin2x * (-ln_hat_z * n_hat_y + ln_hat_y * n_hat_z)
653 + nrSpin2y * (ln_hat_z * n_hat_x - ln_hat_x * n_hat_z)
654 + nrSpin2z * (-ln_hat_y * n_hat_x + ln_hat_x * n_hat_y)
655 )
657 S2z = nrSpin2x * ln_hat_x + nrSpin2y * ln_hat_y + nrSpin2z * ln_hat_z
659 S1 = [S1x, S1y, S1z]
660 S2 = [S2x, S2y, S2z]
662 return S1, S2
665def get_nr_to_lal_rotation_angles(
666 h5_file, sim_metadata, inclination, phi_ref=0, f_ref=None, t_ref=None, tol=1e-6
667):
668 r"""Get the angular coordinates :math:`\theta, \phi`
669 and the rotation angle :math:`\alpha` from the H5 file
671 Parameters
672 ----------
673 h5_file : file object
674 The waveform h5 file handle.
676 inclination : float
677 The inclination angle.
678 phi_ref : float
679 The orbital phase at reference time.
680 Fref, t_ref : float, optional
681 The reference orbital frequency or time
683 sim_metadata : dict
684 The sim_metadata of the waveform file.
685 tol : float
686 The tolerance to use to allow floating point
687 representation errors.
689 Returns
690 -------
691 angles : dict
692 The angular corrdinates Theta, Psi, and the rotation angle Alpha.
693 If available, this also contains the reference time and frequency.
695 Notes
696 -----
698 Variable definitions.
700 theta : Returned inclination angle of source in NR coordinates.
701 psi : Returned azimuth angle of source in NR coordinates.
702 alpha: Returned polarisation angle.
703 h5_file: h5py object of the NR HDF5 file.
704 inclination: inclination of source in LAL source frame.
705 phi_ref: Orbital reference phase.
706 t_ref : Reference time. -1 or None indicates it was not found in the sim_metadata.
707 f_ref: Reference frequency.
709 The reference epoch is defined close to the beginning of the simulation.
710 """
712 # Compute the angles necessary to rotate from the intrinsic NR source frame
713 # into the LAL frame. See DCC-T1600045 for details.
715 # Following section IV of DCC-T1600045
716 # Step 1: Define Phi = phiref
717 orb_phase = phi_ref
719 ##########################################
720 # Step 2: Compute Zref
721 # 2.1 : Check if interpolation is required in IntReq
722 # 2.2 : Get/ Compute the basis vectors of the LAL
723 # frame.
724 # 2.2.1 : If IntReq=yes, given the reference time, interpolate and get
725 # the required basis vectors in the LAL source frame.
726 # 2.2.2 : If no, then check for the default values of the
727 # LAL frame basis vectors in the h5_file
728 # 2.2.3 : If the h5_file does not contain the required default
729 # vectors, then raise an error.
730 # 2.3 : Carryout vector math to get Zref.
731 # 2.1: Compute LN_hat from file. LN_hat = direction of orbital ang. mom.
732 # 2.2: Compute n_hat from file. n_hat = direction from object 2 to object 1
734 ###########################################
735 # Cases
736 ###########################################
737 # req_def_attrs = ["LNhatx", "LNhaty", "LNhatz", "nhatx", "nhaty", "nhatz"]
739 # SXS default attributes required
740 # for computing the LAL source frame.
741 req_def_attrs_sxs = [
742 "reference_time",
743 "reference_mass1",
744 "reference_mass2",
745 "reference_orbital_frequency",
746 "reference_position1",
747 "reference_position2",
748 ]
750 # Check if interpolation is necessary
751 # if t_ref is supplied
753 interp = False
755 if not t_ref:
756 if not f_ref:
757 # No interpolation required. Use available reference values.
758 interp = False
760 else:
761 # Try to get the reference time from orbital frequency
762 try:
763 t_ref = get_ref_time_from_ref_freq(h5_file, f_ref)
765 # Check if interpolation is required
766 interp, avail_ref_time = check_interp_req(h5_file, t_ref)
767 except Exception as excep:
768 print(
769 f"Could not obtain reference time from given reference frequency {f_ref}.",
770 excep,
771 )
772 print("Choosing available reference time")
773 interp = False
774 else:
775 interp, avail_ref_time = check_interp_req(h5_file, t_ref)
777 if interp is False:
778 # Then load default values from the NR data
779 # at hard coded reference time.
781 # Get the available reference time for book keeping
782 t_ref = get_ref_time_from_metadata(sim_metadata)
784 # Check for LAL frame in simulation metadata.
785 # RIT and GT qualify this.
786 # Default attributes in case of no interpolation
787 ref_check_def_h5, absent_attrs_h5 = check_nr_attrs(h5_file)
789 if ref_check_def_h5 is False:
790 # Then the LAL source frame information is not present in the H5 file.
791 # Then this could be SXS or GT data. The LAL source frame need to be computed from
792 # the H5 File or simulation metadata.
794 # Check if LAL source frame info is present in the simulation metadata.
795 ref_check_def_meta, absent_attrs_meta = check_nr_attrs(sim_metadata)
797 if ref_check_def_meta is False:
798 # Then this is SXS data.
800 # Check for raw information in metadata to compute the LAL source frame.
801 ref_check_def_meta_sxs, absent_attrs_meta_sxs = check_nr_attrs(
802 sim_metadata, req_def_attrs_sxs
803 )
805 if ref_check_def_meta_sxs is True:
806 # Compute the LAL source frame from simulation metadata
807 ref_params = compute_lal_source_frame_from_sxs_metadata(
808 sim_metadata
809 )
810 else:
811 raise Exception(
812 "Insufficient information to compute the LAL source frame."
813 f"\n Missing information is {absent_attrs_meta_sxs}."
814 )
815 else:
816 # LAL source frame is present in the simulation metadata
817 ref_params = get_ref_vals(sim_metadata)
818 else:
819 ref_params = get_ref_vals(h5_file)
821 elif interp is True:
822 # Experimental; This assumes all the required atributes needed
823 # to compute the LAL source frame at the given reference time
824 # are present in the H5file only.
826 # Attributes required for interpolation.
827 req_ts_attrs = [
828 "LNhatx-vs-time",
829 "LNhaty-vs-time",
830 "LNhatz-vs-time",
831 "position1x-vs-time",
832 "position1y-vs-time",
833 "position1z-vs-time",
834 "position2x-vs-time",
835 "position2y-vs-time",
836 "position2z-vs-time",
837 ]
839 # Check if time series data of required reference data is present
840 ref_check_interp_req, absent_interp_attrs = check_nr_attrs(
841 h5_file, req_ts_attrs
842 )
844 if ref_check_interp_req is False:
845 raise KeyError(
846 "Insufficient information to compute the LAL source frame at given reference time."
847 f"Missing information is {absent_interp_attrs}."
848 )
849 else:
850 ref_params = compute_lal_source_frame_by_interp(
851 h5_file, req_ts_attrs, t_ref
852 )
854 # Warning 1
855 # Implement this Warning
856 # XLAL_CHECK( ref_time!=XLAL_FAILURE, XLAL_FAILURE, "Error computing reference time.
857 # Try setting fRef equal to the f_low given by the NR simulation or to a value <=0 to deactivate
858 # fRef for a non-precessing simulation.\n")
860 # Get the LAL source frame vectors
861 ln_hat_x = ref_params["LNhatx"]
862 ln_hat_y = ref_params["LNhaty"]
863 ln_hat_z = ref_params["LNhatz"]
865 n_hat_x = ref_params["nhatx"]
866 n_hat_y = ref_params["nhaty"]
867 n_hat_z = ref_params["nhatz"]
869 ln_hat = np.array([ln_hat_x, ln_hat_y, ln_hat_z])
870 n_hat = np.array([n_hat_x, n_hat_y, n_hat_z])
872 # 2.3: Carryout vector math to get Zref in the lal wave frame
873 corb_phase = np.cos(orb_phase)
874 sorb_phase = np.sin(orb_phase)
875 sinclination = np.sin(inclination)
876 cinclination = np.cos(inclination)
878 ln_cross_n = np.cross(ln_hat, n_hat)
879 ln_cross_n_x, ln_cross_n_y, ln_cross_n_z = ln_cross_n
881 z_wave_x = sinclination * (sorb_phase * n_hat_x + corb_phase * ln_cross_n_x)
882 z_wave_y = sinclination * (sorb_phase * n_hat_y + corb_phase * ln_cross_n_y)
883 z_wave_z = sinclination * (sorb_phase * n_hat_z + corb_phase * ln_cross_n_z)
885 z_wave_x += cinclination * ln_hat_x
886 z_wave_y += cinclination * ln_hat_y
887 z_wave_z += cinclination * ln_hat_z
889 z_wave = np.array([z_wave_x, z_wave_y, z_wave_z])
890 z_wave = z_wave / np.linalg.norm(z_wave)
892 #################################################################
893 # Step 3.1: Extract theta and psi from Z in the lal wave frame
894 # NOTE: Theta can only run between 0 and pi, so no problem with arccos here
895 theta = np.arccos(z_wave_z)
897 # Degenerate if Z_wave[2] == 1. In this case just choose psi randomly,
898 # the choice will be cancelled out by alpha correction (I hope!)
900 # If theta is very close to the poles
901 # return a random value
902 if abs(z_wave_z - 1.0) < tol:
903 psi = 0.5
905 else:
906 # psi can run between 0 and 2pi, but only one solution works for x and y */
907 # Possible numerical issues if z_wave_x = sin(theta) */
908 if abs(z_wave_x / np.sin(theta)) > 1.0:
909 if abs(z_wave_x / np.sin(theta)) < (1 + 10 * tol):
910 # LAL tol retained.
911 if (z_wave_x * np.sin(theta)) < 0.0:
912 psi = np.pi
914 else:
915 psi = 0.0
917 else:
918 print(f"z_wave_x = {z_wave_x}")
919 print(f"sin(theta) = {np.sin(theta)}")
920 raise ValueError(
921 "Z_x cannot be bigger than sin(theta). Please contact the developers."
922 )
924 else:
925 psi = np.arccos(z_wave_x / np.sin(theta))
927 y_val = np.sin(psi) * np.sin(theta)
929 # If z_wave[1] is negative, flip psi so that sin(psi) goes negative
930 # while preserving cos(psi) */
931 if z_wave_y < 0.0:
932 psi = 2 * np.pi - psi
933 y_val = np.sin(psi) * np.sin(theta)
935 if abs(y_val - z_wave_y) > (5e3 * tol):
936 # LAL tol retained.
937 print(f"orb_phase = {orb_phase}")
938 print(
939 f"y_val = {y_val}, z_wave_y = {z_wave_y}, abs(y_val - z_wave_y) = {abs(y_val - z_wave_y)}"
940 )
941 raise ValueError("Math consistency failure! Please contact the developers.")
943 # 3.2: Compute the vectors theta_hat and psi_hat
944 # stheta = np.sin(theta)
945 # ctheta = np.cos(theta)
947 spsi = np.sin(psi)
948 cpsi = np.cos(psi)
950 # theta_hat_x = cpsi * ctheta
951 # theta_hat_y = spsi * ctheta
952 # theta_hat_z = -stheta
953 # theta_hat = np.array([theta_hat_x, theta_hat_y, theta_hat_z])
955 psi_hat_x = -spsi
956 psi_hat_y = cpsi
957 psi_hat_z = 0.0
958 psi_hat = np.array([psi_hat_x, psi_hat_y, psi_hat_z])
960 # Step 4: Compute sin(alpha) and cos(alpha)
961 # Rotation angles on the tangent plane
962 # due to spin weight.
964 # n_dot_theta = np.dot(n_hat, theta_hat)
965 # ln_cross_n_dot_theta = np.dot(ln_cross_n, theta_hat)
967 n_dot_psi = np.dot(n_hat, psi_hat)
968 ln_cross_n_dot_psi = np.dot(ln_cross_n, psi_hat)
970 # salpha = corb_phase * n_dot_theta - sorb_phase * ln_cross_n_dot_theta
971 calpha = corb_phase * n_dot_psi - sorb_phase * ln_cross_n_dot_psi
973 if abs(calpha) > 1:
974 calpha_err = abs(calpha) - 1
975 if calpha_err < tol:
976 # This tol could have been much smaller.
977 # Just resuing the default for now.
978 print(
979 f"Correcting the polarization angle for finite precision error {calpha_err}"
980 )
981 calpha = calpha / abs(calpha)
982 else:
983 raise ValueError(
984 "Seems like something is wrong with the polarization angle. Please contact the developers!"
985 )
987 alpha = np.arccos(calpha)
989 angles = {
990 "theta": theta,
991 "psi": psi,
992 "alpha": alpha,
993 "t_ref": t_ref,
994 "f_ref": f_ref,
995 }
997 return angles