1#include <lal/DopplerScan.h>
2#include <lal/LALInitBarycenter.h>
3#include <lal/UserInput.h>
4#include <lal/LALString.h>
5#include <gsl/gsl_sf_trig.h>
6#include <gsl/gsl_rng.h>
7#include "../TwoSpectSpecFunc.h"
29int main(
int argc,
char *argv[] )
38 if ( strcmp(
"H1", uvar.IFO->data[ii] ) == 0 ) {
40 }
else if ( strcmp(
"L1", uvar.IFO->data[ii] ) == 0 ) {
42 }
else if ( strcmp(
"V1", uvar.IFO->data[ii] ) == 0 ) {
44 }
else if ( strcmp(
"H2", uvar.IFO->data[ii] ) == 0 ) {
46 }
else if ( strcmp(
"H2r", uvar.IFO->data[ii] ) == 0 ) {
54 XLAL_ERROR(
XLAL_EINVAL,
"Not using valid interferometer! Expected 'H1', 'H2', 'H2r' (rotated H2), 'L1', or 'V1' not %s.\n", uvar.IFO->data[ii] );
75 gsl_rng_set( rng, 0 );
78 XLAL_CHECK( ( OUTPUT = fopen( uvar.outfilename,
"w" ) ) != NULL,
XLAL_EIO,
"Output file %s could not be opened\n", uvar.outfilename );
80 for (
INT4 n = 0;
n < uvar.skylocations;
n++ ) {
96 cosi0 = 2.0 * gsl_rng_uniform( rng ) - 1.0;
101 psi0 =
LAL_PI * gsl_rng_uniform( rng );
111 REAL8 frequency0 =
frequency + ( gsl_rng_uniform( rng ) - 0.5 ) / uvar.Tsft;
114 REAL4 Fplus0 = multiAMcoefficients->
data[0]->
a->
data[ii] * cos( 2.0 * psi0 ) + multiAMcoefficients->
data[0]->
b->
data[ii] * sin( 2.0 * psi0 );
115 REAL4 Fcross0 = multiAMcoefficients->
data[0]->
b->
data[ii] * cos( 2.0 * psi0 ) - multiAMcoefficients->
data[0]->
a->
data[ii] * sin( 2.0 * psi0 );
116 REAL4 Fplus1 = multiAMcoefficients->
data[1]->
a->
data[ii] * cos( 2.0 * psi0 ) + multiAMcoefficients->
data[1]->
b->
data[ii] * sin( 2.0 * psi0 );
117 REAL4 Fcross1 = multiAMcoefficients->
data[1]->
b->
data[ii] * cos( 2.0 * psi0 ) - multiAMcoefficients->
data[1]->
a->
data[ii] * sin( 2.0 * psi0 );
118 COMPLEX16 RatioTerm0 =
crect( 0.5 * Fplus1 * ( 1.0 + cosi0 * cosi0 ), Fcross1 * cosi0 ) /
crect( 0.5 * Fplus0 * ( 1.0 + cosi0 * cosi0 ), Fcross0 * cosi0 );
120 REAL4 detPhaseArg = 0.0, detPhaseMag = 0.0;
122 for (
INT4 jj = 0; jj < 16 && !loopbroken; jj++ ) {
124 Fplus0 = multiAMcoefficients->
data[0]->
a->
data[ii] * cos( 2.0 * psi ) + multiAMcoefficients->
data[0]->
b->
data[ii] * sin( 2.0 * psi );
125 Fcross0 = multiAMcoefficients->
data[0]->
b->
data[ii] * cos( 2.0 * psi ) - multiAMcoefficients->
data[0]->
a->
data[ii] * sin( 2.0 * psi );
126 Fplus1 = multiAMcoefficients->
data[1]->
a->
data[ii] * cos( 2.0 * psi ) + multiAMcoefficients->
data[1]->
b->
data[ii] * sin( 2.0 * psi );
127 Fcross1 = multiAMcoefficients->
data[1]->
b->
data[ii] * cos( 2.0 * psi ) - multiAMcoefficients->
data[1]->
a->
data[ii] * sin( 2.0 * psi );
128 for (
INT4 kk = 0; kk < 21 && !loopbroken; kk++ ) {
130 if ( !uvar.unrestrictedCosi ) {
139 if ( cabs( complexdenominator ) > 1.0e-6 ) {
140 COMPLEX16 complexval = complexnumerator / complexdenominator;
141 detPhaseMag +=
fmin( cabs( complexval ), 10.0 );
142 detPhaseArg += gsl_sf_angle_restrict_pos( carg( complexval ) );
150 detPhaseMag /= 336.0;
151 detPhaseArg /= 336.0;
159 REAL8 tau = timediff1 - timediff0;
167 REAL8 realSigBinDiff = round( delta0_0 ) - round( delta1_0 );
168 delta1_0 += realSigBinDiff;
171 REAL8 estSigBinDiff = round( delta0 ) - round( delta1 );
172 delta1 += estSigBinDiff;
174 if ( !uvar.rectWindow ) {
175 if ( fabsf( (
REAL4 )delta1_0 ) < (
REAL4 )1.0e-6 ) {
176 if ( fabsf( (
REAL4 )delta0_0 ) < (
REAL4 )1.0e-6 ) {
177 dirichlet0 =
crect( 1.0, 0.0 );
178 }
else if ( fabsf( (
REAL4 )( delta0_0 * delta0_0 - 1.0 ) ) < (
REAL4 )1.0e-6 ) {
179 dirichlet0 =
crect( -2.0, 0.0 );
180 }
else if ( fabsf( (
REAL4 )( delta0_0 - roundf( delta0_0 ) ) ) < (
REAL4 )1.0e-6 ) {
181 dirichlet0 =
crect( 0.0, 0.0 );
184 dirichlet0 = -
LAL_PI *
crectf( cos(
LAL_PI * delta0_0 ), -sin(
LAL_PI * delta0_0 ) ) * delta0_0 * ( delta0_0 * delta0_0 - 1.0 ) / sin(
LAL_PI * delta0_0 );
186 }
else if ( fabsf( (
REAL4 )( delta1_0 * delta1_0 - 1.0 ) ) < (
REAL4 )1.0e-6 ) {
187 if ( fabsf( (
REAL4 )delta0_0 ) < (
REAL4 )1.0e-6 ) {
188 dirichlet0 =
crect( -0.5, 0.0 );
189 }
else if ( fabsf( (
REAL4 )( delta0_0 * delta0_0 - 1.0 ) ) < (
REAL4 )1.0e-6 ) {
190 dirichlet0 =
crect( 1.0, 0.0 );
191 }
else if ( fabsf( (
REAL4 )( delta0_0 - roundf( delta0_0 ) ) ) < (
REAL4 )1.0e-6 ) {
192 dirichlet0 =
crect( 0.0, 0.0 );
195 dirichlet0 = -
LAL_PI_2 *
crectf( -cos(
LAL_PI * delta0_0 ), sin(
LAL_PI * delta0_0 ) ) * delta0_0 * ( delta0_0 * delta0_0 - 1.0 ) / sin(
LAL_PI * delta0_0 );
197 }
else if ( fabsf( (
REAL4 )delta0_0 ) < (
REAL4 )1.0e-6 ) {
198 dirichlet0 = -
LAL_1_PI *
crectf( cos(
LAL_PI * delta1_0 ), sin(
LAL_PI * delta1_0 ) ) * sin(
LAL_PI * delta1_0 ) / ( delta1_0 * ( delta1_0 * delta1_0 - 1.0 ) );
199 }
else if ( fabsf( (
REAL4 )( delta0_0 * delta0_0 - 1.0 ) ) < (
REAL4 )1.0e-6 ) {
200 dirichlet0 =
LAL_2_PI *
crectf( cos(
LAL_PI * delta1_0 ), sin(
LAL_PI * delta1_0 ) ) * sin(
LAL_PI * delta1_0 ) / ( delta1_0 * ( delta1_0 * delta1_0 - 1.0 ) );
201 }
else if ( fabsf( (
REAL4 )( delta0_0 - roundf( delta0_0 ) ) ) < (
REAL4 )1.0e-6 ) {
202 dirichlet0 =
crect( 0.0, 0.0 );
205 dirichlet0 = sin(
LAL_PI * delta1_0 ) / sin(
LAL_PI * delta0_0 ) * ( delta0_0 * ( delta0_0 * delta0_0 - 1.0 ) ) / ( delta1_0 * ( delta1_0 * delta1_0 - 1.0 ) ) *
cpolarf( 1.0,
LAL_PI * ( delta1_0 - delta0_0 ) );
208 if ( fabsf( (
REAL4 )delta1_0 ) < (
REAL4 )1.0e-6 ) {
209 if ( fabsf( (
REAL4 )delta0_0 ) < (
REAL4 )1.0e-6 ) {
210 dirichlet0 =
crect( 1.0, 0.0 );
211 }
else if ( fabsf( (
REAL4 )( delta0_0 - roundf( delta0_0 ) ) ) < (
REAL4 )1.0e-6 ) {
212 dirichlet0 =
crect( 0.0, 0.0 );
217 }
else if ( fabsf( (
REAL4 )delta0_0 ) < (
REAL4 )1.0e-6 ) {
219 }
else if ( fabsf( (
REAL4 )( delta0_0 - roundf( delta0_0 ) ) ) < (
REAL4 )1.0e-6 ) {
220 dirichlet0 =
crect( 0.0, 0.0 );
228 if ( !uvar.rectWindow ) {
229 if ( fabsf( (
REAL4 )delta1 ) < (
REAL4 )1.0e-6 ) {
230 if ( fabsf( (
REAL4 )delta0 ) < (
REAL4 )1.0e-6 ) {
231 dirichlet =
crect( 1.0, 0.0 );
232 }
else if ( fabsf( (
REAL4 )( delta0 * delta0 - 1.0 ) ) < (
REAL4 )1.0e-6 ) {
233 dirichlet =
crect( -2.0, 0.0 );
234 }
else if ( fabsf( (
REAL4 )( delta0 - roundf( delta0 ) ) ) < (
REAL4 )1.0e-6 ) {
235 dirichlet =
crect( 0.0, 0.0 );
239 }
else if ( fabsf( (
REAL4 )( delta1 * delta1 - 1.0 ) ) < (
REAL4 )1.0e-6 ) {
240 if ( fabsf( (
REAL4 )delta0 ) < (
REAL4 )1.0e-6 ) {
241 dirichlet =
crect( -0.5, 0.0 );
242 }
else if ( fabsf( (
REAL4 )( delta0 * delta0 - 1.0 ) ) < (
REAL4 )1.0e-6 ) {
243 dirichlet =
crect( 1.0, 0.0 );
244 }
else if ( fabsf( (
REAL4 )( delta0 - roundf( delta0 ) ) ) < (
REAL4 )1.0e-6 ) {
245 dirichlet =
crect( 0.0, 0.0 );
249 }
else if ( fabsf( (
REAL4 )delta0 ) < (
REAL4 )1.0e-6 ) {
251 }
else if ( fabsf( (
REAL4 )( delta0 * delta0 - 1.0 ) ) < (
REAL4 )1.0e-6 ) {
253 }
else if ( fabsf( (
REAL4 )( delta0 - roundf( delta0 ) ) ) < (
REAL4 )1.0e-6 ) {
254 dirichlet =
crect( 0.0, 0.0 );
256 dirichlet = sin(
LAL_PI * delta1 ) / sin(
LAL_PI * delta0 ) * ( delta0 * ( delta0 * delta0 - 1.0 ) ) / ( delta1 * ( delta1 * delta1 - 1.0 ) ) *
cpolarf( 1.0,
LAL_PI * ( delta1 - delta0 ) );
259 if ( fabsf( (
REAL4 )delta1 ) < (
REAL4 )1.0e-6 ) {
260 if ( fabsf( (
REAL4 )delta0 ) < (
REAL4 )1.0e-6 ) {
261 dirichlet =
crect( 1.0, 0.0 );
262 }
else if ( fabsf( (
REAL4 )( delta0 - roundf( delta0 ) ) ) < (
REAL4 )1.0e-6 ) {
263 dirichlet =
crect( 0.0, 0.0 );
267 }
else if ( fabsf( (
REAL4 )delta0 ) < (
REAL4 )1.0e-6 ) {
269 }
else if ( fabsf( (
REAL4 )( delta0 - roundf( delta0 ) ) ) < (
REAL4 )1.0e-6 ) {
270 dirichlet =
crect( 0.0, 0.0 );
275 dirichlet =
cpolar( 1.0, carg( dirichlet ) );
276 if ( fabs( floor( delta0 ) - floor( delta1 ) ) >= 1.0 ) {
280 COMPLEX16 realRatio = RatioTerm0 * phaseshift0 * conj( dirichlet0 );
281 COMPLEX16 estRatio = RatioTerm * phaseshift * conj( dirichlet );
283 fprintf( OUTPUT,
"%g %g %g %g\n", cabs( realRatio ), gsl_sf_angle_restrict_pos( carg( realRatio ) ), cabs( estRatio ), gsl_sf_angle_restrict_pos( carg( estRatio ) ) );
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
static double double delta
int main(int argc, char *argv[])
INT4 InitUserVars(UserVariables_t *uvar, int argc, char *argv[])
LALDetector * XLALCreateDetector(LALDetector *detector, const LALFrDetector *frDetector, LALDetectorType type)
const LALDetector lalCachedDetectors[LAL_NUM_DETECTORS]
void XLALDestroyMultiDetectorStateSeries(MultiDetectorStateSeries *mdetStates)
Helper function to get rid of a multi-IFO DetectorStateSeries Note, this is "NULL-robust" in the sens...
MultiDetectorStateSeries * XLALGetMultiDetectorStates(const MultiLIGOTimeGPSVector *multiTS, const MultiLALDetector *multiIFO, const EphemerisData *edat, REAL8 tOffset)
Get the detector-time series for the given MultiLIGOTimeGPSVector.
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
void XLALDestroyMultiAMCoeffs(MultiAMCoeffs *multiAMcoef)
Destroy a MultiAMCoeffs structure.
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
#define XLAL_INIT_DECL(var,...)
void * XLALMalloc(size_t n)
char char * XLALStringDuplicate(const char *s)
MultiLIGOTimeGPSVector * XLALMakeMultiTimestamps(LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap, UINT4 numDet)
Same as XLALMakeTimestamps() just for several detectors, additionally specify the number of detectors...
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
MultiSSBtimes * XLALGetMultiSSBtimes(const MultiDetectorStateSeries *multiDetStates, SkyPosition skypos, LIGOTimeGPS refTime, SSBprecision precision)
Multi-IFO version of XLALGetSSBtimes().
void XLALDestroyMultiSSBtimes(MultiSSBtimes *multiSSB)
Destroy a MultiSSBtimes structure.
@ SSBPREC_RELATIVISTICOPT
optimized relativistic, numerically equivalent to SSBPREC_RELATIVISTIC, but faster
COORDINATESYSTEM_EQUATORIAL
#define XLAL_CHECK(assertion,...)
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
LALDetector detectors[LAL_NUM_DETECTORS]
REAL4Vector * b
(weighted) per-SFT antenna-pattern function
REAL4Vector * a
(weighted) per-SFT antenna-pattern function
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
AMCoeffs ** data
noise-weighted AM-coeffs , and
Multi-IFO time-series of DetectorStates.
array of detectors definitions 'LALDetector'
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
Multi-IFO container for SSB timings.
SSBtimes ** data
array of SSBtimes (pointers)
REAL8Vector * Tdot
dT/dt : time-derivative of SSB-time wrt local time for SFT-alpha
REAL8Vector * DeltaT
Time-difference of SFT-alpha - tau0 in SSB-frame.
CHAR * ephemSun
Sun ephemeris file to use.
CHAR * ephemEarth
Earth ephemeris file to use.
REAL8 SFToverlap
overlap SFTs by this many seconds
REAL8 Tsft
SFT time baseline Tsft.