FindChirpIMRSimulation.c

Go to the documentation of this file.

/*
*  Copyright (C) 2007 Duncan Brown, Eirini Messaritaki, Gareth Jones, Jolien Creighton, Patrick Brady, Reinhard Prix, Anand Sengupta, Stephen Fairhurst, Craig Robinson , Thomas Cokelaer
*
*  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., 51 Franklin Street, Fifth Floor, Boston,
*  MA  02110-1301  USA
*/
 
/*-----------------------------------------------------------------------
 *
 * File Name: FindChirpIMRSimulation.c
 *
 * Author: L. M. Goggin
 *
 *-----------------------------------------------------------------------
 */
/** \cond LALINSPIRAL */
 
#include <lal/Units.h>
#include <lal/Date.h>
#include <lal/AVFactories.h>
#include <lal/VectorOps.h>
#include <lal/SeqFactories.h>
#include <lal/DetectorSite.h>
#include <lal/GenerateInspiral.h>
#include <lal/GeneratePPNInspiral.h>
#include <lal/SimulateCoherentGW.h>
#include <Inject.h>
#include <FindChirpIMRSimulation.h>
#include <lal/LIGOLwXML.h>
#include <lal/LIGOMetadataTables.h>
#include <lal/LIGOMetadataInspiralUtils.h>
#include <lal/LALInspiralBank.h>
#include <lal/FindChirp.h>
#include <lal/LALStdlib.h>
#include <lal/LALInspiralBank.h>
#include <lal/GenerateInspiral.h>
#include <lal/GenerateInspRing.h>
#include <lal/TimeSeries.h>
 
 
/**
 * \brief Provides an interface between code build from \ref lalinspiral_findchirp and
 * various simulation packages for injecting chirps into data.
 * \author Brown, D. A. and Creighton, T. D
 *
 * Injects the signals described
 * in the linked list of \c SimInspiralTable structures \c events
 * into the data \c chan. The response function \c resp should
 * contain the response function to use when injecting the signals into the data.
 */
void
LALFindChirpInjectIMR (
    LALStatus                  *status,
    REAL4TimeSeries            *chan,
    SimInspiralTable           *events,
    SimRingdownTable           *ringdownevents,
    COMPLEX8FrequencySeries    *resp,
    INT4                        injectSignalType
    )
 
{
  UINT4                 k;
  DetectorResponse      detector;
  SimInspiralTable     *thisEvent = NULL;
  SimRingdownTable     *thisRingdownEvent = NULL;
  PPNParamStruc         ppnParams;
  CoherentGW            waveform, *wfm;
  INT8                  waveformStartTime;
  REAL4TimeSeries       signalvec;
  COMPLEX8Vector       *unity = NULL;
  CHAR                  warnMsg[512];
#if 0
  UINT4 n;
  UINT4 i;
#endif
 
  INITSTATUS(status);
  ATTATCHSTATUSPTR( status );
 
  ASSERT( chan, status,
      FINDCHIRPH_ENULL, FINDCHIRPH_MSGENULL );
  ASSERT( chan->data, status,
      FINDCHIRPH_ENULL, FINDCHIRPH_MSGENULL );
  ASSERT( chan->data->data, status,
      FINDCHIRPH_ENULL, FINDCHIRPH_MSGENULL );
 
  ASSERT( events, status,
      FINDCHIRPH_ENULL, FINDCHIRPH_MSGENULL );
 
  ASSERT( resp, status,
      FINDCHIRPH_ENULL, FINDCHIRPH_MSGENULL );
  ASSERT( resp->data, status,
      FINDCHIRPH_ENULL, FINDCHIRPH_MSGENULL );
  ASSERT( resp->data->data, status,
      FINDCHIRPH_ENULL, FINDCHIRPH_MSGENULL );
 
 
  /*
   *
   * set up structures and parameters needed
   *
   */
 
  /* fixed waveform injection parameters */
  memset( &ppnParams, 0, sizeof(PPNParamStruc) );
  ppnParams.deltaT   = chan->deltaT;
  ppnParams.lengthIn = 0;
  ppnParams.ppn      = NULL;
 
 
  /*
   *
   * compute the transfer function from the given response function
   *
   */
 
 
  /* allocate memory and copy the parameters describing the freq series */
  memset( &detector, 0, sizeof( DetectorResponse ) );
  detector.transfer = (COMPLEX8FrequencySeries *)
    LALCalloc( 1, sizeof(COMPLEX8FrequencySeries) );
  if ( ! detector.transfer )
  {
    ABORT( status, FINDCHIRPH_EALOC, FINDCHIRPH_MSGEALOC );
  }
  memcpy( &(detector.transfer->epoch), &(resp->epoch),
      sizeof(LIGOTimeGPS) );
  detector.transfer->f0 = resp->f0;
  detector.transfer->deltaF = resp->deltaF;
 
  detector.site = (LALDetector *) LALMalloc( sizeof(LALDetector) );
  /* set the detector site */
  switch ( chan->name[0] )
  {
    case 'H':
      *(detector.site) = lalCachedDetectors[LALDetectorIndexLHODIFF];
      LALWarning( status, "computing waveform for Hanford." );
      break;
    case 'L':
      *(detector.site) = lalCachedDetectors[LALDetectorIndexLLODIFF];
      LALWarning( status, "computing waveform for Livingston." );
      break;
    case 'G':
      *(detector.site) = lalCachedDetectors[LALDetectorIndexGEO600DIFF];
      LALWarning( status, "computing waveform for GEO600." );
      break;
    case 'T':
      *(detector.site) = lalCachedDetectors[LALDetectorIndexTAMA300DIFF];
      LALWarning( status, "computing waveform for TAMA300." );
      break;
    case 'V':
      *(detector.site) = lalCachedDetectors[LALDetectorIndexVIRGODIFF];
      LALWarning( status, "computing waveform for Virgo." );
      break;
    default:
      LALFree( detector.site );
      detector.site = NULL;
      LALWarning( status, "Unknown detector site, computing plus mode "
          "waveform with no time delay" );
      break;
  }
 
  /* set up units for the transfer function */
  if (XLALUnitDivide( &(detector.transfer->sampleUnits),
                      &lalADCCountUnit, &lalStrainUnit ) == NULL) {
    ABORTXLAL(status);
  }
 
  /* invert the response function to get the transfer function */
  LALCCreateVector( status->statusPtr, &( detector.transfer->data ),
      resp->data->length );
  CHECKSTATUSPTR( status );
 
  LALCCreateVector( status->statusPtr, &unity, resp->data->length );
  CHECKSTATUSPTR( status );
  for ( k = 0; k < resp->data->length; ++k )
  {
    unity->data[k] = 1.0;
  }
 
  LALCCVectorDivide( status->statusPtr, detector.transfer->data, unity,
      resp->data );
  CHECKSTATUSPTR( status );
 
  LALCDestroyVector( status->statusPtr, &unity );
  CHECKSTATUSPTR( status );
 
  thisRingdownEvent = ringdownevents;
 
  /*
   *
   * loop over the signals and inject them into the time series
   *
   */
 
 
  for ( thisEvent = events; thisEvent; thisEvent = thisEvent->next)
  {
    /*
     *
     * generate waveform and inject it into the data
     *
     */
 
 
    /* clear the waveform structure */
    memset( &waveform, 0, sizeof(CoherentGW) );
 
    LALGenerateInspiral(status->statusPtr, &waveform, thisEvent, &ppnParams);
    CHECKSTATUSPTR( status );
 
    /* add the ringdown */
    wfm = XLALGenerateInspRing( &waveform, thisEvent, thisRingdownEvent,
       injectSignalType );
 
    if ( !wfm )
    {
      fprintf( stderr, "Failed to generate the waveform \n" );
      if (xlalErrno == XLAL_EFAILED)
      {
        fprintf( stderr, "Too much merger\n");
        XLALDestroyREAL4TimeSeries( chan );
        xlalErrno = XLAL_SUCCESS;
        return;
      }
      else exit ( 1 );
    }
 
 
    waveform = *wfm;
 
    LALInfo( status, ppnParams.termDescription );
 
    if ( thisEvent->geocent_end_time.gpsSeconds )
    {
      /* get the gps start time of the signal to inject */
      waveformStartTime = XLALGPSToINT8NS( &(thisEvent->geocent_end_time) );
      waveformStartTime -= (INT8) ( 1000000000.0 * ppnParams.tc );
    }
    else
    {
      LALInfo( status, "Waveform start time is zero: injecting waveform "
          "into center of data segment" );
 
      /* center the waveform in the data segment */
      waveformStartTime = XLALGPSToINT8NS( &(chan->epoch) );
 
      waveformStartTime += (INT8) ( 1000000000.0 *
          ((REAL8) (chan->data->length - ppnParams.length) / 2.0) * chan->deltaT
          );
    }
 
    snprintf( warnMsg, XLAL_NUM_ELEM(warnMsg),
        "Injected waveform timing:\n"
        "thisEvent->geocent_end_time.gpsSeconds = %d\n"
        "thisEvent->geocent_end_time.gpsNanoSeconds = %d\n"
        "ppnParams.tc = %e\n"
        "waveformStartTime = %" LAL_INT8_FORMAT "\n",
        thisEvent->geocent_end_time.gpsSeconds,
        thisEvent->geocent_end_time.gpsNanoSeconds,
        ppnParams.tc,
        waveformStartTime );
    LALInfo( status, warnMsg );
 
    /* clear the signal structure */
    memset( &signalvec, 0, sizeof(REAL4TimeSeries) );
 
    /* set the start times for injection */
    XLALINT8NSToGPS( &(waveform.a->epoch), waveformStartTime );
    memcpy( &(waveform.f->epoch), &(waveform.a->epoch),
        sizeof(LIGOTimeGPS) );
    memcpy( &(waveform.phi->epoch), &(waveform.a->epoch),
        sizeof(LIGOTimeGPS) );
 
    /* set the start time of the signal vector to the start time of the chan */
    signalvec.epoch = chan->epoch;
 
    /* set the parameters for the signal time series */
    signalvec.deltaT = chan->deltaT;
    if ( ( signalvec.f0 = chan->f0 ) != 0 )
    {
      ABORT( status, FINDCHIRPH_EHETR, FINDCHIRPH_MSGEHETR );
    }
    signalvec.sampleUnits = lalADCCountUnit;
 
    /* simulate the detectors response to the inspiral */
    LALSCreateVector( status->statusPtr, &(signalvec.data), chan->data->length );
    CHECKSTATUSPTR( status );
 
    LALSimulateCoherentGW( status->statusPtr,
        &signalvec, &waveform, &detector );
    CHECKSTATUSPTR( status );
 
/* *****************************************************************************/
#if 0
    FILE *fp;
    char fname[512];
    UINT4 jj, kplus, kcross;
    snprintf( fname, XLAL_NUM_ELEM(fname),
        "waveform-%d-%d-%s.txt",
        thisEvent->geocent_end_time.gpsSeconds,
        thisEvent->geocent_end_time.gpsNanoSeconds,
        thisEvent->waveform );
    fp = fopen( fname, "w" );
    for( jj = 0, kplus = 0, kcross = 1; jj < waveform.phi->data->length;
        ++jj, kplus += 2, kcross +=2 )
    {
      fprintf(fp, "%d %e %e %le %e\n", jj,
          waveform.a->data->data[kplus],
          waveform.a->data->data[kcross],
          waveform.phi->data->data[jj],
          waveform.f->data->data[jj]);
    }
    fclose( fp );
#endif
 
/* ********************************************************************************/
#if 0
    FILE *fp;
    char fname[512];
    UINT4 jj;
    snprintf( fname, XLAL_NUM_ELEM(fname),
        "waveform-%d-%d-%s.txt",
        thisEvent->geocent_end_time.gpsSeconds,
        thisEvent->geocent_end_time.gpsNanoSeconds,
        thisEvent->waveform );
    fp = fopen( fname, "w" );
    for( jj = 0; jj < signalvec.data->length; ++jj )
    {
      fprintf(fp, "%d %e %e \n", jj, signalvec.data->data[jj]);
    }
    fclose( fp );
#endif
 
/* ********************************************************************************/
 
 
    /* inject the signal into the data channel */
    LALSSInjectTimeSeries( status->statusPtr, chan, &signalvec );
    CHECKSTATUSPTR( status );
 
    /* allocate and go to next SimRingdownTable */
    if( thisEvent->next )
      thisRingdownEvent = thisRingdownEvent->next = (SimRingdownTable *)
      calloc( 1, sizeof(SimRingdownTable) );
    else
      thisRingdownEvent->next = NULL;
 
    /* destroy the signal */
    LALSDestroyVector( status->statusPtr, &(signalvec.data) );
    CHECKSTATUSPTR( status );
 
    LALSDestroyVectorSequence( status->statusPtr, &(waveform.a->data) );
    CHECKSTATUSPTR( status );
 
    LALSDestroyVector( status->statusPtr, &(waveform.f->data) );
    CHECKSTATUSPTR( status );
 
    LALDDestroyVector( status->statusPtr, &(waveform.phi->data) );
    CHECKSTATUSPTR( status );
 
    if ( waveform.shift )
    {
      LALSDestroyVector( status->statusPtr, &(waveform.shift->data) );
      CHECKSTATUSPTR( status );
    }
 
    LALFree( waveform.a );
    LALFree( waveform.f );
    LALFree( waveform.phi );
    if ( waveform.shift )
    {
      LALFree( waveform.shift );
    }
 
  }
 
  LALCDestroyVector( status->statusPtr, &( detector.transfer->data ) );
  CHECKSTATUSPTR( status );
 
  if ( detector.site ) LALFree( detector.site );
  LALFree( detector.transfer );
 
  DETATCHSTATUSPTR( status );
  RETURN( status );
}
 
/** \endcond */