LALDemodTest.c 38.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
*  Copyright (C) 2007 Jolien Creighton, Maria Alessandra Papa, Reinhard Prix, Steve Berukoff, Teviet Creighton
*
*  This program is free software; you can redistribute it and/or modify
*  it under the terms of the GNU General Public License as published by
*  the Free Software Foundation; either version 2 of the License, or
*  (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with with program; see the file COPYING. If not, write to the
*  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
*  MA  02111-1307  USA
*/

Reinhard Prix's avatar
Reinhard Prix committed
20
/**
21 22 23 24 25 26
 * \author Berukoff, S.J., Papa, M.A.,
 * \file
 * \ingroup pulsarTODO
 *
 * \brief Performs required tests of LALDemod().
 *
27 28
 * ### Usage ###
 *
29 30 31 32
 * \code
 * LALDemodTest -i <input data file> [-d <gap>] [-n] [-o]
 * \endcode
 *
33
 * ### Description ###
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
 *
 * This routine performs tests on the routine <tt>LALDemod()</tt>.
 * Options:
 * <ul>
 * <li> <tt>-i</tt> -- the input data file (default is 'in.data'; an example is included, format below)</li>
 * <li> <tt>-n</tt> -- add zero-mean Gaussian noise to the signal</li>
 * <li> <tt>-d \<gap\></tt> -- simulate gaps in the data.  The number <tt>\<gaps\></tt> refers to the integral number of SFT timescales between adjacent timestamps.</li>
 * <li> <tt>-o</tt> -- print out result data files</li>
 * </ul>
 *
 * Structure:
 * In more detail, let us begin with a discussion of the structure of the test
 * code, which is composed of several modules.
 * <ul>
 * <li> The first module reads in data from an input parameter data file.
 * The parameters must be listed in the input file in the following order,
 * with the corresponding format:
 * \code
 * total observation time -- float
 * coherent search time -- float
 * factor by which to modify SFT timescale -- float
 * Cross amplitude -- float
 * Plus amplitude -- float
 * DeFT frequency band (centered by default around f0) -- float
 * f0, intrinsic frequency of the signal at the beginning of the
 * observation -- float
 * maximum order of signal spindown parameters -- int
 * signal spindown parameter 1 --  scientific notation
 * signal spindown parameter 2 --  scientific notation
 * signal spindown parameter 3 --  scientific notation
 * signal spindown parameter 4 --  scientific notation
 * signal spindown parameter 5 --  scientific notation
 * signal source right ascension (alpha) -- float (value in DEGREES)
 * signal source declination (delta) -- float (value in DEGREES)
 * maximum order of template spindown parameters -- int
 * template spindown parameter 1 --  scientific notation (NOT
 * template spindown parameter 2 --  scientific notation (NOT scaled by f0)
 * template spindown parameter 3 --  scientific notation (NOT scaled by f0)
 * template spindown parameter 4 --  scientific notation (NOT scaled by f0)
 * template spindown parameter 5 --  scientific notation (NOT scaled by f0)
 * template source right ascension (alpha) -- float (value in DEGREES)
 * template source declination (delta) -- float (value in DEGREES)
 * \endcode
 * Note: Above, the *signal* spindown parameters are scaled by the intrinsic frequency, while the
 * template* spindown parameters are not.  This is due to the difference in definitions between the
 * PulsarSimulateCoherentGW() package, which generates the signal, and this package.
 *
 * </li><li> The next module in the test code, which is optionally executed with
 * the '<tt>-n</tt>' switch, creates noise using LALs
 * <tt>LALNormalDeviates()</tt> routine.  By design, the noise is created in
 * single
 * precision, and is zero-mean and Gaussian.  This noise is added,
 * datum-by-datum, to the time series created in
 * the next module, after the amplitude of the time series has been
 * changed by a
 * factor of \c SNR, which is specified in the input data file.
 *
 * </li><li> The next module to be invoked creates a time series, according to
 * the standard model for pulsars with spindown.  This is done by using the
 * <tt>LALGenerateTaylorCW()</tt> and <tt>LALPulsarSimulateCoherentGW()</tt> functions.  This time series
 * undergoes an FFT, and this transformed data then constitutes the SFT data to be
 * input to the demodulation code.  The fake signal data is characterized
 * by an intrinsic frequency at the beginning of the observation plus some
 * other source parameters.  The DeFT is produced in a band \c f0Band
 * (as specified as an input parameter) and centered at this frequency.
 * The width of the band (plus some extra width of \f$2\cdot 10^{-4}f0 \f$ Hz) determines the
 * sampling frequency of the time series (Nyquist theorem).  In practice this
 * would be the inverse FFT of a data set that has been band-passed around
 * \c f0 and then appropriately down-sampled (e.g. with a lock-in).  The
 * normalization rule for FFT data is the following: if sinusoidal data over a
 * time \f$T\f$ and with amplitude \f$A\f$ is FFT-ed, the sum of the square amplitude of
 * the output of the FFT (power) is equal to \f${A^2 T}\f$.  Thus, the power peak at
 * the sinusoids frequency should be expected to be \f$\sim\f$ \f$\frac{A^2}{2} T\f$,
 * within a factor of 2.  The same normalization rule applies to the DeFT data.
 * Thus by piecing together \f$N\f$ SFTs we expect a DeFT power peak \f$\sim N\f$ higher
 * than that of the SFTs - at least in the case of perfect signal-template match.
 *
 * Let us now spend a few words on the choice of the SFT time baseline.  Given an
 * intrinsic search frequency one can compute the longest time baseline which is
 * still compatible with the requirement that the instantaneous signal frequency
 * during such time baseline does not shift by more than a frequency bin.  This
 * is the default choice for the SFT length, having assumed that the modulation
 * is due to the spin of the Earth and having taken a simple epicyclic model to
 * evaluate the magnitude of this effect.  It is possible to choose a different
 * time baseline by specifying a value for
 * the variable \c gap other than 1.  Note that the SFT time baseline is
 * approximated to the nearest value such that the number of SFT samples is a
 * power of two.  This is also well documented in the code.
 *
 * The set of SFTs does not necessarily come from contiguous data sets: a set of
 * time stamps is created that defines the time of the first sample of each SFT
 * data chunk.  The timestamps which are required in many parts of the code are
 * generated in a small subroutine <tt>times2()</tt>.  This routine takes as input
 * the SFT timescale \c tSFT, the number of SFTs which will be created,
 * \c mObsSFT, and a switch which lets the code know whether to make even
 * timestamps, or timestamps with gaps (see below for more on this).  The
 * subroutine then writes the times to the \c LIGOTimeGPS vector containing
 * the timestamps for the entire test code, and returns this vector.  Note that
 * each datum of the  \c LIGOTimeGPS vector is comprised of two fields; if
 * accessing the \f$i^{th}\f$ datum, the seconds part of the timestamp vector
 * \c ts is <tt>ts[i].gpsSeconds</tt> and the nanoseconds part is
 * <tt>ts[i].gpsNanoSeconds</tt>.  These are the fields which are written in this
 * <tt>times()</tt>.
 *
 * As an important side note, let us discuss the effect that a vector of
 * timestamps with gaps has on the resulting transformed data.  Since each of the
 * timestamps refers to the first datum of each SFT, the existence of the gaps
 * means that instead of transforming a continuous set of data, we are reduced to
 * transforming a piecewise continuous set.  Since we can envision these gaps as
 * simply replacing real data with zeros, we correspondingly should see a power
 * loss in the resulting FFTs signal bins and a broadening of the power spectrum.
 * Since real detectors will clearly have gaps in the data, this effect is
 * obviously something we seek to minimize or eliminate if possible.  This work
 * continues to be under development.
 *
 * The total observation time determines how many SFTs and how many DeFTs are
 * created.  The actual time baseline for both the DeFTs and the total
 * observation time might differ from the ones defined in the input file, the
 * reason being that they are rounded to the nearest multiple of the SFT time
 * baseline.
 *
 * Note that use is made of the <tt>LALBarycenter()</tt> routine (see
 * \ref LALBarycenter.h), which (among other things) provides,  at any given
 * time, the actual instantaneous position and velocity of a  detector at any
 * specified location of the Earth with respect to the SSB.
 *
 * </li><li> Following the creation of a short chunk of time series data, an FFT is
 * performed with the internal FFTW routines.  This outputs a frequency domain
 * chunk which is placed into the \c SFTData array of structures.  This will
 * contain all of the SFT data we need to demodulate, and in the future, will be
 * the storage area for the real data.
 *
 * </li><li> The next module begins the demodulation process.  First, the parameters
 * for the demodulation routine are assigned from values previously calculated in
 * the test code.  Similarly, parameters for the <tt>LALComputeSky()</tt> routine are
 * assigned.  This routine computes the coefficients \f$A_{s\alpha}\f$ and
 * \f$B_{s\alpha}\f$ (see \ref ComputeSky.h) of the spindown parameters
 * for the phase model weve assumed.  These coefficients are used within the
 * <tt>LALDemod()</tt> routine itself. Since they only depend on the template sky
 * position, in a search over many different spin-down parameters they are
 * reused, thus one needs compute them only once.  Then, the <tt>LALComputeAM()</tt>
 * routine is called, to calculate the amplitude modulation filter information.  Finally, at last, the
 * demodulation routine itself is called, and, if the command line option
 * '<tt>-o</tt>' is used,  output are several data files containing demodulated
 * data (these are by default named '<tt>xhat_#</tt>').  These output files have two columns, one
 * for the value of the periodogram and one for the frequency.
 *
 * </li></ul>
 *
183 184 185 186
 * ### Exit codes ###
 *
 *
 * ### Uses ###
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
 *
 * \code
 * lalDebugLevel
 * LALMalloc()
 * LALSCreateVector()
 * LALCreateRandomParams()
 * LALNormalDeviates()
 * LALDestroyRandomParams()
 * LALSDestroyVector()
 * LALCCreateVector()
 * LALCreateForwardRealFFTPlan()
 * LALREAL4VectorFFT()
 * LALCDestroyVector()
 * LALDestroyRealFFTPlan()
 * LALGenerateTaylorCW()
 * LALPulsarSimulateCoherentGW()
 * LALComputeSky()
 * LALFree()
 * LALDemod()
 * LALBarycenter()
 * LALComputeAM()
 * \endcode
 *
210 211
 * ### Notes ###
 *
212 213 214 215 216 217
 * The implementation of the code here is intended to give a general outline of
 * what the demodulation code needs to work.  Most of this test function performs
 * steps (e.g., noise, time- and frequency-series generation) that will be already
 * present in the data.
 *
 */
Adam Mercer's avatar
Adam Mercer committed
218

jolien's avatar
jolien committed
219 220 221 222 223
#ifndef LALDEMODTEST_C
#define LALDEMODTEST_C
#endif

/* Usage */
224
#define USAGE "Usage: %s [-i basicInputsFile] [-n] [-d] [-o]\n"
jolien's avatar
jolien committed
225 226 227 228 229 230


/* Error macro, taken from ResampleTest.c in LAL's pulsar/test package */
#define ERROR( code, msg, statement )                                \
if ( lalDebugLevel & LALERROR )                                      \
{                                                                    \
231
  XLALPrintError( "Error[0] %d: program %s, file %s, line %d, %s\n"   \
jolien's avatar
jolien committed
232
		 "        %s %s\n", (code), *argv, __FILE__,		\
233
		 __LINE__, "$Id$", statement ? statement : "",\
jolien's avatar
jolien committed
234 235 236 237
		 (msg) );                                            \
}                                                                    \
else (void)(0)

238 239 240
#include <lal/LALDemod.h>
#include <lal/LALInitBarycenter.h>
#include <lal/SeqFactories.h>
sberukoff's avatar
sberukoff committed
241
#include <lal/FileIO.h>
242 243 244
#include <unistd.h>
#include <sys/types.h>

sberukoff's avatar
sberukoff committed
245 246 247 248

static void TimeToFloat(REAL8 *f, LIGOTimeGPS *tgps);

static void FloatToTime(LIGOTimeGPS *tgps, REAL8 *f);
jolien's avatar
jolien committed
249

250 251 252
static void times(REAL8 , INT4, LIGOTimeGPS *, INT4 );

static void times2(REAL8 tSFT, INT4 howMany, LIGOTimeGPS **ts, INT4 **sftPerCoh, INT4 sw, INT4 mCohSFT);
253

jolien's avatar
jolien committed
254

255
static PulsarCoherentGW emptySignal;
256

jolien's avatar
jolien committed
257 258
int main(int argc, char **argv)
{
sberukoff's avatar
sberukoff committed
259
  static LALStatus status;
Adam Mercer's avatar
Adam Mercer committed
260

sberukoff's avatar
sberukoff committed
261 262 263 264
  /***** VARIABLE DECLARATION *****/

  ParameterSet *signalParams;
  ParameterSet *templateParams;
265
  char *basicInputsFile;
sberukoff's avatar
sberukoff committed
266 267 268 269 270
  FILE *bif;
  REAL8 tObs, tCoh, tSFT;
  REAL8 oneOverSqrtTSFT;
  REAL8 aCross, aPlus, SNR;
  REAL8 f0;
Adam Mercer's avatar
Adam Mercer committed
271

sberukoff's avatar
sberukoff committed
272 273 274 275
  INT4 mCohSFT, mObsCoh, mObsSFT;
  REAL8 dfSFT, dt;
  INT4 if0Min, if0Max, ifMin, ifMax;
  REAL8 f0Min, f0Max, fMin, f0Band, fWing;
papa's avatar
papa committed
276
  INT4 nDeltaF,ntermsdivbytwo;
277
  LALFstat Fstat;
sberukoff's avatar
sberukoff committed
278 279 280 281 282 283 284

  LIGOTimeGPS *timeStamps;

  REAL4Vector *temp = NULL;
  REAL4Vector *noise = NULL;
  static RandomParams *params;
  INT4 seed=0;
Adam Mercer's avatar
Adam Mercer committed
285

sberukoff's avatar
sberukoff committed
286
  INT4 a,i,k;
Adam Mercer's avatar
Adam Mercer committed
287

sberukoff's avatar
sberukoff committed
288 289 290
  FFT **SFTData;
  RealFFTPlan *pfwd = NULL;
  COMPLEX8Vector *fvec = NULL;
Adam Mercer's avatar
Adam Mercer committed
291

sberukoff's avatar
sberukoff committed
292 293 294
  DemodPar *demParams;
  CSParams *csParams;
  INT4 iSkyCoh;
Adam Mercer's avatar
Adam Mercer committed
295

sberukoff's avatar
sberukoff committed
296
  DeFTPeriodogram **xHat;
Adam Mercer's avatar
Adam Mercer committed
297

sberukoff's avatar
sberukoff committed
298
  const CHAR *noi=NULL;
299
  INT4 deletions=0;
sberukoff's avatar
sberukoff committed
300
  const CHAR *output=NULL;
Adam Mercer's avatar
Adam Mercer committed
301

sberukoff's avatar
sberukoff committed
302 303 304 305 306 307 308 309 310 311 312
  REAL8 factor;
  INT2 arg;

  /* file name if needed for SFT output files */
  CHAR *sftoutname=NULL;

  /* Quantities for use in LALBarycenter package */
  BarycenterInput baryinput;
  EmissionTime emit;
  EarthState earth;
  LALDetector cachedDetector;
Adam Mercer's avatar
Adam Mercer committed
313

sberukoff's avatar
sberukoff committed
314
  EphemerisData *edat=NULL;
315 316
  char earthEphemeris[]="earth98.dat";
  char sunEphemeris[]="sun98.dat";
sberukoff's avatar
sberukoff committed
317 318 319 320 321 322

  REAL4TimeSeries *timeSeries = NULL;

  /*
   *  Variables for AM correction
   */
323
  PulsarCoherentGW cgwOutput = emptySignal;
sberukoff's avatar
sberukoff committed
324
  TaylorCWParamStruc genTayParams;
325
  PulsarDetectorResponse cwDetector;
sberukoff's avatar
sberukoff committed
326 327
  AMCoeffsParams *amParams;
  AMCoeffs amc;
328 329
  INT4 *sftPerCoh;
  INT4 totalSFT=0;
Adam Mercer's avatar
Adam Mercer committed
330

sberukoff's avatar
sberukoff committed
331 332
#define DEBUG 1
#if (0)
Adam Mercer's avatar
Adam Mercer committed
333 334 335 336 337

  INT2 tmp=0;
  REAL8 dtSFT, dfCoh;
  REAL8 fMax;
  INT4 nSFT;
jolien's avatar
jolien committed
338

sberukoff's avatar
sberukoff committed
339
#endif
Adam Mercer's avatar
Adam Mercer committed
340

sberukoff's avatar
sberukoff committed
341
  /***** END VARIABLE DECLARATION *****/
Adam Mercer's avatar
Adam Mercer committed
342 343 344


  /***** PARSE COMMAND LINE OPTIONS *****/
sberukoff's avatar
sberukoff committed
345
  basicInputsFile=(char *)LALMalloc(50*sizeof(char));
346
  snprintf(basicInputsFile, 9, "in.data\0");
sberukoff's avatar
sberukoff committed
347 348
  arg=1;
  while(arg<argc) {
Adam Mercer's avatar
Adam Mercer committed
349

sberukoff's avatar
sberukoff committed
350 351 352
    /* the input file */
    if(!strcmp(argv[arg],"-i"))
      {
353
	strcpy(basicInputsFile, argv[++arg]);
sberukoff's avatar
sberukoff committed
354
	arg++;
355
	if(fopen(basicInputsFile, "r")==NULL)
sberukoff's avatar
sberukoff committed
356 357
	  {
	    ERROR(LALDEMODH_ENOFILE, LALDEMODH_MSGENOFILE, 0);
358
	    XLALPrintError(USAGE, *argv);
sberukoff's avatar
sberukoff committed
359 360 361
	    return LALDEMODH_ENOFILE;
	  }
      }
Adam Mercer's avatar
Adam Mercer committed
362

sberukoff's avatar
sberukoff committed
363 364 365 366 367 368 369 370
    /* turn noise on? if string is not NULL, it will be */
    else if(!strcmp(argv[arg],"-n"))
      {
	noi="n";
	arg++;
      }

    /* make SFT file output (such that the driver code can read in as input data) */
Adam Mercer's avatar
Adam Mercer committed
371
    else if (!strcmp(argv[arg],"-m"))
sberukoff's avatar
sberukoff committed
372 373 374
      {
      sftoutname=argv[++arg];
      arg++;
Adam Mercer's avatar
Adam Mercer committed
375
      }
sberukoff's avatar
sberukoff committed
376 377 378
    /* turn timestamp deletions on? if string is not NULL, it will be */
    else if(!strcmp(argv[arg],"-d"))
      {
379
	deletions=1;
sberukoff's avatar
sberukoff committed
380 381
	arg++;
      }
Adam Mercer's avatar
Adam Mercer committed
382

sberukoff's avatar
sberukoff committed
383 384 385 386 387 388
    /* turn on output? if string is not NULL, it will be */
    else if(!strcmp(argv[arg],"-o"))
      {
	output="o";
	arg++;
      }
Adam Mercer's avatar
Adam Mercer committed
389

sberukoff's avatar
sberukoff committed
390 391 392
    /* default: no input file specified */
    else if(basicInputsFile==NULL)
      {
393
	bif=fopen("in.data", "r");
sberukoff's avatar
sberukoff committed
394
      }
Adam Mercer's avatar
Adam Mercer committed
395

sberukoff's avatar
sberukoff committed
396 397 398
    /* erroneous command line argument */
    else if(basicInputsFile!=NULL && arg<argc)
      {
Adam Mercer's avatar
Adam Mercer committed
399
	ERROR(LALDEMODH_EBADARG, LALDEMODH_MSGEBADARG, 0);
400
	XLALPrintError(USAGE, *argv);
sberukoff's avatar
sberukoff committed
401 402 403 404
	arg=argc;
	return LALDEMODH_EBADARG;
      }
  }
Adam Mercer's avatar
Adam Mercer committed
405

sberukoff's avatar
sberukoff committed
406
  /***** END COMMAND LINE PARSING *****/
jolien's avatar
jolien committed
407

Adam Mercer's avatar
Adam Mercer committed
408

sberukoff's avatar
sberukoff committed
409
  /***** INITIALIZATION OF SIGNAL AND TEMPLATE *****/
jolien's avatar
jolien committed
410

sberukoff's avatar
sberukoff committed
411 412 413 414 415
  /* Allocate space for signal parameters */
  signalParams=LALMalloc(sizeof(ParameterSet));
  signalParams->spind=LALMalloc(sizeof(Spindown));
  signalParams->skyP=LALMalloc(2*sizeof(SkyPos));
  signalParams->spind->spParams=LALMalloc(5*sizeof(REAL8));
jolien's avatar
jolien committed
416

sberukoff's avatar
sberukoff committed
417 418 419 420 421
  /* Allocate space for template parameters */
  templateParams=LALMalloc(sizeof(ParameterSet));
  templateParams->spind=LALMalloc(sizeof(Spindown));
  templateParams->skyP=LALMalloc(2*sizeof(SkyPos));
  templateParams->spind->spParams=LALMalloc(5*sizeof(REAL8));
jolien's avatar
jolien committed
422

sberukoff's avatar
sberukoff committed
423
  /***** END INITIALIZATION *****/
Adam Mercer's avatar
Adam Mercer committed
424 425


sberukoff's avatar
sberukoff committed
426
  /***** GET INPUTS FROM FILES *****/
jolien's avatar
jolien committed
427

428
  bif=fopen(basicInputsFile, "r");
Adam Mercer's avatar
Adam Mercer committed
429

sberukoff's avatar
sberukoff committed
430 431 432
  fscanf(bif, "%lf\n%lf\n%lf\n%lf\n%lf\n%lf\n%lf\n%lf\n%d\n%le\n%le\n%le\n%le\n%le\n%lf\n"
	 " %lf\n%d\n%le\n%le\n%le\n%le\n%le\n%lf\n%lf\n",
	 &tObs, &tCoh, &factor, &aCross, &aPlus, &SNR, &f0Band, &f0,
Adam Mercer's avatar
Adam Mercer committed
433 434
	 &signalParams->spind->m,
	 &signalParams->spind->spParams[0], &signalParams->spind->spParams[1],
sberukoff's avatar
sberukoff committed
435
	 &signalParams->spind->spParams[2], &signalParams->spind->spParams[3],
Adam Mercer's avatar
Adam Mercer committed
436 437
	 &signalParams->spind->spParams[4],
	 &signalParams->skyP->alpha, &signalParams->skyP->delta,
sberukoff's avatar
sberukoff committed
438 439 440
	 &templateParams->spind->m,
	 &templateParams->spind->spParams[0], &templateParams->spind->spParams[1],
	 &templateParams->spind->spParams[2], &templateParams->spind->spParams[3],
Adam Mercer's avatar
Adam Mercer committed
441
	 &templateParams->spind->spParams[4],
sberukoff's avatar
sberukoff committed
442
	 &templateParams->skyP->alpha, &templateParams->skyP->delta);
Adam Mercer's avatar
Adam Mercer committed
443

444
  fclose(bif);
jolien's avatar
jolien committed
445

sberukoff's avatar
sberukoff committed
446
  /***** END FILE INPUT *****/
Adam Mercer's avatar
Adam Mercer committed
447

jolien's avatar
jolien committed
448

sberukoff's avatar
sberukoff committed
449 450
  /***** CALCULATE USEFUL QUANTITIES *****/

Adam Mercer's avatar
Adam Mercer committed
451 452 453 454
  /*	2^n gives the number of samples of the SFT. This will be a signal in a 	*/
  /*	band fSample/2 with fSample=2*(2.e-4 f0+f0Band). This makes sure that	*/
  /*	that the sampling frequency enables us to produce the DeFT frequency 	*/
  /*	band that we want (f0Band) and that there's enought wings of SFT data	*/
sberukoff's avatar
sberukoff committed
455
  /*	(4.e-4*f0) to produce such DeFT band. 					*/
Adam Mercer's avatar
Adam Mercer committed
456 457 458
  /* 	The main formula used here is that the gives the maximum allowed	*/
  /*	time-baseline (Tmax) such that a signal of frequency f0 Doppler 	*/
  /*	modulated due to Earth spin is seen as monochromatic during such	*/
sberukoff's avatar
sberukoff committed
459 460
  /*	observation time: 							*/
  /* 	  	Tmax < 9.6e4/sqrt(f0). 						*/
Adam Mercer's avatar
Adam Mercer committed
461 462 463
  /*	We have thus taken Tmax=9.4e4/sqrt(f0). A different criterion can be	*/
  /*	chosen by setting the variable "factor" to a suitable value. We have	*/
  /*	rounded the resulting number of samples to the nearest smallest power 	*/
sberukoff's avatar
sberukoff committed
464 465 466 467
  /*	of two. 							       	*/

  {
    REAL8 tempf0Band;
papa's avatar
papa committed
468 469

    ntermsdivbytwo=64;
sberukoff's avatar
sberukoff committed
470 471 472 473 474
    /* compute size of SFT wings */
    fWing = 2.0e-4 * f0;
    /* Adjust search band to include wings */
    f0Band += 2.0*fWing;
    tempf0Band = f0Band;
475
    tSFT=1.0e3;/*9.6e4/sqrt(f0);*/
papa's avatar
papa committed
476
    nDeltaF = 2*ntermsdivbytwo+(INT4)(ceil(f0Band*tSFT));
477
      /* tSFT = (REAL8)nDeltaF/f0Band;*/
sberukoff's avatar
sberukoff committed
478 479 480 481 482 483 484 485 486 487 488 489
    oneOverSqrtTSFT = 1.0/sqrt(tSFT);
    /* Number of data in time series for 1 SFT */
    dfSFT=1.0/tSFT;
    /* Get rid of the wings */
    f0Band -= 2.0*fWing;
    /* the index of the right side of the band, NO WINGS */
    if0Max=ceil((f0+f0Band/2.0)/dfSFT);
    /* the index of the left side of the band, NO WINGS */
    if0Min=floor((f0-f0Band/2.0)/dfSFT);
    /* frequency of the right side of the band, NO WINGS */
    f0Max=dfSFT*if0Max;
    /* frequency of the left side of the band, NO WINGS */
Adam Mercer's avatar
Adam Mercer committed
490 491
    f0Min=dfSFT*if0Min;

sberukoff's avatar
sberukoff committed
492 493
    f0Band = tempf0Band;
  }
papa's avatar
papa committed
494
  /* the hard-wired 64 is for the necessary NTERMS_COH_DIV_TWO*/
sberukoff's avatar
sberukoff committed
495
  /* index of the left side of the band, WITH WINGS */
papa's avatar
papa committed
496 497

  ifMin=floor(f0/dfSFT)-nDeltaF/2;
sberukoff's avatar
sberukoff committed
498
  /* indexof the right side of the band, WITH WINGS */
papa's avatar
papa committed
499
  ifMax=ifMin+nDeltaF;
sberukoff's avatar
sberukoff committed
500
  /* frequency of the left side of the band, WITH WINGS */
Adam Mercer's avatar
Adam Mercer committed
501
  fMin=dfSFT*ifMin;
papa's avatar
papa committed
502
  nDeltaF=2*ntermsdivbytwo+(INT4)(ceil(f0Band/dfSFT));
Adam Mercer's avatar
Adam Mercer committed
503
  /* Number of SFTs which make one DeFT */
sberukoff's avatar
sberukoff committed
504 505 506 507 508 509 510 511 512 513
  mCohSFT=ceil(tCoh/tSFT);
  if (mCohSFT==0) mCohSFT=1;
  /* Coherent search time baseline */
  tCoh=tSFT*mCohSFT;
  /* Number of coherent timescales in total obs. time */
  mObsCoh=floor(tObs/tCoh);
  if (mObsCoh ==0) mObsCoh=1;
  /* Total observation time */
  tObs=tCoh*mObsCoh;
  /* Number of SFTs needed during observation time */
Adam Mercer's avatar
Adam Mercer committed
514
  mObsSFT=mCohSFT*mObsCoh;
sberukoff's avatar
sberukoff committed
515

516 517
  printf("%d %d %d\n",mObsSFT,mObsCoh,mCohSFT);

sberukoff's avatar
sberukoff committed
518 519 520 521 522 523
  /* convert input angles from degrees to radians */
  signalParams->skyP->alpha=signalParams->skyP->alpha*LAL_TWOPI/360.0;
  signalParams->skyP->delta=signalParams->skyP->delta*LAL_TWOPI/360.0;
  templateParams->skyP->alpha=templateParams->skyP->alpha*LAL_TWOPI/360.0;
  templateParams->skyP->delta=templateParams->skyP->delta*LAL_TWOPI/360.0;

524 525 526
  /* Convert signal spindown to conform with GenerateTaylorCW() definition */
  for(i=0;i<signalParams->spind->m;i++){signalParams->spind->spParams[i] /= f0;}

Adam Mercer's avatar
Adam Mercer committed
527
  /***** END USEFUL QUANTITIES *****/
sberukoff's avatar
sberukoff committed
528 529 530


  /***** CALL ROUTINE TO GENERATE TIMESTAMPS *****/
Adam Mercer's avatar
Adam Mercer committed
531

532 533 534 535 536 537
  /*times(tSFT, mObsSFT, timeStamps, deletions);*/
  i=0;
  times2(tSFT, mObsCoh, &timeStamps, &sftPerCoh, deletions, mCohSFT);
  /*
    In order to tell many of the loops in this code how many times to iterate,
     we set a variable which is equal to the last entry in the sftPerCoh array.
Adam Mercer's avatar
Adam Mercer committed
538 539
     Note, that for no gaps, this number is simply mObsSFT; with gaps, it's
     something else.  Regardless, this is the value that the subroutines know
540 541 542 543
     as mObsSFT.
  */
  totalSFT = sftPerCoh[mObsCoh];
    /***** END TIMESTAMPS *****/
jolien's avatar
jolien committed
544 545


sberukoff's avatar
sberukoff committed
546
  /***** CREATE NOISE *****/
jolien's avatar
jolien committed
547

sberukoff's avatar
sberukoff committed
548 549 550 551 552
  if(noi!=NULL)
    {
      LALCreateVector(&status, &noise, (UINT4)nDeltaF*mObsSFT);
      LALCreateVector(&status, &temp, (UINT4)nDeltaF*mObsSFT);
      LALCreateRandomParams(&status, &params, seed);
Adam Mercer's avatar
Adam Mercer committed
553 554 555 556

      /* fill temp vector with normal deviates*/
      LALNormalDeviates(&status, temp, params);

sberukoff's avatar
sberukoff committed
557 558 559
      /* rewrite into COMPLEX8 VECTOR *noise */
      i=0;
      while(i<nDeltaF*mObsSFT)
jolien's avatar
jolien committed
560
	{
sberukoff's avatar
sberukoff committed
561 562
	  noise->data[i]=temp->data[i];
	  i++;
Adam Mercer's avatar
Adam Mercer committed
563
	}
jolien's avatar
jolien committed
564

sberukoff's avatar
sberukoff committed
565 566 567 568 569 570 571 572 573
      /* Destroy structures that were used here, if we can */
      LALDestroyRandomParams(&status, &params);
      LALSDestroyVector(&status, &temp);
    }

  /***** END CREATE NOISE *****/

  /* Quantities computed for barycentering */
  edat=(EphemerisData *)LALMalloc(sizeof(EphemerisData));
574 575
  (*edat).ephiles.earthEphemeris = earthEphemeris;
  (*edat).ephiles.sunEphemeris = sunEphemeris;
sberukoff's avatar
sberukoff committed
576

Adam Mercer's avatar
Adam Mercer committed
577
  /* Read in ephemerides */
sberukoff's avatar
sberukoff committed
578
  LALInitBarycenter(&status, edat);
Adam Mercer's avatar
Adam Mercer committed
579 580
  /*Getting detector coords from DetectorSite module of tools package */

sberukoff's avatar
sberukoff committed
581 582 583
  /* Cached options are:
     LALDetectorIndexLHODIFF, LALDetectorIndexLLODIFF,
     LALDetectorIndexVIRGODIFF, LALDetectorIndexGEO600DIFF,
Adam Mercer's avatar
Adam Mercer committed
584
     LALDetectorIndexTAMA300DIFF,LALDetectorIndexCIT40DIFF
sberukoff's avatar
sberukoff committed
585 586 587 588 589 590 591 592 593 594
  */
  cachedDetector = lalCachedDetectors[LALDetectorIndexGEO600DIFF];
  baryinput.site.location[0]=cachedDetector.location[0]/LAL_C_SI;
  baryinput.site.location[1]=cachedDetector.location[1]/LAL_C_SI;
  baryinput.site.location[2]=cachedDetector.location[2]/LAL_C_SI;
  baryinput.alpha=signalParams->skyP->alpha;
  baryinput.delta=signalParams->skyP->delta;
  baryinput.dInv=0.e0;

  /***** CREATE SIGNAL *****/
Adam Mercer's avatar
Adam Mercer committed
595 596

  /* Set up parameter structure */
sberukoff's avatar
sberukoff committed
597
  /* Source position */
598 599
  genTayParams.position.latitude  = signalParams->skyP->alpha;
  genTayParams.position.longitude = signalParams->skyP->delta;
sberukoff's avatar
sberukoff committed
600 601 602 603 604
  genTayParams.position.system = COORDINATESYSTEM_EQUATORIAL;
  /* Source polarization angle */
  /* Note, that the way we compute the statistic, we don't care what this value is. */
  genTayParams.psi = 0.0;
  /* Source initial phase */
Adam Mercer's avatar
Adam Mercer committed
605
  genTayParams.phi0 = 0.0;
sberukoff's avatar
sberukoff committed
606 607 608 609 610 611 612
  /* Source polarization amplitudes */
  genTayParams.aPlus = aPlus;
  genTayParams.aCross = aCross;
  /* Source intrinsic frequency */
  genTayParams.f0 = f0;
  /* Source spindown parameters: create vector and assign values */
  genTayParams.f = NULL;
Adam Mercer's avatar
Adam Mercer committed
613 614
  LALDCreateVector(&status, &(genTayParams.f), signalParams->spind->m);
  for(i=0;i<signalParams->spind->m;i++){genTayParams.f->data[i] = signalParams->spind->spParams[i];}
sberukoff's avatar
sberukoff committed
615 616 617 618 619
  /* Time resolution for source */
  /* Note, that this needs to be sampled only very coarsely! */
  genTayParams.deltaT = 100.0;
  /* Length should include fudge factor to allow for barycentring */
  genTayParams.length = (tObs+1200.0)/genTayParams.deltaT;
620
  memset(&cwDetector, 0, sizeof(PulsarDetectorResponse));
sberukoff's avatar
sberukoff committed
621 622 623
  /* The ephemerides */
  cwDetector.ephemerides = edat;
  /* Specifying the detector site (set above) */
Adam Mercer's avatar
Adam Mercer committed
624 625
  cwDetector.site = &cachedDetector;
  /* The transfer function.
sberukoff's avatar
sberukoff committed
626 627 628 629 630 631 632 633 634 635 636
   * Note, this xfer function has only two points */
  cwDetector.transfer = (COMPLEX8FrequencySeries *)LALMalloc(sizeof(COMPLEX8FrequencySeries));
  memset(cwDetector.transfer, 0, sizeof(COMPLEX8FrequencySeries));
  cwDetector.transfer->epoch = timeStamps[0];
  cwDetector.transfer->f0 = 0.0;
  cwDetector.transfer->deltaF = 16384.0;
  cwDetector.transfer->data = NULL;
  LALCCreateVector(&status, &(cwDetector.transfer->data), 2);
  /* Allocating space for the eventual time series */
  timeSeries = (REAL4TimeSeries *)LALMalloc(sizeof(REAL4TimeSeries));
  timeSeries->data = (REAL4Vector *)LALMalloc(sizeof(REAL4Vector));
Adam Mercer's avatar
Adam Mercer committed
637

sberukoff's avatar
sberukoff committed
638
  /* The length of the TS will be plenty long; below, we have it set so that
Adam Mercer's avatar
Adam Mercer committed
639 640
   * it is sampled at just better than twice the Nyquist frequency, then
   * increased that to make it a power of two (to facilitate a quick FFT)
sberukoff's avatar
sberukoff committed
641 642 643 644 645 646
   */
  {
    INT4 lenPwr;
    lenPwr = ceil(log(tSFT*2.01*f0)/log(2.0));
    timeSeries->data->length = pow(2.0,lenPwr);
  }
Adam Mercer's avatar
Adam Mercer committed
647

sberukoff's avatar
sberukoff committed
648 649 650 651 652 653 654 655 656
  timeSeries->data->data = (REAL4 *)LALMalloc(timeSeries->data->length*sizeof(REAL4));
  dt = timeSeries->deltaT = tSFT/timeSeries->data->length;

  /* unit response function */
  cwDetector.transfer->data->data[0].re = 1.0;
  cwDetector.transfer->data->data[1].re = 1.0;
  cwDetector.transfer->data->data[0].im = 0.0;
  cwDetector.transfer->data->data[1].im = 0.0;
  /* again, the note about the barycentring */
Adam Mercer's avatar
Adam Mercer committed
657 658
  genTayParams.epoch.gpsSeconds = timeStamps[0].gpsSeconds - 600;
  genTayParams.epoch.gpsNanoSeconds = timeStamps[0].gpsNanoSeconds;
659
  memset(&cgwOutput, 0, sizeof(PulsarCoherentGW));
sberukoff's avatar
sberukoff committed
660

Adam Mercer's avatar
Adam Mercer committed
661
  /*
sberukoff's avatar
sberukoff committed
662
   *
Adam Mercer's avatar
Adam Mercer committed
663
   * OKAY, GENERATE THE SIGNAL @ THE SOURCE
sberukoff's avatar
sberukoff committed
664 665
   *
   */
Adam Mercer's avatar
Adam Mercer committed
666
  LALGenerateTaylorCW(&status, &cgwOutput, &genTayParams);
sberukoff's avatar
sberukoff committed
667 668 669

  {
    INT4 len,len2;
Adam Mercer's avatar
Adam Mercer committed
670

sberukoff's avatar
sberukoff committed
671 672
    len=timeSeries->data->length;
    len2=len/2+1;
Adam Mercer's avatar
Adam Mercer committed
673

sberukoff's avatar
sberukoff committed
674 675
    /* Create vector to hold frequency series */
    LALCCreateVector(&status, &fvec, (UINT4)len2);
Adam Mercer's avatar
Adam Mercer committed
676

sberukoff's avatar
sberukoff committed
677 678
    /* Compute measured plan for FFTW */
    LALCreateForwardRealFFTPlan(&status, &pfwd, (UINT4)len, 0);
Adam Mercer's avatar
Adam Mercer committed
679

sberukoff's avatar
sberukoff committed
680
    /* Allocate memory for the SFTData structure */
Adam Mercer's avatar
Adam Mercer committed
681
    /*
sberukoff's avatar
sberukoff committed
682 683 684 685
     * Note that the length allocated for the data
     * array is the width of the band we're interested in,
     * plus f0*4E-4 for the 'wings', plus a bit of wiggle room.
     */
686 687
    SFTData=(FFT **)LALMalloc(totalSFT*sizeof(FFT *));
    for(i=0;i<totalSFT;i++){
sberukoff's avatar
sberukoff committed
688 689 690 691 692 693 694 695 696 697 698 699
      SFTData[i]=(FFT *)LALMalloc(sizeof(FFT));
      SFTData[i]->fft=(COMPLEX8FrequencySeries *)
	LALMalloc(sizeof(COMPLEX8FrequencySeries));
      SFTData[i]->fft->data=(COMPLEX8Vector *)
	LALMalloc(sizeof(COMPLEX8Vector));
      SFTData[i]->fft->data->data=(COMPLEX8 *)
	LALMalloc((nDeltaF+1)*sizeof(COMPLEX8));
    }

    /* Lots of debugging stuff.  If you want to use these, go ahead, but
     * be warned that these files may be HUGE (10s-100s of GB).
     */
jolien's avatar
jolien committed
700

sberukoff's avatar
sberukoff committed
701
    {
Adam Mercer's avatar
Adam Mercer committed
702 703 704 705
      REAL4Vector *tempTS = NULL;

      LALSCreateVector(&status, &tempTS, (UINT4)len);

706
      for(a=0;a<totalSFT;a++)
jolien's avatar
jolien committed
707
	{
sberukoff's avatar
sberukoff committed
708
	  REAL4Vector *ts = timeSeries->data;
Adam Mercer's avatar
Adam Mercer committed
709
	  /*
sberukoff's avatar
sberukoff committed
710
	   *  Note that this epoch is different from the epoch of GenTay.
Adam Mercer's avatar
Adam Mercer committed
711
	   *  The difference is seconds, a little bit more than the
sberukoff's avatar
sberukoff committed
712 713 714 715
	   *  maximum propagation delay of a signal between the Earth
	   *  and the solar system barycentre.
	   */
	  timeSeries->epoch = timeStamps[a];
Adam Mercer's avatar
Adam Mercer committed
716
	  /*
sberukoff's avatar
sberukoff committed
717 718 719 720
	   *
	   *  OKAY, CONVERT SOURCE SIGNAL TO WHAT SOMEBODY SEES AT THE DETECTOR OUTPUT
	   *
	   */
721
	  LALPulsarSimulateCoherentGW(&status, timeSeries, &cgwOutput, &cwDetector);
Adam Mercer's avatar
Adam Mercer committed
722

sberukoff's avatar
sberukoff committed
723 724 725 726 727
	  /* Write time series of correct size to temp Array, for FFT */
	  for(i=0;i<len;i++)
	    {
	      tempTS->data[i] = ts->data[i];
	    }
Adam Mercer's avatar
Adam Mercer committed
728

sberukoff's avatar
sberukoff committed
729 730
	  /* Perform FFTW-LAL Fast Fourier Transform */
	  LALForwardRealFFT(&status, fvec, tempTS, pfwd);
Adam Mercer's avatar
Adam Mercer committed
731 732

	  {
sberukoff's avatar
sberukoff committed
733
	    INT4 fL = ifMin;
Adam Mercer's avatar
Adam Mercer committed
734
	    INT4 cnt = 0;
sberukoff's avatar
sberukoff committed
735 736 737 738
	    while(fL < ifMin+nDeltaF+1)
	      {
		COMPLEX8 *tempSFT=SFTData[a]->fft->data->data;
		COMPLEX8 *fTemp=fvec->data;
Adam Mercer's avatar
Adam Mercer committed
739

sberukoff's avatar
sberukoff committed
740 741 742
		/* Also normalise */
		tempSFT[cnt].re = fTemp[fL].re * oneOverSqrtTSFT;
		tempSFT[cnt].im = fTemp[fL].im * oneOverSqrtTSFT;
743

sberukoff's avatar
sberukoff committed
744
		cnt++; fL++;
Adam Mercer's avatar
Adam Mercer committed
745

sberukoff's avatar
sberukoff committed
746 747 748 749
	      }

	  }
	  /* assign particulars to each SFT */
Adam Mercer's avatar
Adam Mercer committed
750
	  SFTData[a]->fft->data->length = nDeltaF+1;
sberukoff's avatar
sberukoff committed
751 752 753
	  SFTData[a]->fft->epoch = timeStamps[a];
	  SFTData[a]->fft->f0 = f0Min; /* this is the frequency of the first freq in the band */
	  SFTData[a]->fft->deltaF = dfSFT;
754
	  printf("Created SFT number %d of %d\n.",a, mObsSFT);
Adam Mercer's avatar
Adam Mercer committed
755

jolien's avatar
jolien committed
756
	}
sberukoff's avatar
sberukoff committed
757 758
      LALSDestroyVector(&status, &tempTS);
    }
Adam Mercer's avatar
Adam Mercer committed
759

sberukoff's avatar
sberukoff committed
760 761
    /*
     * Note, we have to destroy the memory that GenTay allocates.
Adam Mercer's avatar
Adam Mercer committed
762
     * This is currently (04.02) not documented!
sberukoff's avatar
sberukoff committed
763
     */
Adam Mercer's avatar
Adam Mercer committed
764

sberukoff's avatar
sberukoff committed
765 766
    if(&(cgwOutput.a) !=NULL) {
      LALSDestroyVectorSequence(&status, &(cgwOutput.a->data));
Adam Mercer's avatar
Adam Mercer committed
767
      LALFree(cgwOutput.a);
sberukoff's avatar
sberukoff committed
768 769 770 771 772 773 774 775 776 777
    }
    if(&(cgwOutput.f) !=NULL) {
      LALSDestroyVector(&status, &(cgwOutput.f->data));
      LALFree(cgwOutput.f);
    }
    if(&(cgwOutput.phi) !=NULL) {
      LALDDestroyVector(&status, &(cgwOutput.phi->data));
      LALFree(cgwOutput.phi);
    }

Adam Mercer's avatar
Adam Mercer committed
778
    if(noi!=NULL) {LALDestroyVector(&status, &noise);}
sberukoff's avatar
sberukoff committed
779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802
    LALCDestroyVector(&status, &fvec);
    LALFree(timeSeries->data->data);
    LALFree(timeSeries->data);
    LALFree(timeSeries);
    LALDestroyRealFFTPlan(&status, &pfwd);
  }
  LALDDestroyVector(&status,&(genTayParams.f));
  LALCDestroyVector(&status, &(cwDetector.transfer->data));
  LALFree(cwDetector.transfer);

  /***** END CREATE SIGNAL *****/


 /* BEGIN AMPLITUDE MODULATION */

  /* Allocate space for amParams stucture */
  /* Here, amParams->das is the Detector and Source info */
  amParams = (AMCoeffsParams *)LALMalloc(sizeof(AMCoeffsParams));
  amParams->das = (LALDetAndSource *)LALMalloc(sizeof(LALDetAndSource));
  amParams->das->pSource = (LALSource *)LALMalloc(sizeof(LALSource));
  /* Fill up AMCoeffsParams structure */
  amParams->baryinput = &baryinput;
  amParams->earth = &earth;
  amParams->edat = edat;
Adam Mercer's avatar
Adam Mercer committed
803
  amParams->das->pDetector = &cachedDetector;
804 805
  amParams->das->pSource->equatorialCoords.latitude = signalParams->skyP->alpha;
  amParams->das->pSource->equatorialCoords.longitude = signalParams->skyP->delta;
sberukoff's avatar
sberukoff committed
806 807 808
  amParams->das->pSource->orientation = 0.0;
  amParams->das->pSource->equatorialCoords.system = COORDINATESYSTEM_EQUATORIAL;
  amParams->polAngle = genTayParams.psi;
809
  amParams->mObsSFT = totalSFT;
sberukoff's avatar
sberukoff committed
810 811 812
 /* Allocate space for AMCoeffs */
  amc.a = NULL;
  amc.b = NULL;
813 814
  LALSCreateVector(&status, &(amc.a), (UINT4)totalSFT);
  LALSCreateVector(&status, &(amc.b), (UINT4)totalSFT);
Adam Mercer's avatar
Adam Mercer committed
815 816 817

 /*
  * Compute timestamps for middle of each ts chunk
sberukoff's avatar
sberukoff committed
818
  * Note, we have decided to use the values of 'a' and 'b'
Adam Mercer's avatar
Adam Mercer committed
819
  * that correspond to the midpoint of the timeSeries for
sberukoff's avatar
sberukoff committed
820 821 822
  * which the SFT denotes.  In practice, this can be a flaw
  * since 'a' is a sinusoid with T=45ksec and T_b=90ksec; both
  * of these have amplitude about 0.4.  Thus, for 'a', on a timescale
Adam Mercer's avatar
Adam Mercer committed
823
  * of 10ksec, the value changes significantly.
sberukoff's avatar
sberukoff committed
824 825
  */
  {
826 827 828 829 830 831 832 833 834
    LIGOTimeGPS *midTS;
	REAL8 t;

	midTS = (LIGOTimeGPS *)LALCalloc(totalSFT,sizeof(LIGOTimeGPS));
    /* compute timestamps of end and beg */
    amParams->ts0 = timeStamps[0];
    TimeToFloat(&t, &(timeStamps[0]));
    t+=tObs;
    FloatToTime(&(amParams->tsEnd), &t);
Adam Mercer's avatar
Adam Mercer committed
835

836
    for(k=0; k<totalSFT; k++)
sberukoff's avatar
sberukoff committed
837 838
      {
	REAL8 teemp=0.0;
Adam Mercer's avatar
Adam Mercer committed
839

sberukoff's avatar
sberukoff committed
840
	TimeToFloat(&teemp, &(timeStamps[k]));
841
	teemp += 0.5*tSFT;
sberukoff's avatar
sberukoff committed
842 843 844 845
	FloatToTime(&(midTS[k]), &teemp);
      }
    /* Compute the AM coefficients */
    LALComputeAM(&status, &amc, midTS, amParams);
846
    LALFree(midTS);
sberukoff's avatar
sberukoff committed
847 848 849
  }

  /***** DEMODULATE SIGNAL *****/
Adam Mercer's avatar
Adam Mercer committed
850

sberukoff's avatar
sberukoff committed
851 852
 /* Allocate space and set quantity values for demodulation parameters */
  demParams=(DemodPar *)LALMalloc(sizeof(DemodPar));
Adam Mercer's avatar
Adam Mercer committed
853
  demParams->skyConst=(REAL8 *)LALMalloc((2*templateParams->spind->m *
854
					  (totalSFT+1)+2*totalSFT+3)*sizeof(REAL8));
sberukoff's avatar
sberukoff committed
855 856 857 858 859 860 861 862
  demParams->spinDwn=(REAL8 *)LALMalloc(templateParams->spind->m*sizeof(REAL8));
  demParams->if0Max=if0Max;
  demParams->if0Min=if0Min;
  demParams->mCohSFT=mCohSFT;
  demParams->mObsCoh=mObsCoh;
  demParams->ifMin=ifMin;
  demParams->spinDwnOrder=templateParams->spind->m;
  demParams->amcoe = &amc;
863 864 865 866 867
  demParams->sftPerCoh = sftPerCoh;
  for(i=0;i<signalParams->spind->m;i++)
    {
      demParams->spinDwn[i]=templateParams->spind->spParams[i];
    }
Adam Mercer's avatar
Adam Mercer committed
868

869
  /* Allocate space and set quantities for call to LALComputeSky() */
870 871 872
  csParams=(CSParams *)LALMalloc(sizeof(CSParams));
  csParams->skyPos=(REAL8 *)LALMalloc(2*sizeof(REAL8));
  csParams->skyPos[0]=templateParams->skyP->alpha;
sberukoff's avatar
sberukoff committed
873 874 875
 csParams->skyPos[1]=templateParams->skyP->delta;
 csParams->tGPS=timeStamps;
 csParams->spinDwnOrder=templateParams->spind->m;
876
 csParams->mObsSFT=totalSFT;
sberukoff's avatar
sberukoff committed
877
 csParams->tSFT=tSFT;
Adam Mercer's avatar
Adam Mercer committed
878
 csParams->edat=edat;
sberukoff's avatar
sberukoff committed
879 880 881 882
 csParams->emit=&emit;
 csParams->earth=&earth;
 csParams->baryinput=&baryinput;

883
 iSkyCoh=0;
Adam Mercer's avatar
Adam Mercer committed
884

sberukoff's avatar
sberukoff committed
885
  /* Call COMPUTESKY() */
886
  LALComputeSky(&status, demParams->skyConst, iSkyCoh, csParams);
Adam Mercer's avatar
Adam Mercer committed
887

sberukoff's avatar
sberukoff committed
888 889 890
  /* Deallocate space for ComputeSky parameters */
  LALFree(csParams->skyPos);
  LALFree(csParams);
Adam Mercer's avatar
Adam Mercer committed
891

sberukoff's avatar
sberukoff committed
892 893 894 895 896
  /* Allocate memory for demodulated data */
  xHat=(DeFTPeriodogram **)LALMalloc(mObsCoh*sizeof(DeFTPeriodogram *));
  for(i=0; i<mObsCoh; i++)
    {
      xHat[i]=(DeFTPeriodogram *)LALMalloc(sizeof(DeFTPeriodogram));
Adam Mercer's avatar
Adam Mercer committed
897
      xHat[i]->fft=(REAL8FrequencySeries *)
sberukoff's avatar
sberukoff committed
898
	LALMalloc(sizeof(REAL8FrequencySeries));
899
      xHat[i]->fA=(COMPLEX16FrequencySeries *)LALMalloc(sizeof(COMPLEX16FrequencySeries));
Adam Mercer's avatar
Adam Mercer committed
900
      xHat[i]->fB=(COMPLEX16FrequencySeries *)LALMalloc(sizeof(COMPLEX16FrequencySeries));
sberukoff's avatar
sberukoff committed
901 902 903
      xHat[i]->fft->data=(REAL8Vector *)LALMalloc(sizeof(REAL8Vector));
      xHat[i]->fft->data->data=(REAL8 *)LALMalloc((UINT4)((if0Max-if0Min+1)*mCohSFT)*sizeof(REAL8));
      xHat[i]->fft->data->length=(UINT4)((if0Max-if0Min+1)*mCohSFT);
904
      xHat[i]->fA->data=(COMPLEX16Vector *)LALMalloc(sizeof(COMPLEX16Vector));
Adam Mercer's avatar
Adam Mercer committed
905
      xHat[i]->fA->data->data=(COMPLEX16 *)LALMalloc((UINT4)((if0Max-if0Min+1)*mCohSFT)*sizeof(COMPLEX16));
906 907
      xHat[i]->fA->data->length=(UINT4)((if0Max-if0Min+1)*mCohSFT);
      xHat[i]->fB->data=(COMPLEX16Vector *)LALMalloc(sizeof(COMPLEX16Vector));
Adam Mercer's avatar
Adam Mercer committed
908
      xHat[i]->fB->data->data=(COMPLEX16 *)LALMalloc((UINT4)((if0Max-if0Min+1)*mCohSFT)*sizeof(COMPLEX16));
909
      xHat[i]->fB->data->length=(UINT4)((if0Max-if0Min+1)*mCohSFT);
sberukoff's avatar
sberukoff committed
910
    }
Adam Mercer's avatar
Adam Mercer committed
911

sberukoff's avatar
sberukoff committed
912
  for(k=0; k<mObsCoh; k++)
sberukoff's avatar
sberukoff committed
913
    {
sberukoff's avatar
sberukoff committed
914
      demParams->iCoh=k;
Adam Mercer's avatar
Adam Mercer committed
915

sberukoff's avatar
sberukoff committed
916 917 918
      /**************************/
      /*       DEMODULATE       */
      /**************************/
919 920 921
      Fstat.F = xHat+k;
      /*      LALDemod(&status, *(xHat+k), SFTData, demParams); */
      LALDemod(&status, &Fstat, SFTData, demParams);
922

923
    }
Adam Mercer's avatar
Adam Mercer committed
924

925
/***** END DEMODULATION *****/
Adam Mercer's avatar
Adam Mercer committed
926 927


sberukoff's avatar
sberukoff committed
928 929 930 931 932 933 934 935
  /***** DEALLOCATION *****/

  /* Deallocate AM  */
  LALSDestroyVector(&status, &(amc.a));
  LALSDestroyVector(&status, &(amc.b));
  LALFree(amParams->das->pSource);
  LALFree(amParams->das);
  LALFree(amParams);
Adam Mercer's avatar
Adam Mercer committed
936

sberukoff's avatar
sberukoff committed
937
  /* Deallocate SFTData structure, since we don't need it anymore */
938
  for(i=0;i<totalSFT;i++)
sberukoff's avatar
sberukoff committed
939 940 941 942 943 944 945
    {
      LALFree(SFTData[i]->fft->data->data);
      LALFree(SFTData[i]->fft->data);
      LALFree(SFTData[i]->fft);
      LALFree(SFTData[i]);
    }
  LALFree(SFTData);
Adam Mercer's avatar
Adam Mercer committed
946

sberukoff's avatar
sberukoff committed
947 948 949 950 951 952
  /* Deallocate memory used by demodulated data */
  for(i=0; i<mObsCoh; i++)
    {
      LALFree(xHat[i]->fft->data->data);
      LALFree(xHat[i]->fft->data);
      LALFree(xHat[i]->fft);
953 954 955 956 957 958
      LALFree(xHat[i]->fA->data->data);
      LALFree(xHat[i]->fA->data);
      LALFree(xHat[i]->fA);
      LALFree(xHat[i]->fB->data->data);
      LALFree(xHat[i]->fB->data);
      LALFree(xHat[i]->fB);
sberukoff's avatar
sberukoff committed
959 960 961
      LALFree(xHat[i]);
    }
  LALFree(xHat);
Adam Mercer's avatar
Adam Mercer committed
962

sberukoff's avatar
sberukoff committed
963 964 965 966 967 968 969 970 971 972 973 974 975 976 977
  /* Deallocate template and signal params */
  LALFree(templateParams->spind->spParams);
  LALFree(templateParams->spind);
  LALFree(templateParams->skyP);
  LALFree(templateParams);

  LALFree(signalParams->spind->spParams);
  LALFree(signalParams->spind);
  LALFree(signalParams->skyP);
  LALFree(signalParams);

  /* Deallocate demodulation params */
  LALFree(demParams->skyConst);
  LALFree(demParams->spinDwn);
  LALFree(demParams);
Adam Mercer's avatar
Adam Mercer committed
978

sberukoff's avatar
sberukoff committed
979 980 981
  LALFree(basicInputsFile);
  /* Anything else */

982 983 984 985
  /* Deallocate timestamps */
  LALFree(timeStamps);
  LALFree(sftPerCoh);

sberukoff's avatar
sberukoff committed
986 987 988
  LALFree(edat->ephemS);
  LALFree(edat->ephemE);
  LALFree(edat);
Adam Mercer's avatar
Adam Mercer committed
989
  LALCheckMemoryLeaks();
sberukoff's avatar
sberukoff committed
990
  return 0;
jolien's avatar
jolien committed
991 992 993
}


sberukoff's avatar
sberukoff committed
994
/***** This is the routine which computes the timestamps *****/