LALSimIMRPhenom.c   LALSimIMRPhenom.c 
skipping to change at line 446 skipping to change at line 446
return phenParams; return phenParams;
} }
/*********************************************************************/ /*********************************************************************/
/* Compute phenomenological parameters for aligned-spin binaries */ /* Compute phenomenological parameters for aligned-spin binaries */
/* Ref. Eq.(2) and Table I of http://arxiv.org/pdf/0909.2867 */ /* Ref. Eq.(2) and Table I of http://arxiv.org/pdf/0909.2867 */
/* */ /* */
/* Takes solar masses. Populates and returns a new BBHPhenomParams */ /* Takes solar masses. Populates and returns a new BBHPhenomParams */
/* structure. */ /* structure. */
/*********************************************************************/ /*********************************************************************/
static BBHPhenomParams *ComputeIMRPhenomBParams(const REAL8 m1, const REAL8 const m2, const REAL8 chi) { static BBHPhenomParams *ComputeIMRPhenomBParams(const REAL8 m1, const REAL8 m2, const REAL8 chi) {
/* calculate the total mass and symmetric mass ratio */ /* calculate the total mass and symmetric mass ratio */
const REAL8 totalMass = m1 + m2; const REAL8 totalMass = m1 + m2;
const REAL8 eta = m1 * m2 / (totalMass * totalMass); const REAL8 eta = m1 * m2 / (totalMass * totalMass);
const REAL8 piM = totalMass * LAL_PI * LAL_MTSUN_SI; const REAL8 piM = totalMass * LAL_PI * LAL_MTSUN_SI;
/* spinning phenomenological waveforms */ /* spinning phenomenological waveforms */
const REAL8 etap2 = eta*eta; const REAL8 etap2 = eta*eta;
const REAL8 chip2 = chi*chi; const REAL8 chip2 = chi*chi;
const REAL8 etap3 = etap2*eta; const REAL8 etap3 = etap2*eta;
const REAL8 etap2chi = etap2*chi; const REAL8 etap2chi = etap2*chi;
skipping to change at line 667 skipping to change at line 667
const REAL8 phi0, /**< phase at peak */ const REAL8 phi0, /**< phase at peak */
const REAL8 deltaF, /**< frequency resolution */ const REAL8 deltaF, /**< frequency resolution */
const REAL8 m1, /**< mass of companion 1 [solar mass es] */ const REAL8 m1, /**< mass of companion 1 [solar mass es] */
const REAL8 m2, /**< mass of companion 2 [solar mass es] */ const REAL8 m2, /**< mass of companion 2 [solar mass es] */
const REAL8 chi, /**< mass-weighted aligned-spin para meter */ const REAL8 chi, /**< mass-weighted aligned-spin para meter */
const REAL8 f_min, /**< start frequency */ const REAL8 f_min, /**< start frequency */
const REAL8 f_max, /**< end frequency; if 0 */ const REAL8 f_max, /**< end frequency; if 0 */
const REAL8 distance, /**< distance of source */ const REAL8 distance, /**< distance of source */
const BBHPhenomParams *params /**< from ComputeIMRPhenomBParams */ const BBHPhenomParams *params /**< from ComputeIMRPhenomBParams */
) { ) {
const LIGOTimeGPS ligotimegps_zero = {0, 0}; static LIGOTimeGPS ligotimegps_zero = {0, 0};
size_t i; size_t i;
const REAL8 fMerg = params->fMerger; const REAL8 fMerg = params->fMerger;
const REAL8 fRing = params->fRing; const REAL8 fRing = params->fRing;
const REAL8 sigma = params->sigma; const REAL8 sigma = params->sigma;
const REAL8 totalMass = m1 + m2; const REAL8 totalMass = m1 + m2;
const REAL8 eta = m1 * m2 / (totalMass * totalMass); const REAL8 eta = m1 * m2 / (totalMass * totalMass);
const REAL8 piM = LAL_PI * totalMass * LAL_MTSUN_SI;
/* compute the amplitude pre-factor */ /* compute the amplitude pre-factor */
REAL8 amp0 = pow(LAL_MTSUN_SI*totalMass, 5./6.) * pow(fMerg, -7./6.) REAL8 amp0 = pow(LAL_MTSUN_SI*totalMass, 5./6.) * pow(fMerg, -7./6.)
/ pow(LAL_PI, 2./3.) * sqrt(5. * eta / 24.) / (distance / LAL_C_SI); / pow(LAL_PI, 2./3.) * sqrt(5. * eta / 24.) / (distance / LAL_C_SI);
/***********************************************************************/ /***********************************************************************/
/* these are the parameters required for the "new" phenomenological IMR /* these are the parameters required for the "new" phenomenological IMR
* waveforms*/ * waveforms*/
/***********************************************************************/ /***********************************************************************/
/* PN corrections to the frequency domain amplitude of the (2,2) mode */ /* PN corrections to the frequency domain amplitude of the (2,2) mode */
const REAL8 alpha2 = -323./224. + 451.*eta/168.; const REAL8 alpha2 = -323./224. + 451.*eta/168.;
const REAL8 alpha3 = (27./8. - 11.*eta/6.)*chi; const REAL8 alpha3 = (27./8. - 11.*eta/6.)*chi;
/* leading order power law of the merger amplitude */ /* leading order power law of the merger amplitude */
const REAL8 mergPower = -2./3.; static REAL8 mergPower = -2./3.;
/* spin-dependent corrections to the merger amplitude */ /* spin-dependent corrections to the merger amplitude */
const REAL8 epsilon_1 = 1.4547*chi - 1.8897; const REAL8 epsilon_1 = 1.4547*chi - 1.8897;
const REAL8 epsilon_2 = -1.8153*chi + 1.6557; const REAL8 epsilon_2 = -1.8153*chi + 1.6557;
/* normalisation constant of the inspiral amplitude */ /* normalisation constant of the inspiral amplitude */
REAL8 vMerg = cbrt(LAL_PI * totalMass * LAL_MTSUN_SI * fMerg); const REAL8 vMerg = cbrt(piM * fMerg);
REAL8 vRing = cbrt(LAL_PI * totalMass * LAL_MTSUN_SI * fRing); const REAL8 vRing = cbrt(piM * fRing);
REAL8 w1 = (1. + alpha2 * vMerg * vMerg + alpha3 * vMerg * vMerg * vMerg) / (1. + epsilon_1 * vMerg + epsilon_2 * vMerg * vMerg); REAL8 w1 = (1. + alpha2 * vMerg * vMerg + alpha3 * piM * fMerg) / (1. + e psilon_1 * vMerg + epsilon_2 * vMerg * vMerg);
REAL8 w2 = w1 * (LAL_PI * sigma / 2.) * pow(fRing / fMerg, mergPower) REAL8 w2 = w1 * (LAL_PI * sigma / 2.) * pow(fRing / fMerg, mergPower)
* (1. + epsilon_1 * vRing + epsilon_2 * vRing * vRing); * (1. + epsilon_1 * vRing + epsilon_2 * vRing * vRing);
/* allocate htilde */ /* allocate htilde */
size_t n = NextPow2(f_max / deltaF) + 1; size_t n = NextPow2(f_max / deltaF) + 1;
*htilde = XLALCreateCOMPLEX16FrequencySeries("htilde: FD waveform", &ligo timegps_zero, 0.0, deltaF, &lalStrainUnit, n); *htilde = XLALCreateCOMPLEX16FrequencySeries("htilde: FD waveform", &ligo timegps_zero, 0.0, deltaF, &lalStrainUnit, n);
memset((*htilde)->data->data, 0, n * sizeof(COMPLEX16)); memset((*htilde)->data->data, 0, n * sizeof(COMPLEX16));
XLALUnitDivide(&((*htilde)->sampleUnits), &((*htilde)->sampleUnits), &lal SecondUnit); XLALUnitDivide(&((*htilde)->sampleUnits), &((*htilde)->sampleUnits), &lal SecondUnit);
if (!(*htilde)) XLAL_ERROR(XLAL_EFUNC); if (!(*htilde)) XLAL_ERROR(XLAL_EFUNC);
/* now generate the waveform at all frequency bins except DC and Nyquist /* now generate the waveform */
*/ size_t ind_max = (size_t) (f_max / deltaF);
for (i=1; i < n - 1; i++) { for (i = (size_t) (f_min / deltaF); i < ind_max; i++) {
REAL8 ampEff, psiEff; REAL8 ampEff, psiEff;
REAL8 v, v2, v3, v4, v5, v6, v7, v8; REAL8 v, v2, v3, v4, v5, v6, v7, v8;
/* Fourier frequency corresponding to this bin */ /* Fourier frequency corresponding to this bin */
REAL8 f = i * deltaF; REAL8 f = i * deltaF;
REAL8 fNorm = f / fMerg;
/* PN expansion parameter */ /* PN expansion parameter */
v = cbrt(LAL_PI * totalMass * LAL_MTSUN_SI * f); v = cbrt(piM * f);
v2 = v*v; v3 = v2*v; v4 = v2*v2; v5 = v4*v; v6 = v3*v3; v7 = v6*v, v8 = v7*v; v2 = v*v; v3 = v2*v; v4 = v2*v2; v5 = v4*v; v6 = v3*v3; v7 = v6*v, v8 = v7*v;
/* compute the amplitude */ /* compute the amplitude */
if ((f < f_min) || (f > f_max)) if (f <= fMerg)
continue; ampEff = pow(f / fMerg, -7./6.)*(1. + alpha2 * v2 + alpha3 * v3);
else if (f <= fMerg)
ampEff = pow(fNorm, -7./6.)*(1. + alpha2 * v2 + alpha3 * v3);
else if ((f > fMerg) & (f <= fRing))
ampEff = w1 * pow(fNorm, mergPower) * (1. + epsilon_1 * v + epsilon_2
* v2);
else if (f > fRing) else if (f > fRing)
ampEff = w2 * LorentzianFn(f, fRing, sigma); ampEff = w2 * LorentzianFn(f, fRing, sigma);
else { else /* fMerg < f <= fRing */
XLALDestroyCOMPLEX16FrequencySeries(*htilde); ampEff = w1 * pow(f / fMerg, mergPower) * (1. + epsilon_1 * v + epsil
*htilde = NULL; on_2 * v2);
XLAL_ERROR(XLAL_EDOM);
}
/* now compute the phase */ /* now compute the phase */
psiEff = -phi0 /* phi is flipped relative to IMRPhenomA */ psiEff = -phi0 /* phi is flipped relative to IMRPhenomA */
+ 3./(128.*eta*v5)*(1 + params->psi2*v2 + 3./(128.*eta*v5)*(1 + params->psi2*v2
+ params->psi3*v3 + params->psi4*v4 + params->psi3*v3 + params->psi4*v4
+ params->psi5*v5 + params->psi6*v6 + params->psi5*v5 + params->psi6*v6
+ params->psi7*v7 + params->psi8*v8); + params->psi7*v7 + params->psi8*v8);
/* generate the waveform */ /* generate the waveform */
((*htilde)->data->data)[i].re = amp0 * ampEff * cos(psiEff); ((*htilde)->data->data)[i].re = amp0 * ampEff * cos(psiEff);
skipping to change at line 903 skipping to change at line 897
peak_amp_sq = amp_sq; peak_amp_sq = amp_sq;
} }
} }
return peak_ind; return peak_ind;
} }
static int apply_phase_shift(const REAL8TimeSeries *hp, const REAL8TimeSeri es *hc, const REAL8 shift) { static int apply_phase_shift(const REAL8TimeSeries *hp, const REAL8TimeSeri es *hc, const REAL8 shift) {
REAL8 *hpdata = hp->data->data; REAL8 *hpdata = hp->data->data;
REAL8 *hcdata = hc->data->data; REAL8 *hcdata = hc->data->data;
size_t k = hp->data->length; size_t k = hp->data->length;
const double cs = cos(shift);
const double ss = sin(shift);
for (;k--;) { for (;k--;) {
const REAL8 temp_hpdata = hpdata[k] * cos(shift) - hcdata[k] * sin( const REAL8 temp_hpdata = hpdata[k] * cs - hcdata[k] * ss;
shift); hcdata[k] = hpdata[k] * ss + hcdata[k] * cs;
hcdata[k] = hpdata[k] * sin(shift) + hcdata[k] * cos(shift);
hpdata[k] = temp_hpdata; hpdata[k] = temp_hpdata;
} }
return 0; return 0;
} }
static int apply_inclination(const REAL8TimeSeries *hp, const REAL8TimeSeri es *hc, const REAL8 inclination) { static int apply_inclination(const REAL8TimeSeries *hp, const REAL8TimeSeri es *hc, const REAL8 inclination) {
REAL8 inclFacPlus, inclFacCross; REAL8 inclFacPlus, inclFacCross;
REAL8 *hpdata = hp->data->data; REAL8 *hpdata = hp->data->data;
REAL8 *hcdata = hc->data->data; REAL8 *hcdata = hc->data->data;
size_t k = hp->data->length; size_t k = hp->data->length;
 End of changes. 13 change blocks. 
26 lines changed or deleted 20 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/