LALSimIMRPSpinInspiralRD.c   LALSimIMRPSpinInspiralRD.c 
skipping to change at line 52 skipping to change at line 52
#endif #endif
#define minIntLen 8 #define minIntLen 8
/* use error codes above 1024 to avoid conflicts with GSL */ /* use error codes above 1024 to avoid conflicts with GSL */
#define LALPSIRDPN_TEST_ENERGY 1025 #define LALPSIRDPN_TEST_ENERGY 1025
#define LALPSIRDPN_TEST_OMEGADOT 1026 #define LALPSIRDPN_TEST_OMEGADOT 1026
#define LALPSIRDPN_TEST_OMEGANAN 1028 #define LALPSIRDPN_TEST_OMEGANAN 1028
#define LALPSIRDPN_TEST_OMEGAMATCH 1029 #define LALPSIRDPN_TEST_OMEGAMATCH 1029
#define LALPSIRDPN_TEST_OMEGANONPOS 1031 #define LALPSIRDPN_TEST_OMEGANONPOS 1031
#define LALPSIRDPN_TEST_OMEGACUT 1032
typedef struct LALPSpinInspiralRDstructparams { typedef struct LALPSpinInspiralRDstructparams {
REAL8 dt; REAL8 dt;
REAL8 eta; ///< symmetric mass ratio REAL8 eta; ///< symmetric mass ratio
REAL8 dm; ///< \f$m_1-m_2\f$ REAL8 dm; ///< \f$m_1-m_2\f$
REAL8 m1m2; ///< \f$m_1/m_2\f$ REAL8 m1m2; ///< \f$m_1/m_2\f$
REAL8 m2m1; ///< \f$m_2/m_1\f$ REAL8 m2m1; ///< \f$m_2/m_1\f$
REAL8 m1m; REAL8 m1m;
REAL8 m2m; REAL8 m2m;
REAL8 m1msq; REAL8 m1msq;
REAL8 m2msq; REAL8 m2msq;
REAL8 m; REAL8 m;
REAL8 wdotorb[8]; ///< Coefficients of the analytic PN expansio n of \f$ \dot\omega_orb\f $ REAL8 wdotorb[8]; ///< Coefficients of the analytic PN expansio n of \f$ \dot\omega_orb\f$
REAL8 wdotorblog; ///< Log coefficient of the PN expansion of o f \f$\dot\omega_orb\f$ REAL8 wdotorblog; ///< Log coefficient of the PN expansion of o f \f$\dot\omega_orb\f$
REAL8 wdotspin15S1LNh; REAL8 wdotspin15S1LNh;
REAL8 wdotspin15S2LNh; REAL8 wdotspin15S2LNh;
REAL8 wdotspin20S1S2; REAL8 wdotspin20S1S2;
REAL8 wdotspin20S1S1; ///< Coeff. of the \f$s_1s_1\f$ cntrb. to \f$ \dot\omega\f$ REAL8 wdotspin20S1S1; ///< Coeff. of the \f$s_1s_1\f$ cntrb. to \f$ \dot\omega\f$
REAL8 wdotspin20S1S2LNh; REAL8 wdotspin20S1S2LNh;
REAL8 wdotspin25S1LNh; REAL8 wdotspin25S1LNh;
REAL8 wdotspin25S2LNh; ///< Coeff. of the \f$s_2\cdot \hat L_N\f$ cn trb. to \f$\dot\omega\f$ REAL8 wdotspin25S2LNh; ///< Coeff. of the \f$s_2\cdot \hat L_N\f$ cn trb. to \f$\dot\omega\f$
REAL8 wdotspin30S1LNh; REAL8 wdotspin30S1LNh;
REAL8 wdotspin30S2LNh; REAL8 wdotspin30S2LNh;
skipping to change at line 94 skipping to change at line 93
REAL8 epnspin15S1dotLNh; ///< Coeff. of the \f$S_1\cdot L\f$ term in en ergy REAL8 epnspin15S1dotLNh; ///< Coeff. of the \f$S_1\cdot L\f$ term in en ergy
REAL8 epnspin15S2dotLNh; ///< Coeff. of the \f$S_2\cdot L\f$ term in en ergy REAL8 epnspin15S2dotLNh; ///< Coeff. of the \f$S_2\cdot L\f$ term in en ergy
REAL8 epnspin20S1S2; ///< Coeff. of the \f$S_1\cdot S_2\f$ term in energy REAL8 epnspin20S1S2; ///< Coeff. of the \f$S_1\cdot S_2\f$ term in energy
REAL8 epnspin20S1S2dotLNh; ///< Coeff. of the \f$S_{1,2}\cdot L\f$ term i n energy REAL8 epnspin20S1S2dotLNh; ///< Coeff. of the \f$S_{1,2}\cdot L\f$ term i n energy
REAL8 epnspin20S1S1; ///< Coeff. of the \f$S_1\cdot S_1\f$ term in energy REAL8 epnspin20S1S1; ///< Coeff. of the \f$S_1\cdot S_1\f$ term in energy
REAL8 epnspin20S1S1dotLNh; REAL8 epnspin20S1S1dotLNh;
REAL8 epnspin20S2S2; ///< Coeff. of the \f$S_2\cdot S_2\f$ term in energy REAL8 epnspin20S2S2; ///< Coeff. of the \f$S_2\cdot S_2\f$ term in energy
REAL8 epnspin20S2S2dotLNh; REAL8 epnspin20S2S2dotLNh;
REAL8 epnspin25S1dotLNh; REAL8 epnspin25S1dotLNh;
REAL8 epnspin25S2dotLNh; REAL8 epnspin25S2dotLNh;
REAL8 OmCutoff;
REAL8 lengths; REAL8 lengths;
REAL8 omOffset; REAL8 omOffset;
REAL8 polarization; REAL8 polarization;
int length; int length;
UINT4 inspiralOnly; UINT4 inspiralOnly;
} LALPSpinInspiralRDparams; } LALPSpinInspiralRDparams;
static REAL8 OmMatch(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL 8 S2S2) { static REAL8 OmMatch(REAL8 LNhS1, REAL8 LNhS2, REAL8 S1S1, REAL8 S1S2, REAL 8 S2S2) {
const REAL8 omM = 0.0555; const REAL8 omM = 0.0555;
skipping to change at line 278 skipping to change at line 276
/** /**
* Convenience function to set up LALPSpinInspiralRDparams struct * Convenience function to set up LALPSpinInspiralRDparams struct
*/ */
static int XLALPSpinInspiralRDparamsSetup( static int XLALPSpinInspiralRDparamsSetup(
LALPSpinInspiralRDparams *mparams, /** Output: RDparams structure */ LALPSpinInspiralRDparams *mparams, /** Output: RDparams structure */
UINT4 inspiralOnly, /** Only generate inspiral */ UINT4 inspiralOnly, /** Only generate inspiral */
REAL8 deltaT, /** sampling interval */ REAL8 deltaT, /** sampling interval */
REAL8 fLow, /** Starting frequency */ REAL8 fLow, /** Starting frequency */
REAL8 fCutoff, /** CHECKME: Cutoff frequency? */
REAL8 m1, /** Mass 1 */ REAL8 m1, /** Mass 1 */
REAL8 m2, /** Mass 2 */ REAL8 m2, /** Mass 2 */
LALSimInspiralInteraction interaction, /** Spin interaction */ LALSimInspiralInteraction interaction, /** Spin interaction */
UINT4 order /** twice PN Order in Phase */ UINT4 order /** twice PN Order in Phase */
) )
{ {
REAL8 totalMass = m1+m2; REAL8 totalMass = m1+m2;
REAL8 eta = m1*m2/(totalMass * totalMass); REAL8 eta = m1*m2/(totalMass * totalMass);
REAL8 chirpMass = pow(m1*m2,0.6)/pow(totalMass,0.2); REAL8 chirpMass = pow(m1*m2,0.6)/pow(totalMass,0.2);
REAL8 ST[9]; /* SpinTaylor terms */ REAL8 ST[9]; /* SpinTaylor terms */
XLALSimInspiralSpinTaylorCoeffs(ST,eta); XLALSimInspiralSpinTaylorCoeffs(ST,eta);
mparams->inspiralOnly = inspiralOnly; mparams->inspiralOnly = inspiralOnly;
mparams->dt = deltaT; mparams->dt = deltaT;
mparams->OmCutoff = fCutoff*totalMass * LAL_MTSUN_SI * (REAL8) LAL_PI ;
mparams->lengths = (5.0 / 256.0) / LAL_PI * pow(LAL_PI * chirpMass * LAL_MTSUN_SI * fLow,-5.0 / 3.0) / fLow; mparams->lengths = (5.0 / 256.0) / LAL_PI * pow(LAL_PI * chirpMass * LAL_MTSUN_SI * fLow,-5.0 / 3.0) / fLow;
mparams->omOffset = 0.006; mparams->omOffset = 0.006;
/* setup coefficients for PN equations */ /* setup coefficients for PN equations */
mparams->m = totalMass; mparams->m = totalMass;
mparams->m2m1 = m2 / m1; mparams->m2m1 = m2 / m1;
mparams->m1m2 = m1 / m2; mparams->m1m2 = m1 / m2;
mparams->m1m = m1 / totalMass; mparams->m1m = m1 / totalMass;
mparams->m2m = m2 / totalMass; mparams->m2m = m2 / totalMass;
mparams->m1msq = mparams->m1m * mparams->m1m; mparams->m1msq = mparams->m1m * mparams->m1m;
skipping to change at line 754 skipping to change at line 750
return LALPSIRDPN_TEST_OMEGANONPOS; return LALPSIRDPN_TEST_OMEGANONPOS;
} }
else if (dvalues[1] < 0.0) { else if (dvalues[1] < 0.0) {
/* omegadot < 0 */ /* omegadot < 0 */
return LALPSIRDPN_TEST_OMEGADOT; return LALPSIRDPN_TEST_OMEGADOT;
} }
else if (isnan(omega)) { else if (isnan(omega)) {
/* omega is nan */ /* omega is nan */
return LALPSIRDPN_TEST_OMEGANAN; return LALPSIRDPN_TEST_OMEGANAN;
} }
else if ((params->inspiralOnly==1)&&(omega>params->OmCutoff)) {
return LALPSIRDPN_TEST_OMEGACUT;
}
else if ((params->inspiralOnly!=1)&&(omega>omegaMatch)) { else if ((params->inspiralOnly!=1)&&(omega>omegaMatch)) {
return LALPSIRDPN_TEST_OMEGAMATCH; return LALPSIRDPN_TEST_OMEGAMATCH;
} }
else else
return GSL_SUCCESS; return GSL_SUCCESS;
} }
static int XLALSpinInspiralFillH2Modes( static int XLALSpinInspiralFillH2Modes(
REAL8Vector* h2P2, REAL8Vector* h2P2,
REAL8Vector* h2M2, REAL8Vector* h2M2,
skipping to change at line 1113 skipping to change at line 1106
REAL8Vector* h4M2, REAL8Vector* h4M2,
REAL8Vector* h4P1, REAL8Vector* h4P1,
REAL8Vector* h4M1, REAL8Vector* h4M1,
REAL8Vector* h40, REAL8Vector* h40,
REAL8Vector* freq, REAL8Vector* freq,
REAL8Vector* phase, REAL8Vector* phase,
LALPSpinInspiralPhenPars *phenPars LALPSpinInspiralPhenPars *phenPars
) )
{ {
UINT4 j; INT4 j;
INT4 k; INT4 k;
INT4 kend=0; INT4 kMatch=0;
INT4 jend=0; INT4 jMatch=0;
INT4 Npoints=10;
INT4 intlen; INT4 intlen;
UINT4 intreturn; INT4 intreturn;
LALSpinInspiralAngle trigAngle; LALSpinInspiralAngle trigAngle;
REAL8Array *yout; REAL8Array *yout;
ark4GSLIntegrator *integrator; ark4GSLIntegrator *integrator;
REAL8 Psi; REAL8 Psi;
REAL8 alpha=0.; REAL8 alpha=0.;
REAL8 alphaold; REAL8 alphaold;
REAL8 v,v2; REAL8 v,v2;
skipping to change at line 1144 skipping to change at line 1138
REAL8 LNhxy; REAL8 LNhxy;
REAL8 LNhS1; REAL8 LNhS1;
REAL8 LNhS2; REAL8 LNhS2;
REAL8 S1S1; REAL8 S1S1;
REAL8 S1S2; REAL8 S1S2;
REAL8 S2S2; REAL8 S2S2;
REAL8 omegaMatch; REAL8 omegaMatch;
REAL8 c1,c2; REAL8 c1,c2;
INT4 errcode;
REAL8 *yin = (REAL8 *) LALMalloc(sizeof(REAL8) * neqs); REAL8 *yin = (REAL8 *) LALMalloc(sizeof(REAL8) * neqs);
/* allocate the integrator */ /* allocate the integrator */
integrator = XLALAdaptiveRungeKutta4Init(neqs,XLALSpinInspiralDerivatives ,XLALSpinInspiralTest,1.0e-6,1.0e-6); integrator = XLALAdaptiveRungeKutta4Init(neqs,XLALSpinInspiralDerivatives ,XLALSpinInspiralTest,1.0e-6,1.0e-6);
if (!integrator) { if (!integrator) {
fprintf(stderr,"**** LALPSpinInspiralRD ERROR ****: Cannot allocate ada ptive integrator.\n"); fprintf(stderr,"**** LALPSpinInspiralRD ERROR ****: Cannot allocate ada ptive integrator.\n");
if (XLALClearErrno() == XLAL_ENOMEM) if (XLALClearErrno() == XLAL_ENOMEM)
XLAL_ERROR( XLAL_ENOMEM ); XLAL_ERROR( XLAL_ENOMEM );
else else
XLAL_ERROR( XLAL_EDOM ); XLAL_ERROR( XLAL_EDOM );
} }
/* stop the integration only when the test is true */ /* stop the integration only when the test is true */
integrator->stopontestonly = 1; integrator->stopontestonly = 1;
/* run the integration; note: time is measured in units of total mass */ /* run the integration; note: time is measured in units of total mass */
Mass = mparams->m * LAL_MTSUN_SI; Mass = mparams->m * LAL_MTSUN_SI;
dt = mparams->dt; dt = mparams->dt;
for (j=0; j<neqs; j++) yin[j]=yinit[j]; for (UINT4 jeq=0; jeq<neqs; jeq++) yin[jeq]=yinit[jeq];
REAL8 S1x0=yinit[5]; REAL8 S1x0=yinit[5];
REAL8 S1y0=yinit[6]; REAL8 S1y0=yinit[6];
REAL8 S1z0=yinit[7]; REAL8 S1z0=yinit[7];
REAL8 S2x0=yinit[8]; REAL8 S2x0=yinit[8];
REAL8 S2y0=yinit[9]; REAL8 S2y0=yinit[9];
REAL8 S2z0=yinit[10]; REAL8 S2z0=yinit[10];
intlen = XLALAdaptiveRungeKutta4(integrator,(void *)mparams,yin,0.0,mpara ms->lengths/Mass,dt/Mass,&yout); intlen = XLALAdaptiveRungeKutta4Hermite(integrator,(void *)mparams,yin,0. 0,mparams->lengths/Mass,dt/Mass,&yout);
intreturn = integrator->returncode; intreturn = integrator->returncode;
XLALAdaptiveRungeKutta4Free(integrator); XLALAdaptiveRungeKutta4Free(integrator);
if (intlen == XLAL_FAILURE) if (intlen == XLAL_FAILURE)
{ {
XLALPrintError("Error in Adaptive Integrator\n"); XLALPrintError("Error in Adaptive Integrator\n");
XLAL_ERROR(XLAL_EFUNC); XLAL_ERROR(XLAL_EFUNC);
} }
/* End integration*/ /* End integration*/
skipping to change at line 1200 skipping to change at line 1192
if (XLALClearErrno() == XLAL_ENOMEM) { if (XLALClearErrno() == XLAL_ENOMEM) {
XLAL_ERROR( XLAL_ENOMEM); XLAL_ERROR( XLAL_ENOMEM);
} else { } else {
fprintf(stderr,"**** LALPSpinInspiralRD ERROR ****: integration faile d with errorcode %d, integration length %d\n",intreturn,intlen); fprintf(stderr,"**** LALPSpinInspiralRD ERROR ****: integration faile d with errorcode %d, integration length %d\n",intreturn,intlen);
XLAL_ERROR( XLAL_EFAILED); XLAL_ERROR( XLAL_EFAILED);
} }
} }
/* if we have enough space, compute the waveform components; otherwise ab ort */ /* if we have enough space, compute the waveform components; otherwise ab ort */
if ( intlen >= mparams->length ) { if ( intlen >= mparams->length ) {
fprintf(stderr, "%i\n", intlen);
fprintf(stderr, "%i\n", mparams->length);
fprintf(stderr,"**** LALPSpinInspiralRD ERROR ****: no space to write i n waveforms: %d vs. %d\n",intlen,mparams->length); fprintf(stderr,"**** LALPSpinInspiralRD ERROR ****: no space to write i n waveforms: %d vs. %d\n",intlen,mparams->length);
XLALPrintError("**** LALPSpinInspiralRD ERROR ****: error generating se
cond derivatives\n");
XLALPrintError(" m: : %12.5f %12.5f\n",m
params->m1m*mparams->m,mparams->m2m*mparams->m);
XLALPrintError(" S1: : %12.5f %12.5f %12
.5f\n",S1x0,S1y0,S1z0);
XLALPrintError(" S2: : %12.5f %12.5f %12
.5f\n",S2x0,S2y0,S2z0);
XLAL_ERROR(XLAL_ESIZE); XLAL_ERROR(XLAL_ESIZE);
} }
if ( intlen < minIntLen ) { if ( intlen < minIntLen ) {
fprintf(stderr,"**** LALPSpinInspiralRD ERROR ****: incorrect integrati on with length %d\n",intlen); fprintf(stderr,"**** LALPSpinInspiralRD ERROR ****: incorrect integrati on with length %d\n",intlen);
XLAL_ERROR(XLAL_ESIZE); XLAL_ERROR(XLAL_ESIZE);
} }
/* End of integration checks*/ /* End of integration checks*/
REAL8 *Phi = &yout->data[1*intlen]; REAL8 *Phi = &yout->data[1*intlen];
skipping to change at line 1229 skipping to change at line 1223
REAL8 *S1z = &yout->data[8*intlen]; REAL8 *S1z = &yout->data[8*intlen];
REAL8 *S2x = &yout->data[9*intlen]; REAL8 *S2x = &yout->data[9*intlen];
REAL8 *S2y = &yout->data[10*intlen]; REAL8 *S2y = &yout->data[10*intlen];
REAL8 *S2z = &yout->data[11*intlen]; REAL8 *S2z = &yout->data[11*intlen];
REAL8 *energy = &yout->data[12*intlen]; REAL8 *energy = &yout->data[12*intlen];
mparams->polarization=2.*atan2(LNhy[1],LNhx[1])-atan2(LNhy[2],LNhx[2]); mparams->polarization=2.*atan2(LNhy[1],LNhx[1])-atan2(LNhy[2],LNhx[2]);
if (mparams->inspiralOnly!=1) { if (mparams->inspiralOnly!=1) {
j=intlen; INT4 errcode=XLAL_SUCCESS;
jend=0;
INT4 Npoints=10; j=intlen;
do { do {
j--; j--;
LNhS1=(LNhx[j]*S1x[j]+LNhy[j]*S1y[j]+LNhz[j]*S1z[j])/mparams->m1msq; LNhS1=(LNhx[j]*S1x[j]+LNhy[j]*S1y[j]+LNhz[j]*S1z[j])/mparams->m1msq;
LNhS2=(LNhx[j]*S2x[j]+LNhy[j]*S2y[j]+LNhz[j]*S2z[j])/mparams->m2msq; LNhS2=(LNhx[j]*S2x[j]+LNhy[j]*S2y[j]+LNhz[j]*S2z[j])/mparams->m2msq;
S1S1=(S1x[j]*S1x[j]+S1y[j]*S1y[j]+S1z[j]*S1z[j])/mparams->m1msq/mpara ms->m1msq; S1S1=(S1x[j]*S1x[j]+S1y[j]*S1y[j]+S1z[j]*S1z[j])/mparams->m1msq/mpara ms->m1msq;
S1S2=(S1x[j]*S2x[j]+S1y[j]*S2y[j]+S1z[j]*S2z[j])/mparams->m1msq/mpara ms->m2msq; S1S2=(S1x[j]*S2x[j]+S1y[j]*S2y[j]+S1z[j]*S2z[j])/mparams->m1msq/mpara ms->m2msq;
S2S2=(S2x[j]*S2x[j]+S2y[j]*S2y[j]+S2z[j]*S2z[j])/mparams->m2msq/mpara ms->m2msq; S2S2=(S2x[j]*S2x[j]+S2y[j]*S2y[j]+S2z[j]*S2z[j])/mparams->m2msq/mpara ms->m2msq;
omegaMatch=OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2); omegaMatch=OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2);
if (omegaMatch>omega[j]) { if (omegaMatch>omega[j]) {
if (omega[j-1]<omega[j]) jend=j; if (omega[j-1]<omega[j]) jMatch=j;
// The numerical integrator sometimes stops and stores twice the last // The numerical integrator sometimes stops and stores twice the last
// omega value, this 'if' instruction avoids keeping two identical // omega value, this 'if' instruction avoids keeping two identical
// values of omega at the end of the integra tion. // values of omega at the end of the integra tion.
} }
} while ((j>0)&&(jend==0)); } while ((j>0)&&(jMatch==0));
if (omegaMatch<omega[jend]) { if (omegaMatch<omega[jMatch]) {
fprintf(stderr,"*** LALPSpinInspiralRD ERROR ***: Impossible to attac h phenom. part\n"); fprintf(stderr,"*** LALPSpinInspiralRD ERROR ***: Impossible to attac h phenom. part\n");
XLAL_ERROR(XLAL_EFAILED); XLAL_ERROR(XLAL_EFAILED);
} }
kend=Npoints-1; // Data structure are copied into Npoints-long
if (omega[jend+1]>omega[jend]) { // REAL8Array for interpolation and derivative computation
jend++; if ( ((INT4)Npoints) > intlen) Npoints = intlen;
kend--;
} if ( (omega[jMatch+1]>omega[jMatch]) && ((jMatch+1)<intlen) )
kMatch=Npoints-2;
else
kMatch=Npoints-1;
//We keep until the point where omega > omegaMatch for better derivativ e //We keep until the point where omega > omegaMatch for better derivativ e
// computation, but do the matching at the last point at which //computation, but do the matching at the last point at which
// omega < omegaMatch // omega < omegaMatch
if (Npoints > jend) Npoints = jend+1;
REAL8Vector *omega_s = XLALCreateREAL8Vector(Npoints); REAL8Vector *omega_s = XLALCreateREAL8Vector(Npoints);
REAL8Vector *LNhx_s = XLALCreateREAL8Vector(Npoints); REAL8Vector *LNhx_s = XLALCreateREAL8Vector(Npoints);
REAL8Vector *LNhy_s = XLALCreateREAL8Vector(Npoints); REAL8Vector *LNhy_s = XLALCreateREAL8Vector(Npoints);
REAL8Vector *LNhz_s = XLALCreateREAL8Vector(Npoints); REAL8Vector *LNhz_s = XLALCreateREAL8Vector(Npoints);
REAL8Vector *alpha_s = XLALCreateREAL8Vector(Npoints); REAL8Vector *alpha_s = XLALCreateREAL8Vector(Npoints);
REAL8Vector *domega = XLALCreateREAL8Vector(Npoints); REAL8Vector *domega = XLALCreateREAL8Vector(Npoints);
REAL8Vector *dLNhx = XLALCreateREAL8Vector(Npoints); REAL8Vector *dLNhx = XLALCreateREAL8Vector(Npoints);
REAL8Vector *dLNhy = XLALCreateREAL8Vector(Npoints); REAL8Vector *dLNhy = XLALCreateREAL8Vector(Npoints);
REAL8Vector *dLNhz = XLALCreateREAL8Vector(Npoints); REAL8Vector *dLNhz = XLALCreateREAL8Vector(Npoints);
REAL8Vector *diota = XLALCreateREAL8Vector(Npoints); REAL8Vector *diota = XLALCreateREAL8Vector(Npoints);
REAL8Vector *dalpha = XLALCreateREAL8Vector(Npoints); REAL8Vector *dalpha = XLALCreateREAL8Vector(Npoints);
REAL8Vector *ddomega = XLALCreateREAL8Vector(Npoints); REAL8Vector *ddomega = XLALCreateREAL8Vector(Npoints);
REAL8Vector *ddiota = XLALCreateREAL8Vector(Npoints); REAL8Vector *ddiota = XLALCreateREAL8Vector(Npoints);
REAL8Vector *ddalpha = XLALCreateREAL8Vector(Npoints); REAL8Vector *ddalpha = XLALCreateREAL8Vector(Npoints);
for (k=0;k<Npoints;k++) { for (k=0;k<Npoints;k++) {
j=k+jend-Npoints+1; j=k+jMatch-kMatch;
omega_s->data[k] = omega[j]; omega_s->data[k] = omega[j];
LNhx_s->data[k] = LNhx[j]; LNhx_s->data[k] = LNhx[j];
LNhy_s->data[k] = LNhy[j]; LNhy_s->data[k] = LNhy[j];
LNhz_s->data[k] = LNhz[j]; LNhz_s->data[k] = LNhz[j];
} }
errcode = XLALGenerateWaveDerivative(domega,omega_s,dt); errcode = XLALGenerateWaveDerivative(domega,omega_s,dt);
errcode += XLALGenerateWaveDerivative(dLNhx,LNhx_s,dt); errcode += XLALGenerateWaveDerivative(dLNhx,LNhx_s,dt);
errcode += XLALGenerateWaveDerivative(dLNhy,LNhy_s,dt); errcode += XLALGenerateWaveDerivative(dLNhy,LNhy_s,dt);
errcode += XLALGenerateWaveDerivative(dLNhz,LNhz_s,dt); errcode += XLALGenerateWaveDerivative(dLNhz,LNhz_s,dt);
skipping to change at line 1319 skipping to change at line 1313
} }
errcode = XLALGenerateWaveDerivative(ddiota,diota,dt); errcode = XLALGenerateWaveDerivative(ddiota,diota,dt);
errcode += XLALGenerateWaveDerivative(ddalpha,dalpha,dt); errcode += XLALGenerateWaveDerivative(ddalpha,dalpha,dt);
errcode += XLALGenerateWaveDerivative(ddomega,domega,dt); errcode += XLALGenerateWaveDerivative(ddomega,domega,dt);
if (errcode != 0) { if (errcode != 0) {
XLALPrintError("**** LALPSpinInspiralRD ERROR ****: error generating second derivatives\n"); XLALPrintError("**** LALPSpinInspiralRD ERROR ****: error generating second derivatives\n");
XLALPrintError(" m: : %12.5f %12.5f\n" ,mparams->m1m*mparams->m,mparams->m2m*mparams->m); XLALPrintError(" m: : %12.5f %12.5f\n" ,mparams->m1m*mparams->m,mparams->m2m*mparams->m);
XLALPrintError(" S1: : %12.5f %12.5f % 12.5f\n",S1x0,S1y0,S1z0); XLALPrintError(" S1: : %12.5f %12.5f % 12.5f\n",S1x0,S1y0,S1z0);
XLALPrintError(" S2: : %12.5f %12.5f % 12.5f\n",S2x0,S2y0,S2z0); XLALPrintError(" S2: : %12.5f %12.5f % 12.5f\n",S2x0,S2y0,S2z0);
XLALPrintError(" omM %12.5f om[%d] %12.5f\n",omegaMatch,jend,om ega); XLALPrintError(" omM %12.5f om[%d] %12.5f\n",omegaMatch,jMatch, omega);
XLAL_ERROR(XLAL_EFAILED); XLAL_ERROR(XLAL_EFAILED);
} }
if (ddomega->data[kend]<0.) { if (ddomega->data[kMatch]<0.) {
fprintf(stdout,"*** LALPSpinInspiralRD WARNING: the attach of the phe nom. phase has been shifted back: m1 %12.6f m2 %12.6f\n",mparams->m1m*mpar ams->m,mparams->m2m*mparams->m); fprintf(stdout,"*** LALPSpinInspiralRD WARNING: the attach of the phe nom. phase has been shifted back: m1 %12.6f m2 %12.6f\n",mparams->m1m*mpar ams->m,mparams->m2m*mparams->m);
fprintf(stdout," Integration returned %d\n 1025: Energy increases\ fprintf(stdout," Integration returned %d\n 1025: Energy increases\
n 1026: Omegadot -ve\n 1028: Omega NAN\n 1029: Omega > Omegamatch\n n 1026: Omegadot -ve\n 1028: Omega NAN\n 1029: Omega > Omegamatch\n
1031: Omega -ve\n 1032: Omega > OmegaCut %12.6e\n",intreturn,mparams->Om 1031: Omega -ve\n",intreturn);
Cutoff); while ((kMatch>0)&&(ddomega->data[kMatch]<0.)) {
while ((kend>0)&&(ddomega->data[kend]<0.)) { kMatch--;
kend--; jMatch--;
jend--;
} }
} }
phenPars->intreturn = intreturn; phenPars->intreturn = intreturn;
phenPars->energy = energy[jend]; phenPars->energy = energy[jMatch];
phenPars->omega = omega_s->data[kend]; phenPars->omega = omega_s->data[kMatch];
phenPars->domega = domega->data[kend]; phenPars->domega = domega->data[kMatch];
phenPars->ddomega = ddomega->data[kend]; phenPars->ddomega = ddomega->data[kMatch];
phenPars->diota = diota->data[kend]; phenPars->diota = diota->data[kMatch];
phenPars->ddiota = ddiota->data[kend]; phenPars->ddiota = ddiota->data[kMatch];
phenPars->dalpha = dalpha->data[kend]; phenPars->dalpha = dalpha->data[kMatch];
phenPars->ddalpha = ddalpha->data[kend]; phenPars->ddalpha = ddalpha->data[kMatch];
phenPars->countback = jend; phenPars->countback = jMatch;
phenPars->Psi = Phi[jend]; phenPars->Psi = Phi[jMatch];
phenPars->endtime = ((REAL8) jend)*dt; phenPars->endtime = ((REAL8) jMatch)*dt;
phenPars->ci = LNhz[jend]; phenPars->ci = LNhz[jMatch];
phenPars->LNhS1 = LNhS1; phenPars->LNhS1 = LNhS1;
phenPars->LNhS2 = LNhS2; phenPars->LNhS2 = LNhS2;
phenPars->S1S2 = S1S2; phenPars->S1S2 = S1S2;
phenPars->S1S1 = S1S1; phenPars->S1S1 = S1S1;
phenPars->S2S2 = S2S2; phenPars->S2S2 = S2S2;
XLALDestroyREAL8Vector(omega_s); XLALDestroyREAL8Vector(omega_s);
XLALDestroyREAL8Vector(LNhx_s); XLALDestroyREAL8Vector(LNhx_s);
XLALDestroyREAL8Vector(LNhy_s); XLALDestroyREAL8Vector(LNhy_s);
XLALDestroyREAL8Vector(LNhz_s); XLALDestroyREAL8Vector(LNhz_s);
skipping to change at line 1367 skipping to change at line 1361
XLALDestroyREAL8Vector(dLNhy); XLALDestroyREAL8Vector(dLNhy);
XLALDestroyREAL8Vector(dLNhz); XLALDestroyREAL8Vector(dLNhz);
XLALDestroyREAL8Vector(diota); XLALDestroyREAL8Vector(diota);
XLALDestroyREAL8Vector(dalpha); XLALDestroyREAL8Vector(dalpha);
XLALDestroyREAL8Vector(domega); XLALDestroyREAL8Vector(domega);
XLALDestroyREAL8Vector(ddomega); XLALDestroyREAL8Vector(ddomega);
XLALDestroyREAL8Vector(ddiota); XLALDestroyREAL8Vector(ddiota);
XLALDestroyREAL8Vector(ddalpha); XLALDestroyREAL8Vector(ddalpha);
} }
else { else {
jend=intlen-1; jMatch=intlen-1;
phenPars->intreturn = intreturn; phenPars->intreturn = intreturn;
phenPars->energy = 0.; phenPars->energy = 0.;
phenPars->omega = 0.; phenPars->omega = 0.;
phenPars->domega = 0.; phenPars->domega = 0.;
phenPars->ddomega = 0.; phenPars->ddomega = 0.;
phenPars->diota = 0.; phenPars->diota = 0.;
phenPars->ddiota = 0.; phenPars->ddiota = 0.;
phenPars->dalpha = 0.; phenPars->dalpha = 0.;
phenPars->ddalpha = 0.; phenPars->ddalpha = 0.;
phenPars->countback = intlen-1; phenPars->countback = intlen-1;
skipping to change at line 1413 skipping to change at line 1407
else { else {
if ((S1x[0]*S1x[0]+S1y[0]*S1y[0]+S2x[0]*S2x[0]+S2y[0]*S2y[0])>0.) { if ((S1x[0]*S1x[0]+S1y[0]*S1y[0]+S2x[0]*S2x[0]+S2y[0]*S2y[0])>0.) {
c1=0.75+mparams->eta/2-0.75*mparams->dm; c1=0.75+mparams->eta/2-0.75*mparams->dm;
c2=0.75+mparams->eta/2+0.75*mparams->dm; c2=0.75+mparams->eta/2+0.75*mparams->dm;
alpha=atan2(-c1*S1x[0]-c2*S2x[0],c1*S1y[0]+c2*S2y[0]); alpha=atan2(-c1*S1x[0]-c2*S2x[0],c1*S1y[0]+c2*S2y[0]);
} }
else else
alpha=0.; alpha=0.;
} }
for (j=0;j<=(UINT4)jend;j++) { for (j=0;j<=jMatch;j++) {
freq->data[j]=omega[j]; freq->data[j]=omega[j];
v=cbrt(omega[j]); v=cbrt(omega[j]);
v2=v*v; v2=v*v;
// amp22= -2.0 * params->mu * LAL_MRSUN_SI/(params->distance) * sqrt( 1 6.*LAL_PI/5.)*v2; // amp22= -2.0 * params->mu * LAL_MRSUN_SI/(params->distance) * sqrt( 1 6.*LAL_PI/5.)*v2;
// amp20 = amp22*sqrt(3/2) // amp20 = amp22*sqrt(3/2)
// Y22 \pm Y2-2= sqrt(5/PI) ((1+cos^2 t)/4, (cos t)/2) // Y22 \pm Y2-2= sqrt(5/PI) ((1+cos^2 t)/4, (cos t)/2)
// Y21 \pm Y2-1= sqrt(5/PI) ((sin t)/2, (sin 2t)/4) // Y21 \pm Y2-1= sqrt(5/PI) ((sin t)/2, (sin 2t)/4)
// Y20 = sqrt(15/2 PI) (sin^2 t)/4 // Y20 = sqrt(15/2 PI) (sin^2 t)/4
skipping to change at line 1459 skipping to change at line 1453
trigAngle.c8i2 = trigAngle.c4i2 * trigAngle.c4i2; trigAngle.c8i2 = trigAngle.c4i2 * trigAngle.c4i2;
trigAngle.s8i2 = trigAngle.s4i2 * trigAngle.s4i2; trigAngle.s8i2 = trigAngle.s4i2 * trigAngle.s4i2;
//alphaoold = alphaold; //alphaoold = alphaold;
alphaold = alpha; alphaold = alpha;
if ((LNhy[j]*LNhy[j]+LNhx[j]*LNhx[j])>0.) { if ((LNhy[j]*LNhy[j]+LNhx[j]*LNhx[j])>0.) {
alpha = atan2(LNhy[j], LNhx[j]); alpha = atan2(LNhy[j], LNhx[j]);
} }
else alpha = alphaold; else alpha = alphaold;
errcode = XLALSpinInspiralFillH2Modes(h2P2,h2M2,h2P1,h2M1,h20,j,amp22, v,mparams->eta,mparams->dm,Psi,alpha,&trigAngle); XLALSpinInspiralFillH2Modes(h2P2,h2M2,h2P1,h2M1,h20,j,amp22,v,mparams-> eta,mparams->dm,Psi,alpha,&trigAngle);
/*if (j>2) { /*if (j>2) {
if ((alphaold*alphaoold)<0.) { if ((alphaold*alphaoold)<0.) {
if ( fabs(cos(2.*(Phi[j-1]+alphaold))-cos(2.*(Phi[j-2]+alph aoold)))>0.2) { if ( fabs(cos(2.*(Phi[j-1]+alphaold))-cos(2.*(Phi[j-2]+alph aoold)))>0.2) {
fprintf(stdout,"*** LALPSpinInspiralRD WARNING ***: Possibl e problem with coordinate singularity:\n Step %d LNhy: %12.6e LNhx: %12.6e Psi+alpha: %12.6e alpha %12.6e\n Step %d LNhy: %12.6e LNhx: %12.6e Psi +alpha: %12.6e alpha %12.6e\n Step %d LNhy: %12.6e LNhx: %12.6e Psi+alp ha: %12.6e alpha %12.6e\n",j,LNhy[j],LNhx[j],(Phi[j]+alpha)/LAL_PI,alpha/L AL_PI,j-1,LNhy[j-1],LNhx[j-1],(Phi[j-1]+alphaold)/LAL_PI,alphaold/LAL_PI,j- 2,LNhy[j-2],LNhx[j-2],(Phi[j-2]+alphaoold)/LAL_PI,alphaoold/LAL_PI); fprintf(stdout,"*** LALPSpinInspiralRD WARNING ***: Possibl e problem with coordinate singularity:\n Step %d LNhy: %12.6e LNhx: %12.6e Psi+alpha: %12.6e alpha %12.6e\n Step %d LNhy: %12.6e LNhx: %12.6e Psi +alpha: %12.6e alpha %12.6e\n Step %d LNhy: %12.6e LNhx: %12.6e Psi+alp ha: %12.6e alpha %12.6e\n",j,LNhy[j],LNhx[j],(Phi[j]+alpha)/LAL_PI,alpha/L AL_PI,j-1,LNhy[j-1],LNhx[j-1],(Phi[j-1]+alphaold)/LAL_PI,alphaold/LAL_PI,j- 2,LNhy[j-2],LNhx[j-2],(Phi[j-2]+alphaoold)/LAL_PI,alphaoold/LAL_PI);
fprintf(stdout," m: (%12.6e,%12.6e)\n", mparams- >m1m*mparams->m, mparams->m2m*mparams->m); fprintf(stdout," m: (%12.6e,%12.6e)\n", mparams- >m1m*mparams->m, mparams->m2m*mparams->m);
fprintf(stdout," S1: (%9.6f,%9.6f,%9.6f)\n",yini t[5]/mparams->m1msq,yinit[6]/mparams->m1msq,yinit[7]/mparams->m1msq); fprintf(stdout," S1: (%9.6f,%9.6f,%9.6f)\n",yini t[5]/mparams->m1msq,yinit[6]/mparams->m1msq,yinit[7]/mparams->m1msq);
fprintf(stdout," S2: (%9.6f,%9.6f,%9.6f)\n",yini t[8]/mparams->m2msq,yinit[9]/mparams->m2msq,yinit[10]/mparams->m2msq); fprintf(stdout," S2: (%9.6f,%9.6f,%9.6f)\n",yini t[8]/mparams->m2msq,yinit[9]/mparams->m2msq,yinit[10]/mparams->m2msq);
} }
} }
}*/ }*/
errcode += XLALSpinInspiralFillH3Modes(h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h3 0,j,amp33,v,mparams->eta,mparams->dm,Psi,alpha,&trigAngle); XLALSpinInspiralFillH3Modes(h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h30,j,amp33,v ,mparams->eta,mparams->dm,Psi,alpha,&trigAngle);
errcode += XLALSpinInspiralFillH4Modes(h4P4,h4M4,h4P3,h4M3,h4P2,h4M2,h4 P1,h4M1,h40,j,amp44,v,mparams->eta,mparams->dm,Psi,alpha,&trigAngle); XLALSpinInspiralFillH4Modes(h4P4,h4M4,h4P3,h4M3,h4P2,h4M2,h4P1,h4M1,h40 ,j,amp44,v,mparams->eta,mparams->dm,Psi,alpha,&trigAngle);
} }
if (yin) LALFree(yin); if (yin) LALFree(yin);
if (yout) XLALDestroyREAL8Array(yout); if (yout) XLALDestroyREAL8Array(yout);
phenPars->alpha=alpha; phenPars->alpha=alpha;
return errcode; return XLAL_SUCCESS;
} /* End of the inspiral part created via the adaptive integration method * / } /* End of the inspiral part created via the adaptive integration method * /
int XLALSimIMRPSpinFinalMassSpin( int XLALSimIMRPSpinFinalMassSpin(
REAL8 *finalMass, REAL8 *finalMass,
REAL8 *finalSpin, REAL8 *finalSpin,
REAL8 m1, REAL8 m1,
REAL8 m2, REAL8 m2,
REAL8 s1x, REAL8 s1x,
REAL8 s1y, REAL8 s1y,
skipping to change at line 2102 skipping to change at line 2096
REAL8 f_min, /**< start frequency */ REAL8 f_min, /**< start frequency */
REAL8 r, /**< distance of source */ REAL8 r, /**< distance of source */
REAL8 iota, /**< inclination of source (rad) */ REAL8 iota, /**< inclination of source (rad) */
REAL8 s1x, /**< x-component of dimensionless spin for object 1 */ REAL8 s1x, /**< x-component of dimensionless spin for object 1 */
REAL8 s1y, /**< y-component of dimensionless spin for object 1 */ REAL8 s1y, /**< y-component of dimensionless spin for object 1 */
REAL8 s1z, /**< z-component of dimensionless spin for object 1 */ REAL8 s1z, /**< z-component of dimensionless spin for object 1 */
REAL8 s2x, /**< x-component of dimensionless spin for object 2 */ REAL8 s2x, /**< x-component of dimensionless spin for object 2 */
REAL8 s2y, /**< y-component of dimensionless spin for object 2 */ REAL8 s2y, /**< y-component of dimensionless spin for object 2 */
REAL8 s2z, /**< z-component of dimensionless spin for object 2 */ REAL8 s2z, /**< z-component of dimensionless spin for object 2 */
int phaseO, /**< twice post-Newtonian phase order */ int phaseO, /**< twice post-Newtonian phase order */
InputAxis axisChoice /**< Choice of axis for input spin params * InputAxis axisChoice, /**< Choice of axis for input spin params *
/ /
int inspiralOnly /**< 0 generate RD, 1 generate inspiralOnly
*/
) )
{ {
LIGOTimeGPS t_start = LIGOTIMEGPSZERO; LIGOTimeGPS t_start = LIGOTIMEGPSZERO;
UINT4 length; // signal vector length UINT4 length; // signal vector length
REAL8 tn = XLALSimInspiralTaylorLength(deltaT, m1, m2, f_min, phaseO); REAL8 tn = XLALSimInspiralTaylorLength(deltaT, m1, m2, f_min, phaseO);
REAL8 x = 1.1 * (tn + 1. ) / deltaT; REAL8 x = 1.1 * (tn + 1. ) / deltaT;
length = ceil(log10(x)/log10(2.)); length = ceil(log10(x)/log10(2.));
length = pow(2, length); length = pow(2, length);
skipping to change at line 2202 skipping to change at line 2197
/* The number of Ring Down modes is hard-coded here */ /* The number of Ring Down modes is hard-coded here */
const UINT4 nmodes=2; const UINT4 nmodes=2;
/* Nmodes should be restricted to either 1 or 2*/ /* Nmodes should be restricted to either 1 or 2*/
UINT4 errcode; UINT4 errcode;
REAL8 finalMass,finalSpin; REAL8 finalMass,finalSpin;
REAL8 energy=0.; REAL8 energy=0.;
REAL8 omegaMatch; REAL8 omegaMatch;
REAL8 frOmRD,omegaRD; REAL8 frOmRD,omegaRD;
REAL8 fCutoff = 0.0; /* Set only for inspiral only waveforms */
REAL8 mass1 = m1 / LAL_MSUN_SI; /* mass of object 1 in solar masses */ REAL8 mass1 = m1 / LAL_MSUN_SI; /* mass of object 1 in solar masses */
REAL8 mass2 = m2 / LAL_MSUN_SI; /* mass of object 2 in solar masses */ REAL8 mass2 = m2 / LAL_MSUN_SI; /* mass of object 2 in solar masses */
REAL8 totalMass = mass1 + mass2; REAL8 totalMass = mass1 + mass2;
REAL8 eta = mass1 * mass2 / (totalMass * totalMass); REAL8 eta = mass1 * mass2 / (totalMass * totalMass);
REAL8 mass = (m1 + m2) * LAL_G_SI / pow(LAL_C_SI, 3.0); /* convert m from kilograms to seconds */ REAL8 mass = (m1 + m2) * LAL_G_SI / pow(LAL_C_SI, 3.0); /* convert m from kilograms to seconds */
REAL8 mu = eta * totalMass; /* Reduced mass in solar masses */ REAL8 mu = eta * totalMass; /* Reduced mass in solar masses */
unitHz = mass * (REAL8) LAL_PI; unitHz = mass * (REAL8) LAL_PI;
/* /*
if ((signalvec2)||(hh)) if ((signalvec2)||(hh))
params->nStartPad = 0;*/ /* must be zero for templates and injection */ params->nStartPad = 0;*/ /* must be zero for templates and injection */
/* -- length in seconds from Newtonian formula; */ /* -- length in seconds from Newtonian formula; */
dt = deltaT; dt = deltaT;
/* setup coefficients for PN equations */ /* setup coefficients for PN equations */
if(XLALPSpinInspiralRDparamsSetup(&mparams, /** Output: RDparams structur e */ if(XLALPSpinInspiralRDparamsSetup(&mparams, /** Output: RDparams structur e */
0, /** Do not Only generate inspiral */ inspiralOnly, /** Do not Only generate ins piral */
deltaT, /** sampling interval */ deltaT, /** sampling interval */
f_min, /** Starting frequency */ f_min, /** Starting frequency */
fCutoff, /** CHECKME: Cutoff frequenc mass1, /** Mass 1 */
y? */ mass2, /** Mass 2 */
mass1, /** Mass 1 */
mass2, /** Mass 2 */
LAL_SIM_INSPIRAL_INTERACTION_ALL_SPIN, /** Spin interaction */ LAL_SIM_INSPIRAL_INTERACTION_ALL_SPIN, /** Spin interaction */
phaseO /** PN Order in Phase */ )) phaseO /** PN Order in Phase */ ))
XLAL_ERROR(XLAL_EFUNC); XLAL_ERROR(XLAL_EFUNC);
/* Check that initial frequency is smaller than omegamatch ~ xxyy for m=1 00 Msun */ /* Check that initial frequency is smaller than omegamatch ~ xxyy for m=1 00 Msun */
initphi = phi_start; initphi = phi_start;
initomega = f_min*unitHz; initomega = f_min*unitHz;
/* Check that initial frequency is smaller than omegamatch ~ xxyy for m=1 00 Msun */ /* Check that initial frequency is smaller than omegamatch ~ xxyy for m=1 00 Msun */
skipping to change at line 2252 skipping to change at line 2245
omegaMatch = OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2); omegaMatch = OmMatch(LNhS1,LNhS2,S1S1,S1S2,S2S2);
if ( initomega > omegaMatch ) { if ( initomega > omegaMatch ) {
/*if ((params->spin1[0]==params->spin1[1])&&(params->spin1[1]==params-> spin2[0])&&(params->spin2[0]==params->spin2[1])&&(params->spin2[1]==0.)) { /*if ((params->spin1[0]==params->spin1[1])&&(params->spin1[1]==params-> spin2[0])&&(params->spin2[0]==params->spin2[1])&&(params->spin2[1]==0.)) {
//Beware, this correspond to a shift of the initial phase! //Beware, this correspond to a shift of the initial phase!
initomega = 0.95*omegaMatch; initomega = 0.95*omegaMatch;
fprintf(stdout,"*** LALPSpinInspiralRD WARNING ***: Initial frequency reset from %12.6e to %12.6e Hz, m:(%12.4e,%12.4e)\n",params->fLower,initom ega/unitHz,params->mass1,params->mass2); fprintf(stdout,"*** LALPSpinInspiralRD WARNING ***: Initial frequency reset from %12.6e to %12.6e Hz, m:(%12.4e,%12.4e)\n",params->fLower,initom ega/unitHz,params->mass1,params->mass2);
}*/ }*/
/*else {*/ /*else {*/
XLALPrintError("**** LALPSpinInspiralRD ERROR ****: Initial frequency t XLALPrintError("**** LALPSpinInspiralRD ERROR ****: the product of init
oo high: %11.5e for omM ~ %11.5e and m:(%8.3f, %8.3f)\n",f_min,omegaMatch/u ial frequency times masses is too high: %11.5e for omM ~ %11.5e\n",f_min*ma
nitHz,mass1,mass2); ss*LAL_PI,omegaMatch);
XLALPrintError(" please consider dec
reasing initial freq %8.3f Hz or m:(%8.3f, %8.3f) Msun\n",f_min,mass1,mass2
);
XLAL_ERROR(XLAL_EFAILED); XLAL_ERROR(XLAL_EFAILED);
/*}*/ /*}*/
} }
/* Here we use the following convention: /* Here we use the following convention:
the coordinates of the spin vectors spin1,2 and the inclination the coordinates of the spin vectors spin1,2 and the inclination
variable refers to different physical parameters according to the valu e of variable refers to different physical parameters according to the valu e of
axisChoice: axisChoice:
* OrbitalL: inclination denotes the angle between the view direction * OrbitalL: inclination denotes the angle between the view direction
skipping to change at line 2532 skipping to change at line 2526
errcode = XLALSpinInspiralAdaptiveEngine(neqs,yinit,amp22ini,&mparams,h2P 2,h2M2,h2P1,h2M1,h20,h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h30,h4P4,h4M4,h4P3,h4M3, h4P2,h4M2,h4P1,h4M1,h40,fap,phap,&phenPars); errcode = XLALSpinInspiralAdaptiveEngine(neqs,yinit,amp22ini,&mparams,h2P 2,h2M2,h2P1,h2M1,h20,h3P3,h3M3,h3P2,h3M2,h3P1,h3M1,h30,h4P4,h4M4,h4P3,h4M3, h4P2,h4M2,h4P1,h4M1,h40,fap,phap,&phenPars);
if(errcode) XLAL_ERROR(XLAL_EFUNC); if(errcode) XLAL_ERROR(XLAL_EFUNC);
intreturn=phenPars.intreturn; intreturn=phenPars.intreturn;
count= phenPars.countback; count= phenPars.countback;
/* report on abnormal termination: /* report on abnormal termination:
Termination is fine if omegamatch is passed or if energy starts Termination is fine if omegamatch is passed or if energy starts
increasing */ increasing */
if ( (intreturn!=LALPSIRDPN_TEST_OMEGACUT) && (intreturn != LALPSIRDPN_TE ST_OMEGAMATCH) && (intreturn != LALPSIRDPN_TEST_ENERGY) ) if ( (intreturn != LALPSIRDPN_TEST_OMEGAMATCH) && (intreturn != LALPSIRDP N_TEST_ENERGY) )
{ {
XLALPrintWarning("** LALPSpinInspiralRD WARNING **: integration termi nated with code %d.\n",intreturn); XLALPrintWarning("** LALPSpinInspiralRD WARNING **: integration termi nated with code %d.\n",intreturn);
fprintf(stderr," 1025: Energy increases\n 1026: Omegadot -ve\n 102 8: Omega NAN\n 1029: Omega > Omegamatch\n 1031: Omega -ve\n 1032: Omega > OmegaCut %12.6e\n",mparams.OmCutoff); fprintf(stderr," 1025: Energy increases\n 1026: Omegadot -ve\n 102 8: Omega NAN\n 1029: Omega > Omegamatch\n 1031: Omega -ve\n");
XLALPrintWarning(" Waveform parameters were m1 = %14.6e, m2 = %14.6e , inc = %10.6f,\n", mass1, mass2, iota); XLALPrintWarning(" Waveform parameters were m1 = %14.6e, m2 = %14.6e , inc = %10.6f,\n", mass1, mass2, iota);
XLALPrintWarning(" S1 = (%10.6f,%10.6f,%10. 6f)\n", s1x, s1y, s1z); XLALPrintWarning(" S1 = (%10.6f,%10.6f,%10. 6f)\n", s1x, s1y, s1z);
XLALPrintWarning(" S2 = (%10.6f,%10.6f,%10. 6f)\n", s2x, s2y, s2z); XLALPrintWarning(" S2 = (%10.6f,%10.6f,%10. 6f)\n", s2x, s2y, s2z);
} }
if (intreturn==LALPSIRDPN_TEST_OMEGAMATCH) { if (intreturn==LALPSIRDPN_TEST_OMEGAMATCH) {
tim = t0 = phenPars.endtime; tim = t0 = phenPars.endtime;
tAs = t0 + 2. * phenPars.domega / phenPars.ddomega; tAs = t0 + 2. * phenPars.domega / phenPars.ddomega;
om1 = phenPars.domega * tAs * (1. - t0 / tAs) * (1. - t0 / tAs); om1 = phenPars.domega * tAs * (1. - t0 / tAs) * (1. - t0 / tAs);
 End of changes. 40 change blocks. 
73 lines changed or deleted 72 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/