22#include <lal/LALInspiral.h>
23#include <lal/FindRoot.h>
24#include <lal/SeqFactories.h>
25#include <lal/NRWaveInject.h>
26#include <lal/RealFFT.h>
27#include <lal/TimeFreqFFT.h>
28#include <lal/TimeSeries.h>
29#include <lal/SeqFactories.h>
30#include <lal/VectorOps.h>
31#include <lal/BBHPhenomCoeffs.h>
32#include <lal/AVFactories.h>
37#define UNUSED __attribute__ ((unused))
42typedef struct tagBBHPhenomParams{
110 if (!signalvec || !(signalvec->
data) || !
params) {
114 if ((signalvec->
length < 2)
138 if (!signalvec || !(signalvec->
data) || !
params) {
142 if ((signalvec->
length < 2)
145 || (
params->spin2[0] != 0) || (
params->spin2[1] != 0)) {
164 if (!signalvec1 || !signalvec2 ||
165 !(signalvec1->
data) || !(signalvec2->
data)) {
191 if (!signalvec1 || !signalvec2
192 || !(signalvec1->
data) || !(signalvec2->
data)
220 if (!signalvec1 || !(signalvec1->
data)) {
224 if (signalvec1->
length <= 2) {
244 if (!signalvec1 || !(signalvec1->
data)
245 || !signalvec2 || !(signalvec2->
data)
250 if ((signalvec1->
length <= 2) || (signalvec2->
length <= 2)) {
271 REAL8 dt, cosI, fLower, peakAmp, fCut, fRes, f, totalMass, softWin, z1, z2;
272 REAL8 fLowerOrig, eta, tau0, winFLo, sigLo, sigHi, tF0, expectedAmplRatio;
273 REAL8 phaseShift, sig1, sig2, startPhaseOrig, phiC;
275 UINT4 i, j, k,
l,
n, peakAmpIdx, sigLength, iLower;
276 REAL4Vector *signalFD1 = NULL, *signalFD2 = NULL, *signalTD1 = NULL, *signalTD2 = NULL, *
a=NULL, *fVec=NULL;
278 REAL4FFTPlan *revPlan = NULL;
290 || (
params->spin2[0] != 0) || (
params->spin2[1] != 0)) {
299 switch (
params->approximant) {
302 if ((
params->spin1[2] != 0) || (
params->spin2[2] != 0)) {
324 fLowerOrig =
params->fLower;
325 startPhaseOrig =
params->startPhase;
332 fCut = (1.025)*phenParams.
fCut;
335 if (fLower > fLowerOrig) fLower = fLowerOrig;
336 if (fLower < 0.5) fLower = 0.5;
337 if (fCut >
params->tSampling/2.-100.) fCut =
params->tSampling/2.-100.;
346 n = pow(2., ceil(log2(tau0*
params->tSampling)));
351 if (!signalFD1 || !signalFD2 || !signalTD1 || !signalTD2) {
360 switch (
params->approximant) {
385 winFLo = (fLowerOrig + fLower)/2.;
390 signalFD1->
data[0] = 0.;
391 signalFD2->data[0] = 0.;
392 for (k = 1; k <=
n/2; k++) {
394 softWin = (1+tanh((4.*(f-winFLo)/sigLo)))*(1-tanh(4.*(f-fCut)/sigHi))/4.;
395 signalFD1->
data[k] *= softWin;
396 signalFD1->
data[
n-k] *= softWin;
397 signalFD2->data[k] *= softWin;
398 signalFD2->data[
n-k] *= softWin;
418 for (
i = 0;
i <
n;
i++) {
419 signalTD1->data[
i] *=
params->tSampling/
n;
420 signalTD2->data[
i] *=
params->tSampling/
n;
426 for (
i=0;
i< windowLength;
i++){
427 signalTD1->data[
n-
i-1] *=
i/windowLength;
428 signalTD2->data[
n-
i-1] *=
i/windowLength;
442 a->data[
i] = sqrt(pow(signalTD1->data[
i],2.) + pow(signalTD2->data[
i],2.));
443 phi->
data[
i] = -atan2(signalTD2->data[
i], signalTD1->data[
i]);
446 if (
a->data[
i] > peakAmp) {
447 peakAmp =
a->data[
i];
466 for (
i=0;
i< fVec->length;
i++) {
467 if (
a->data[
i] < expectedAmplRatio*peakAmp) {
470 if ((iLower == 0) && (fVec->data[
i] >= fLowerOrig)) iLower =
i;
477 params->fLower = fLowerOrig;
479 params->startPhase = startPhaseOrig;
483 phaseShift = phi->
data[peakAmpIdx] +
params->startPhase;
484 phiC = phi->
data[peakAmpIdx] -
params->startPhase;
485 for (
i=iLower;
i < fVec->length;
i++) {
486 sig1 = signalTD1->data[
i]*cos(phaseShift) - signalTD2->data[
i]*sin(phaseShift);
487 sig2 = signalTD1->data[
i]*sin(phaseShift) + signalTD2->data[
i]*cos(phaseShift);
488 signalTD1->data[
i] = sig1;
489 signalTD2->data[
i] = sig2;
490 phi->
data[
i] -= phiC;
498 if (signalvec1) sigLength = signalvec1->
length;
499 else if (signalvec2) sigLength = signalvec2->
length;
500 else if (freqVec) sigLength = freqVec->
length;
501 else if (phiVec) sigLength = phiVec->
length;
502 else if (aVec) sigLength = aVec->
length / 2;
503 else if (h) sigLength = h->
length / 2;
506 z1 = -0.5*(1. + cosI*cosI);
509 for (
i=iLower;
i < fVec->length;
i++) {
514 if (signalvec1) signalvec1->
data[j] = signalTD1->data[
i];
515 if (signalvec2) signalvec2->
data[j] = signalTD2->data[
i];
516 if (freqVec) freqVec->
data[j] = fVec->data[
i];
517 if (phiVec) phiVec->
data[j] = phi->
data[
i];
519 aVec->
data[k] =
a->data[
i];
523 h->
data[k] = z1 * signalTD1->data[
i];
524 h->
data[
l] = z2 * signalTD2->data[
i];
563 if (!
params || !waveform) {
574 if (waveform->
a || waveform->
h || waveform->
f
575 || waveform->
phi || waveform->
shift) {
582 snprintf(message, 256,
"WARNING: Amp Order has been reset to %d",
params->ampOrder);
597 memset(
a->data, 0, 2 * paramsInit.
nbins *
sizeof(
REAL4));
600 if (!h || !
a || !ff || !phi) {
616 if (phi->
data[
i] != 0.0)
break;
625 snprintf(message, 256,
"fFinal = %f",
params->fFinal);
627 s = 0.5 * (phi->
data[count-1] - phi->
data[0]);
628 snprintf(message, 256,
"cycles = %f",
s/3.14159);
631 snprintf(message, 256,
"The waveform has only %f cycles; we don't "
632 "keep waveform with less than 2 cycles.", (
double)
s/ (
double)
LAL_PI );
638 phiC = phi->
data[count-1] ;
639 for (
i=0;
i<count;
i++) {
646 if (!(waveform->
h) || !(waveform->
a)) {
657 if (!(waveform->
h->
data) || !(waveform->
a->
data) || !(waveform->
f) || !(waveform->
phi)) {
678 waveform->
psi = ppnParams->
psi;
684 ppnParams->
tc = (double)(count-1) /
params->tSampling;
685 ppnParams->
length = count;
725 REAL8 totalMass, piM, eta, fMerg_a, fMerg_b, fMerg_c, fRing_a, fRing_b, etap2;
726 REAL8 fRing_c, sigma_a, sigma_b, sigma_c, fCut_a, fCut_b, fCut_c;
727 REAL8 psi0_a, psi0_b, psi0_c, psi2_a, psi2_b, psi2_c, psi3_a, psi3_b, psi3_c;
728 REAL8 psi4_a, psi4_b, psi4_c, psi6_a, psi6_b, psi6_c, psi7_a, psi7_b, psi7_c;
786 phenParams->
fCut = (fCut_a*etap2 + fCut_b*eta + fCut_c)/piM;
787 phenParams->
fMerger = (fMerg_a*etap2 + fMerg_b*eta + fMerg_c)/piM;
788 phenParams->
fRing = (fRing_a*etap2 + fRing_b*eta + fRing_c)/piM;
789 phenParams->
sigma = (sigma_a*etap2 + sigma_b*eta + sigma_c)/piM;
791 phenParams->
psi0 = (psi0_a*etap2 + psi0_b*eta + psi0_c)/(eta*pow(piM, 5./3.));
792 phenParams->
psi1 = 0.;
793 phenParams->
psi2 = (psi2_a*etap2 + psi2_b*eta + psi2_c)/(eta*pow(piM, 3./3.));
794 phenParams->
psi3 = (psi3_a*etap2 + psi3_b*eta + psi3_c)/(eta*pow(piM, 2./3.));
795 phenParams->
psi4 = (psi4_a*etap2 + psi4_b*eta + psi4_c)/(eta*pow(piM, 1./3.));
796 phenParams->
psi5 = 0.;
797 phenParams->
psi6 = (psi6_a*etap2 + psi6_b*eta + psi6_c)/(eta*pow(piM, -1./3.));
798 phenParams->
psi7 = (psi7_a*etap2 + psi7_b*eta + psi7_c)/(eta*pow(piM, -2./3.));
814 REAL8 etap2, chip2, etap3, etap2chi, etachip2, etachi;
816 if (!
params || !phenParams)
return;
822 delta = sqrt(1.-4.*eta);
831 etap2chi = etap2*chi;
832 etachip2 = eta*chip2;
835 phenParams->
psi0 = 3./(128.*eta);
837 phenParams->
psi2 = 3715./756. +
838 -9.2091e+02*eta + 4.9213e+02*etachi + 1.3503e+02*etachip2 +
839 6.7419e+03*etap2 + -1.0534e+03*etap2chi +
842 phenParams->
psi3 = -16.*
LAL_PI + 113.*chi/3. +
843 1.7022e+04*eta + -9.5659e+03*etachi + -2.1821e+03*etachip2 +
844 -1.2137e+05*etap2 + 2.0752e+04*etap2chi +
847 phenParams->
psi4 = 15293365./508032. - 405.*chip2/8. +
848 -1.2544e+05*eta + 7.5066e+04*etachi + 1.3382e+04*etachip2 +
849 8.7354e+05*etap2 + -1.6573e+05*etap2chi +
852 phenParams->
psi6 = -8.8977e+05*eta + 6.3102e+05*etachi + 5.0676e+04*etachip2 +
853 5.9808e+06*etap2 + -1.4148e+06*etap2chi +
856 phenParams->
psi7 = 8.6960e+05*eta + -6.7098e+05*etachi + -3.0082e+04*etachip2 +
857 -5.8379e+06*etap2 + 1.5145e+06*etap2chi +
860 phenParams->
psi8 = -3.6600e+05*eta + 3.0670e+05*etachi + 6.3176e+02*etachip2 +
861 2.4265e+06*etap2 + -7.2180e+05*etap2chi +
864 phenParams->
fMerger = 1. - 4.4547*pow(1.-chi,0.217) + 3.521*pow(1.-chi,0.26) +
865 6.4365e-01*eta + 8.2696e-01*etachi + -2.7063e-01*etachip2 +
866 -5.8218e-02*etap2 + -3.9346e+00*etap2chi +
869 phenParams->
fRing = (1. - 0.63*pow(1.-chi,0.3))/2. +
870 1.4690e-01*eta + -1.2281e-01*etachi + -2.6091e-02*etachip2 +
871 -2.4900e-02*etap2 + 1.7013e-01*etap2chi +
874 phenParams->
sigma = (1. - 0.63*pow(1.-chi,0.3))*pow(1.-chi,0.45)/4. +
875 -4.0979e-01*eta + -3.5226e-02*etachi + 1.0082e-01*etachip2 +
876 1.8286e+00*etap2 + -2.0169e-02*etap2chi +
879 phenParams->
fCut = 3.2361e-01 + 4.8935e-02*chi + 1.3463e-02*chip2 +
880 -1.3313e-01*eta + -8.1719e-02*etachi + 1.4512e-01*etachip2 +
881 -2.7140e-01*etap2 + 1.2788e-01*etap2chi +
884 phenParams->
fCut /= piM;
886 phenParams->
fRing /= piM;
887 phenParams->
sigma /= piM;
889 phenParams->
psi1 = 0.;
890 phenParams->
psi5 = 0.;
902 REAL8 df, shft, phi, amp0, ampEff=0, psiEff, fMerg, fNorm;
903 REAL8 f, fRing, sigma, totalMass, eta;
917 totalMass = insp_template->
mass1 + insp_template->
mass2;
918 eta = insp_template->
mass1 * insp_template->
mass2 / pow(totalMass, 2.);
928 *(signalvec->
data+0) = 0.;
929 *(signalvec->
data+
n/2) = 0.;
934 for (
i=1;
i<nby2;
i++) {
943 if ((f < insp_template->fLower) || (f >
params->fCut))
946 ampEff = amp0*pow(fNorm, -7./6.);
947 else if ((f > fMerg) & (f <= fRing))
948 ampEff = amp0*pow(fNorm, -2./3.);
949 else if (f > fRing) {
951 ampEff *= amp0*
LAL_PI_2*pow(fRing/fMerg,-2./3.)*sigma;
955 psiEff = shft*f + phi
956 +
params->psi0*pow(f,-5./3.)
957 +
params->psi1*pow(f,-4./3.)
958 +
params->psi2*pow(f,-3./3.)
959 +
params->psi3*pow(f,-2./3.)
960 +
params->psi4*pow(f,-1./3.)
962 +
params->psi6*pow(f,1./3.)
963 +
params->psi7*pow(f,2./3.);
966 *(signalvec->
data+
i) = (
REAL4) (ampEff * cos(psiEff));
967 *(signalvec->
data+j) = (
REAL4) (ampEff * sin(psiEff));
982 REAL8 df, shft, phi, amp0, ampEff=0, psiEff, fMerg, fNorm;
983 REAL8 f, fRing, sigma, totalMass, eta;
985 REAL8 v, alpha2, alpha3, w1, vMerg, v2, v3, v4, v5, v6, v7, v8;
986 REAL8 epsilon_1, epsilon_2,
w2, vRing, chi, mergPower,
delta;
1005 totalMass = insp_template->
mass1 + insp_template->
mass2;
1006 eta = insp_template->
eta = insp_template->
mass1 * insp_template->
mass2 / pow(totalMass, 2.);
1007 delta = sqrt(1.-4.*eta);
1018 *(signalvec->
data+0) = 0.;
1019 *(signalvec->
data+
n/2) = 0.;
1027 alpha2 = -323./224. + 451.*insp_template->
eta/168.;
1028 alpha3 = (27./8. - 11.*insp_template->
eta/6.)*chi;
1034 epsilon_1 = 1.4547*chi - 1.8897;
1035 epsilon_2 = -1.8153*chi + 1.6557;
1041 w1 = 1. + alpha2*pow(vMerg,2.) + alpha3*pow(vMerg,3.);
1042 w1 = w1/(1. + epsilon_1*vMerg + epsilon_2*vMerg*vMerg);
1043 w2 = w1*(
LAL_PI*sigma/2.)*pow(fRing/fMerg, mergPower)*(1. + epsilon_1*vRing
1044 + epsilon_2*vRing*vRing);
1050 for (
i=1;
i<nby2;
i++) {
1061 v2 = v*v; v3 = v2*v; v4 = v2*v2; v5 = v4*v; v6 = v3*v3; v7 = v6*v, v8 = v7*v;
1064 if ((f < insp_template->fLower) || (f >
params->fCut))
1066 else if (f <= fMerg)
1067 ampEff = pow(fNorm, -7./6.)*(1. + alpha2*v2 + alpha3*v3);
1068 else if ((f > fMerg) & (f <= fRing))
1069 ampEff = w1*pow(fNorm, mergPower)*(1. + epsilon_1*v + epsilon_2*v2);
1074 psiEff = shft*f + phi
1075 + 3./(128.*eta*v5)*(1 +
params->psi2*v2
1081 *(signalvec->
data+
i) = (
REAL4) (amp0 * ampEff * cos(psiEff));
1082 *(signalvec->
data+j) = (
REAL4) (-amp0 * ampEff * sin(psiEff));
1106 hpDot->
data[0] = 0.0;
1107 hpDot->
data[len-1] = 0.0;
1108 hcDot->data[0] = 0.0;
1109 hcDot->data[len-1] = 0.0;
1110 for( k = 1; k < len-1; k++) {
1112 hcDot->data[k] = 1./(2.*
dt)*(hc->
data[k+1]-hc->
data[k-1]);
1117 for( k = 0; k < len; k++) {
1118 fOfT = -hcDot->data[k] * hp->
data[k] +hpDot->
data[k] * hc->
data[k];
1120 fOfT /= (pow(hp->
data[k],2.) + pow(hc->
data[k], 2.));
1136 return sigma / (2 *
LAL_PI * ((freq - fRing)*(freq - fRing)
1137 + sigma*sigma / 4.0));
1144 REAL8 UNUSED deltaT) {
1152 currentFreq = freq->
data[kMid];
1157 if (currentFreq > cutFreq && k > 0) {
1158 while (currentFreq > cutFreq)
1159 currentFreq = freq->
data[k--];
1163 while (currentFreq < cutFreq && k < len)
1164 currentFreq = freq->
data[k++];
1169 for (k = 0; k <
k0; k++)
1187 switch (
params->approximant) {
1209 switch (
params->approximant) {
int XLALInspiralInit(InspiralTemplate *params, InspiralInit *paramsInit)
static double double delta
#define ABORT(statusptr, code, mesg)
#define CHECKSTATUSPTR(statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
void XLALDestroyREAL4VectorSequence(REAL4VectorSequence *vecseq)
REAL4VectorSequence * XLALCreateREAL4VectorSequence(UINT4 length, UINT4 veclen)
#define BBHPHENOMCOEFFSH_PSI4_X
#define BBHPHENOMCOEFFSH_PSI3_Y
#define BBHPHENOMCOEFFSH_PSI2_Z
#define BBHPHENOMCOEFFSH_PSI6_Z
#define BBHPHENOMCOEFFSH_FRING_C
#define BBHPHENOMCOEFFSH_PSI2_Y
#define BBHPHENOMCOEFFSH_FRING_B
#define BBHPHENOMCOEFFSH_SIGMA_A
#define BBHPHENOMCOEFFSH_FRING_A
#define BBHPHENOMCOEFFSH_PSI0_X
#define BBHPHENOMCOEFFSH_PSI3_X
#define BBHPHENOMCOEFFSH_PSI3_Z
#define BBHPHENOMCOEFFSH_PSI6_X
#define BBHPHENOMCOEFFSH_PSI6_Y
#define BBHPHENOMCOEFFSH_FMERG_B
#define BBHPHENOMCOEFFSH_FCUT_B
#define BBHPHENOMCOEFFSH_PSI2_X
#define BBHPHENOMCOEFFSH_PSI7_X
#define BBHPHENOMCOEFFSH_PSI0_Z
#define BBHPHENOMCOEFFSH_PSI4_Y
#define BBHPHENOMCOEFFSH_SIGMA_B
#define BBHPHENOMCOEFFSH_SIGMA_C
#define BBHPHENOMCOEFFSH_FMERG_A
#define BBHPHENOMCOEFFSH_PSI4_Z
#define BBHPHENOMCOEFFSH_PSI7_Y
#define BBHPHENOMCOEFFSH_FMERG_C
#define BBHPHENOMCOEFFSH_FCUT_C
#define BBHPHENOMCOEFFSH_FCUT_A
#define BBHPHENOMCOEFFSH_PSI7_Z
#define BBHPHENOMCOEFFSH_PSI0_Y
#define GENERATEPPNINSPIRALH_EFSTOP
Reached requested termination frequency.
#define LALINSPIRALH_ENULL
Arguments contained an unexpected null pointer.
#define LALINSPIRALH_ESWITCH
Unknown case in switch.
void * XLALMalloc(size_t n)
REAL4FFTPlan * XLALCreateReverseREAL4FFTPlan(UINT4 size, int measurelvl)
void XLALDestroyREAL4FFTPlan(REAL4FFTPlan *plan)
int XLALREAL4VectorFFT(REAL4Vector *_LAL_RESTRICT_ output, const REAL4Vector *_LAL_RESTRICT_ input, const REAL4FFTPlan *plan)
REAL4TimeSeries * XLALCreateREAL4TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
void XLALDestroyREAL8TimeSeries(REAL8TimeSeries *series)
void XLALDestroyREAL4TimeSeries(REAL4TimeSeries *series)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
const LALUnit lalStrainUnit
const LALUnit lalHertzUnit
const LALUnit lalDimensionlessUnit
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
void XLALDestroyREAL4Vector(REAL4Vector *vector)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
int XLALREAL8VectorUnwrapAngle(REAL8Vector *out, const REAL8Vector *in)
void XLALError(const char *func, const char *file, int line, int errnum)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_PRINT_DEPRECATION_WARNING(replacement)
REAL4TimeVectorSeries * a
REAL4TimeVectorSeries * h
The inspiral waveform parameter structure containing information about the waveform to be generated.
REAL8 startTime
starting time of the waveform (in sec); if different from zero, the waveform will start with an insta...
INT4 nStartPad
Number of leading elements in the signal generation to be set to zero (input) If template is requeste...
REAL8 eta
symmetric mass ratio (input/output)
REAL8 spin2[3]
Spin vector of the secondary.
REAL8 mass1
Mass of the primary in solar mass (input/output)
REAL8 distance
Distance to the binary in seconds.
REAL8 startPhase
starting phase of the waveform in radians (input)
REAL8 spin1[3]
Spin vector of the primary.
REAL8 mass2
Mass of the secondary in solar mass (mass1 need not be larger than mass2 (input/output)
REAL8 tSampling
Sampling rate in Hz (input)
This structure stores the parameters for constructing a restricted post-Newtonian waveform.
const CHAR * termDescription
The termination code description (above)
INT4 termCode
The termination condition (above) that stopped computation of the waveform.
REAL4 dfdt
The maximum value of encountered over any timestep used in generating the waveform.
REAL4 phi
The phase at coalescence (or arbitrary reference phase for a post -Newtonian approximation),...
REAL4 psi
polarization angle (radians)
REAL4 fStart
The actual starting frequency of the waveform, in Hz (normally close but not identical to fStartIn)
SkyPosition position
location of source on sky
UINT4 length
The length of the generated waveform.
REAL4 fStop
The frequency at the termination of the waveform, in Hz.
REAL8 tc
The time from the start of the waveform to coalescence (in the point-mass approximation),...
REAL4 fStartIn
The requested starting frequency of the waveform, in Hz.
REAL8 deltaT
The requested sampling interval of the waveform, in s.
REAL4VectorSequence * data