Commit f8ecd1d7 authored by Maria Haney's avatar Maria Haney

Merge branch 'PN_modes_master_rebase' into 'master'

Precessing PN modes

See merge request lscsoft/lalsuite!1010
parents cac7e583 89440468
......@@ -50,6 +50,7 @@ extern "C" {
* @defgroup LALSimIMRPSpinInspiralRD_c LALSimIMRPSpinInspiralRD.c
* @defgroup LALSimIMRTidal_c LALSimIMRLackeyTidal2013.c
* @defgroup LALSimPrecessingNRSur_c LALSimIMRPrecessingNRSur.c
* @defgroup LALSimIMRNRWaveforms_c LALSimIMRNRWaveforms.c
* @}
*
* @addtogroup LALSimIMR_h
......@@ -106,7 +107,7 @@ int XLALSimIMRPhenomPCalculateModelParametersFromSourceFrame(REAL8 *chi1_l, REAL
int XLALSimIMREOBNRv2DominantMode(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, const REAL8 phiC, const REAL8 deltaT, const REAL8 m1SI, const REAL8 m2SI, const REAL8 fLower, const REAL8 distance, const REAL8 inclination);
int XLALSimIMREOBNRv2AllModes(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, const REAL8 phiC, const REAL8 deltaT, const REAL8 m1SI, const REAL8 m2SI, const REAL8 fLower, const REAL8 distance, const REAL8 inclination);
SphHarmTimeSeries *XLALSimIMREOBNRv2Modes(const REAL8 phiRef, const REAL8 deltaT, const REAL8 m1, const REAL8 m2, const REAL8 fLower, const REAL8 distance);
SphHarmTimeSeries *XLALSimIMREOBNRv2Modes(const REAL8 deltaT, const REAL8 m1, const REAL8 m2, const REAL8 fLower, const REAL8 distance);
/* in module LALSimIMRSpinAlignedEOB.c */
......@@ -342,6 +343,9 @@ int XLALSimInspiralNRWaveformGetSpinsFromHDF5File(
const char *NRDataFile /**< Location of NR HDF file */
);
/* The following XLALSimInspiralNRWaveformGetHplusHcross() generates polarizations
* reading directly the NR files and does not return l,m modes.
*/
int XLALSimInspiralNRWaveformGetHplusHcross(
REAL8TimeSeries **hplus, /**< OUTPUT h_+ vector */
REAL8TimeSeries **hcross, /**< OUTPUT h_x vector */
......@@ -363,6 +367,25 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
LALValue* ModeArray /**< Container for the ell and m modes to generate. To generate all available modes pass NULL */
);
/* The following XLALSimInspiralNRWaveformGetHlms() reads NR file to output l,m modes.
*/
INT4 XLALSimInspiralNRWaveformGetHlms(SphHarmTimeSeries **hlms, /**< OUTPUT */
REAL8 deltaT, /**< sampling interval (s) */
REAL8 m1, /**< mass of companion 1 (kg) */
REAL8 m2, /**< mass of companion 2 (kg) */
REAL8 r, /**< distance of source (m) */
REAL8 fStart, /**< start GW frequency (Hz) */
REAL8 fRef, /**< reference GW frequency (Hz) */
REAL8 s1x, /**< initial value of S1x */
REAL8 s1y, /**< initial value of S1y */
REAL8 s1z, /**< initial value of S1z */
REAL8 s2x, /**< initial value of S2x */
REAL8 s2y, /**< initial value of S2y */
REAL8 s2z, /**< initial value of S2z */
const char *NRDataFile, /**< Location of NR HDF file */
LALValue* ModeArray /**< Container for the ell and m modes to generate. To generate all available modes pass NULL */
);
/* in module LALSimIMRPrecessingNRSur.c */
int XLALSimInspiralPrecessingNRSurPolarizations(
......@@ -399,7 +422,7 @@ SphHarmTimeSeries *XLALSimInspiralPrecessingNRSurModes(
REAL8 S2z, /**< z-component of the dimensionless spin of object 2 */
REAL8 fMin, /**< start GW frequency (Hz) */
REAL8 fRef, /**< reference GW frequency (Hz) */
REAL8 distnace, /**< distance of source (m) */
REAL8 distance, /**< distance of source (m) */
LALDict* LALparams, /**< Dict with extra parameters */
Approximant approximant /**< approximant (NRSur7dq2 or NRSur7dq4) */
);
......@@ -434,7 +457,6 @@ int XLALPrecessingNRSurDynamics(
Approximant approximant /**< approximant (NRSur7dq2 or NRSur7dq4). */
);
/* in module LALSimNRTunedTides.c */
double XLALSimNRTunedTidesComputeKappa2T(
REAL8 m1_SI, /**< Mass of companion 1 (kg) */
......
......@@ -1748,7 +1748,6 @@ XLALSimIMREOBNRv2AllModes(
* SWSH modes in a SphHarmTimeSeries struct.
*/
SphHarmTimeSeries *XLALSimIMREOBNRv2Modes(
const REAL8 phiRef, /**< Orbital phase at coalescence (radians) */
const REAL8 deltaT, /**< Sampling interval (s) */
const REAL8 m1, /**< First component mass (kg) */
const REAL8 m2, /**< Second component mass (kg) */
......@@ -1757,7 +1756,7 @@ SphHarmTimeSeries *XLALSimIMREOBNRv2Modes(
)
{
SphHarmTimeSeries *hlms = NULL;
if ( XLALSimIMREOBNRv2Generator(NULL, NULL, &hlms, phiRef, deltaT, m1, m2,
if ( XLALSimIMREOBNRv2Generator(NULL, NULL, &hlms, 0., deltaT, m1, m2,
fLower, distance, 0., 1) == XLAL_FAILURE )
{
XLAL_ERROR_NULL( XLAL_EFUNC );
......
/*
* Copyright (C) 2015 Ian Harry, Patricia Schmidt
* Copyright (C) 2015 Ian Harry, Patricia Schmidt, Riccardo Sturani
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -52,6 +52,7 @@
#include <lal/LALSimIMR.h>
#include <lal/LALConfig.h>
#include <lal/SphericalHarmonics.h>
#include <lal/LALSimInspiralPrecess.h>
#include <lal/H5FileIO.h>
......@@ -624,14 +625,14 @@ int XLALSimInspiralNRWaveformGetSpinsFromHDF5File(
/* Everything needs to be declared as unused in case HDF is not enabled. */
int XLALSimInspiralNRWaveformGetHplusHcross(
UNUSED REAL8TimeSeries **hplus, /**< Output h_+ vector */
UNUSED REAL8TimeSeries **hcross, /**< Output h_x vector */
UNUSED REAL8 phiRef, /**< orbital phase at reference pt. */
UNUSED REAL8 inclination, /**< inclination angle */
UNUSED static INT4 XLALSimIMRNRWaveformGetModes(
UNUSED SphHarmTimeSeries **Hlms, /**< OUTPUT */
UNUSED LIGOTimeGPS *epoch, /**< OUTPUT */
UNUSED UINT4 *length, /**< OUTPUT */
UNUSED REAL8 deltaT, /**< sampling interval (s) */
UNUSED REAL8 m1, /**< mass of companion 1 (kg) */
UNUSED REAL8 m2, /**< mass of companion 2 (kg) */
UNUSED REAL8 m1, /**< mass of companion 1 (solar units) */
UNUSED REAL8 m2, /**< mass of companion 2 (solar units) */
UNUSED REAL8 r, /**< distance of source (m) */
UNUSED REAL8 fStart, /**< start GW frequency (Hz) */
UNUSED REAL8 fRef, /**< reference GW frequency (Hz) */
......@@ -641,13 +642,13 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
UNUSED REAL8 s2x, /**< initial value of S2x */
UNUSED REAL8 s2y, /**< initial value of S2y */
UNUSED REAL8 s2z, /**< initial value of S2z */
UNUSED const char *NRDataFile, /**< Location of NR HDF file */
UNUSED LALH5File* file, /**< pointer to location of NR HDF file */
UNUSED LALValue* ModeArray /**< Container for the ell and m modes to generate. To generate all available modes pass NULL */
)
{
#ifndef LAL_HDF5_ENABLED
XLAL_ERROR(XLAL_EFAILED, "HDF5 support not enabled");
#else
#ifndef LAL_HDF5_ENABLED
XLAL_ERROR_NULL(XLAL_FAILURE, "HDF5 support not enabled");
#else
/* Declarations */
UINT4 curr_idx, nr_file_format;
INT4 model, modem;
......@@ -655,34 +656,20 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
REAL8 nrEta;
REAL8 S1x, S1y, S1z, S2x, S2y, S2z;
REAL8 Mflower, time_start_M, time_start_s, time_end_M, time_end_s;
REAL8 est_start_time, curr_h_real, curr_h_imag;
REAL8 theta, psi, calpha, salpha;
REAL8 est_start_time;
REAL8 distance_scale_fac;
COMPLEX16 curr_ylm;
REAL8TimeSeries *hplus_corr;
REAL8TimeSeries *hcross_corr;
COMPLEX16TimeSeries *hlm;
SphHarmTimeSeries *hlms_tmp=NULL;
/* These keys follow a strict formulation and cannot be longer than 11
* characters */
char amp_key[30];
char phase_key[30];
gsl_vector *tmpVector=NULL;
LALH5File *file, *group;
LALH5File *group;
LIGOTimeGPS tmpEpoch = LIGOTIMEGPSZERO;
REAL8Vector *curr_amp, *curr_phase;
/* Use solar masses for units. NR files will use
* solar masses as well, so easier for that conversion
*/
m1 = m1 / LAL_MSUN_SI;
m2 = m2 / LAL_MSUN_SI;
file = XLALH5FileOpen(NRDataFile, "r");
if (file == NULL)
{
XLAL_ERROR(XLAL_EIO, "NR SIMULATION DATA FILE %s NOT FOUND.\n", NRDataFile);
}
/* Sanity checks on physical parameters passed to waveform
* generator to guarantee consistency with NR data file.
*/
......@@ -801,31 +788,12 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
}
array_length = (UINT4)(ceil( (time_end_s - time_start_s) / deltaT));
/* Compute correct angles for hplus and hcross following LAL convention. */
theta = psi = calpha = salpha = 0.;
XLALSimInspiralNRWaveformGetRotationAnglesFromH5File(&theta, &psi, &calpha,
&salpha, file, inclination, phiRef, fRef*(m1+m2));
*length=array_length;
/* Create the return time series, use arbitrary epoch here. We set this
* properly later. */
XLALGPSAdd(&tmpEpoch, time_start_s);
*hplus = XLALCreateREAL8TimeSeries("H_PLUS", &tmpEpoch, 0.0, deltaT,
&lalStrainUnit, array_length );
*hcross = XLALCreateREAL8TimeSeries("H_CROSS", &tmpEpoch, 0.0, deltaT,
&lalStrainUnit, array_length );
hplus_corr = XLALCreateREAL8TimeSeries("H_PLUS", &tmpEpoch, 0.0, deltaT,
&lalStrainUnit, array_length );
hcross_corr = XLALCreateREAL8TimeSeries("H_CROSS", &tmpEpoch, 0.0, deltaT,
&lalStrainUnit, array_length );
for (curr_idx = 0; curr_idx < array_length; curr_idx++)
{
hplus_corr->data->data[curr_idx] = 0.0;
hcross_corr->data->data[curr_idx] = 0.0;
}
*epoch=tmpEpoch;
/* Create the distance scale factor */
distance_scale_fac = (m1 + m2) * LAL_MRSUN_SI / r;
......@@ -833,6 +801,7 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
/* Generate the waveform */
/* NOTE: We assume that for a given ell mode, all m modes are present */
INT4 NRLmax;
XLALH5FileQueryScalarAttributeValue(&NRLmax, file, "Lmax");
INT4 modearray_needs_destroying=0;
......@@ -846,6 +815,8 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
}
}
/* else Use the ModeArray given */
hlm=XLALCreateCOMPLEX16TimeSeries("hlm",&tmpEpoch,0.0,deltaT,&lalStrainUnit,array_length);
memset(hlm->data->data, 0, array_length * sizeof(COMPLEX16));
for (model=2; model < (NRLmax + 1) ; model++)
{
......@@ -861,7 +832,6 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
}
XLAL_PRINT_INFO("generating model = %i modem = %i\n", model, modem);
snprintf(amp_key, sizeof(amp_key), "amp_l%d_m%d", model, modem);
snprintf(phase_key, sizeof(phase_key), "phase_l%d_m%d", model, modem);
......@@ -881,31 +851,132 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
XLALSimInspiralNRWaveformGetDataFromHDF5File(&curr_phase, file, (m1 + m2),
time_start_s, array_length, deltaT, phase_key);
curr_ylm = XLALSpinWeightedSphericalHarmonic(theta, psi, -2,
model, modem);
for (curr_idx = 0; curr_idx < array_length; curr_idx++)
{
curr_h_real = curr_amp->data[curr_idx]
* cos(curr_phase->data[curr_idx]) * distance_scale_fac;
curr_h_imag = curr_amp->data[curr_idx]
* sin(curr_phase->data[curr_idx]) * distance_scale_fac;
hlm->data->data[curr_idx]= (curr_amp->data[curr_idx]*cos(curr_phase->data[curr_idx]) + I*curr_amp->data[curr_idx]*sin(curr_phase->data[curr_idx]) ) * distance_scale_fac;
}
/* Note that the hlm built here do not respect the LAL convention
* the function XLALSimIMRNRWaveformGetHlm will perform the pi/2 rotation
* necessary to bring the in the right frame.
* See https://dcc.ligo.org/LIGO-T1900080 for details.
*/
hlms_tmp=XLALSphHarmTimeSeriesAddMode(hlms_tmp,hlm,model,modem);
XLALDestroyREAL8Vector(curr_amp);
XLALDestroyREAL8Vector(curr_phase);
}
}
XLALDestroyCOMPLEX16TimeSeries(hlm);
if (modearray_needs_destroying)
XLALDestroyValue(ModeArray);
hplus_corr->data->data[curr_idx] = hplus_corr->data->data[curr_idx]
+ curr_h_real * creal(curr_ylm) - curr_h_imag * cimag(curr_ylm);
*Hlms=hlms_tmp;
hcross_corr->data->data[curr_idx] = hcross_corr->data->data[curr_idx]
- curr_h_real * cimag(curr_ylm) - curr_h_imag * creal(curr_ylm);
return XLAL_SUCCESS;
#endif
}
}
/* Everything needs to be declared as unused in case HDF is not enabled. */
INT4 XLALSimInspiralNRWaveformGetHplusHcross(
UNUSED REAL8TimeSeries **hplus, /**< Output h_+ vector */
UNUSED REAL8TimeSeries **hcross, /**< Output h_x vector */
UNUSED REAL8 phiRef, /**< orbital phase at reference pt. */
UNUSED REAL8 inclination, /**< inclination angle */
UNUSED REAL8 deltaT, /**< sampling interval (s) */
UNUSED REAL8 m1_SI, /**< mass of companion 1 (kg) */
UNUSED REAL8 m2_SI, /**< mass of companion 2 (kg) */
UNUSED REAL8 r, /**< distance of source (m) */
UNUSED REAL8 fStart, /**< start GW frequency (Hz) */
UNUSED REAL8 fRef, /**< reference GW frequency (Hz) */
UNUSED REAL8 s1x, /**< initial value of S1x */
UNUSED REAL8 s1y, /**< initial value of S1y */
UNUSED REAL8 s1z, /**< initial value of S1z */
UNUSED REAL8 s2x, /**< initial value of S2x */
UNUSED REAL8 s2y, /**< initial value of S2y */
UNUSED REAL8 s2z, /**< initial value of S2z */
UNUSED const char *NRDataFile, /**< Location of NR HDF file */
UNUSED LALValue* ModeArray /**< Container for the ell and m modes to generate. To generate all available modes pass NULL */
)
{
#ifndef LAL_HDF5_ENABLED
XLAL_ERROR(XLAL_EFAILED, "HDF5 support not enabled");
#else
/* Declarations */
UINT4 curr_idx, array_length;
INT4 model, modem, NRLmax;
REAL8 m1,m2;
REAL8 theta, psi, calpha, salpha;
COMPLEX16 curr_ylm, tmp;
COMPLEX16TimeSeries *curr_hlm=NULL;
REAL8TimeSeries *hplus_corr;
REAL8TimeSeries *hcross_corr;
XLALDestroyREAL8Vector(curr_amp);
XLALDestroyREAL8Vector(curr_phase);
/* These keys follow a strict formulation and cannot be longer than 11
* characters */
LALH5File *file=NULL;
LIGOTimeGPS tmpEpoch = LIGOTIMEGPSZERO;
SphHarmTimeSeries *tmp_Hlms=NULL;
}
/* Use solar masses for units. NR files will use
* solar masses as well, so easier for that conversion
*/
m1 = m1_SI / LAL_MSUN_SI;
m2 = m2_SI / LAL_MSUN_SI;
file = XLALH5FileOpen(NRDataFile, "r");
if (file == NULL)
{
XLAL_ERROR(XLAL_EIO, "NR SIMULATION DATA FILE %s NOT FOUND.\n", NRDataFile);
}
INT4 err_code = XLALSimIMRNRWaveformGetModes(&tmp_Hlms,&tmpEpoch,&array_length, \
deltaT, m1, m2, r, fStart, fRef, \
s1x,s1y,s1z,s2x,s2y,s2z, \
file, ModeArray);
if (err_code!=XLAL_SUCCESS)
XLAL_ERROR(XLAL_FAILURE);
/* Compute correct angles for hplus and hcross following LAL convention. */
theta = psi = calpha = salpha = 0.;
XLALSimInspiralNRWaveformGetRotationAnglesFromH5File(&theta, &psi, &calpha,
&salpha, file, inclination, phiRef, fRef*(m1+m2));
XLALH5FileClose(file);
*hplus = XLALCreateREAL8TimeSeries("H_PLUS", &tmpEpoch, 0.0, deltaT,
&lalStrainUnit, array_length );
*hcross = XLALCreateREAL8TimeSeries("H_CROSS", &tmpEpoch, 0.0, deltaT,
&lalStrainUnit, array_length );
hplus_corr = XLALCreateREAL8TimeSeries("H_PLUS", &tmpEpoch, 0.0, deltaT,
&lalStrainUnit, array_length );
hcross_corr = XLALCreateREAL8TimeSeries("H_CROSS", &tmpEpoch, 0.0, deltaT,
&lalStrainUnit, array_length );
for (curr_idx = 0; curr_idx < array_length; curr_idx++)
{
hplus_corr->data->data[curr_idx] = 0.0;
hcross_corr->data->data[curr_idx] = 0.0;
}
NRLmax=XLALSphHarmTimeSeriesGetMaxL(tmp_Hlms);
for (model=2; model < (NRLmax + 1) ; model++)
{
for (modem=-model; modem < (model+1); modem++)
{
curr_ylm = XLALSpinWeightedSphericalHarmonic(theta, psi, -2, model, modem);
curr_hlm = XLALSphHarmTimeSeriesGetMode(tmp_Hlms, model, modem);
if (curr_hlm) {
for (curr_idx = 0; curr_idx < array_length; curr_idx++)
{
tmp=curr_hlm->data->data[curr_idx]*curr_ylm;
hplus_corr->data->data[curr_idx] += creal(tmp);
hcross_corr->data->data[curr_idx] -= cimag(tmp);
}
}
}
}
XLALDestroySphHarmTimeSeries(tmp_Hlms);
/* Correct for the "alpha" angle as given in T1500606 or arxiv:1703.01076
* to translate from the NR wave frame to LAL wave-frame.
*
......@@ -925,10 +996,108 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
XLALDestroyREAL8TimeSeries(hplus_corr);
XLALDestroyREAL8TimeSeries(hcross_corr);
return XLAL_SUCCESS;
#endif
}
/* Everything needs to be declared as unused in case HDF is not enabled. */
INT4 XLALSimInspiralNRWaveformGetHlms(UNUSED SphHarmTimeSeries **Hlms, /**< OUTPUT */
UNUSED REAL8 deltaT, /**< sampling interval (s) */
UNUSED REAL8 m1_SI, /**< mass of companion 1 (kg) */
UNUSED REAL8 m2_SI, /**< mass of companion 2 (kg) */
UNUSED REAL8 r, /**< distance of source (m) */
UNUSED REAL8 fStart, /**< start GW frequency (Hz) */
UNUSED REAL8 fRef, /**< reference GW frequency (Hz) */
UNUSED REAL8 s1x, /**< initial value of S1x */
UNUSED REAL8 s1y, /**< initial value of S1y */
UNUSED REAL8 s1z, /**< initial value of S1z */
UNUSED REAL8 s2x, /**< initial value of S2x */
UNUSED REAL8 s2y, /**< initial value of S2y */
UNUSED REAL8 s2z, /**< initial value of S2z */
UNUSED const char *NRDataFile, /**< Location of NR HDF file */
UNUSED LALValue* ModeArray /**< Container for the ell and m modes to generate. To generate all available modes pass NULL */
)
{
#ifndef LAL_HDF5_ENABLED
XLAL_ERROR_NULL(XLAL_FAILURE, "HDF5 support not enabled");
#else
/* Declarations */
UINT4 curr_idx, array_length;
INT4 model, modem, NRLmax;
REAL8 m1,m2;
REAL8 theta, psi, calpha, salpha;
COMPLEX16TimeSeries *tmp_hlm=NULL;
SphHarmTimeSeries *tmp_Hlms=NULL;
LALH5File *file;
LIGOTimeGPS tmpEpoch = LIGOTIMEGPSZERO;
/* Use solar masses for units. NR files will use
* solar masses as well, so easier for that conversion
*/
m1 = m1_SI / LAL_MSUN_SI;
m2 = m2_SI / LAL_MSUN_SI;
file = XLALH5FileOpen(NRDataFile, "r");
if (file == NULL)
{
XLAL_ERROR(XLAL_EIO, "NR SIMULATION DATA FILE %s NOT FOUND.\n", NRDataFile);
}
INT4 err_code = XLALSimIMRNRWaveformGetModes(&tmp_Hlms,&tmpEpoch,&array_length, \
deltaT, m1, m2, r, fStart, fRef, \
s1x,s1y,s1z,s2x,s2y,s2z, \
file, ModeArray);
if (err_code!=XLAL_SUCCESS)
XLAL_ERROR(XLAL_FAILURE);
/* Compute correct angles for hplus and hcross following LAL convention. */
psi = calpha = salpha = 0.;
XLALSimInspiralNRWaveformGetRotationAnglesFromH5File(&theta, &psi, &calpha,
&salpha, file, 0., 0., fRef*(m1+m2));
XLALH5FileClose(file);
if (modearray_needs_destroying)
XLALDestroyValue(ModeArray);
NRLmax= XLALSphHarmTimeSeriesGetMaxL(tmp_Hlms);
COMPLEX16 facm=-1;
for (model=2; model < (NRLmax + 1) ; model++)
{
for (modem=-model; modem < (model+1); modem++)
{
tmp_hlm=XLALSphHarmTimeSeriesGetMode(tmp_Hlms,model,modem);
if (tmp_hlm) {
for (curr_idx = 0; curr_idx < array_length; curr_idx++)
tmp_hlm->data->data[curr_idx]*=facm;
*Hlms=XLALSphHarmTimeSeriesAddMode(*Hlms,tmp_hlm,model,modem);
}
facm*=I;
}
facm=1./facm;
}
XLALDestroySphHarmTimeSeries(tmp_Hlms);
/* Correct for the "alpha" angle as given in T1500606 or arxiv:1703.01076
* to translate from the NR wave frame to LAL wave-frame.
*/
REAL8TimeSeries *alphats = XLALCreateREAL8TimeSeries("alpha", &tmpEpoch, 0.0, deltaT, &lalDimensionlessUnit, array_length );
REAL8TimeSeries *thetats = XLALCreateREAL8TimeSeries("theta", &tmpEpoch, 0.0, deltaT, &lalDimensionlessUnit, array_length );
REAL8TimeSeries *psits = XLALCreateREAL8TimeSeries("tpsi", &tmpEpoch, 0.0, deltaT, &lalDimensionlessUnit, array_length );
REAL8 alpha=atan2(salpha,calpha);
for (curr_idx = 0; curr_idx < array_length; curr_idx++) {
alphats->data->data[curr_idx] = alpha;
thetats->data->data[curr_idx] =-theta;
psits->data->data[curr_idx] =-psi;
}
err_code=XLALSimInspiralPrecessionRotateModes(*Hlms,alphats,thetats,psits);
if (err_code!=XLAL_SUCCESS)
XLAL_ERROR(XLAL_EINVAL);
XLALDestroyREAL8TimeSeries(alphats);
XLALDestroyREAL8TimeSeries(thetats);
XLALDestroyREAL8TimeSeries(psits);
return XLAL_SUCCESS;
#endif
#endif
}
......@@ -40,6 +40,8 @@
#include <lal/BandPassTimeSeries.h>
#include <lal/Units.h>
#include <lal/LALSimBlackHoleRingdown.h>
#include <lal/LALSimInspiralPrecess.h>
#include <lal/LALSimInspiralWaveformParams.h>
#include "LALSimInspiralPNCoefficients.c"
#include "check_series_macros.h"
......@@ -596,17 +598,10 @@ int XLALSimInspiralChooseTDWaveform(
E1y = 1.;
E1z = 0.;
polariz+=LAL_PI/2.;
/* Maximum PN amplitude order for precessing waveforms is
* MAX_PRECESSING_AMP_PN_ORDER */
amplitudeO = amplitudeO <= MAX_PRECESSING_AMP_PN_ORDER ?
amplitudeO : MAX_PRECESSING_AMP_PN_ORDER;
/* Call the waveform driver routine */
ret = XLALSimInspiralSpinTaylorT5(hplus, hcross, phiRef, v0, deltaT,
ret = XLALSimInspiralSpinTaylorT5(hplus, hcross, phiRef, deltaT,
m1, m2, f_min, f_ref, distance, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, lambda1, lambda2,
quadparam1, quadparam2,
LALparams,
phaseO, amplitudeO);
LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, LALparams);
break;
// need to make a consistent choice for SpinTaylorT4 and PSpinInspiralRD waveform inputs
......@@ -627,16 +622,10 @@ int XLALSimInspiralChooseTDWaveform(
E1y = 1.;
E1z = 0.;
polariz+=LAL_PI/2.;
/* Maximum PN amplitude order for precessing waveforms is
* MAX_PRECESSING_AMP_PN_ORDER */
amplitudeO = amplitudeO <= MAX_PRECESSING_AMP_PN_ORDER ?
amplitudeO : MAX_PRECESSING_AMP_PN_ORDER;
/* Call the waveform driver routine */
ret = XLALSimInspiralSpinTaylorT4(hplus, hcross, phiRef, v0, deltaT,
ret = XLALSimInspiralSpinTaylorT4(hplus, hcross, phiRef, deltaT,
m1, m2, f_min, f_ref, distance, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, lambda1, lambda2,
quadparam1, quadparam2, LALparams,
phaseO, amplitudeO);
LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, LALparams);
break;
case SpinTaylorT1:
......@@ -652,15 +641,10 @@ int XLALSimInspiralChooseTDWaveform(
E1y = 1.;
E1z = 0.;
polariz+=LAL_PI/2.;
/* Maximum PN amplitude order for precessing waveforms is
* MAX_PRECESSING_AMP_PN_ORDER */
amplitudeO = amplitudeO <= MAX_PRECESSING_AMP_PN_ORDER ?
amplitudeO : MAX_PRECESSING_AMP_PN_ORDER;
/* Call the waveform driver routine */
ret = XLALSimInspiralSpinTaylorT1(hplus, hcross, phiRef, v0, deltaT,
m1, m2, f_min, f_ref, distance, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, lambda1, lambda2,
quadparam1, quadparam2, LALparams,
phaseO, amplitudeO);
ret = XLALSimInspiralSpinTaylorT1(hplus, hcross, phiRef, deltaT,
m1, m2, f_min, f_ref, distance, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z,
LNhatx, LNhaty, LNhatz, E1x, E1y, E1z, LALparams);
break;
case SpinDominatedWf:
......@@ -1024,7 +1008,6 @@ int XLALSimInspiralChooseTDWaveform(
S1z, S2z, LALparams);
break;
default:
XLALPrintError("TD version of approximant not implemented in lalsimulation\n");
XLAL_ERROR(XLAL_EINVAL);
......@@ -2534,12 +2517,18 @@ int XLALSimInspiralChooseWaveform(
* Interface to compute a set of -2 spin-weighted spherical harmonic modes
* for a binary inspiral for a given waveform approximant.
* PN Approximants (TaylorT1 - T4), EOBNRv2 (EOBNRv2HM), NRSur7dq2, NRSur7dq4
* and NRHybSur3dq8 are implemented.
* NRHybSur3dq8 and spin-precessing SpintaylorT1, T5, T4 are implemented.
*
* The EOBNRv2 model returns the (2,2), (2,1), (3,3), (4,4), and (5,5) modes.
* Note that the inclination parameter is not passed to create hlm modes,
* hence to recover the correct h+,x one has to combine the hlm modes with
* Euler angles alpha=0, iota=inclination, psi=0,Pi/2 (according to the approximat) i.e.
* (h+ + I hx) (psi,iota,alpha)= e^(-2Ialpha) Sum_{l,m} Y_lm(-iota,-psi) h_lm,
* or equivalently rotate h_lm -> h'_lm=DWigner(psi,iota,alpha) h_lm
* and then obtain
* (h+ + I hx) = Sum_{l,m} Y_lm(0,0) h'_lm,
*/
SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
REAL8 phiRef, /**< reference orbital phase (rad) */
REAL8 deltaT, /**< sampling interval (s) */
REAL8 m1, /**< mass of companion 1 (kg) */
REAL8 m2, /**< mass of companion 2 (kg) */
......@@ -2559,6 +2548,7 @@ SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
{
REAL8 v0 = 1.;
SphHarmTimeSeries *hlm = NULL;
INT4 errCode=0;
/* SEOBNR flag for precessing model version. 3 for SEOBNRv3, 300 for SEOBNRv3_opt, 401 for SEOBNRv4P, 402 for SEOBNRv4PHM */
UINT4 PrecEOBversion;
......@@ -2599,6 +2589,7 @@ SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
REAL8 lambda2 = XLALSimInspiralWaveformParamsLookupTidalLambda2(LALpars);
int amplitudeO = XLALSimInspiralWaveformParamsLookupPNAmplitudeOrder(LALpars);
int phaseO =XLALSimInspiralWaveformParamsLookupPNPhaseOrder(LALpars);
UINT4 l;
switch (approximant)
{
......@@ -2611,7 +2602,7 @@ SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALpars) )
ABORT_NONDEFAULT_MODES_CHOICE_NULL(LALpars);
/* Call the waveform driver routine */
hlm = XLALSimInspiralTaylorT1PNModes(phiRef, v0,
hlm = XLALSimInspiralTaylorT1PNModes(v0,
deltaT, m1, m2, f_min, f_ref, r, lambda1, lambda2,
XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALpars), amplitudeO,
phaseO, lmax);
......@@ -2625,7 +2616,7 @@ SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALpars) )
ABORT_NONDEFAULT_MODES_CHOICE_NULL(LALpars);
/* Call the waveform driver routine */
hlm = XLALSimInspiralTaylorT2PNModes(phiRef, v0,
hlm = XLALSimInspiralTaylorT2PNModes(v0,
deltaT, m1, m2, f_min, f_ref, r, lambda1, lambda2,
XLALSimInspiralWaveformParamsLookupPNTidalOrder(LALpars), amplitudeO,
phaseO, lmax);
......@@ -2639,7 +2630,7 @@ SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALpars) )
ABORT_NONDEFAULT_MODES_CHOICE_NULL(LALpars);
/* Call the waveform driver routine */
hlm = XLALSimInspiralTaylorT3PNModes(phiRef, v0,
hlm = XLALSimInspiralTaylorT3PNModes(v0,
deltaT, m1, m2, f_min, f_ref, r, lambda1, lambda2,
0, amplitudeO, phaseO, lmax);
break;
......@@ -2652,7 +2643,7 @@ SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALpars) )
ABORT_NONDEFAULT_MODES_CHOICE_NULL(LALpars);
/* Call the waveform driver routine */
hlm = XLALSimInspiralTaylorT4PNModes(phiRef, v0,
hlm = XLALSimInspiralTaylorT4PNModes(v0,
deltaT, m1, m2, f_min, f_ref, r, lambda1, lambda2,
0, amplitudeO,
phaseO, lmax);
......@@ -2667,9 +2658,9 @@ SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
if( !XLALSimInspiralWaveformParamsModesChoiceIsDefault(LALpars) )
ABORT_NONDEFAULT_MODES_CHOICE_NULL(LALpars);
/* Call the waveform driver routine */
hlm = XLALSimIMREOBNRv2Modes(phiRef, deltaT, m1, m2, f_min, r);
hlm = XLALSimIMREOBNRv2Modes(deltaT, m1, m2, f_min, r);
// EOB driver only outputs modes with m>0, add m<0 modes by symmetry
size_t l, j;
size_t j;
int m;
for( l=2; l <= XLALSphHarmTimeSeriesGetMaxL( hlm ); l++ ) {
for( m=-l; m<0; m++){
......@@ -2764,11 +2755,66 @@ SphHarmTimeSeries *XLALSimInspiralChooseTDModes(
hlm = XLALSimIMRSpinPrecEOBModes(deltaT, m1, m2, f_min, r, spin1, spin2, PrecEOBversion,LALpars);
break;
case SpinTaylorT1:
case SpinTaylorT2: