Commit aff93c45 authored by LING SUN's avatar LING SUN Committed by Karl Wette
Browse files

Allow CW codes to accept arbitrary aPlus/aCross amplitudes

parent 5991066a
......@@ -83,6 +83,7 @@
#define TRUE (1==1)
#define FALSE (1==0)
#define SQ(x) ( (x) * (x) )
/*----- SWITCHES -----*/
/*----- Macros -----*/
......@@ -2098,15 +2099,40 @@ write_PulsarCandidate_to_fp ( FILE *fp, const PulsarCandidate *pulsarParams, co
fprintf (fp, "\n");
/* Amplitude parameters with error-estimates */
fprintf (fp, "h0 = % .6g;\n", pulsarParams->Amp.h0 );
fprintf (fp, "dh0 = % .6g;\n", pulsarParams->dAmp.h0 );
fprintf (fp, "cosi = % .6g;\n", pulsarParams->Amp.cosi );
fprintf (fp, "dcosi = % .6g;\n", pulsarParams->dAmp.cosi );
fprintf (fp, "aPlus = % .6g;\n", pulsarParams->Amp.aPlus );
fprintf (fp, "daPlus = % .6g;\n", pulsarParams->dAmp.aPlus );
fprintf (fp, "aCross = % .6g;\n", pulsarParams->Amp.aCross );
fprintf (fp, "daCross = % .6g;\n", pulsarParams->dAmp.aCross );
fprintf (fp, "phi0 = % .6g;\n", pulsarParams->Amp.phi0 );
fprintf (fp, "dphi0 = % .6g;\n", pulsarParams->dAmp.phi0 );
fprintf (fp, "psi = % .6g;\n", pulsarParams->Amp.psi );
fprintf (fp, "dpsi = % .6g;\n", pulsarParams->dAmp.psi );
/* To allow arbitrary aPlus/aCross amplitudes, we have switched to output-variables A^\nu = (aPlus, aCross, phi0, psi)
* and hence using a different Jacobian matrix. The dh0 and dcosi recovered below are only valid when GW is emitted at
* twice the spin frequency, and are derived from dA+, dAx. Since both the current and original calculations ignore
* off-diagonal items in computing the errors, the dh0 and dcosi below are expected to be slightly different from the
* original estimate when the output-variables are B^\nu = (h0, cosi, phi0, psi). Since they are only used for sanity
* checks, it should not impact any usage. */
if ( fabs(pulsarParams->Amp.aPlus) >= fabs(pulsarParams->Amp.aCross) ){ /* Assume GW at twice the spin frequency only */
REAL8 h0, dh0, cosi, dcosi;
h0 = pulsarParams->Amp.aPlus + sqrt( SQ(pulsarParams->Amp.aPlus) - SQ(pulsarParams->Amp.aCross) );
if (h0 > 0)
cosi = pulsarParams->Amp.aCross/h0;
else
cosi = 0;
dh0 = (pulsarParams->Amp.aPlus + pulsarParams->dAmp.aPlus) + sqrt( SQ(pulsarParams->Amp.aPlus + pulsarParams->dAmp.aPlus) - SQ(fabs(pulsarParams->Amp.aCross) + pulsarParams->dAmp.aCross) ) - h0;
if ( (h0+dh0) > 0 )
dcosi = fabs( (fabs(pulsarParams->Amp.aCross) + pulsarParams->dAmp.aCross) / (h0 + dh0) - fabs(cosi) );
else
dcosi = 0;
fprintf (fp, "h0 = % .6g;\n", h0 );
fprintf (fp, "dh0 = % .6g;\n", dh0 );
fprintf (fp, "cosi = % .6g;\n", cosi );
fprintf (fp, "dcosi = % .6g;\n", dcosi );
}
fprintf (fp, "\n");
/* Doppler parameters */
......
......@@ -325,13 +325,13 @@ InitPFS ( ConfigVariables *cfg, UserInput_t *uvar )
/* ----- internally we always use h0, cosi */
if ( have_h0 )
{
cfg->pap.h0=uvar->h0;
cfg->pap.cosi=uvar->cosi;
cfg->pap.aPlus = 0.5 * uvar->h0 * (1.0 + SQ(uvar->cosi));
cfg->pap.aCross = uvar->h0 * uvar->cosi;
}
else
{
cfg->pap.h0 = uvar->aPlus + sqrt( SQ( uvar->aPlus ) - SQ( uvar->aCross ) );
cfg->pap.cosi= uvar->aCross / cfg->pap.h0;
cfg->pap.aPlus = uvar->aPlus;
cfg->pap.aCross = uvar->aCross;
}
cfg->pap.psi=uvar->psi;
cfg->pap.phi0=uvar->phi0;
......
......@@ -482,6 +482,7 @@ XLALsynthesizeSignals ( gsl_matrix **A_Mu_i, /**< [OUT] list of numDraws 4D lin
REAL8 phi0Max = phi0Min + AmpRange.phi0Band;
REAL8 SNR = AmpRange.SNR;
REAL8 res_rho2;
REAL8 h0, cosi;
int gslstat;
......@@ -536,8 +537,10 @@ XLALsynthesizeSignals ( gsl_matrix **A_Mu_i, /**< [OUT] list of numDraws 4D lin
{
PulsarAmplitudeParams Amp;
Amp.h0 = gsl_ran_flat ( rng, h0NatMin, h0NatMax );
Amp.cosi = gsl_ran_flat ( rng, cosiMin, cosiMax );
h0 = gsl_ran_flat ( rng, h0NatMin, h0NatMax );
cosi = gsl_ran_flat ( rng, cosiMin, cosiMax );
Amp.aPlus = 0.5 * h0 * (1.0 + SQ(cosi));
Amp.aCross = h0 * cosi;
Amp.psi = gsl_ran_flat ( rng, psiMin, psiMax );
Amp.phi0 = gsl_ran_flat ( rng, phi0Min, phi0Max );
......@@ -580,7 +583,8 @@ XLALsynthesizeSignals ( gsl_matrix **A_Mu_i, /**< [OUT] list of numDraws 4D lin
/* if specified SNR: rescale signal to this SNR */
if ( SNR > 0 ) {
REAL8 rescale_h0 = SNR / sqrt ( res_rho2 );
Amp.h0 *= rescale_h0;
Amp.aPlus *= rescale_h0;
Amp.aCross *= rescale_h0;
res_rho2 = SQ(SNR);
gsl_vector_scale ( &A_Mu.vector, rescale_h0);
gsl_vector_scale ( &s_mu.vector, rescale_h0);
......@@ -588,8 +592,8 @@ XLALsynthesizeSignals ( gsl_matrix **A_Mu_i, /**< [OUT] list of numDraws 4D lin
gsl_vector_set ( *rho2_i, row, res_rho2 );
gsl_matrix_set ( *Amp_i, row, 0, Amp.h0 );
gsl_matrix_set ( *Amp_i, row, 1, Amp.cosi );
gsl_matrix_set ( *Amp_i, row, 0, h0 );
gsl_matrix_set ( *Amp_i, row, 1, cosi );
gsl_matrix_set ( *Amp_i, row, 2, Amp.psi );
gsl_matrix_set ( *Amp_i, row, 3, Amp.phi0 );
......
......@@ -255,12 +255,8 @@ main(int argc, char *argv[])
params.pulsar.position.system = COORDINATESYSTEM_EQUATORIAL;
params.pulsar.position.longitude = GV.pulsar.Doppler.Alpha;
params.pulsar.position.latitude = GV.pulsar.Doppler.Delta;
{
REAL8 h0 = GV.pulsar.Amp.h0;
REAL8 cosi = GV.pulsar.Amp.cosi;
params.pulsar.aPlus = 0.5 * h0 * ( 1.0 + SQ(cosi) );
params.pulsar.aCross = h0 * cosi;
}
params.pulsar.aPlus = GV.pulsar.Amp.aPlus;
params.pulsar.aCross = GV.pulsar.Amp.aCross;
params.pulsar.phi0 = GV.pulsar.Amp.phi0;
params.pulsar.psi = GV.pulsar.Amp.psi;
......@@ -659,26 +655,21 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
XLAL_ERROR ( XLAL_EINVAL, "Need BOTH --aPlus and --aCross !\n\n");
}
if ( have_h0 && have_cosi )
{
cfg->pulsar.Amp.h0 = uvar->h0;
cfg->pulsar.Amp.cosi = uvar->cosi;
{ /* translate {h_0, cosi} into A_{+,x} */
/* assume at 2f */
REAL8 h0 = uvar->h0;
REAL8 cosi = uvar->cosi;
cfg->pulsar.Amp.aPlus = 0.5 * h0 * (1.0 + SQ(cosi));
cfg->pulsar.Amp.aCross = h0 * cosi;
}
else if ( have_aPlus && have_aCross )
{ /* translate A_{+,x} into {h_0, cosi} */
REAL8 disc;
if ( fabs(uvar->aCross) > uvar->aPlus ) {
XLAL_ERROR ( XLAL_EINVAL, "Invalid input parameters: |aCross| = %g must be <= than aPlus = %g.\n", fabs(uvar->aCross), uvar->aPlus );
}
disc = sqrt ( SQ(uvar->aPlus) - SQ(uvar->aCross) );
cfg->pulsar.Amp.h0 = uvar->aPlus + disc;
if ( cfg->pulsar.Amp.h0 > 0 )
cfg->pulsar.Amp.cosi = uvar->aCross / cfg->pulsar.Amp.h0; // avoid division by 0!
else
cfg->pulsar.Amp.cosi = 0;
{
cfg->pulsar.Amp.aPlus = uvar->aPlus;
cfg->pulsar.Amp.aCross = uvar->aCross;
}
else {
cfg->pulsar.Amp.h0 = 0.0;
cfg->pulsar.Amp.cosi = 0.0;
cfg->pulsar.Amp.aPlus = 0.0;
cfg->pulsar.Amp.aCross = 0.0;
}
cfg->pulsar.Amp.phi0 = uvar->phi0;
cfg->pulsar.Amp.psi = uvar->psi;
......@@ -1251,8 +1242,8 @@ XLALInitUserVars ( UserVariables_t *uvar, int argc, char *argv[] )
XLALRegisterUvarMember( Alpha, RAJ, 0, OPTIONAL, "Sky: equatorial J2000 right ascension (in radians or hours:minutes:seconds)");
XLALRegisterUvarMember( Delta, DECJ, 0, OPTIONAL, "Sky: equatorial J2000 declination (in radians or degrees:minutes:seconds)");
XLALRegisterUvarMember( h0, REAL8, 0, OPTIONAL, "Overall signal-amplitude h0");
XLALRegisterUvarMember( cosi, REAL8, 0, OPTIONAL, "cos(iota) of inclination-angle iota");
XLALRegisterUvarMember( h0, REAL8, 0, OPTIONAL, "Overall signal-amplitude h0 (for emission at twice spin frequency only)");
XLALRegisterUvarMember( cosi, REAL8, 0, OPTIONAL, "cos(iota) of inclination-angle iota (for emission at twice spin frequency only)");
XLALRegisterUvarMember( aPlus, REAL8, 0, OPTIONAL, "ALTERNATIVE to {--h0,--cosi}: A_+ amplitude");
XLALRegisterUvarMember( aCross, REAL8, 0, OPTIONAL, "ALTERNATIVE to {--h0,--cosi}: A_x amplitude");
......
......@@ -450,8 +450,9 @@ XLALInitCode ( ConfigVariables *cfg, const UserVariables_t *uvar, const char *ap
}
/* ----- get parameter-space point from user-input) */
cfg->signalParams.Amp.h0 = uvar->h0;
cfg->signalParams.Amp.cosi = uvar->cosi;
cfg->signalParams.Amp.aPlus = 0.5 * uvar->h0 * (1.0 + SQ(uvar->cosi));
cfg->signalParams.Amp.aCross = uvar->h0 * uvar->cosi;
cfg->signalParams.Amp.psi = uvar->psi;
cfg->signalParams.Amp.phi0 = uvar->phi0;
......@@ -585,11 +586,21 @@ XLALOutputDopplerMetric ( FILE *fp, const DopplerPhaseMetric *Pmetric, const Dop
}
fprintf ( fp, "%%%% DetectorMotionType = '%s'\n", XLALDetectorMotionName(meta->detMotionType) );
fprintf ( fp, "h0 = %g;\ncosi = %g;\npsi = %g;\nphi0 = %g;\n", Amp->h0, Amp->cosi, Amp->psi, Amp->phi0 );
fprintf ( fp, "aPlus = %g;\naCross = %g;\npsi = %g;\nphi0 = %g;\n", Amp->aPlus, Amp->aCross, Amp->psi, Amp->phi0 );
if ( fabs(Amp->aPlus) >= fabs(Amp->aCross) ){ /* Assume GW at twice the spin frequency only */
REAL8 h0, cosi;
h0 = Amp->aPlus + sqrt( SQ(Amp->aPlus) - SQ(Amp->aCross) );
if (h0 > 0)
cosi = Amp->aCross/h0;
else
cosi = 0;
fprintf ( fp, "h0 = %g;\ncosi = %g;\n", h0, cosi);
}
fprintf ( fp, "%%%% DopplerPoint = {\n");
fprintf ( fp, "refTime = %.1f;\n", XLALGPSGetREAL8 ( &doppler->refTime ) );
fprintf ( fp, "Alpha = %f;\nDelta = %f;\n", doppler->Alpha, doppler->Delta );
fprintf ( fp, "fkdot = [%f, %g, %g, %g ];\n", doppler->fkdot[0], doppler->fkdot[1], doppler->fkdot[2], doppler->fkdot[3] );
if ( doppler->asini > 0 )
{
fprintf ( fp, "%%%% orbit = { \n");
......
......@@ -283,8 +283,8 @@ REAL4VectorAligned * coherentlyAddSFTs(const MultiSFTVector *multiSFTvector, con
PulsarParams XLAL_INIT_DECL(pulsarParams);
LIGOTimeGPS SSBepoch = multissb->data[0]->refTime, spinEpoch;
XLALGPSAdd(&SSBepoch, multissb->data[0]->DeltaT->data[0]-params->SFToverlap);
pulsarParams.Amp.h0 = 0.0; //arbitrary, we only care about frequency
pulsarParams.Amp.cosi = 0.0; //arbitrary, we only care about frequency
pulsarParams.Amp.aPlus = 0.0; //arbitrary, we only care about frequency
pulsarParams.Amp.aCross = 0.0; //arbitrary, we only care about frequency
pulsarParams.Amp.psi = 0.0; //arbitrary, we only care about frequency
pulsarParams.Amp.phi0 = 0.0; //arbitrary, we only care about frequency
if (NSparams->assumeNSrefTime == NULL) {
......@@ -637,13 +637,10 @@ REAL4TimeSeries * computeNSfreqTS(const PulsarParams *pulsarParams, LIGOTimeGPS
{
XLAL_CHECK_NULL( pulsarParams != NULL && duration > 0, XLAL_EINVAL );
REAL8 h0 = pulsarParams->Amp.h0;
REAL8 cosi = pulsarParams->Amp.cosi;
SpinOrbitCWParamStruc XLAL_INIT_DECL(sourceParams);
sourceParams.psi = pulsarParams->Amp.psi;
sourceParams.aPlus = 0.5*h0*(1.0+cosi*cosi);
sourceParams.aCross = h0*cosi;
sourceParams.aPlus = pulsarParams->Amp.aPlus;
sourceParams.aCross = pulsarParams->Amp.aCross;
sourceParams.phi0 = pulsarParams->Amp.phi0;
sourceParams.f0 = pulsarParams->Doppler.fkdot[0];
sourceParams.position.latitude = pulsarParams->Doppler.Delta;
......@@ -1328,7 +1325,10 @@ MultiSFTVector * generateSFTdata(UserInput_t *uvar, const MultiLALDetector *dete
XLAL_CHECK_NULL ( (marginalizedSignalData = XLALCreateREAL8Vector(multiSFTvector->data[0]->data->data->length)) != NULL, XLAL_EFUNC );
memset(marginalizedSignalData->data, 0, sizeof(REAL8)*marginalizedSignalData->length);
for (UINT4 jj=0; jj<300; jj++) {
oneSignal->data[0].Amp.cosi = 2.0*gsl_rng_uniform(rng) - 1.0;
REAL8 cosi = 2.0*gsl_rng_uniform(rng) - 1.0;
REAL8 h0 = 1.0;
oneSignal->data[0].Amp.aPlus = 0.5 * h0 * (1.0 + cosi * cosi);
oneSignal->data[0].Amp.aCross = h0 * cosi;
oneSignal->data[0].Amp.psi = LAL_TWOPI*gsl_rng_uniform(rng);
oneSignal->data[0].Amp.phi0 = LAL_TWOPI*gsl_rng_uniform(rng);
oneSignal->data[0].Doppler.argp = LAL_TWOPI*gsl_rng_uniform(rng);
......
......@@ -393,10 +393,8 @@ XLALGenerateCWSignalTS ( const PulsarParams *pulsarParams, ///< input CW pulsar-
XLAL_CHECK_NULL ( fHet >= 0, XLAL_EDOM );
// translate amplitude params
REAL8 h0 = pulsarParams->Amp.h0;
REAL8 cosi = pulsarParams->Amp.cosi;
REAL8 aPlus = 0.5 * h0 * ( 1.0 + SQ(cosi) );
REAL8 aCross = h0 * cosi;
REAL8 aPlus = pulsarParams->Amp.aPlus;
REAL8 aCross = pulsarParams->Amp.aCross;
// translate 'modern' fkdot into 'old-style' spindown-vector
UINT4 s_max;
for ( s_max = PULSAR_MAX_SPINS-1; s_max > 0; s_max -- )
......@@ -825,21 +823,13 @@ XLALReadPulsarParams ( PulsarParams *pulsarParams, ///< [out] pulsar parameters
// {h0,cosi} or {aPlus, aCross} mutually exclusive sets
XLAL_CHECK ( ! ( have_h0 && have_aPlus ), XLAL_EINVAL );
if ( have_aPlus ) /* translate A_{+,x} into {h_0, cosi} */
if ( have_h0 ) /* translate {h_0, cosi} into A_{+,x}*/
{
XLAL_CHECK ( fabs ( aCross ) <= aPlus, XLAL_EDOM, "ERROR: |aCross| (= %g) must be <= aPlus (= %g).\n", fabs(aCross), aPlus );
REAL8 disc = sqrt ( SQ(aPlus) - SQ(aCross) );
h0 = aPlus + disc;
if ( h0 > 0 ) {
cosi = aCross / h0; // avoid division by 0!
}
} //if {aPlus, aCross}
XLAL_CHECK ( h0 >= 0, XLAL_EDOM );
pulsarParams->Amp.h0 = h0;
XLAL_CHECK ( (cosi >= -1) && (cosi <= 1), XLAL_EDOM );
pulsarParams->Amp.cosi= cosi;
XLAL_CHECK ( h0 >= 0, XLAL_EDOM );
XLAL_CHECK ( (cosi >= -1) && (cosi <= 1), XLAL_EDOM );
pulsarParams->Amp.aPlus = 0.5 * h0 * (1.0 + SQ(cosi));
pulsarParams->Amp.aCross = h0 * cosi;
} //if {h0, cosi}
// ----- psi
REAL8 psi = 0; BOOLEAN have_psi;
......
......@@ -67,7 +67,6 @@ XLALEstimatePulsarAmplitudeParams ( PulsarCandidate *pulsarParams, ///< [in,out]
REAL8 beta;
REAL8 phi0, psi;
REAL8 b1, b2, b3;
REAL8 h0, cosi;
REAL8 cosphi0, sinphi0, cos2psi, sin2psi;
......@@ -188,42 +187,33 @@ XLALEstimatePulsarAmplitudeParams ( PulsarCandidate *pulsarParams, ///< [in,out]
__func__, tolerance );
}
// translate A_{+,x} into {h_0, cosi}
h0 = aPlus + sqrt ( disc ); // not yet normalized !
cosi = aCross / h0;
// ========== Estimate the errors ==========
// ----- compute derivatives \partial A^\mu / \partial B^\nu, where
// we consider the output-variables B^\nu = (h0, cosi, phi0, psi)
// ----- compute derivatives \partial A^\mu / \partial A^\nu, where
// we consider the output-variables A^\nu = (aPlus, aCross, phi0, psi)
// where aPlus = 0.5 * h0 * (1 + cosi^2) and aCross = h0 * cosi
{ // Ahat^mu is defined as A^mu with the replacements: A_+ --> A_x, and A_x --> h0
REAL8 A1hat = aCross * cosphi0 * cos2psi - h0 * sinphi0 * sin2psi;
REAL8 A2hat = aCross * cosphi0 * sin2psi + h0 * sinphi0 * cos2psi;
REAL8 A3hat = - aCross * sinphi0 * cos2psi - h0 * cosphi0 * sin2psi;
REAL8 A4hat = - aCross * sinphi0 * sin2psi + h0 * cosphi0 * cos2psi;
{ // Eqn (114) in T0900149 v5
// ----- A1 = aPlus * cosphi0 * cos2psi - aCross * sinphi0 * sin2psi; -----
gsl_matrix_set (Jh_Mu_nu, 0, 0, A1h / h0 ); /* dA1/h0 */
gsl_matrix_set (Jh_Mu_nu, 0, 1, A1hat ); /* dA1/dcosi */
gsl_matrix_set (Jh_Mu_nu, 0, 0, cosphi0 * cos2psi ); /* dA1/daPlus */
gsl_matrix_set (Jh_Mu_nu, 0, 1, - sinphi0 * sin2psi ); /* dA1/daCross */
gsl_matrix_set (Jh_Mu_nu, 0, 2, A3h ); /* dA1/dphi0 */
gsl_matrix_set (Jh_Mu_nu, 0, 3, - 2.0 * A2h ); /* dA1/dpsi */
// ----- A2 = aPlus * cosphi0 * sin2psi + aCross * sinphi0 * cos2psi; -----
gsl_matrix_set (Jh_Mu_nu, 1, 0, A2h / h0 ); /* dA2/h0 */
gsl_matrix_set (Jh_Mu_nu, 1, 1, A2hat ); /* dA2/dcosi */
gsl_matrix_set (Jh_Mu_nu, 1, 0, cosphi0 * sin2psi ); /* dA2/daPlus */
gsl_matrix_set (Jh_Mu_nu, 1, 1, sinphi0 * cos2psi ); /* dA2/daCross */
gsl_matrix_set (Jh_Mu_nu, 1, 2, A4h ); /* dA2/dphi0 */
gsl_matrix_set (Jh_Mu_nu, 1, 3, 2.0 * A1h ); /* dA2/dpsi */
// ----- A3 = - aPlus * sinphi0 * cos2psi - aCross * cosphi0 * sin2psi; -----
gsl_matrix_set (Jh_Mu_nu, 2, 0, A3h / h0 ); /* dA3/h0 */
gsl_matrix_set (Jh_Mu_nu, 2, 1, A3hat ); /* dA3/dcosi */
gsl_matrix_set (Jh_Mu_nu, 2, 0, - sinphi0 * cos2psi ); /* dA3/daPlus */
gsl_matrix_set (Jh_Mu_nu, 2, 1, - cosphi0 * sin2psi ); /* dA3/daCross */
gsl_matrix_set (Jh_Mu_nu, 2, 2, - A1h ); /* dA3/dphi0 */
gsl_matrix_set (Jh_Mu_nu, 2, 3, - 2.0 * A4h ); /* dA3/dpsi */
// ----- A4 = - aPlus * sinphi0 * sin2psi + aCross * cosphi0 * cos2psi; -----
gsl_matrix_set (Jh_Mu_nu, 3, 0, A4h / h0 ); /* dA4/h0 */
gsl_matrix_set (Jh_Mu_nu, 3, 1, A4hat ); /* dA4/dcosi */
gsl_matrix_set (Jh_Mu_nu, 3, 0, - sinphi0 * sin2psi ); /* dA4/daPlus */
gsl_matrix_set (Jh_Mu_nu, 3, 1, cosphi0 * cos2psi ); /* dA4/daCross */
gsl_matrix_set (Jh_Mu_nu, 3, 2, - A2h ); /* dA4/dphi0 */
gsl_matrix_set (Jh_Mu_nu, 3, 3, 2.0 * A3h ); /* dA4/dpsi */
}
......@@ -258,14 +248,14 @@ XLALEstimatePulsarAmplitudeParams ( PulsarCandidate *pulsarParams, ///< [in,out]
XLAL_CHECK( XLALExtrapolatePulsarPhase( &phi0, pulsarParams->Doppler.fkdot, phi0, dtau ) == XLAL_SUCCESS, XLAL_EFUNC );
// fill candidate-struct with the obtained signal-parameters and error-estimations
pulsarParams->Amp.h0 = normAmu * h0;
pulsarParams->Amp.cosi = cosi;
pulsarParams->Amp.aPlus = normAmu * aPlus;
pulsarParams->Amp.aCross = normAmu * aCross;
pulsarParams->Amp.phi0 = phi0;
pulsarParams->Amp.psi = psi;
// read out principal estimation-errors from diagonal elements of inverse Fisher-matrix
pulsarParams->dAmp.h0 = normAmu * sqrt( gsl_matrix_get (Jh_Mu_nu, 0, 0 ) );
pulsarParams->dAmp.cosi = sqrt( gsl_matrix_get (Jh_Mu_nu, 1, 1 ) );
pulsarParams->dAmp.aPlus = normAmu * sqrt( gsl_matrix_get (Jh_Mu_nu, 0, 0 ) );
pulsarParams->dAmp.aCross = normAmu * sqrt( gsl_matrix_get (Jh_Mu_nu, 1, 1 ) );
pulsarParams->dAmp.phi0 = sqrt( gsl_matrix_get (Jh_Mu_nu, 2, 2 ) );
pulsarParams->dAmp.psi = sqrt( gsl_matrix_get (Jh_Mu_nu, 3, 3 ) );
// return also the full Amplitude-Fisher matrix: invert Jh_Mu_nu
......@@ -296,8 +286,8 @@ XLALAmplitudeParams2Vect ( PulsarAmplitudeVect A_Mu, ///< [out] Canonical ampli
)
{
REAL8 aPlus = 0.5 * Amp.h0 * ( 1.0 + SQ(Amp.cosi) );
REAL8 aCross = Amp.h0 * Amp.cosi;
REAL8 aPlus = Amp.aPlus;
REAL8 aCross = Amp.aCross;
REAL8 cos2psi = cos ( 2.0 * Amp.psi );
REAL8 sin2psi = sin ( 2.0 * Amp.psi );
REAL8 cosphi0 = cos ( Amp.phi0 );
......@@ -324,7 +314,7 @@ XLALAmplitudeVect2Params ( PulsarAmplitudeParams *Amp, ///< [out] Physical ampl
)
{
REAL8 h0Ret, cosiRet, psiRet, phi0Ret;
REAL8 psiRet, phi0Ret;
REAL8 A1, A2, A3, A4, Asq, Da, disc;
REAL8 Ap2, Ac2, aPlus, aCross;
......@@ -365,9 +355,6 @@ XLALAmplitudeVect2Params ( PulsarAmplitudeParams *Amp, ///< [out] Physical ampl
phi0Ret += LAL_PI;
}
h0Ret = aPlus + sqrt ( disc );
cosiRet = aCross / h0Ret;
// make unique by fixing the gauge to be psi in [-pi/4, pi/4], phi0 in [0, 2*pi]
while ( psiRet > LAL_PI_4 ) {
psiRet -= LAL_PI_2;
......@@ -386,8 +373,8 @@ XLALAmplitudeVect2Params ( PulsarAmplitudeParams *Amp, ///< [out] Physical ampl
}
// Return final answer
Amp->h0 = h0Ret;
Amp->cosi = cosiRet;
Amp->aPlus = aPlus;
Amp->aCross = aCross;
Amp->psi = psiRet;
Amp->phi0 = phi0Ret;
......@@ -404,25 +391,25 @@ XLALAmplitudeVect2Params ( PulsarAmplitudeParams *Amp, ///< [out] Physical ampl
* E[2F] = 4 + SNR^2 for a template perfectly matching the signal, ie SNR^2 = (signal|signal).
*/
REAL8
XLALComputeOptimalSNR2FromMmunu ( const PulsarAmplitudeParams pap, /**< [in] PulsarAmplitudeParameter {h0, cosi, psi, phi0} */
XLALComputeOptimalSNR2FromMmunu ( const PulsarAmplitudeParams pap, /**< [in] PulsarAmplitudeParameter {aPlus, aCross, psi, phi0} */
const AntennaPatternMatrix Mmunu /**< [in] Antenna-pattern Matrix */
)
{
// Calculate the SNR using the expressions from the 'CFSv2 notes' T0900149-v5
// https://dcc.ligo.org/LIGO-T0900149-v5/public, namely Eqs.(77)-(80), keeping in mind
// that Mmunu.{Ad,Bd,Cd} = Nsft * {A,B,C}, and Tdata = Nsft * Tsft
REAL8 cos2psi = cos ( 2.0 * pap.psi );
REAL8 sin2psi = sin ( 2.0 * pap.psi );
REAL8 cos2psiSq = SQ( cos2psi );
REAL8 sin2psiSq = SQ( sin2psi );
REAL8 etaSq = SQ(pap.cosi);
REAL8 al1 = 0.25 * SQ( 1.0 + etaSq ) * cos2psiSq + etaSq * sin2psiSq; // Eq.(78)
REAL8 al2 = 0.25 * SQ( 1.0 + etaSq ) * sin2psiSq + etaSq * cos2psiSq; // Eq.(79)
REAL8 al3 = 0.25 * SQ( 1.0 - etaSq ) * sin2psi * cos2psi; // Eq.(80)
REAL8 h0sq_al1 = SQ( pap.aPlus ) * cos2psiSq + SQ( pap.aCross ) * sin2psiSq; // Eq.(78) generalised to (aPlus,aCross)
REAL8 h0sq_al2 = SQ( pap.aPlus ) * sin2psiSq + SQ( pap.aCross ) * cos2psiSq; // Eq.(79) generalised to (aPlus,aCross)
REAL8 h0sq_al3 = ( SQ( pap.aPlus ) - SQ( pap.aCross ) ) * sin2psi * cos2psi; // Eq.(80) generalised to (aPlus,aCross)
// SNR^2: Eq.(77)
REAL8 rho2 = SQ(pap.h0) * ( al1 * Mmunu.Ad + al2 * Mmunu.Bd + 2.0 * al3 * Mmunu.Cd ) * Mmunu.Sinv_Tsft;
REAL8 rho2 = ( h0sq_al1 * Mmunu.Ad + h0sq_al2 * Mmunu.Bd + 2.0 * h0sq_al3 * Mmunu.Cd ) * Mmunu.Sinv_Tsft;
return rho2;
......
......@@ -56,10 +56,10 @@ extern "C" {
/** Type containing the JKS 'amplitude parameters' {h0, cosi, phi0, psi} */
typedef struct tagPulsarAmplitudeParams {
REAL8 h0; /**< overall signal amplitude */
REAL8 cosi; /**< cos(iota) of inclination angle iota of spin-axis wrt line of sight */
REAL8 psi; /**< polarization angle psi */
REAL8 phi0; /**< initial signal-phase (at some reference time) */
REAL8 aPlus; /**< Signal amplitude (plus) */
REAL8 aCross; /**< Signal amplitude (cross) */
} PulsarAmplitudeParams;
/** Struct for 'canonical' coordinates in amplitude-params space A^mu = {A1, A2, A3, A4} */
......
......@@ -537,6 +537,8 @@ XLALSynthesizeTransientAtoms ( InjParams_t *injParamsOut, /**< [out] return su
const MultiNoiseWeights *multiNoiseWeights /** [in] per-detector noise weights SX^-1/S^-1, no per-SFT variation (can be NULL for unit weights) */
)
{
REAL8 h0, cosi;
/* check input */
XLAL_CHECK_NULL ( rng && multiAMBuffer && multiDetStates, XLAL_EINVAL, "Invalid NULL input!\n");
XLAL_CHECK_NULL ( !multiNoiseWeights || multiNoiseWeights->data, XLAL_EINVAL, "Invalid NULL input for multiNoiseWeights->data!\n" );
......@@ -580,13 +582,15 @@ XLALSynthesizeTransientAtoms ( InjParams_t *injParamsOut, /**< [out] return su
/* ----- draw amplitude vector A^mu from given ranges in {h0, cosi, psi, phi0} */
PulsarAmplitudeParams Amp;
if ( AmpPrior.fixedSNR > 0 ) /* special treatment of fixed-SNR: use h0=1, later rescale signal */
Amp.h0 = 1.0;
h0 = 1.0;
else if ( AmpPrior.fixedSNR == 0 )/* same as setting h0 = 0 */
Amp.h0 = 0;
h0 = 0;
else /* otherwise, draw from h0-prior */
Amp.h0 = XLALDrawFromPDF1D ( AmpPrior.pdf_h0Nat, rng );
h0 = XLALDrawFromPDF1D ( AmpPrior.pdf_h0Nat, rng );
Amp.cosi = XLALDrawFromPDF1D ( AmpPrior.pdf_cosi, rng );
cosi = XLALDrawFromPDF1D ( AmpPrior.pdf_cosi, rng );
Amp.aPlus = 0.5 * h0 * (1.0 + SQ(cosi));
Amp.aCross = h0 * cosi;
Amp.psi = XLALDrawFromPDF1D ( AmpPrior.pdf_psi, rng );
Amp.phi0 = XLALDrawFromPDF1D ( AmpPrior.pdf_phi0, rng );
......@@ -637,7 +641,8 @@ XLALSynthesizeTransientAtoms ( InjParams_t *injParamsOut, /**< [out] return su
XLAL_ERROR_NULL ( XLAL_EFUNC );
}
Amp.h0 *= rescale; /* rescale amplitude-params for consistency */
Amp.aPlus *= rescale; /* rescale amplitude-params for consistency */
Amp.aCross *= rescale;
UINT4 i; for (i=0; i < 4; i ++) A_Mu[i] *= rescale;
rho2 *= SQ(rescale); /* rescale reported optimal SNR */
......@@ -768,7 +773,7 @@ write_InjParams_to_fp ( FILE * fp, /**< [in] file-pointer to output file */
ret = fprintf ( fp, " %5.3f %6.3f %6.3f %7.3g %6.3f %6.3f %6.3f % 9.3g % 9.3g % 9.3g % 9.3g %10.5g %10.5g %10.5g %10.5g %9.3f %9.3f %1d %9.3g%s\n",
par->skypos.longitude, par->skypos.latitude, /* skypos */
par->SNR, /* SNR */
par->ampParams.h0, par->ampParams.cosi, par->ampParams.psi, par->ampParams.phi0, /* amplitude params {h0,cosi,psi,phi0}*/
par->ampParams.aPlus, par->ampParams.aCross, par->ampParams.psi, par->ampParams.phi0, /* amplitude params {aPlus,aCross,psi,phi0}*/
par->ampVect[0], par->ampVect[1], par->ampVect[2], par->ampVect[3], /* ampltiude vector A^mu */
par->multiAM.Mmunu.Ad, par->multiAM.Mmunu.Bd, par->multiAM.Mmunu.Cd, par->multiAM.Mmunu.Dd, /* antenna-pattern matrix components */
t0_d, tau_d, par->transientWindow.type, /* transient-window params */
......
......@@ -1514,7 +1514,12 @@ XLALComputeDopplerFstatMetric ( const DopplerMetricParams *metricParams, /**<
}
/* ----- compute the F-metric gF_ij and related matrices ---------- */
REAL8 cosi = metricParams->signalParams.Amp.cosi;
REAL8 h0 = metricParams->signalParams.Amp.aPlus + sqrt ( POW2(metricParams->signalParams.Amp.aPlus) - POW2(metricParams->signalParams.Amp.aCross) );
REAL8 cosi;
if (h0 > 0)
cosi = metricParams->signalParams.Amp.aCross/h0;
else
cosi = 0;
REAL8 psi = metricParams->signalParams.Amp.psi;
if ( (metric_k = XLALComputeFmetricFromAtoms ( atoms, cosi, psi)) == NULL ) {
......
......@@ -94,12 +94,14 @@ main ( int argc, char *argv[] )
// ----- CW sources to injet ----------
REAL8 Freq = 100.0;
REAL8 h0 = 1.0;
REAL8 cosi = 0.5;
PulsarParamsVector *injectSources;
XLAL_CHECK ( (injectSources = XLALCreatePulsarParamsVector(1)) != NULL, XLAL_EFUNC );
injectSources->data[0].Amp.h0 = 1;
injectSources->data[0].Amp.cosi = 0.5;
injectSources->data[0].Amp.aPlus = 0.5 * h0 * (1.0 + cosi * cosi);
injectSources->data[0].Amp.aCross = h0 * cosi;
injectSources->data[0].Amp.psi = 0.1;
injectSources->data[0].Amp.phi0 = 1.2;
......
......@@ -163,7 +163,8 @@ test_XLALComputeDopplerMetrics ( void )
// prepare metric parameters for modern XLALComputeDopplerFstatMetric() and mid-old XLALOldDopplerFstatMetric()
const DopplerCoordinateSystem coordSys = { 4, { DOPPLERCOORD_FREQ, DOPPLERCOORD_ALPHA, DOPPLERCOORD_DELTA, DOPPLERCOORD_F1DOT } };
const PulsarAmplitudeParams Amp = { 0.03, -0.3, 0.5, 0.0 }; // h0, cosi, psi, phi0
//const PulsarAmplitudeParams Amp = { 0.03, -0.3, 0.5, 0.0 }; // h0, cosi, psi, phi0
const PulsarAmplitudeParams Amp = { 0.5, 0.0, 0.01635, -0.009 }; // psi, phi0, aPlus, aCross
const PulsarDopplerParams dop = {
.refTime = startTimeGPS,
.Alpha = Alpha,
......
......@@ -851,11 +851,10 @@ InitCode ( ConfigVariables *cfg,
XLAL_CHECK ( cfg->multiAMcoe != NULL, XLAL_EFUNC );
/* ----- compute amplitude-factors alpha1, alpha2, alpha3 ----- */
REAL8 cosi = metricParams->signalParams.Amp.cosi;
REAL8 psi = metricParams->signalParams.Amp.psi;
REAL8 Aplus = 0.5 * ( 1.0 + SQ(cosi) );
REAL8 Across = cosi;
REAL8 Aplus = metricParams->signalParams.Amp.aPlus;
REAL8 <