GeneratePulsarSignal.c

Go to the documentation of this file.

/*
 * Copyright (C) 2004, 2005 Reinhard Prix
 * Copyright (C) 2004, 2005 Greg Mendell
 *
 *  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
 */

/* NOTES: */
/* 07/14/04 gam; add functions LALComputeSkyAndZeroPsiAMResponse and LALFastGeneratePulsarSFTs */
/* 10/08/04 gam; fix indexing into trig lookup tables (LUTs) by having table go from -2*pi to 2*pi */
/* 10/12/04 gam; When computing fCross and fPlus need to use 2.0*psi. */
/* 09/07/05 gam; Add Dterms parameter to LALFastGeneratePulsarSFTs; use this to fill in SFT bins with fake data as per LALDemod else fill in bin with zero */

/**
 * \author Reinhard Prix, Greg Mendell
 * \date 2005
 * \file 
 * \ingroup moduleGeneratePulsarSignal
 *
 * \brief Module to generate simulated pulsar-signals.
 *
 * $Id: GeneratePulsarSignal.c,v 1.57 2008/08/14 17:04:11 evang Exp $
 *
 */
#include <math.h>

#include <lal/AVFactories.h>
#include <lal/SeqFactories.h>
#include <lal/RealFFT.h>
#include <lal/SFTutils.h>
#include <lal/Window.h>

#include <lal/GeneratePulsarSignal.h>

/*----------------------------------------------------------------------*/
/* Internal helper functions */
static int check_timestamp_bounds (const LIGOTimeGPSVector *timestamps, LIGOTimeGPS t0, LIGOTimeGPS t1);
static void checkNoiseSFTs (LALStatus *, const SFTVector *sfts, REAL8 f0, REAL8 f1, REAL8 deltaF);
static void correct_phase (LALStatus *, SFTtype *sft, LIGOTimeGPS tHeterodyne);

/*----------------------------------------------------------------------*/

NRCSID( GENERATEPULSARSIGNALC, "$Id: GeneratePulsarSignal.c,v 1.57 2008/08/14 17:04:11 evang Exp $");

extern INT4 lalDebugLevel;

static REAL8 eps = 1.e-14;      /* maximal REAL8 roundoff-error (used for determining if some number is an INT) */

/* ----- DEFINES ----- */
#define GPS2REAL(gps) ((gps).gpsSeconds + 1e-9 * (gps).gpsNanoSeconds )

/*---------- Global variables ----------*/
/* empty init-structs for the types defined in here */
static LALStatus emptyStatus;   
static SpinOrbitCWParamStruc emptyCWParams;
static CoherentGW emptySignal;

const PulsarSignalParams empty_PulsarSignalParams;
const SFTParams empty_SFTParams;
const SFTandSignalParams empty_SFTandSignalParams;


/** Generate a time-series at the detector for a given pulsar.
 */
void
LALGeneratePulsarSignal (LALStatus *status, 
                         REAL4TimeSeries **signal,         /**< output time-series */
                         const PulsarSignalParams *params) /**< input params */
{ 
  SpinOrbitCWParamStruc sourceParams = emptyCWParams; 
  CoherentGW sourceSignal = emptySignal;
  DetectorResponse detector;
  REAL8 SSBduration;
  LIGOTimeGPS t0, t1, tmpTime;
  REAL4TimeSeries *output;
  UINT4 i;

  INITSTATUS( status, "LALGeneratePulsarSignal", GENERATEPULSARSIGNALC );
  ATTATCHSTATUSPTR (status);

  ASSERT (signal != NULL, status, GENERATEPULSARSIGNALH_ENULL, GENERATEPULSARSIGNALH_MSGENULL);
  ASSERT (*signal == NULL, status,  GENERATEPULSARSIGNALH_ENONULL,  GENERATEPULSARSIGNALH_MSGENONULL);

  /*----------------------------------------------------------------------
   *
   * First call GenerateSpinOrbitCW() to generate the source-signal
   *
   *----------------------------------------------------------------------*/
  sourceParams.psi = params->pulsar.psi;
  sourceParams.aPlus = params->pulsar.aPlus;
  sourceParams.aCross = params->pulsar.aCross;
  sourceParams.phi0 = params->pulsar.phi0;
  sourceParams.f0 = params->pulsar.f0;
  /* set source position: make sure it's "normalized", i.e. [0<=alpha<2pi]x[-pi/2<=delta<=pi/2] */
  TRY( LALNormalizeSkyPosition(status->statusPtr, &(sourceParams.position), &(params->pulsar.position)), status);

  /* if pulsar is in binary-orbit, set binary parameters */
  if (params->orbit)
    {
      /*------------------------------------------------------------ */
      /* temporary fix for comparison with Chris' code */
      /*
        TRY (LALConvertGPS2SSB (status->statusPtr, &tmpTime, params->orbit->orbitEpoch, params), status);
        sourceParams.orbitEpoch = tmpTime;
      */
      sourceParams.orbitEpoch =  params->orbit->tp;
      sourceParams.omega = params->orbit->argp;
      /* ------- here we do conversion to Teviets preferred variables -------*/
      sourceParams.rPeriNorm = params->orbit->asini*(1.0 - params->orbit->ecc);
      sourceParams.oneMinusEcc = 1.0 - params->orbit->ecc;
      sourceParams.angularSpeed = (LAL_TWOPI/params->orbit->period)*sqrt((1.0 +params->orbit->ecc)/pow((1.0 - params->orbit->ecc),3.0));
    }
  else
    sourceParams.rPeriNorm = 0.0;               /* this defines an isolated pulsar */

  if ( params->pulsar.refTime.gpsSeconds != 0)
    sourceParams.spinEpoch = params->pulsar.refTime;   /* pulsar reference-time in SSB frame (TDB) */
  else  /* if not given: use startTime converted to SSB as tRef ! */
    {
      TRY ( LALConvertGPS2SSB(status->statusPtr, &tmpTime, params->startTimeGPS, params), status);
      sourceParams.spinEpoch = tmpTime;
    }

  /* sampling-timestep and length for source-parameters */
  /* in seconds; hardcoded; was 60s in makefakedata_v2,
   * but for fast binaries (e.g. SCO-X1) we need faster sampling 
   * This does not seem to affect performance a lot (~4% in makefakedata),
   * but we'll nevertheless make this sampling faster for binaries and slower
   * for isolated pulsars */
  if (params->orbit)
    sourceParams.deltaT = 5;    /* for binaries */
  else
    sourceParams.deltaT = 60;   /* for isolated pulsars */

  /* start-time in SSB time */
  TRY (LALConvertGPS2SSB(status->statusPtr, &t0, params->startTimeGPS, params), status);
  t0.gpsSeconds -= (UINT4)sourceParams.deltaT; /* start one time-step earlier to be safe */

  /* end time in SSB */
  t1 = params->startTimeGPS;
  TRY ( LALAddFloatToGPS(status->statusPtr, &t1, &t1, params->duration), status);
  TRY (LALConvertGPS2SSB(status->statusPtr, &t1, t1, params), status);   /* convert time to SSB */

  /* get duration of source-signal */
  TRY (LALDeltaFloatGPS(status->statusPtr, &SSBduration, &t1, &t0), status);
  SSBduration += 2.0 * sourceParams.deltaT; /* add two time-steps to be safe */

  sourceParams.epoch = t0; 
  sourceParams.length = (UINT4) ceil( SSBduration / sourceParams.deltaT );

  /* we use frequency-spindowns, but GenerateSpinOrbitCW wants f_k = fkdot / (f0 * k!) */
  sourceParams.f = NULL;
  if (params->pulsar.spindown)
    {
      UINT4 kFact = 1;
      TRY ( LALDCreateVector(status->statusPtr, &(sourceParams.f), params->pulsar.spindown->length), status);
      for (i=0; i < sourceParams.f->length; i++ )
        {
          sourceParams.f->data[i] = params->pulsar.spindown->data[i] / (kFact * params->pulsar.f0);
          kFact *= i + 2;
        }
    }

  /* finally, call the function to generate the source waveform */
  TRY ( LALGenerateSpinOrbitCW(status->statusPtr, &sourceSignal, &sourceParams), status);

  /* free spindown-vector right away, so we don't forget */
  if (sourceParams.f) {
    TRY ( LALDDestroyVector(status->statusPtr, &(sourceParams.f) ), status);
  }

  /* check that sampling interval was short enough */
  if ( sourceParams.dfdt > 2.0 )  /* taken from makefakedata_v2 */
    {
      LALPrintError ("GenerateSpinOrbitCW() returned df*dt = %f > 2.0", sourceParams.dfdt);
      ABORT (status, GENERATEPULSARSIGNALH_ESAMPLING, GENERATEPULSARSIGNALH_MSGESAMPLING);
    }

  /*----------------------------------------------------------------------
   *
   * Now call the function to translate the source-signal into a (heterodyned)
   * signal at the detector 
   *
   *----------------------------------------------------------------------*/
  /* first set up the detector-response */
  detector.transfer = params->transfer;
  detector.site = params->site;
  detector.ephemerides = params->ephemerides;

  /* *contrary* to makefakedata_v2, we use the GPS start-time of the timeseries as
   * the heterodyning epoch, in order to make sure that the global phase of the 
   * SFTs it correct: this is necessary to allow correct parameter-estimation on the 
   * pulsar signal-phase in heterodyned SFTs 
   */
  detector.heterodyneEpoch = params->startTimeGPS;

  /* ok, we  need to prepare the output time-series */
  if ( (output = LALCalloc (1, sizeof (*output) )) == NULL) {
    ABORT (status,  GENERATEPULSARSIGNALH_EMEM,  GENERATEPULSARSIGNALH_MSGEMEM);
  }

  /* NOTE: a timeseries of length N*dT has no timestep at N*dT !! (convention) */
  LALCreateVector(status->statusPtr, &(output->data), (UINT4) ceil( params->samplingRate * params->duration) );
  BEGINFAIL(status) {
    LALFree (output);
  } ENDFAIL(status);

  output->deltaT = 1.0 / params->samplingRate;
  output->f0 = params->fHeterodyne;
  output->epoch = params->startTimeGPS;

  
  TRY ( LALSimulateCoherentGW(status->statusPtr, output, &sourceSignal, &detector ), status );

  /* set 'name'-field of timeseries to contain the right "channel prefix" for the detector */
  {
    CHAR *name = XLALGetChannelPrefix ( params->site->frDetector.name );
    if ( name == NULL ) {
      ABORT (status, GENERATEPULSARSIGNALH_EDETECTOR, GENERATEPULSARSIGNALH_MSGEDETECTOR );
    }
    strcpy ( output->name, name );
    LALFree ( name );
  }
                               
  /*----------------------------------------------------------------------*/
  /* Free all allocated memory that is not returned */
  TRY (LALSDestroyVectorSequence( status->statusPtr, &(sourceSignal.a->data)), status);
  LALFree( sourceSignal.a );
  TRY (LALSDestroyVector(status->statusPtr, &(sourceSignal.f->data ) ), status);
  LALFree(sourceSignal.f);
  TRY (LALDDestroyVector(status->statusPtr, &(sourceSignal.phi->data )), status);
  LALFree(sourceSignal.phi);

  *signal = output;

  DETATCHSTATUSPTR (status);
  RETURN (status);

} /* LALGeneratePulsarSignal() */



/** Turn the given time-series into (v2-)SFTs and add noise if given.
 */
void
LALSignalToSFTs (LALStatus *status, 
                 SFTVector **outputSFTs,        /**< [out] SFT-vector */
                 const REAL4TimeSeries *signal, /**< input time-series */
                 const SFTParams *params)       /**< params for output-SFTs */
{
  UINT4 numSFTs;                        /* number of SFTs */
  UINT4 numTimesteps;                   /* number of time-samples in an Tsft */
  UINT4 numBins;                        /* number of frequency-bins in SFT */
  UINT4 iSFT;
  UINT4 idatabin;
  REAL8 Band, f0, deltaF, dt;
  LIGOTimeGPSVector *timestamps = NULL;
  REAL4Vector timeStretch = {0,0};
  LIGOTimeGPS tStart;                   /* start time of input time-series */
  LIGOTimeGPS tLast;                    /* start-time of last _possible_ SFT */
  LIGOTimeGPS tmpTime;
  REAL8 duration, delay;
  UINT4 index0n;                        /* first frequency-bin to use from noise-SFT */
  SFTtype *thisSFT, *thisNoiseSFT;      /* SFT-pointers */
  SFTVector *sftvect = NULL;            /* return value. For better readability */
  UINT4 j;
  RealFFTPlan *pfwd = NULL;             /* FFTW plan */
  LIGOTimeGPS tPrev;                    /* real timstamp of previous SFT */
  UINT4 totalIndex;                     /* timestep-index to start next FFT from */
  INT4 relIndexShift;                   /* relative index-shift from previous SFT (number of timesteps) */

  INITSTATUS( status, "LALSignalToSFTs", GENERATEPULSARSIGNALC);
  ATTATCHSTATUSPTR( status );
  
  ASSERT (outputSFTs != NULL, status, GENERATEPULSARSIGNALH_ENULL, GENERATEPULSARSIGNALH_MSGENULL);
  ASSERT (*outputSFTs == NULL, status,  GENERATEPULSARSIGNALH_ENONULL,  GENERATEPULSARSIGNALH_MSGENONULL);
  ASSERT (signal != NULL, status, GENERATEPULSARSIGNALH_ENULL, GENERATEPULSARSIGNALH_MSGENULL);
  ASSERT (params != NULL, status, GENERATEPULSARSIGNALH_ENULL, GENERATEPULSARSIGNALH_MSGENULL);

  /* UPGRADING switch: complain loudly if user didn't set 'make_v2SFTs' to encourage upgrading */
  if (  params-> make_v2SFTs != 1 )
    {
      fprintf (stderr, "\n********************************************************************************\n\n");
      fprintf (stderr, "     WARNING: LALSignalToSFTs() now returns properly *v2-normalized* SFTs\n");
      fprintf (stderr, "    (see  http://www.lsc-group.phys.uwm.edu/lal/slug/nightly/doxygen/html/group__SFTfileIO.html \n");
      fprintf (stderr, "    for more details on what that means).\n\n");
      fprintf (stderr, "    Please adapt your code correspondingly and set SFTParams.make_v2SFTs=1 to acknowledge this.\n");
      fprintf (stderr, "    This parameter and warning will be removed once the transition is complete.\n\n");
      fprintf (stderr, "********************************************************************************\n\n");
    }

  f0 = signal->f0;                              /* lowest frequency */
  dt = signal->deltaT;          /* timeseries timestep */
  Band = 1.0 / (2.0 * dt);              /* NOTE: frequency-band is determined by sampling-rate! */
  deltaF = 1.0 / params->Tsft;                  /* frequency-resolution */

  /* if noiseSFTs are given: check they are consistent with signal! */
  if (params->noiseSFTs) {
    TRY (checkNoiseSFTs(status->statusPtr, params->noiseSFTs, f0, f0 + Band, deltaF), status);
  }
    
  /* make sure that number of timesamples/SFT is an integer (up to possible rounding errors) */
  {
    REAL8 timesteps = params->Tsft / dt;                /* this is a float!*/
    numTimesteps = (UINT4) (timesteps + 0.5);           /* round to closest int */
    ASSERT ( fabs(timesteps - numTimesteps)/timesteps < eps, status, 
             GENERATEPULSARSIGNALH_EINCONSBAND, GENERATEPULSARSIGNALH_MSGEINCONSBAND);
  }
  /* Prepare FFT: compute plan for FFTW */
  TRY (LALCreateForwardRealFFTPlan(status->statusPtr, &pfwd, numTimesteps, 0), status);         

  /* get some info about time-series */
  tStart = signal->epoch;                                       /* start-time of time-series */

  /* get last possible start-time for an SFT */
  duration =  (UINT4) (1.0* signal->data->length * dt + 0.5); /* total duration rounded to seconds */
  TRY ( LALAddFloatToGPS(status->statusPtr, &tLast, &tStart, duration - params->Tsft), status);

  /* for simplicity we _always_ work with timestamps. 
   * Therefore, we have to generate them now if none have been provided by the user. */
  if (params->timestamps == NULL) 
    {
      LIGOTimeGPS lastTs;
      TRY(LALMakeTimestamps(status->statusPtr, &timestamps, tStart, duration, params->Tsft ), status);
      /* see if the last timestamp is valid (can fit a full SFT in?), if not, drop it */
      lastTs = timestamps->data[timestamps->length-1];
      if ( GPS2REAL(lastTs) > GPS2REAL(tLast) )
        timestamps->length --;
    }
  else  /* if given, use those, and check they are valid */
    {
      timestamps = params->timestamps;
    }
  
  /* check that all timestamps lie within [tStart, tLast] */
  if ( check_timestamp_bounds(timestamps, tStart, tLast) != 0) {
    ABORT (status, GENERATEPULSARSIGNALH_ETIMEBOUND, GENERATEPULSARSIGNALH_MSGETIMEBOUND);
  }

  /* check that we have the right number of noise-SFTs */
  if (  params->noiseSFTs && ( params->noiseSFTs->length != timestamps->length)  ) {
    ABORT ( status, GENERATEPULSARSIGNALH_ENUMSFTS,  GENERATEPULSARSIGNALH_MSGENUMSFTS );
  }

  /* prepare SFT-vector for return */
  numSFTs = timestamps->length;                 /* number of SFTs to produce */
  numBins = (UINT4)(numTimesteps/2) + 1;                /* number of frequency-bins per SFT */

  LALCreateSFTVector (status->statusPtr, &sftvect, numSFTs, numBins);
  BEGINFAIL (status) {
    if (params->timestamps == NULL)
      LALDestroyTimestampVector(status->statusPtr, &timestamps);
  } ENDFAIL (status);

  tPrev = tStart;       /* initialize */
  totalIndex = 0;       /* start from first timestep by default */

  /* main loop: apply FFT the requested time-stretches */
  for (iSFT = 0; iSFT < numSFTs; iSFT++)
    {
      REAL8 realDelay;

      thisSFT = &(sftvect->data[iSFT]); /* point to current SFT-slot */

      /* find the start-bin for this SFT in the time-series */
      TRY ( LALDeltaFloatGPS(status->statusPtr, &delay, &(timestamps->data[iSFT]), &tPrev), status);
      /* round properly: picks *closest* timestep (==> "nudging") !!  */
      relIndexShift = (INT4) (delay / signal->deltaT + 0.5);    
      totalIndex += relIndexShift;

      timeStretch.length = numTimesteps;
      timeStretch.data = signal->data->data + totalIndex; /* point to the right sample-bin */

      /* fill the header of the i'th output SFT */
      realDelay = (REAL4)(relIndexShift * signal->deltaT);  /* avoid rounding-errors*/
      TRY (LALAddFloatToGPS(status->statusPtr, &tmpTime,&tPrev, realDelay),status);

      strcpy ( thisSFT->name, signal->name );
      /* set the ACTUAL timestamp! (can be different from requested one ==> "nudging") */
      thisSFT->epoch = tmpTime;                 
      thisSFT->f0 = signal->f0;                 /* minimum frequency */
      thisSFT->deltaF = 1.0 / params->Tsft;     /* frequency-spacing */

      tPrev = tmpTime;                          /* prepare next loop */

      /* ok, issue at least a warning if we have "nudged" an SFT-timestamp */
      if (lalDebugLevel > 0)
        {
          REAL8 diff;
          TRY (LALDeltaFloatGPS(status->statusPtr, &diff, &(timestamps->data[iSFT]),&tmpTime),status);
          if (diff != 0) 
            {
              LALPrintError ("Warning: timestamp %d had to be 'nudged' by %e s to fit"
                             "with time-series\n", iSFT, diff);
              /* double check if magnitude of nudging seems reasonable .. */
              if ( fabs(diff) >= signal->deltaT ) 
                {
                  LALPrintError ("WARNING: nudged by more than deltaT=%e... "
                                 "this sounds wrong! (We better stop)\n");
                  ABORT (status, GENERATEPULSARSIGNALH_ENULL, GENERATEPULSARSIGNALH_MSGENULL );
                }
            } /* if nudging */
        } /* if lalDebugLevel */

      /* Now window the current time series stretch, if necessary */
      if ( params->window ) {
          const float A = 1.0 / sqrt(params->window->sumofsquares / params->window->data->length);
          for( idatabin = 0; idatabin < timeStretch.length; idatabin++ ) {
              timeStretch.data[idatabin] *= A * params->window->data->data[idatabin];
          }
      }

      /* the central step: FFT the ith time-stretch into an SFT-slot */
      LALForwardRealFFT(status->statusPtr, thisSFT->data, &timeStretch, pfwd);
      BEGINFAIL(status) {
        LALDestroySFTVector(status->statusPtr, &sftvect);
      } ENDFAIL(status);


      /* normalize DFT-data to conform to v2 ( ie. COMPLEX8FrequencySeries ) specification ==> multiply DFT by dt */
      {
        UINT4 i;
        COMPLEX8 *data = &( thisSFT->data->data[0] );
        for ( i=0; i < numBins ; i ++ )
        {
          data->re *= dt;
          data->im *= dt;
          data ++;
        } /* for i < numBins */
      } /* normalize data */
      
      /* correct heterodyning-phase, IF NECESSARY */
      if ( ( (INT4)signal->f0 != signal->f0  )
           || (signal->epoch.gpsNanoSeconds != 0)
           || (thisSFT->epoch.gpsNanoSeconds != 0) )
        {
          /* theterodyne = signal->epoch!*/
          correct_phase(status->statusPtr, thisSFT, signal->epoch);
          BEGINFAIL (status) {
            LALDestroySFTVector(status->statusPtr, &sftvect);
          } ENDFAIL (status);
        } /* if phase-correction necessary */

      /* Now add the noise-SFTs if given */
      if (params->noiseSFTs)
        {
          COMPLEX8 *data, *noise;
          thisNoiseSFT = &(params->noiseSFTs->data[iSFT]);
          index0n = (UINT4) ((thisSFT->f0 - thisNoiseSFT->f0) / thisSFT->deltaF);

          data = &(thisSFT->data->data[0]);
          noise = &(thisNoiseSFT->data->data[index0n]);
          for (j=0; j < numBins; j++)
            {
              data->re += noise->re;
              data->im += noise->im;
              data++;
              noise++;
            } /* for j < numBins */
        
        } /* if noiseSFTs */

    } /* for iSFT < numSFTs */ 

  /* free stuff */
  LALDestroyRealFFTPlan(status->statusPtr, &pfwd);

  /* did we get timestamps or did we make them? */
  if (params->timestamps == NULL)
    {
      LALFree(timestamps->data);
      LALFree(timestamps);
      timestamps = NULL;
    }

  *outputSFTs = sftvect;

  DETATCHSTATUSPTR( status );
  RETURN (status);

} /* LALSignalToSFTs() */



/* 07/14/04 gam */
/**
 * Wrapper for LALComputeSky() and  LALComputeDetAMResponse() that finds the sky
 * constants and \f$F_+\f$ and \f$F_\times\f$ for use with LALFastGeneratePulsarSFTs().
 *  Uses the output of LALComputeSkyAndZeroPsiAMResponse() and the same inputs as
 * LALGeneratePulsarSignal() and LALSignalToSFTs().
 * This function used LALComputeSkyBinary() if params->pSigParams->orbit is not
 * NULL, else it uses LALComputeSky() to find the skyConsts.
 * NOTE THAT THIS FUNCTION COMPUTES \f$F_+\f$ and \f$F_x\f$ for ZERO Psi!!!
 * LALFastGeneratePulsarSFTs() used these to find \f$F_+\f$ and \f$F_x\f$ for NONZERO Psi.
 */
void
LALComputeSkyAndZeroPsiAMResponse (LALStatus *status,
                                   SkyConstAndZeroPsiAMResponse *output,        /**< output */
                                   const SFTandSignalParams *params)            /**< params */
{
  INT4 i;
  INT4 numSFTs;                      /* number of SFTs */
  BarycenterInput baryinput;         /* Stores detector location and other barycentering data */  
  CSParams *csParams   = NULL;       /* ComputeSky parameters */
  CSBParams *csbParams = NULL;       /* ComputeSkyBinary parameters */
  SkyPosition tmp;    
  EarthState earth;
  EmissionTime emit;
  LALDetAMResponse response;  /* output of LALComputeDetAMResponse */
  LALDetAndSource      *das;  /* input for LALComputeDetAMResponse */
  LALGPSandAcc   timeAndAcc;  /* input for LALComputeDetAMResponse */
  REAL8 halfTsft;             /* half the time of one SFT */
  LIGOTimeGPS midTS;          /* midpoint time for an SFT */

  INITSTATUS( status, "LALComputeSkyAndZeroPsiAMResponse", GENERATEPULSARSIGNALC);
  ATTATCHSTATUSPTR( status );
  
  numSFTs = params->pSFTParams->timestamps->length; /* number of SFTs */
  halfTsft = 0.5*params->pSFTParams->Tsft;          /* half the time of one SFT */

  /* setup baryinput for LALComputeSky */
  baryinput.site = *(params->pSigParams->site);
  /* account for a quirk in LALBarycenter(): -> see documentation of type BarycenterInput */
  baryinput.site.location[0] /= LAL_C_SI;
  baryinput.site.location[1] /= LAL_C_SI;
  baryinput.site.location[2] /= LAL_C_SI;
  if (params->pSigParams->pulsar.position.system != COORDINATESYSTEM_EQUATORIAL) {
      ABORT (status, GENERATEPULSARSIGNALH_EBADCOORDS, GENERATEPULSARSIGNALH_MSGEBADCOORDS);
  }
  TRY( LALNormalizeSkyPosition (status->statusPtr, &tmp, &(params->pSigParams->pulsar.position)), status);
  baryinput.alpha = tmp.longitude;
  baryinput.delta = tmp.latitude;
  baryinput.dInv = 0.e0;      /* following makefakedata_v2 */
  
  if (params->pSigParams->orbit) {
    /* LALComputeSkyBinary parameters */
    csbParams=(CSBParams *)LALMalloc(sizeof(CSBParams));
    csbParams->skyPos=(REAL8 *)LALMalloc(2*sizeof(REAL8));
    if (params->pSigParams->pulsar.spindown) {
       csbParams->spinDwnOrder=params->pSigParams->pulsar.spindown->length;
    } else {
       csbParams->spinDwnOrder=0;
    }
    csbParams->mObsSFT=numSFTs;
    csbParams->tSFT=params->pSFTParams->Tsft;
    csbParams->tGPS=params->pSFTParams->timestamps->data;
    csbParams->skyPos[0]=params->pSigParams->pulsar.position.longitude;
    csbParams->skyPos[1]=params->pSigParams->pulsar.position.latitude;
    csbParams->OrbitalEccentricity = params->pSigParams->orbit->ecc; /* Orbital eccentricy */
    csbParams->ArgPeriapse = params->pSigParams->orbit->argp;       /* argument of periapsis (radians) */
    csbParams->TperiapseSSB = params->pSigParams->orbit->tp; /* time of periapsis passage (in SSB) */
    /* compute semi-major axis and orbital period */
    csbParams->SemiMajorAxis = params->pSigParams->orbit->asini;
    csbParams->OrbitalPeriod = params->pSigParams->orbit->period;
    csbParams->baryinput=&baryinput;
    csbParams->emit = &emit;
    csbParams->earth = &earth;
    csbParams->edat=params->pSigParams->ephemerides;

    /* Call LALComputeSkyBinary */
    TRY ( LALComputeSkyBinary (status->statusPtr, output->skyConst, 0, csbParams), status);
    LALFree(csbParams->skyPos);
    LALFree(csbParams);
  } else {
    /* LALComputeSky parameters */
    csParams=(CSParams *)LALMalloc(sizeof(CSParams));
    csParams->tGPS=params->pSFTParams->timestamps->data;
    csParams->skyPos=(REAL8 *)LALMalloc(2*sizeof(REAL8));
    csParams->mObsSFT=numSFTs;
    csParams->tSFT=params->pSFTParams->Tsft;
    csParams->edat=params->pSigParams->ephemerides;
    csParams->baryinput=&baryinput;
    if (params->pSigParams->pulsar.spindown) {
       csParams->spinDwnOrder=params->pSigParams->pulsar.spindown->length;
    } else {
       csParams->spinDwnOrder=0;
    }
    csParams->skyPos[0]=params->pSigParams->pulsar.position.longitude;
    csParams->skyPos[1]=params->pSigParams->pulsar.position.latitude;
    csParams->earth = &earth;
    csParams->emit = &emit;

    /* Call LALComputeSky */
    TRY ( LALComputeSky (status->statusPtr, output->skyConst, 0, csParams), status);
    LALFree(csParams->skyPos);
    LALFree(csParams);
  }
  
  /* Set up das, the Detector and Source info */
  das = (LALDetAndSource *)LALMalloc(sizeof(LALDetAndSource));
  das->pSource = (LALSource *)LALMalloc(sizeof(LALSource));
  das->pDetector = params->pSigParams->site;
  das->pSource->equatorialCoords.latitude = params->pSigParams->pulsar.position.latitude;
  das->pSource->equatorialCoords.longitude = params->pSigParams->pulsar.position.longitude;
  das->pSource->orientation = 0.0;  /* NOTE THIS FUNCTION COMPUTE F_+ and F_x for ZERO Psi!!! */
  das->pSource->equatorialCoords.system = params->pSigParams->pulsar.position.system;
  timeAndAcc.accuracy=LALLEAPSEC_STRICT;
 
  /* loop that calls LALComputeDetAMResponse to find F_+ and F_x at the midpoint of each SFT for ZERO Psi */
  for(i=0; i<numSFTs; i++) {
      /* Find mid point from timestamp, half way through SFT. */
      TRY ( LALAddFloatToGPS (status->statusPtr, &midTS, &(params->pSFTParams->timestamps->data[i]), halfTsft), status);
      timeAndAcc.gps=midTS;
      TRY ( LALComputeDetAMResponse(status->statusPtr, &response, das, &timeAndAcc), status);
      output->fPlusZeroPsi[i] = response.plus;
      output->fCrossZeroPsi[i] = response.cross;
  }
  LALFree(das->pSource);
  LALFree(das);

  DETATCHSTATUSPTR( status );
  RETURN (status);
} /* LALComputeSkyAndZeroPsiAMResponse */

/* 07/14/04 gam */
/** 
 * Fast generation of Fake SFTs for a pure pulsar signal.
 * Uses the output of LALComputeSkyAndZeroPsiAMResponse and the same inputs
 * as LALGeneratePulsarSignal and LALSignalToSFTs.  The fake signal is
 * Taylor expanded to first order about the midpoint time of each SFT.
 * Analytic expressions are used to find each SFT
 */
void 
LALFastGeneratePulsarSFTs (LALStatus *status,
                           SFTVector **outputSFTs,
                           const SkyConstAndZeroPsiAMResponse *input,
                           const SFTandSignalParams *params)
{
  INT4 numSFTs;                 /* number of SFTs */
  REAL4 N;                      /* N = number of time-samples that would have been used to generate SFTs directly */
  INT4 iSFT;                    /* index that gives which SFT in an SFTVector */
  INT4 SFTlen;                  /* number of frequency bins in an SFT */
  REAL8 tSFT, f0, band, f0Signal, deltaF;
  REAL4 fPlus, fCross, psi, phi0Signal;
  /* REAL4 halfAPlus, halfACross, cosPsi, sinPsi; */ /* 10/12/04 gam */
  REAL4 halfAPlus, halfACross, cos2Psi, sin2Psi;
  REAL8 real