diff --git a/gwfast/signal.py b/gwfast/signal.py index 4314f51..4d1131c 100644 --- a/gwfast/signal.py +++ b/gwfast/signal.py @@ -7,22 +7,20 @@ import os import jax - #Enable 64bit on JAX, fundamental #from jax.config import config #config.update("jax_enable_x64", True) jax.config.update("jax_enable_x64", True) #config.update("TF_CPP_MIN_LOG_LEVEL", 0) -os.environ['XLA_PYTHON_CLIENT_PREALLOCATE']='false' -os.environ['XLA_PYTHON_CLIENT_ALLOCATOR']='platform' +os.environ['XLA_PYTHON_CLIENT_PREALLOCATE'] = 'false' +os.environ['XLA_PYTHON_CLIENT_ALLOCATOR'] = 'platform' # We use both the original numpy, denoted as onp, and the JAX implementation of numpy, denoted as np import numpy as onp import jax.numpy as np -from jax.interpreters import xla -from jax import pmap, vmap, jacrev, jit #jacfwd +from jax import pmap, vmap, jacrev, jit import time import h5py import numdifftools as ndt @@ -33,13 +31,12 @@ from gwfast import gwfastGlobals as glob - class GWSignal(object): """ Class to compute the GW signal emitted by a coalescing binary system as seen by a detector on Earth. - + The functions defined within this class allow to get e.g. the amplitude of the signal, its phase, SNR and Fisher matrix elements. - + :param WaveFormModel wf_model: Object containing the waveform model. :param str psd_path: Full path to the file containing the detector's *Power Spectral Density*, PSD, or *Amplitude Spectral Density*, ASD, including the file extension. The file is assumed to have two columns, the first containing the frequencies (in :math:`\\rm Hz`) and the second containing the detector's PSD/ASD at each frequency. :param str detector_shape: The shape of the detector, to be chosen among ``'L'`` for an L-shaped detector (90°-arms) and ``'T'`` for a triangular detector (3 nested detectors with 60°-arms). @@ -56,21 +53,21 @@ class GWSignal(object): :param float DutyFactor: Duty factor of the detector, between 0 and 1, representing the percentage of time the detector (each detector independently in the case of a triangular detector) is supposed to be operational. :param bool, optional compute2arms: Boolean specifying if, in the case of a triangular detector, the computation can be performed only in two of the instruments, using the null-stream to get the signal in the third instrument, speeding up the computation by 1/3. :param bool, optional jitCompileDerivs: Boolean specifying if the derivatives function has to be jit compiled. - + """ ''' Inputs are an object containing the waveform model, the coordinates of the detector (latitude and longitude in deg), its shape (L or T), the angle with respect to East of the bisector of the arms (deg) and its ASD or PSD (given in a .txt file containing two columns: one with the frequencies and one with the ASD or PSD values, remember ASD=sqrt(PSD)) - + ''' - def __init__(self, wf_model, + def __init__(self, wf_model, psd_path=None, detector_shape = 'T', det_lat=40.44, det_long=9.45, - det_xax=0., + det_xax=0., verbose=True, is_ASD=True, useEarthMotion = False, @@ -85,54 +82,54 @@ def __init__(self, wf_model, """ if (detector_shape!='L') and (detector_shape!='T'): raise ValueError('Enter valid detector configuration') - + if psd_path is None: raise ValueError('Enter a valid PSD or ASD path') - + if verbose: if not is_ASD: print('Using PSD from file %s ' %psd_path) else: print('Using ASD from file %s ' %psd_path) - + if (useEarthMotion) and (wf_model.objType == 'BBH') and (verbose): print('WARNING: the motion of Earth gives a negligible contribution for BBH signals, consider switching it off to make the code run faster') if (not useEarthMotion) and (wf_model.objType == 'BNS') and (verbose): print('WARNING: the motion of Earth gives a relevant contribution for BNS signals, consider switching it on') if (not useEarthMotion) and (wf_model.objType == 'NSBH') and (verbose): print('WARNING: the motion of Earth gives a relevant contribution for NSBH signals, consider switching it on') - + self.wf_model = wf_model - + self.psd_base_path = ('/').join(psd_path.split('/')[:-1]) self.psd_file_name = psd_path.split('/')[-1] - + self.verbose = verbose self.detector_shape = detector_shape - + self.det_lat_rad = det_lat*np.pi/180. self.det_long_rad = det_long*np.pi/180. - + self.det_xax_rad = det_xax*np.pi/180. - + self.IntTablePath = IntTablePath #This is the percentage of time each arm of the detector (or the whole detector for an L) is supposed to be operational, between 0 and 1, default is None, resulting in a detector always online self.DutyFactor = DutyFactor - + noise = onp.loadtxt(psd_path, usecols=(0,1)) f = noise[:,0] if is_ASD: S = (noise[:,1])**2 else: S = noise[:,1] - + self.strainFreq = f self.noiseCurve = S - + import scipy.integrate as igt mask = self.strainFreq >= fmin self.strainInteg = igt.cumtrapz(self.strainFreq[mask]**(-7./3.)/S[mask], self.strainFreq[mask], initial = 0) - + self.useEarthMotion = useEarthMotion self.noMotion = noMotion if self.noMotion and self.useEarthMotion: @@ -140,25 +137,25 @@ def __init__(self, wf_model, self.useEarthMotion = False self.fmin = fmin #Hz self.fmax = fmax #Hz or None - + if detector_shape == 'L': self.angbtwArms = 0.5*np.pi elif detector_shape == 'T': self.angbtwArms = np.pi/3. - + self.IntegInterpArr = None self.compute2arms = compute2arms - + onp.random.seed(None) self.seedUse = onp.random.randint(2**32 - 1, size=1) self.jitCompileDerivs = jitCompileDerivs - + if not self.wf_model.is_LAL: self._init_jax() else: self._SignalDerivatives_use = self._SignalDerivatives - - + + def _init_jax(self): """ JAX initialisation method @@ -172,12 +169,12 @@ def _init_jax(self): if self.verbose: print('Jax local device count: %s' %str(jax.local_device_count())) print('Jax device count: %s' %str(jax.device_count())) - + if self.jitCompileDerivs: self._SignalDerivatives_use = jit(self._SignalDerivatives, static_argnames=['use_chi1chi2', 'use_m1m2', 'computeAnalyticalDeriv', 'use_prec_ang', 'computeDerivFinDiff', 'stepNDT', 'methodNDT']) else: self._SignalDerivatives_use = self._SignalDerivatives - + inj_params_init = {'Mc': np.array([77.23905294]), 'Phicoal': np.array([3.28297867]), 'chi1z': np.array([0.2018924]), @@ -207,36 +204,36 @@ def _init_jax(self): self.detector_shape = 'L' # Get a faster Initialization with an L _ = self.SNRInteg(inj_params_init, res=10) _ = self.FisherMatr(inj_params_init, res=10) - + if self.verbose: print('Done.') self.verbose = verboseOr self.detector_shape = detector_shapeOr - + def _clear_cache(self): if self.jitCompileDerivs: print('Clearing cache...') self._SignalDerivatives_use = jax.jit(self._SignalDerivatives, static_argnums=(15,16,17,18,19)) - + def _update_seed(self, seed=None): """ Update the seed for the duty cycle with a random value or a user input value. - + :param int, optional seed: User input value for the seed. - + """ onp.random.seed(None) if seed is None: self.seedUse = onp.random.randint(2**32 - 1, size=1) else: self.seedUse = seed - + def _tabulateIntegrals(self, res=200, store=True, Mcmin=.9, Mcmax=9., etamin=.1): """ Compute the table of integrals to use :py:class:`GWSignal.SNRFastInsp`. - + .. deprecated:: 1.0.0 - + """ def IntegrandC(f, Mc, tcoal, n): t = tcoal - 2.18567 * ((1.21/Mc)**(5./3.)) * ((100/f[:,onp.newaxis])**(8./3.))/(3600.*24) @@ -244,18 +241,18 @@ def IntegrandC(f, Mc, tcoal, n): def IntegrandS(f, Mc, tcoal, n): t = tcoal - 2.18567 * ((1.21/Mc)**(5./3.)) * ((100/f[:,onp.newaxis])**(8./3.))/(3600.*24) return (f[:,onp.newaxis]**(-7./3.))*np.sin(n*2.*np.pi*t) - + Mcgrid = onp.linspace(Mcmin, Mcmax, res) etagrid = onp.linspace(etamin, 0.25, res) tcgrid = onp.linspace(0.,2.*np.pi,res) - + Igrid = onp.zeros((res,res,res,9)) - + if self.verbose: print('Computing table of integrals...\n') - + in_time=time.time() - + for i,Mc in enumerate(Mcgrid): for j,eta in enumerate(etagrid): tmpev = {'Mc':np.array([Mc, ]), 'eta':np.array([eta])} @@ -269,18 +266,18 @@ def IntegrandS(f, Mc, tcoal, n): tmpIntegrandS = IntegrandS(self.strainFreq[mask], Mc, tcgrid, m+1.) Igrid[i,j,:,m] = onp.trapz(tmpIntegrandC/noisegrids.T, fgrids.T, axis=0) Igrid[i,j,:,m+4] = onp.trapz(tmpIntegrandS/noisegrids.T, fgrids.T, axis=0) - + tmpIntegrand = IntegrandC(self.strainFreq[mask], Mc, tcgrid, 0.) Igrid[i,j,:,8] = onp.trapz(tmpIntegrand/noisegrids.T, fgrids.T, axis=0) - + if self.verbose: print('Done in %.2fs \n' %(time.time() - in_time)) - + if store: print('Saving result...') if not os.path.isdir(os.path.join(self.psd_base_path, 'Integral_Tables')): os.mkdir(os.path.join(self.psd_base_path, 'Integral_Tables')) - + with h5py.File(os.path.join(self.psd_base_path, 'Integral_Tables', type(self.wf_model).__name__+str(res)+'.h5'), 'w') as out: out.create_dataset('Mc', data=Mcgrid, compression='gzip', shuffle=True) out.create_dataset('eta', data=etagrid, compression='gzip', shuffle=True) @@ -290,15 +287,15 @@ def IntegrandS(f, Mc, tcoal, n): out.attrs['etamin'] = etamin out.attrs['Mcmin'] = Mcmin out.attrs['Mcmax'] = Mcmax - + return Igrid, Mcgrid, etagrid, tcgrid - + def _make_SNRig_interpolator(self, ): """ Make interpolator of the table of integrals to use :py:class:`GWSignal.SNRFastInsp`. - + .. deprecated:: 1.0.0 - + """ from scipy.interpolate import RegularGridInterpolator if self.IntTablePath is not None: @@ -317,34 +314,34 @@ def _make_SNRig_interpolator(self, ): if self.verbose: print('Attributes of pre-computed integrals: ') print([(k, inp.attrs[k]) for k in inp.attrs.keys()]) - + else: print('Tabulating integrals...') Igrid, Mcs, etas, tcs = self._tabulateIntegrals() - + else: print('Tabulating integrals...') Igrid, Mcs, etas, tcs = self._tabulateIntegrals() - + self.IntegInterpArr = onp.array([]) for i in range(9): # The interpolator contains in the elements from i=0 to 3 the integrals of cos((i+1) Om t)f^{-7/3} # in the elements from 4 to 7 the integrals of the sine (from lower to higher i), and in element 8 the # integral of f^-7/3 alone - + self.IntegInterpArr = onp.append(self.IntegInterpArr,RegularGridInterpolator((Mcs, etas, tcs), Igrid[:,:,:,i])) - - + + def _ra_dec_from_th_phi(self, theta, phi): return utils.ra_dec_from_th_phi_rad(theta, phi) - - + + def _PatternFunction(self, theta, phi, t, psi, rot=0.): """ Compute the value of the so-called pattern functions of the detector for a set of sky coordinates, GW polarisation(s) and time(s). - + For the definition of the pattern functions see `arXiv:gr-qc/9804014 `_ eq. (10)--(13). - + :param array or float theta: The :math:`\\theta` sky position angle(s), in :math:`\\rm rad`. :param array or float phi: The :math:`\phi` sky position angle(s), in :math:`\\rm rad`. :param array or float t: The time(s) given as GMST. @@ -352,11 +349,11 @@ def _PatternFunction(self, theta, phi, t, psi, rot=0.): :param float rot: Further rotation of the interferometer with respect to the :py:data:`self.xax` orientation, in degrees, needed for the triangular geometry. In this case, the three arms will have orientations 1 --> :py:data:`self.xax`, 2 --> :py:data:`self.xax` + 60°, 3 --> :py:data:`self.xax` + 120°. :return: Plus and cross pattern functions of the detector evaluated at the given parameters. :rtype: tuple(array, array) or tuple(float, float) - + """ # See P. Jaranowski, A. Krolak, B. F. Schutz, PRD 58, 063001, eq. (10)--(13) - - + + def afun(ra, dec, t, rot): phir = self.det_long_rad a1 = 0.0625*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*(3.-np.cos(2.*dec))*np.cos(2.*(ra - phir - 2.*np.pi*t)) @@ -365,82 +362,82 @@ def afun(ra, dec, t, rot): a4 = 0.5*np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.sin(2.*dec)*np.sin(ra - phir - 2.*np.pi*t) a5 = 3.*0.25*np.sin(2*(self.det_xax_rad+rot))*(np.cos(self.det_lat_rad)*np.cos(dec))**2. return a1 - a2 + a3 - a4 + a5 - + def bfun(ra, dec, t, rot): phir = self.det_long_rad b1 = np.cos(2*(self.det_xax_rad+rot))*np.sin(self.det_lat_rad)*np.sin(dec)*np.cos(2.*(ra - phir - 2.*np.pi*t)) b2 = 0.25*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*np.sin(dec)*np.sin(2.*(ra - phir - 2.*np.pi*t)) b3 = np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.cos(dec)*np.cos(ra - phir - 2.*np.pi*t) b4 = 0.5*np.sin(2*(self.det_xax_rad+rot))*np.sin(2.*self.det_lat_rad)*np.cos(dec)*np.sin(ra - phir - 2.*np.pi*t) - + return b1 + b2 + b3 + b4 - + rot_rad = rot*np.pi/180. - + ras, decs = self._ra_dec_from_th_phi(theta, phi) afac = afun(ras, decs, t, rot_rad) bfac = bfun(ras, decs, t, rot_rad) - + Fp = np.sin(self.angbtwArms)*(afac*np.cos(2.*psi) + bfac*np.sin(2*psi)) Fc = np.sin(self.angbtwArms)*(bfac*np.cos(2.*psi) - afac*np.sin(2*psi)) - + return Fp, Fc - + def _phiPhase(self, theta, phi, t, iota, psi, Fp=None,Fc=None): #The polarization phase contribution (the change in F+ and Fx with time influences also the phase) - + if (Fp is None) or (Fc is None): Fp, Fc = self._PatternFunction(theta, phi, t, psi) - + phiP = -np.arctan2(np.cos(iota)*Fc,0.5*(1.+((np.cos(iota))**2))*Fp) - + #The contriution to the amplitude is negligible, so we do not compute it - + return phiP - + def _DeltLoc(self, theta, phi, t): """ Compute the time needed to go from Earth center to detector location for a set of sky coordinates and time(s). The result is given in seconds. - + :param array or float theta: The :math:`\\theta` sky position angle(s), in :math:`\\rm rad`. :param array or float phi: The :math:`\phi` sky position angle(s), in :math:`\\rm rad`. :param array or float t: The time(s) given as GMST. - + :return: Time shift(s) to go from Earth center to detector location. :rtype: array or float - + """ # Time needed to go from Earth center to detector location - + ras, decs = self._ra_dec_from_th_phi(theta, phi) - + comp1 = np.cos(decs)*np.cos(ras)*np.cos(self.det_lat_rad)*np.cos(self.det_long_rad + 2.*np.pi*t) comp2 = np.cos(decs)*np.sin(ras)*np.cos(self.det_lat_rad)*np.sin(self.det_long_rad + 2.*np.pi*t) comp3 = np.sin(decs)*np.sin(self.det_lat_rad) # The minus sign arises from the definition of the unit vector pointing to the source Delt = - glob.REarth*(comp1+comp2+comp3)/glob.clight - + return Delt # in seconds - + def GWAmplitudes(self, evParams, f, rot=0.): """ Compute the amplitude of the signal(s) as seen by the detector, as a function of the parameters, at given frequencies. - + :param dict(array, array, ...) evParams: Dictionary containing the parameters of the event(s), as in :py:data:`events`. :param array or float f: The frequency(ies) at which to perform the calculation, in :math:`\\rm Hz`. :param float rot: Further rotation of the interferometer with respect to the :py:data:`self.xax` orientation, in degrees, needed for the triangular geometry. :return: Plus and cross amplitudes at the detector, evaluated at the given parameters and frequency(ies). :rtype: tuple(array, array) or tuple(float, float) - + """ - # evParams are all the parameters characterizing the event(s) under exam. It has to be a dictionary containing the entries: + # evParams are all the parameters characterizing the event(s) under exam. It has to be a dictionary containing the entries: # Mc -> chirp mass (Msun), dL -> luminosity distance (Gpc), theta & phi -> sky position (rad), iota -> inclination angle of orbital angular momentum to l.o.s toward the detector, # psi -> polarisation angle, tcoal -> time of coalescence as GMST (fraction of days), eta -> symmetric mass ratio, Phicoal -> GW frequency at coalescence. # chi1z, chi2z -> dimensionless spin components aligned to orbital angular momentum [-1;1], Lambda1,2 -> tidal parameters of the objects, # f is the frequency (Hz) - + theta, phi, iota, psi, tcoal = evParams['theta'], evParams['phi'], evParams['iota'], evParams['psi'], evParams['tcoal'] - + if self.noMotion: t = 0. t = t + self._DeltLoc(theta, phi, t)/(3600.*24.) @@ -453,7 +450,7 @@ def GWAmplitudes(self, evParams, f, rot=0.): t = t + self._DeltLoc(theta, phi, t)/(3600.*24.) # wfAmpl = self.wf_model.Ampl(f, **evParams) Fp, Fc = self._PatternFunction(theta, phi, t, psi, rot=rot) - + if (self.wf_model.is_HigherModes) or (self.wf_model.is_Precessing): # If the waveform includes higher modes or precessing spins, it is not possible to compute amplitude and phase separately, make all together hp, hc = self.wf_model.hphc(f, **evParams) @@ -462,19 +459,19 @@ def GWAmplitudes(self, evParams, f, rot=0.): wfAmpl = self.wf_model.Ampl(f, **evParams) Ap = wfAmpl*Fp*0.5*(1.+(np.cos(iota))**2) Ac = wfAmpl*Fc*np.cos(iota) - + return Ap, Ac - + def GWPhase(self, evParams, f): """ Compute the complete phase of the signal(s), as a function of the parameters, at given frequencies. - + :param dict(array, array, ...) evParams: Dictionary containing the parameters of the event(s), as in :py:data:`events`. :param array or float f: The frequency(ies) at which to perform the calculation, in :math:`\\rm Hz`. - + :return: Complete signal phase, evaluated at the given parameters and frequency(ies). :rtype: array or float - + """ # Phase of the GW signal tcoal, Phicoal = evParams['tcoal'], evParams['Phicoal'] @@ -485,7 +482,7 @@ def GWPhase(self, evParams, f): def GWstrain(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=0., is_m1m2=False, is_chi1chi2=False, is_prec_ang=False, is_Lam1Lam2=False, return_single_comp=None): """ Compute the full GW strain (complex) as a function of the parameters, at given frequencies. - + :param array or float f: The frequency(ies) at which to perform the calculation, in :math:`\\rm Hz`. :param array or float Mc: The chirp mass(es), :math:`{\cal M}_c`, in units of :math:`\\rm M_{\odot}`. If ``is_m1m2=True`` this is interpreted as the primary mass, :math:`m_1`, in units of :math:`\\rm M_{\odot}`. :param array or float eta: The symmetric mass ratio(s), :math:`\eta`. If ``is_m1m2=True`` this is interpreted as the secondary mass, :math:`m_2`, in units of :math:`\\rm M_{\odot}`. @@ -513,19 +510,19 @@ def GWstrain(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, :param str return_single_comp: String specifying if a single component of the signal should be returned, to be chosen among ``Ap`` and ``Ac``, to return the plus and cross amplitude, :math:`A_+` and :math:`A_{\\times}`, respectively, and ``Psip`` and ``Psic``, to return the plus and cross phase, :math:`\Phi_+` and :math:`\Phi_{\\times}`, respectively. :return: Complete signal strain (complex), evaluated at the given parameters and frequency(ies). :rtype: array or float - + """ # Full GW strain expression (complex) # Here we have the decompressed parameters and we put them back in a dictionary just to have an easier # implementation of the JAX module for derivatives - + if is_m1m2: # Interpret Mc as m1 and eta as m2 McUse, etaUse = utils.Mceta_from_m1m2(Mc, eta) else: McUse = Mc etaUse = eta - + if not self.wf_model.is_Precessing: if is_chi1chi2: # Interpret chiS as chi1z and chiA as chi2z @@ -546,21 +543,21 @@ def GWstrain(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, else: # convert angles and iota iota, chi1xUse, chi1yUse, chi1z, chi2xUse, chi2yUse, chi2z = utils.TransformPrecessing_angles2comp(thetaJN=iota, phiJL=chi1y, theta1=chi1x, theta2=chi2x, phi12=chi2y, chi1=chiS, chi2=chiA, Mc=McUse, eta=etaUse, fRef=self.fmin, phiRef=0.) - + evParams = {'Mc':McUse, 'dL':dL, 'theta':theta, 'phi':phi, 'iota':iota, 'psi':psi, 'tcoal':tcoal, 'eta':etaUse, 'Phicoal':Phicoal, 'chi1z':chi1z, 'chi2z':chi2z, 'chi1x':chi1xUse, 'chi2x':chi2xUse, 'chi1y':chi1yUse, 'chi2y':chi2yUse} - + if self.wf_model.is_tidal: if not is_Lam1Lam2: Lambda1, Lambda2 = utils.Lam12_from_Lamt_delLam(LambdaTilde, deltaLambda, etaUse) else: Lambda1, Lambda2 = LambdaTilde, deltaLambda - + evParams['Lambda1'] = Lambda1 evParams['Lambda2'] = Lambda2 - + if self.wf_model.is_eccentric: evParams['ecc'] = ecc - + if self.useEarthMotion: # Compute Doppler contribution t = tcoal - self.wf_model.tau_star(f, **evParams)/(3600.*24.) @@ -578,109 +575,91 @@ def GWstrain(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, t = tcoal tmpDeltLoc = self._DeltLoc(theta, phi, t) # in seconds t = t + tmpDeltLoc/(3600.*24.) - + phiL = (2.*np.pi*f)*tmpDeltLoc - if (self.wf_model.is_HigherModes) or (self.wf_model.is_Precessing): - # If the waveform includes higher modes or precessing spins, it is not possible to compute amplitude and phase separately, make all together - Fp, Fc = self._PatternFunction(theta, phi, t, psi, rot=rot) - hp, hc = self.wf_model.hphc(f, **evParams) - hp = hp*Fp*np.exp(1j*(phiD + phiL + 2.*np.pi*f*(tcoal*3600.*24.) - Phicoal)) - hc = hc*Fc*np.exp(1j*(phiD + phiL + 2.*np.pi*f*(tcoal*3600.*24.) - Phicoal)) - + needs_HM = (self.wf_model.is_HigherModes) or (self.wf_model.is_Precessing) + is_LAL = self.wf_model.is_LAL + + if not (needs_HM or is_LAL): + # Return the simplest waveform construction + Ap, Ac = self.GWAmplitudes(evParams, f, rot=rot) + Psi = self.GWPhase(evParams, f) + Psi = Psi + phiD + phiL + if return_single_comp is not None: if (return_single_comp == 'Ap'): - return np.abs(hp) + return Ap elif (return_single_comp == 'Ac'): - return np.abs(hc) + return Ac elif (return_single_comp == 'Psip'): - return np.unwrap(np.angle(hp)) + return Psi # np.unwrap(Psi) elif (return_single_comp == 'Psic'): - return np.unwrap(np.angle(hc)) + return Psi + np.pi*0.5 # np.unwrap(Psi + np.pi*0.5) elif (return_single_comp == 'At'): - return np.abs(hp + hc) + return np.abs(Ap + 1j*Ac) elif (return_single_comp == 'Psit'): - return np.unwrap(np.angle(hp + hc)) + return Psi + np.arctan2(np.real(Ac), np.real(Ap)) else: raise ValueError('Single component to return has to be among Ap, Ac, Psip, Psic') else: - return hp + hc - else: - if self.wf_model.is_LAL: - # If the waveform comes from LAL, and does not include HM or precessing spins, it is pointless to perform twice the computation just to add the cos(iota) factors. We thus evaluate hphc once and add them here - Fp, Fc = self._PatternFunction(theta, phi, t, psi, rot=rot) - hp, hc = self.wf_model.hphc(f, **evParams) - hp = hp*Fp*np.exp(1j*(phiD + phiL + 2.*np.pi*f*(tcoal*3600.*24.) - Phicoal))*0.5*(1.+(np.cos(iota))**2) - hc = hc*Fc*np.exp(1j*(phiD + phiL + 2.*np.pi*f*(tcoal*3600.*24.) - Phicoal))*np.cos(iota) - - if return_single_comp is not None: - if (return_single_comp == 'Ap'): - return np.abs(hp) - elif (return_single_comp == 'Ac'): - return np.abs(hc) - elif (return_single_comp == 'Psip'): - return np.unwrap(np.angle(hp)) - elif (return_single_comp == 'Psic'): - return np.unwrap(np.angle(hc)) - elif (return_single_comp == 'At'): - return np.abs(hp + hc) - elif (return_single_comp == 'Psit'): - return np.unwrap(np.angle(hp + hc)) - else: - raise ValueError('Single component to return has to be among Ap, Ac, Psip, Psic') - else: - return hp + hc + return (Ap + 1j*Ac)*np.exp(Psi*1j) + #return np.sqrt(Ap*Ap + Ac*Ac)*np.exp((Psi+phiP)*1j) + + # If the waveform includes higher modes or precessing spins, it is not possible to compute amplitude and phase separately, make all together + Fp, Fc = self._PatternFunction(theta, phi, t, psi, rot=rot) + hp, hc = self.wf_model.hphc(f, **evParams) + hp *= Fp*np.exp(1j*(phiD + phiL + 2.*np.pi*f*(tcoal*3600.*24.) - Phicoal)) + hc *= Fc*np.exp(1j*(phiD + phiL + 2.*np.pi*f*(tcoal*3600.*24.) - Phicoal)) + + if is_LAL: + hp *= 0.5*(1.+(np.cos(iota))**2) + hc *= np.cos(iota) + + if return_single_comp is not None: + if (return_single_comp == 'Ap'): + return np.abs(hp) + elif (return_single_comp == 'Ac'): + return np.abs(hc) + elif (return_single_comp == 'Psip'): + return np.unwrap(np.angle(hp)) + elif (return_single_comp == 'Psic'): + return np.unwrap(np.angle(hc)) + elif (return_single_comp == 'At'): + return np.abs(hp + hc) + elif (return_single_comp == 'Psit'): + return np.unwrap(np.angle(hp + hc)) else: - Ap, Ac = self.GWAmplitudes(evParams, f, rot=rot) - Psi = self.GWPhase(evParams, f) - Psi = Psi + phiD + phiL - - if return_single_comp is not None: - if (return_single_comp == 'Ap'): - return Ap - elif (return_single_comp == 'Ac'): - return Ac - elif (return_single_comp == 'Psip'): - return Psi #np.unwrap(Psi) - elif (return_single_comp == 'Psic'): - return Psi + np.pi*0.5 #np.unwrap(Psi + np.pi*0.5) - elif (return_single_comp == 'At'): - return np.abs(Ap + 1j*Ac) - elif (return_single_comp == 'Psit'): - return Psi + np.arctan2(np.real(Ac),np.real(Ap)) - else: - raise ValueError('Single component to return has to be among Ap, Ac, Psip, Psic') - else: - return (Ap + 1j*Ac)*np.exp(Psi*1j) - #return np.sqrt(Ap*Ap + Ac*Ac)*np.exp((Psi+phiP)*1j) - - + raise ValueError('Single component to return has to be among Ap, Ac, Psip, Psic') + else: + return hp + hc + def SNRInteg(self, evParams, res=1000, return_all=False): """ Compute the *signal-to-noise-ratio*, SNR, as a function of the parameters of the event(s). - + :param dict(array, array, ...) evParams: Dictionary containing the parameters of the event(s), as in :py:data:`events`. :param int res: The resolution of the frequency grid to use. :param bool, optional return_all: Boolean specifying if, in the case of a triangular detector, the SNRs of the individual instruments have to be returned separately. In this case the return type is *list(array, array, array)*. - + :return: SNR(s) as a function of the parameters of the event(s). The shape is :math:`(N_{\\rm events})`. :rtype: 1-D array - + """ # SNR calculation performing the frequency integral for each signal # This is computationally more expensive, but needed for complex waveform models if self.DutyFactor is not None: onp.random.seed(self.seedUse) - + utils.check_evparams(evParams) - + #if not np.isscalar(evParams['Mc']): # SNR = np.zeros(len(np.asarray(evParams['Mc']))) #else: # SNR = 0. - + allSNRsq=[] - + if self.wf_model.is_Precessing: try: _ =evParams['chi1x'] @@ -702,7 +681,7 @@ def SNRInteg(self, evParams, res=1000, return_all=False): evParams['chi2z'] = evParams['chiS'] - evParams['chiA'] except KeyError: raise ValueError('Two among chi1z, chi2z and chiS, chiA have to be provided.') - + if self.wf_model.is_tidal: try: _=evParams['Lambda1'] @@ -711,18 +690,18 @@ def SNRInteg(self, evParams, res=1000, return_all=False): evParams['Lambda1'], evParams['Lambda2'] = utils.Lam12_from_Lamt_delLam(evParams['LambdaTilde'], evParams['deltaLambda'], evParams['eta']) except KeyError: raise ValueError('Two among Lambda1, Lambda2 and LambdaTilde and deltaLambda have to be provided.') - + fcut = self.wf_model.fcut(**evParams) - + if self.fmax is not None: fcut = np.where(fcut > self.fmax, self.fmax, fcut) - + fminarr = np.full(fcut.shape, self.fmin) fgrids = np.geomspace(fminarr,fcut,num=int(res)) # Out of the provided PSD range, we use a constant value of 1, which results in completely negligible contributions strainGrids = np.interp(fgrids, self.strainFreq, self.noiseCurve, left=1., right=1.) - - if self.detector_shape=='L': + + if self.detector_shape=='L': Aps, Acs = self.GWAmplitudes(evParams, fgrids) Atot = Aps*Aps + Acs*Acs SNRsq = np.trapezoid(Atot/strainGrids, fgrids, axis=0) @@ -765,7 +744,7 @@ def SNRInteg(self, evParams, res=1000, return_all=False): allSNRsq.append(tmpSNRsq3) #SNR = np.sqrt(tmpSNRsq1 + tmpSNRsq2 + tmpSNRsq3) allSNRsq = np.array(allSNRsq) - + if return_all: if self.detector_shape=='T': return 2*np.sqrt(allSNRsq) @@ -775,18 +754,17 @@ def SNRInteg(self, evParams, res=1000, return_all=False): return 2*np.sqrt(allSNRsq.sum(axis=0)) else: return np.squeeze(2*np.sqrt(allSNRsq), axis=0) - + # The factor of two arises by cutting the integral from 0 to infinity - - - def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', + + def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', use_m1m2=False, use_chi1chi2=True, use_prec_ang=True, computeDerivFinDiff=False, computeAnalyticalDeriv=True, return_all=False, **kwargs): """ Compute the *Fisher information matrix*, FIM, as a function of the parameters of the event(s). - + :param dict(array, array, ...) evParams: Dictionary containing the parameters of the event(s), as in :py:data:`events`. :param int res: The resolution of the frequency grid to use. :param float df: The spacing of the frequency grid to use, in :math:`\\rm Hz`. Alternative to ``res``. @@ -800,19 +778,19 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', :param kwargs: Optional arguments to be passed to :py:class:`gwfast.signal.GWSignal._SignalDerivatives`, such as ``methodNDT``. :return: FIM(s) as a function of the parameters of the event(s). The shape is :math:`(N_{\\rm parameters}`, :math:`N_{\\rm parameters}`, :math:`N_{\\rm events})`. :rtype: 3-D array - + """ # If use_m1m2=True the Fisher is computed w.r.t. m1 and m2, not Mc and eta # If use_chi1chi2=True the Fisher is computed w.r.t. chi1z and chi2z, not chiS and chiA if self.DutyFactor is not None: onp.random.seed(self.seedUse) - + utils.check_evparams(evParams) - - + + McOr, dL, theta, phi = evParams['Mc'].astype('complex128'), evParams['dL'].astype('complex128'), evParams['theta'].astype('complex128'), evParams['phi'].astype('complex128') iota, psi, tcoal, etaOr, Phicoal = evParams['iota'].astype('complex128'), evParams['psi'].astype('complex128'), evParams['tcoal'].astype('complex128'), evParams['eta'].astype('complex128'), evParams['Phicoal'].astype('complex128') - + if use_m1m2: # In this case Mc represents m1 and eta represents m2 Mc, eta = utils.m1m2_from_Mceta(McOr, etaOr) @@ -830,17 +808,17 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', evParams['chi2z'] = evParams['chiS'] - evParams['chiA'] except KeyError: raise ValueError('Two among chi1z, chi2z and chiS, chiA have to be provided.') - + chi1z, chi2z = evParams['chi1z'].astype('complex128'), evParams['chi2z'].astype('complex128') - + if use_chi1chi2: # In this case chiS represents chi1z and chiA represents chi2z chiS, chiA = chi1z, chi2z else: chiS, chiA = 0.5*(chi1z + chi2z), 0.5*(chi1z - chi2z) - + chi1x, chi2x, chi1y, chi2y = np.zeros(Mc.shape), np.zeros(Mc.shape), np.zeros(Mc.shape), np.zeros(Mc.shape) - + else: try: _=evParams['chi1x'] @@ -857,7 +835,7 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', else: # In this case iota=thetaJN, chi1y=phiJL, chi1x=tilt1, chi2x=tilt2, chi2y=phi12, chiS=chi1, chiA=chi2 iota, chi1y, chi1x, chi2x, chi2y, chiS, chiA = utils.TransformPrecessing_comp2angles(evParams['iota'].astype('complex128'), evParams['chi1x'].astype('complex128'), evParams['chi1y'].astype('complex128'), evParams['chi1z'].astype('complex128'), evParams['chi2x'].astype('complex128'), evParams['chi2y'].astype('complex128'), evParams['chi2z'].astype('complex128'), McOr, etaOr, fRef=self.fmin, phiRef=0.) - + if self.wf_model.is_tidal: try: Lambda1, Lambda2 = evParams['Lambda1'].astype('complex128'), evParams['Lambda2'].astype('complex128') @@ -867,10 +845,10 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', except KeyError: raise ValueError('Two among Lambda1, Lambda2 and LambdaTilde and deltaLambda have to be provided.') LambdaTilde, deltaLambda = utils.Lamt_delLam_from_Lam12(Lambda1, Lambda2, etaOr) - + else: Lambda1, Lambda2, LambdaTilde, deltaLambda = np.zeros(Mc.shape), np.zeros(Mc.shape), np.zeros(Mc.shape), np.zeros(Mc.shape) - + if self.wf_model.is_eccentric: try: ecc = evParams['ecc'].astype('complex128') @@ -878,12 +856,12 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', raise ValueError('Eccentricity has to be provided.') else: ecc = np.zeros(Mc.shape) - + fcut = self.wf_model.fcut(**evParams) - + if self.fmax is not None: fcut = np.where(fcut > self.fmax, self.fmax, fcut) - + fminarr = np.full(fcut.shape, self.fmin) if res is None and df is not None: res = np.floor( np.real((1+(fcut-fminarr)/df))) @@ -894,29 +872,29 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', fgrids = np.linspace(fminarr, fcut, num=int(res)) elif spacing=='geom': fgrids = np.geomspace(fminarr, fcut, num=int(res)) - + # Out of the provided PSD range, we use a constant value of 1, which results in completely negligible conntributions strainGrids = np.interp(fgrids, self.strainFreq, self.noiseCurve, left=1., right=1.) nParams = self.wf_model.nParams tcelem = self.wf_model.ParNums['tcoal'] - + if (self.wf_model.is_LAL) and (not computeDerivFinDiff): computeDerivFinDiff=True if self.verbose: print('Using LAL or TEOBResumS waveforms it is not possible to compute the derivatives using JAX automatic differentiation routines, being the functions written in C. Proceeding using numdifftools for numerical differentiation (finite differences)') - + allFishers=[] - - if self.detector_shape=='L': + + if self.detector_shape=='L': # Compute derivatives FisherDerivs = self._SignalDerivatives_use(fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=0., use_m1m2=use_m1m2, use_chi1chi2=use_chi1chi2, use_prec_ang=use_prec_ang, computeAnalyticalDeriv=computeAnalyticalDeriv, computeDerivFinDiff=computeDerivFinDiff, **kwargs) # Change the units of the tcoal derivative from days to seconds (this improves conditioning) FisherDerivs = onp.array(FisherDerivs) FisherDerivs[tcelem,:,:] /= (3600.*24.) - + FisherIntegrands = (onp.conjugate(FisherDerivs[:,:,onp.newaxis,:])*FisherDerivs.transpose(1,0,2)) - + Fisher = onp.zeros((nParams,nParams,len(Mc))) # This for is unavoidable for alpha in range(nParams): @@ -939,7 +917,7 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', FisherDerivs = onp.array(FisherDerivs) FisherDerivs[tcelem,:,:] /= (3600.*24.) FisherIntegrands = (onp.conjugate(FisherDerivs[:,:,onp.newaxis,:])*FisherDerivs.transpose(1,0,2)) - + tmpFisher = onp.zeros((nParams,nParams,len(Mc))) # This for is unavoidable if self.verbose: @@ -949,7 +927,7 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', for beta in range(alpha,nParams): tmpElem = FisherIntegrands[alpha,:,beta,:].T tmpFisher[alpha,beta, :] = onp.trapz(tmpElem.real/strainGrids.real, fgrids.real, axis=0)*4. - + tmpFisher[beta,alpha, :] = tmpFisher[alpha,beta, :] if self.DutyFactor is not None: excl = onp.random.choice([0,1],len(evParams['Mc']), p=[1.-self.DutyFactor,self.DutyFactor]) @@ -958,7 +936,7 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', #Fisher += tmpFisher else: # The signal in 3 arms sums to zero for geometrical reasons, so we can use this to skip some calculations - + # Compute derivatives FisherDerivs1 = self._SignalDerivatives_use(fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=0., use_m1m2=use_m1m2, use_chi1chi2=use_chi1chi2, use_prec_ang=use_prec_ang, computeAnalyticalDeriv=computeAnalyticalDeriv, computeDerivFinDiff=computeDerivFinDiff, **kwargs) # Change the units of the tcoal derivative from days to seconds (this improves conditioning) @@ -966,7 +944,7 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', FisherDerivs1[tcelem,:,:] /= (3600.*24.) FisherIntegrands = (onp.conjugate(FisherDerivs1[:,:,onp.newaxis,:])*FisherDerivs1.transpose(1,0,2)) - + tmpFisher = onp.zeros((nParams,nParams,len(Mc))) if self.verbose: print('Filling matrix for arm 1...') @@ -975,20 +953,20 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', for beta in range(alpha,nParams): tmpElem = FisherIntegrands[alpha,:,beta,:].T tmpFisher[alpha,beta, :] = onp.trapz(tmpElem.real/strainGrids.real, fgrids.real, axis=0)*4. - + tmpFisher[beta,alpha, :] = tmpFisher[alpha,beta, :] if self.DutyFactor is not None: excl = onp.random.choice([0,1],len(evParams['Mc']), p=[1.-self.DutyFactor,self.DutyFactor]) tmpFisher = tmpFisher*excl #Fisher += tmpFisher allFishers.append(tmpFisher) - - + + FisherDerivs2 = self._SignalDerivatives_use(fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=60., use_m1m2=use_m1m2, use_chi1chi2=use_chi1chi2, use_prec_ang=use_prec_ang, computeAnalyticalDeriv=computeAnalyticalDeriv, computeDerivFinDiff=computeDerivFinDiff, **kwargs) FisherDerivs2 = onp.array(FisherDerivs2) FisherDerivs2[tcelem,:,:] /= (3600.*24.) FisherIntegrands = (onp.conjugate(FisherDerivs2[:,:,onp.newaxis,:])*FisherDerivs2.transpose(1,0,2)) - + tmpFisher = onp.zeros((nParams,nParams,len(Mc))) # This for is unavoidable if self.verbose: @@ -997,18 +975,18 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', for beta in range(alpha,nParams): tmpElem = FisherIntegrands[alpha,:,beta,:].T tmpFisher[alpha,beta, :] = onp.trapz(tmpElem.real/strainGrids.real, fgrids.real, axis=0)*4. - + tmpFisher[beta,alpha, :] = tmpFisher[alpha,beta, :] if self.DutyFactor is not None: excl = onp.random.choice([0,1],len(evParams['Mc']), p=[1.-self.DutyFactor,self.DutyFactor]) tmpFisher = tmpFisher*excl #Fisher += tmpFisher allFishers.append(tmpFisher) - + FisherDerivs3 = - (FisherDerivs1 + FisherDerivs2) - + FisherIntegrands = (onp.conjugate(FisherDerivs3[:,:,onp.newaxis,:])*FisherDerivs3.transpose(1,0,2)) - + tmpFisher = onp.zeros((nParams,nParams,len(Mc))) # This for is unavoidable if self.verbose: @@ -1017,27 +995,25 @@ def FisherMatr(self, evParams, res=1000, df=None, spacing='geom', for beta in range(alpha,nParams): tmpElem = FisherIntegrands[alpha,:,beta,:].T tmpFisher[alpha,beta, :] = onp.trapz(tmpElem.real/strainGrids.real, fgrids.real, axis=0)*4. - + tmpFisher[beta,alpha, :] = tmpFisher[alpha,beta, :] if self.DutyFactor is not None: excl = onp.random.choice([0,1],len(evParams['Mc']), p=[1.-self.DutyFactor,self.DutyFactor]) tmpFisher = tmpFisher*excl #Fisher += tmpFisher allFishers.append(tmpFisher) - + if return_all: return allFishers elif self.detector_shape=='T': return onp.array(allFishers).sum(axis=0) else: return allFishers[0] - - - + def _SignalDerivatives(self, fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=0., use_m1m2=False, use_chi1chi2=True, use_prec_ang=True, computeDerivFinDiff=False, computeAnalyticalDeriv=True, stepNDT=MaxStepGenerator(base_step=1e-5), methodNDT='central', **kwargs): """ Compute the derivatives of the GW strain with respect to the parameters of the event(s) at given frequencies (in :math:`\\rm Hz`). - + :param array or float fgrids: The frequency(ies) at which to perform the calculation, in :math:`\\rm Hz`. :param array or float Mc: The chirp mass(es), :math:`{\cal M}_c`, in units of :math:`\\rm M_{\odot}`. If ``use_m1m2=True`` this is interpreted as the primary mass, :math:`m_1`, in units of :math:`\\rm M_{\odot}`. :param array or float eta: The symmetric mass ratio(s), :math:`\eta`. If ``use_m1m2=True`` this is interpreted as the secondary mass, :math:`m_2`, in units of :math:`\\rm M_{\odot}`. @@ -1068,17 +1044,17 @@ def _SignalDerivatives(self, fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, :param str methodNDT: The method to use in the computation with numerical differentiation (finite differences). This can be ``'central'``, ``'complex'``, ``'multicomplex'``, ``'forward'`` or ``'backward'``. :return: Complete signal strain derivatives (complex), evaluated at the given parameters and frequency(ies). :rtype: array - + """ # `numdifftools.step_generators `_ if self.verbose: print('Computing derivatives...') # Function to compute the derivatives of a GW signal, both with JAX (automatic differentiation) and NumDiffTools (finite differences). It offers the possibility to compute directly the derivative of the complex signal. It is also possible to compute analytically the derivatives w.r.t. dL, theta, phi, psi, tcoal and Phicoal, and also iota in absence of HM or precessing spins. - + if (self.wf_model.is_newtonian): if self.verbose: print('WARNING: In the Newtonian inspiral case the mass ratio and spins do not enter the waveform, and the corresponding Fisher matrix elements vanish, we then discard them.\n') - + if computeAnalyticalDeriv: derivargs = (1) inputNumdL, inputNumiota = 1, 2 @@ -1100,27 +1076,27 @@ def _SignalDerivatives(self, fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, derivargs = (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17) if not self.wf_model.is_tidal: derivargs = derivargs[:-2] - + if self.wf_model.is_eccentric: derivargs = derivargs + (18,) - + nParams = self.wf_model.nParams - + if not computeDerivFinDiff: if self.wf_model.is_holomorphic: GWstrainUse = lambda f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc: self.GWstrain(f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=rot, is_m1m2=use_m1m2, is_chi1chi2=use_chi1chi2, is_prec_ang=use_prec_ang) - + FisherDerivs = np.asarray(vmap(jacrev(GWstrainUse, argnums=derivargs, holomorphic=True))(fgrids.T, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc)) else: # In the non holomorphic case, to improve the accuracy, we compute separately the derivatives of the real and imaginary part of the strain as real functions fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc = np.real(fgrids), np.real(Mc), np.real(eta), np.real(dL), np.real(theta), np.real(phi), np.real(iota), np.real(psi), np.real(tcoal), np.real(Phicoal), np.real(chiS), np.real(chiA), np.real(chi1x), np.real(chi2x), np.real(chi1y), np.real(chi2y), np.real(LambdaTilde), np.real(deltaLambda), np.real(ecc) - + GWstrainUse_real = lambda f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc: np.real(self.GWstrain(f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=rot, is_m1m2=use_m1m2, is_chi1chi2=use_chi1chi2, is_prec_ang=use_prec_ang)) GWstrainUse_imag = lambda f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc: np.imag(self.GWstrain(f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=rot, is_m1m2=use_m1m2, is_chi1chi2=use_chi1chi2, is_prec_ang=use_prec_ang)) - + realDerivs = np.asarray(vmap(jacrev(GWstrainUse_real, argnums=derivargs))(fgrids.T, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc)) imagDerivs = np.asarray(vmap(jacrev(GWstrainUse_imag, argnums=derivargs))(fgrids.T, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc)) - + FisherDerivs = realDerivs + 1j*imagDerivs else: if self.wf_model.is_newtonian: @@ -1208,7 +1184,7 @@ def _SignalDerivatives(self, fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, else: GWstrainUse = lambda pars: self.GWstrain(fgrids, pars[0], pars[1], dL, theta, phi, pars[2], psi, tcoal, Phicoal, pars[3], pars[4], chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, pars[5], rot=rot, is_m1m2=use_m1m2, is_chi1chi2=use_chi1chi2) evpars = [Mc, eta, iota, chiS, chiA, ecc] - + dh = ndt.Jacobian(GWstrainUse, step=stepNDT, method=methodNDT, order=2, n=1) FisherDerivs = np.asarray(dh(evpars)) if len(FisherDerivs.shape) == 2: #len(Mc) == 1: @@ -1221,7 +1197,7 @@ def _SignalDerivatives(self, fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, NAnalyticalDerivs = 7 else: NAnalyticalDerivs = 6 - + dL_deriv, theta_deriv, phi_deriv, iota_deriv, psi_deriv, tc_deriv, Phicoal_deriv = self._AnalyticalDerivatives(fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=rot, use_m1m2=use_m1m2, use_chi1chi2=use_chi1chi2, use_prec_ang=use_prec_ang) if (not self.wf_model.is_HigherModes) and (not self.wf_model.is_Precessing): if not self.wf_model.is_newtonian: @@ -1232,13 +1208,13 @@ def _SignalDerivatives(self, fgrids, Mc, eta, dL, theta, phi, iota, psi, tcoal, else: tmpsplit1, tmpsplit2, tmpsplit3, _ = onp.vsplit(FisherDerivs, onp.array([inputNumdL, inputNumiota, nParams-NAnalyticalDerivs])) FisherDerivs = np.vstack((tmpsplit1, np.asarray(dL_deriv).T[np.newaxis,:], np.asarray(theta_deriv).T[np.newaxis,:], np.asarray(phi_deriv).T[np.newaxis,:], tmpsplit2, np.asarray(psi_deriv).T[np.newaxis,:], np.asarray(tc_deriv).T[np.newaxis,:], np.asarray(Phicoal_deriv).T[np.newaxis,:], tmpsplit3)) - + return FisherDerivs - + def _AnalyticalDerivatives(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, Phicoal, chiS, chiA, chi1x, chi2x, chi1y, chi2y, LambdaTilde, deltaLambda, ecc, rot=0., use_m1m2=False, use_chi1chi2=False, use_prec_ang=False): """ Compute analytical derivatives with respect to ``dL``, ``theta``, ``phi``, ``psi``, ``tcoal``, ``Phicoal`` and ``iota`` (the latter only for the fundamental mode in the non-precessing case). - + :param array or float f: The frequency(ies) at which to perform the calculation, in :math:`\\rm Hz`. :param array or float Mc: The chirp mass(es), :math:`{\cal M}_c`, in units of :math:`\\rm M_{\odot}`. If ``use_m1m2=True`` this is interpreted as the primary mass, :math:`m_1`, in units of :math:`\\rm M_{\odot}`. :param array or float eta: The symmetric mass ratio(s), :math:`\eta`. If ``use_m1m2=True`` this is interpreted as the secondary mass, :math:`m_2`, in units of :math:`\\rm M_{\odot}`. @@ -1264,7 +1240,7 @@ def _AnalyticalDerivatives(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, P :param bool, optional use_prec_ang: Boolean specifying if the ``iota`` input should be interpreted as the inclination angle with respect to total angular momentum, ``chiS`` and ``chiA`` as the primary and secondary spin magnitudes, ``chi1x`` and ``chi2x`` as the primary and secondary spin tilts, ``chi1y`` as the azimuthal angle of orbital angular momentum relative to total angular momentum and ``chi2y`` as the difference in azimuthal angle between spin vectors. :return: Analytical derivatives with respect to ``dL``, ``theta``, ``phi``, ``iota``, ``psi``, ``tcoal`` and ``Phicoal``. If the :py:class:`self.wf_model` is precessing or includes higher order modes the derivative with respect to ``iota`` will be ``None`` :rtype: tuple(array, array, array, array, array, array, array) - + """ # Module to compute analytically the derivatives w.r.t. dL, theta, phi, psi, tcoal, Phicoal and also iota in absence of HM or precessing spins. Each derivative is inserted into its own function with representative name, for ease of check. if use_m1m2: @@ -1273,7 +1249,7 @@ def _AnalyticalDerivatives(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, P else: McUse = Mc etaUse = eta - + if not self.wf_model.is_Precessing: if use_chi1chi2: # Interpret chiS as chi1z and chiA as chi2z @@ -1294,18 +1270,18 @@ def _AnalyticalDerivatives(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, P else: # convert angles and iota iota, chi1xUse, chi1yUse, chi1z, chi2xUse, chi2yUse, chi2z = utils.TransformPrecessing_angles2comp(thetaJN=iota, phiJL=chi1y, theta1=chi1x, theta2=chi2x, phi12=chi2y, chi1=chiS, chi2=chiA, Mc=McUse, eta=etaUse, fRef=self.fmin, phiRef=0.) - + evParams = {'Mc':McUse, 'dL':dL, 'theta':theta, 'phi':phi, 'iota':iota, 'psi':psi, 'tcoal':tcoal, 'eta':etaUse, 'Phicoal':Phicoal, 'chi1z':chi1z, 'chi2z':chi2z, 'chi1x':chi1xUse, 'chi2x':chi2xUse, 'chi1y':chi1yUse, 'chi2y':chi2yUse} - + if self.wf_model.is_tidal: Lambda1, Lambda2 = utils.Lam12_from_Lamt_delLam(LambdaTilde, deltaLambda, etaUse) - + evParams['Lambda1'] = Lambda1 evParams['Lambda2'] = Lambda2 - + if self.wf_model.is_eccentric: evParams['ecc'] = ecc - + if (not self.wf_model.is_HigherModes) and (not self.wf_model.is_Precessing): wfPhiGw = self.wf_model.Phi(f, **evParams) wfAmpl = self.wf_model.Ampl(f, **evParams) @@ -1313,7 +1289,7 @@ def _AnalyticalDerivatives(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, P else: # If the waveform includes higher modes, it is not possible to compute amplitude and phase separately, make all together wfhp, wfhc = self.wf_model.hphc(f, **evParams) - + if self.useEarthMotion: # Compute Doppler contribution tnoloc = tcoal - self.wf_model.tau_star(f, **evParams)/(3600.*24.) @@ -1331,11 +1307,11 @@ def _AnalyticalDerivatives(self, f, Mc, eta, dL, theta, phi, iota, psi, tcoal, P tnoloc = tcoal tmpDeltLoc = self._DeltLoc(theta, phi, tnoloc) # in seconds t = tnoloc + tmpDeltLoc/(3600.*24.) - + phiL = (2.*np.pi*f)*tmpDeltLoc - + rot_rad = rot*np.pi/180. - + def afun(ra, dec, t, rot): phir = self.det_long_rad a1 = 0.0625*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*(3.-np.cos(2.*dec))*np.cos(2.*(ra - phir - 2.*np.pi*t)) @@ -1344,44 +1320,44 @@ def afun(ra, dec, t, rot): a4 = 0.5*np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.sin(2.*dec)*np.sin(ra - phir - 2.*np.pi*t) a5 = 3.*0.25*np.sin(2*(self.det_xax_rad+rot))*(np.cos(self.det_lat_rad)*np.cos(dec))**2. return a1 - a2 + a3 - a4 + a5 - + def bfun(ra, dec, t, rot): phir = self.det_long_rad b1 = np.cos(2*(self.det_xax_rad+rot))*np.sin(self.det_lat_rad)*np.sin(dec)*np.cos(2.*(ra - phir - 2.*np.pi*t)) b2 = 0.25*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*np.sin(dec)*np.sin(2.*(ra - phir - 2.*np.pi*t)) b3 = np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.cos(dec)*np.cos(ra - phir - 2.*np.pi*t) b4 = 0.5*np.sin(2*(self.det_xax_rad+rot))*np.sin(2.*self.det_lat_rad)*np.cos(dec)*np.sin(ra - phir - 2.*np.pi*t) - + return b1 + b2 + b3 + b4 - + ras, decs = self._ra_dec_from_th_phi(theta, phi) afac = afun(ras, decs, t, rot_rad) bfac = bfun(ras, decs, t, rot_rad) - + Fp = np.sin(self.angbtwArms)*(afac*np.cos(2.*psi) + bfac*np.sin(2*psi)) Fc = np.sin(self.angbtwArms)*(bfac*np.cos(2.*psi) - afac*np.sin(2*psi)) - + hp, hc = wfhp*Fp*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)), wfhc*Fc*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) - + def psi_par_deriv(): - + Fp_psider = 2*np.sin(self.angbtwArms)*(-afac*np.sin(2.*psi) + bfac*np.cos(2*psi)) Fc_psider = 2*np.sin(self.angbtwArms)*(-bfac*np.sin(2.*psi) - afac*np.cos(2*psi)) - + return wfhp*Fp_psider*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) + wfhc*Fc_psider*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) - + def phi_par_deriv(): - + def Delt_loc_phider(ra, dec, t): - + comp1 = -np.cos(dec)*np.sin(ra)*np.cos(self.det_lat_rad)*np.cos(self.det_long_rad + 2.*np.pi*t) comp2 = np.cos(dec)*np.cos(ra)*np.cos(self.det_lat_rad)*np.sin(self.det_long_rad + 2.*np.pi*t) - + Delt_phider = - glob.REarth*(comp1+comp2)/glob.clight - + return Delt_phider/(3600.*24.) # in days - + def afun_phider(ra, dec, t, rot): phir = self.det_long_rad a1 = 0.0625*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*(3.-np.cos(2.*dec))*(-2.)*np.sin(2.*(ra - phir - 2.*np.pi*t)) @@ -1390,41 +1366,41 @@ def afun_phider(ra, dec, t, rot): a4 = 0.5*np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.sin(2.*dec)*np.cos(ra - phir - 2.*np.pi*t) return a1 - a2 + a3 - a4 - + def bfun_phider(ra, dec, t, rot): phir = self.det_long_rad - + b1 = np.cos(2*(self.det_xax_rad+rot))*np.sin(self.det_lat_rad)*np.sin(dec)*(-2.)*np.sin(2.*(ra - phir - 2.*np.pi*t)) b2 = 0.25*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*np.sin(dec)*2.*np.cos(2.*(ra - phir - 2.*np.pi*t)) b3 = -np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.cos(dec)*np.sin(ra - phir - 2.*np.pi*t) b4 = 0.5*np.sin(2*(self.det_xax_rad+rot))*np.sin(2.*self.det_lat_rad)*np.cos(dec)*np.cos(ra - phir - 2.*np.pi*t) - + return b1 + b2 + b3 + b4 locDt_phider = Delt_loc_phider(ras, decs, tnoloc) afac_phider = afun_phider(ras, decs, t, rot_rad)*(1.-2.*np.pi*locDt_phider) bfac_phider = bfun_phider(ras, decs, t, rot_rad)*(1.-2.*np.pi*locDt_phider) - + Fp_phider = np.sin(self.angbtwArms)*(afac_phider*np.cos(2.*psi) + bfac_phider*np.sin(2*psi)) Fc_phider = np.sin(self.angbtwArms)*(bfac_phider*np.cos(2.*psi) - afac_phider*np.sin(2*psi)) - + ampP_phider = wfhp*Fp_phider*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) ampC_phider = wfhc*Fc_phider*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) phiD_phideriv = 0. phiL_phideriv = 2.*np.pi*f*locDt_phider*(3600.*24.) - + return ampP_phider + 1j*(phiD_phideriv + phiL_phideriv)*hp + ampC_phider + 1j*(phiD_phideriv + phiL_phideriv)*hc - + def theta_par_deriv(): def Delt_loc_thder(ra, dec, t): - + comp1 = np.sin(dec)*np.cos(ra)*np.cos(self.det_lat_rad)*np.cos(self.det_long_rad + 2.*np.pi*t) comp2 = np.sin(dec)*np.sin(ra)*np.cos(self.det_lat_rad)*np.sin(self.det_long_rad + 2.*np.pi*t) comp3 = -np.cos(dec)*np.sin(self.det_lat_rad) - + Delt_thder = - glob.REarth*(comp1+comp2+comp3)/glob.clight - + return Delt_thder/(3600.*24.) # in days - + def afun_thder(ra, dec, t, rot, loc_thder): phir = self.det_long_rad a1 = 0.0625*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*(-2.*np.sin(2.*dec))*np.cos(2.*(ra - phir - 2.*np.pi*t)) + 0.0625*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad)) *(3.-np.cos(2.*dec))*np.sin(2.*(ra - phir - 2.*np.pi*t))*(4.*np.pi*loc_thder) @@ -1433,41 +1409,41 @@ def afun_thder(ra, dec, t, rot, loc_thder): a4 = -0.5*np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*2.*np.cos(2.*dec)*np.sin(ra - phir - 2.*np.pi*t) - 0.5*np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.sin(2.*dec)*np.cos(ra - phir - 2.*np.pi*t)*(2.*np.pi*loc_thder) a5 = 2.*3.*0.25*np.sin(2*(self.det_xax_rad+rot))*((np.cos(self.det_lat_rad))**2)*np.cos(dec)*np.sin(dec) return a1 - a2 + a3 - a4 + a5 - + def bfun_thder(ra, dec, t, rot, loc_thder): phir = self.det_long_rad b1 = -np.cos(2*(self.det_xax_rad+rot))*np.sin(self.det_lat_rad)*np.cos(dec)*np.cos(2.*(ra - phir - 2.*np.pi*t)) + np.cos(2*(self.det_xax_rad+rot))*np.sin(self.det_lat_rad)*np.sin(dec)*np.sin(2.*(ra - phir - 2.*np.pi*t))*(4.*np.pi*loc_thder) b2 = -0.25*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*np.cos(dec)*np.sin(2.*(ra - phir - 2.*np.pi*t)) - 0.25*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*np.sin(dec)*np.cos(2.*(ra - phir - 2.*np.pi*t))*(4.*np.pi*loc_thder) b3 = np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.sin(dec)*np.cos(ra - phir - 2.*np.pi*t) + np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.cos(dec)*np.sin(ra - phir - 2.*np.pi*t)*(2.*np.pi*loc_thder) b4 = 0.5*np.sin(2*(self.det_xax_rad+rot))*np.sin(2.*self.det_lat_rad)*np.sin(dec)*np.sin(ra - phir - 2.*np.pi*t) - 0.5*np.sin(2*(self.det_xax_rad+rot))*np.sin(2.*self.det_lat_rad)*np.cos(dec)*np.cos(ra - phir - 2.*np.pi*t)*(2.*np.pi*loc_thder) - + return b1 + b2 + b3 + b4 - + locDt_thder = Delt_loc_thder(ras, decs, tnoloc) afac_thder = afun_thder(ras, decs, t, rot_rad, locDt_thder) bfac_thder = bfun_thder(ras, decs, t, rot_rad, locDt_thder) - + Fp_thder = np.sin(self.angbtwArms)*(afac_thder*np.cos(2.*psi) + bfac_thder*np.sin(2*psi)) Fc_thder = np.sin(self.angbtwArms)*(bfac_thder*np.cos(2.*psi) - afac_thder*np.sin(2*psi)) - + ampP_thder = wfhp*Fp_thder*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) ampC_thder = wfhc*Fc_thder*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) phiD_thderiv = 0. phiL_thderiv = 2.*np.pi*f*locDt_thder*(3600.*24.) - + return ampP_thder + 1j*(phiD_thderiv + phiL_thderiv)*hp + ampC_thder + 1j*(phiD_thderiv + phiL_thderiv)*hc - + def tcoal_par_deriv(): - + def Delt_loc_tcder(ra, dec, t): - + comp1 = -np.cos(decs)*np.cos(ras)*np.cos(self.det_lat_rad)*np.sin(self.det_long_rad + 2.*np.pi*t) comp2 = np.cos(decs)*np.sin(ras)*np.cos(self.det_lat_rad)*np.cos(self.det_long_rad + 2.*np.pi*t) - + Delt_tcder = - 2.*np.pi*glob.REarth*(comp1+comp2)/glob.clight - + return Delt_tcder/(3600.*24.) # in days - + def afun_tcder(ra, dec, t, rot): phir = self.det_long_rad a1 = 0.0625*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*(3.-np.cos(2.*dec))*(-2.)*np.sin(2.*(ra - phir - 2.*np.pi*t)) @@ -1476,38 +1452,38 @@ def afun_tcder(ra, dec, t, rot): a4 = 0.5*np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.sin(2.*dec)*np.cos(ra - phir - 2.*np.pi*t) return a1 - a2 + a3 - a4 - + def bfun_tcder(ra, dec, t, rot): phir = self.det_long_rad b1 = np.cos(2*(self.det_xax_rad+rot))*np.sin(self.det_lat_rad)*np.sin(dec)*(-2.)*np.sin(2.*(ra - phir - 2.*np.pi*t)) b2 = 0.25*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2.*self.det_lat_rad))*np.sin(dec)*2.*np.cos(2.*(ra - phir - 2.*np.pi*t)) b3 = -np.cos(2*(self.det_xax_rad+rot))*np.cos(self.det_lat_rad)*np.cos(dec)*np.sin(ra - phir - 2.*np.pi*t) b4 = 0.5*np.sin(2*(self.det_xax_rad+rot))*np.sin(2.*self.det_lat_rad)*np.cos(dec)*np.cos(ra - phir - 2.*np.pi*t) - + return b1 + b2 + b3 + b4 locDt_tcder = Delt_loc_tcder(ras, decs, tnoloc) afac_tcder = -2.*np.pi*afun_tcder(ras, decs, t, rot_rad)*(1.+locDt_tcder) bfac_tcder = -2.*np.pi*bfun_tcder(ras, decs, t, rot_rad)*(1.+locDt_tcder) - + Fp_tcder = np.sin(self.angbtwArms)*(afac_tcder*np.cos(2.*psi) + bfac_tcder*np.sin(2*psi)) Fc_tcder = np.sin(self.angbtwArms)*(bfac_tcder*np.cos(2.*psi) - afac_tcder*np.sin(2*psi)) - + ampP_tcder = wfhp*Fp_tcder*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) ampC_tcder = wfhc*Fc_tcder*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) phiD_tcderiv = 0. phiL_tcderiv = 2.*np.pi*f*locDt_tcder*(3600.*24.) return ampP_tcder + 1j*(phiD_tcderiv + phiL_tcderiv + 2.*np.pi*f*3600.*24.)*hp + ampC_tcder + 1j*(phiD_tcderiv + phiL_tcderiv + 2.*np.pi*f*3600.*24.)*hc - + def iota_par_deriv(): - + if (not self.wf_model.is_HigherModes) and (not self.wf_model.is_Precessing): wfhp_iotader, wfhc_iotader = -wfAmpl*np.exp(-1j*wfPhiGw)*(np.cos(iota)*np.sin(iota)), -1j*wfAmpl*np.exp(-1j*wfPhiGw)*np.sin(iota) return wfhp_iotader*Fp*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) + wfhc_iotader*Fc*np.exp(1j*(2.*np.pi*f*(tcoal*3600.*24.) - Phicoal + phiD + phiL)) else: # This derivative is computed numerically if the waveform contains higher modes return None - + dL_deriv = -(hp+hc)/dL Phicoal_deriv = -1j*(hp+hc) psi_deriv = psi_par_deriv() @@ -1515,125 +1491,125 @@ def iota_par_deriv(): theta_deriv = theta_par_deriv() tc_deriv = tcoal_par_deriv() iota_deriv = iota_par_deriv() - + return dL_deriv, theta_deriv, phi_deriv, iota_deriv, psi_deriv, tc_deriv, Phicoal_deriv - + def optimal_location(self, tcoal, is_tGPS=False): """ Compute the optimal sky position for a signal to be seen by the detector at a given time. - + The computation assumes :math:`\psi = 0`. - + :param float tcoal: The time at which to compute the optimal location, as GMST in days. :param bool, optional is_tGPS: Boolean specifying if the provided time is a GPS time (in seconds) rather than a GMST. - + :return: Optimal :math:`\\theta` and :math:`\phi` sky coordinates, in :math:`\\rm rad`. :rtype: array(float, float) - + """ # Function to compute the optimal theta and phi for a signal to be seen by the detector at a given GMST. The boolean is_tGPS can be used to specify whether the provided time is a GPS time rather than a GMST, so that it will be converted. # For a triangle the best location is the same of an L in the same place, as can be shown by explicit geometrical computation. # Even if considering Earth rotation, the highest SNR will still be obtained if the source is in the optimal location close to the merger. from scipy.optimize import minimize - + if is_tGPS: tc = utils.GPSt_to_LMST(tcoal, lat=0., long=0.) else: tc = tcoal - + def pattern_fixedtpsi(pars, tc=tc): theta, phi = pars Fp, Fc = self._PatternFunction(theta, phi, t=tc, psi=0) return -np.sqrt(Fp**2 + Fc**2) # we actually minimize the pattern function times -1, which is the same as maximizing it return minimize(pattern_fixedtpsi, [1.,1.], bounds=((0.,onp.pi), (0.,2.*onp.pi))).x - + def SNRFastInsp(self, evParams, checkInterp=False): """ Compute the inspiral SNR taking into account Earth rotation, without the need of performing an integral for each event - + .. deprecated:: 1.0.0 Use the standard function :py:class:`GWSignal.SNRInteg`. """ - # This module allows to compute the inspiral SNR taking into account Earth rotation, without the need + # This module allows to compute the inspiral SNR taking into account Earth rotation, without the need # of performing an integral for each event - - Mc, dL, theta, phi, iota, psi, tcoal, eta = evParams['Mc'], evParams['dL'], evParams['theta'], evParams['phi'], evParams['iota'], evParams['psi'], evParams['tcoal'], evParams['eta'] - + + Mc, dL, theta, phi, iota, psi, tcoal, eta = evParams['Mc'], evParams['dL'], evParams['theta'], evParams['phi'], evParams['iota'], evParams['psi'], evParams['tcoal'], evParams['eta'] + ras, decs = self._ra_dec_from_th_phi(theta, phi) - - + + if not np.isscalar(Mc): SNR = np.zeros(Mc.shape) else: SNR = 0 - + # Factor in front of the integral in the inspiral only case fac = np.sqrt(5./6.)/np.pi**(2./3.)*(glob.GMsun_over_c3*Mc)**(5./6.)*glob.clightGpc/dL#*np.exp(-logdL) - + fcut = self.wf_model.fcut(**evParams) if self.fmax is not None: fcut = np.where(fcut > self.fmax, self.fmax, fcut) mask = self.strainFreq >= self.fmin - + def CoeffsRot(ra, dec, psi, rot=0.): rot = rot*np.pi/180. rasDet = ra - self.det_long_rad # Referring to overleaf, I now call VC2 the last vector appearing in the C2 expression, VS2 the one in the S2 expression and so on # e1 is the first element and e2 the second - + VC2e1 = 0.0675*np.cos(2*rasDet)*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2*dec))*(3.-np.cos(2.*self.det_lat_rad)) - 0.25*np.sin(2*rasDet)*np.cos(2*(self.det_xax_rad+rot))*(3.-np.cos(2*dec))*np.sin(self.det_lat_rad) VC2e2 = 0.25*np.sin(2*rasDet)*np.sin(2*(self.det_xax_rad+rot))*np.sin(dec)*(3.-np.cos(2.*self.det_lat_rad)) + np.cos(2*rasDet)*np.cos(2*(self.det_xax_rad+rot))*np.sin(dec)*np.sin(self.det_lat_rad) C2p = (np.cos(2.*psi)*VC2e1 + np.sin(2*psi)*VC2e2)*np.sin(self.angbtwArms) C2c = (-np.sin(2*psi)*VC2e1 + np.cos(2.*psi)*VC2e2)*np.sin(self.angbtwArms) - + VS2e1 = 0.0675*np.sin(2*rasDet)*np.sin(2*(self.det_xax_rad+rot))*(3.-np.cos(2*dec))*(3.-np.cos(2.*self.det_lat_rad)) + 0.25*np.cos(2*rasDet)*np.cos(2*(self.det_xax_rad+rot))*(3.-np.cos(2*dec))*np.sin(self.det_lat_rad) VS2e2 = -0.25*np.cos(2*rasDet)*np.sin(2*(self.det_xax_rad+rot))*np.sin(dec)*(3.-np.cos(2.*self.det_lat_rad)) + np.sin(2*rasDet)*np.cos(2*(self.det_xax_rad+rot))*np.sin(dec)*np.sin(self.det_lat_rad) S2p = (np.cos(2.*psi)*VS2e1 + np.sin(2*psi)*VS2e2)*np.sin(self.angbtwArms) S2c = (-np.sin(2*psi)*VS2e1 + np.cos(2.*psi)*VS2e2)*np.sin(self.angbtwArms) - + VC1e1 = 0.25*np.cos(rasDet)*np.sin(2.*(self.det_xax_rad+rot))*np.sin(2*dec)*np.sin(2.*self.det_lat_rad) - 0.5*np.sin(rasDet)*np.cos(2.*(self.det_xax_rad+rot))*np.sin(2*dec)*np.cos(self.det_lat_rad) VC1e2 = np.cos(rasDet)*np.cos(2*(self.det_xax_rad+rot))*np.cos(dec)*np.cos(self.det_lat_rad) + 0.5*np.sin(rasDet)*np.sin(2.*(self.det_xax_rad+rot))*np.cos(dec)*np.sin(2.*self.det_lat_rad) C1p = (np.cos(2.*psi)*VC1e1 + np.sin(2*psi)*VC1e2)*np.sin(self.angbtwArms) C1c = (-np.sin(2*psi)*VC1e1 + np.cos(2.*psi)*VC1e2)*np.sin(self.angbtwArms) - + VS1e1 = 0.25*np.sin(rasDet)*np.sin(2.*(self.det_xax_rad+rot))*np.sin(2*dec)*np.sin(2.*self.det_lat_rad) + 0.5*np.cos(rasDet)*np.cos(2.*(self.det_xax_rad+rot))*np.sin(2*dec)*np.cos(self.det_lat_rad) VS1e2 = np.sin(rasDet)*np.cos(2*(self.det_xax_rad+rot))*np.cos(dec)*np.cos(self.det_lat_rad) - 0.5*np.cos(rasDet)*np.sin(2.*(self.det_xax_rad+rot))*np.cos(dec)*np.sin(2.*self.det_lat_rad) S1p = (np.cos(2.*psi)*VS1e1 + np.sin(2*psi)*VS1e2)*np.sin(self.angbtwArms) S1c = (-np.sin(2*psi)*VS1e1 + np.cos(2.*psi)*VS1e2)*np.sin(self.angbtwArms) - + C0p = 0.75*np.cos(2.*psi)*np.sin(2.*(self.det_xax_rad+rot))*((np.cos(dec)*np.cos(self.det_lat_rad))**2)*np.sin(self.angbtwArms) C0c = -0.75*np.sin(2.*psi)*np.sin(2.*(self.det_xax_rad+rot))*((np.cos(dec)*np.cos(self.det_lat_rad))**2)*np.sin(self.angbtwArms) - + return np.array([C2p, C2c]), np.array([S2p, S2c]), np.array([C1p, C1c]), np.array([S1p, S1c]), np.array([C0p, C0c]) - + def FpFcsqInt(C2s, S2s, C1s, S1s, C0s, Igs, iota): - Fp4 = 0.5*(C2s[0]**2 - S2s[0]**2)*Igs[3] + C2s[0]*S2s[0]*Igs[7] - + Fp4 = 0.5*(C2s[0]**2 - S2s[0]**2)*Igs[3] + C2s[0]*S2s[0]*Igs[7] + Fp3 = (C2s[0]*C1s[0] - S2s[0]*S1s[0])*Igs[2] + (C2s[0]*S1s[0] + S2s[0]*C1s[0])*Igs[6] - + Fp2 = (0.5*(C1s[0]**2 - S1s[0]**2) + 2.*C2s[0]*C0s[0])*Igs[1] + (2.*C0s[0]*S2s[0] + C1s[0]*S1s[0])*Igs[5] - + Fp1 = (2.*C0s[0]*C1s[0] + C1s[0]*C2s[0] + S2s[0]*S1s[0])*Igs[0] + (2.*C0s[0]*S1s[0] + C1s[0]*S2s[0] - S1s[0]*C2s[0])*Igs[4] - + Fp0 = (C0s[0]**2 + 0.5*(C1s[0]**2 + C2s[0]**2 + S1s[0]**2 + S2s[0]**2))*Igs[8] - + FpsqInt = Fp4 + Fp3 + Fp2 + Fp1 + Fp0 - - Fc4 = 0.5*(C2s[1]**2 - S2s[1]**2)*Igs[3] + C2s[1]*S2s[1]*Igs[7] - + + Fc4 = 0.5*(C2s[1]**2 - S2s[1]**2)*Igs[3] + C2s[1]*S2s[1]*Igs[7] + Fc3 = (C2s[1]*C1s[1] - S2s[1]*S1s[1])*Igs[2] + (C2s[1]*S1s[1] + S2s[1]*C1s[1])*Igs[6] - + Fc2 = (0.5*(C1s[1]**2 - S1s[1]**2) + 2.*C2s[1]*C0s[1])*Igs[1] + (2.*C0s[1]*S2s[1] + C1s[1]*S1s[1])*Igs[5] - + Fc1 = (2.*C0s[1]*C1s[1] + C1s[1]*C2s[1] + S2s[1]*S1s[1])*Igs[0] + (2.*C0s[1]*S1s[1] + C1s[1]*S2s[1] - S1s[1]*C2s[1])*Igs[4] - + Fc0 = (C0s[1]**2 + 0.5*(C1s[1]**2 + C2s[1]**2 + S1s[1]**2 + S2s[1]**2))*Igs[8] - + FcsqInt = Fc4 + Fc3 + Fc2 + Fc1 + Fc0 - + return FpsqInt*(0.5*(1.+(np.cos(iota))**2))**2, FcsqInt*(np.cos(iota))**2 - + if not self.useEarthMotion: t = tcoal - self.wf_model.tau_star(self.fmin, **evParams)/(3600.*24) if self.detector_shape=='L': @@ -1661,11 +1637,11 @@ def IntegrandC(f, Mc, tcoal, n): def IntegrandS(f, Mc, tcoal, n): t = tcoal - 2.18567 * ((1.21/Mc)**(5./3.)) * ((100/f)**(8./3.))/(3600.*24) return (f**(-7./3.))*np.sin(n*2.*np.pi*t) - + fminarr = np.full(fcut.shape, self.fmin) fgrids = np.geomspace(fminarr,fcut,num=int(5000)) strainGrids = np.interp(fgrids, self.strainFreq, self.noiseCurve, left=1., right=1.) - + for m in range(4): tmpIntegrandC = IntegrandC(fgrids, Mc, tcoal, m+1.) tmpIntegrandS = IntegrandS(fgrids, Mc, tcoal, m+1.) @@ -1673,8 +1649,8 @@ def IntegrandS(f, Mc, tcoal, n): Igs[m+4,:] = onp.trapz(tmpIntegrandS/strainGrids, fgrids, axis=0) tmpIntegrand = IntegrandC(fgrids, Mc, tcoal, 0.) Igs[8,:] = onp.trapz(tmpIntegrand/strainGrids, fgrids, axis=0) - - + + if self.detector_shape=='L': C2s, S2s, C1s, S1s, C0s = CoeffsRot(ras, decs, psi, rot=0.) FpsqInt, FcsqInt = FpFcsqInt(C2s, S2s, C1s, S1s, C0s, Igs, iota) @@ -1682,19 +1658,19 @@ def IntegrandS(f, Mc, tcoal, n): SNR = fac * np.sqrt(QsqInt) elif self.detector_shape=='T': for i in range(3): - C2s, S2s, C1s, S1s, C0s = CoeffsRot(ras, decs, psi, rot=i*60.) + C2s, S2s, C1s, S1s, C0s = CoeffsRot(ras, decs, psi, rot=i*60.) FpsqInt, FcsqInt = FpFcsqInt(C2s, S2s, C1s, S1s, C0s, Igs, iota) QsqInt = FpsqInt + FcsqInt tmpSNR = fac * np.sqrt(QsqInt) SNR = SNR + tmpSNR*tmpSNR SNR = np.sqrt(SNR) - + return SNR - + def WFOverlap(self, WF1, WF2, evParams1, evParams2, res=1000, return_separate=False, **kwargs): """ Compute the *overlap* of two waveforms in a single detector on two sets of parameters, for one or multiple events. - + :param WaveFormModel WF1: Object containing the first waveform model to analyse. :param WaveFormModel WF2: Object containing the second waveform model to analyse. :param dict(array, array, ...) evParams1: Dictionary containing the parameters of the event(s) for the first waveform model, as in :py:data:`events`. @@ -1702,14 +1678,14 @@ def WFOverlap(self, WF1, WF2, evParams1, evParams2, res=1000, return_separate=Fa :param int res: The resolution of the frequency grid to use. :param bool, optional return_all: Boolean specifying if, instead of returning the overlap, the function has to return separately product at the numerator of the definition, :math:`(h_1|h_2)`, and the SNRs at the denominator. This is needed to compute the overlap for a detector network. In this case the return type is *tuple(array, array, array)*. :param unused kwargs: Optional arguments. - + :return: Overlap(s) of the two waveforms. The shape is :math:`(N_{\\rm events})`. :rtype: 1-D array - + """ utils.check_evparams(evParams1) utils.check_evparams(evParams2) - + # Checks on imput parameters for waveform 1 if WF1.is_Precessing: @@ -1732,7 +1708,7 @@ def WFOverlap(self, WF1, WF2, evParams1, evParams2, res=1000, return_separate=Fa evParams1['chi2z'] = evParams1['chiS'] - evParams1['chiA'] except KeyError: raise ValueError('Two among chi1z, chi2z and chiS, chiA have to be provided.') - + if WF1.is_tidal: try: _=evParams1['LambdaTilde'] @@ -1743,7 +1719,7 @@ def WFOverlap(self, WF1, WF2, evParams1, evParams2, res=1000, return_separate=Fa raise ValueError('Two among Lambda1, Lambda2 and LambdaTilde and deltaLambda have to be provided.') else: evParams1['LambdaTilde'], evParams1['deltaLambda'] = np.zeros_like(evParams1['Mc']), np.zeros_like(evParams1['Mc']) - + if not WF1.is_eccentric: evParams1['ecc'] = np.zeros_like(evParams1['Mc']) @@ -1769,7 +1745,7 @@ def WFOverlap(self, WF1, WF2, evParams1, evParams2, res=1000, return_separate=Fa evParams2['chi2z'] = evParams2['chiS'] - evParams2['chiA'] except KeyError: raise ValueError('Two among chi1z, chi2z and chiS, chiA have to be provided.') - + if WF2.is_tidal: try: _=evParams2['LambdaTilde'] @@ -1780,7 +1756,7 @@ def WFOverlap(self, WF1, WF2, evParams1, evParams2, res=1000, return_separate=Fa raise ValueError('Two among Lambda1, Lambda2 and LambdaTilde and deltaLambda have to be provided.') else: evParams2['LambdaTilde'], evParams2['deltaLambda'] = np.zeros_like(evParams2['Mc']), np.zeros_like(evParams2['Mc']) - + if not WF2.is_eccentric: evParams2['ecc'] = np.zeros_like(evParams2['Mc'])