static REAL8VectorList *
XLALREAL8VectorListAddEntry (REAL8VectorList *head, const REAL8Vector *entry)
{
  UINT4 dim;
  REAL8VectorList *ptr = NULL;  /* running list-pointer */
  REAL8VectorList *newElement = NULL;   /* new list-element */
  /* check illegal input */
  if ( (head == NULL) || (entry == NULL) )
    return NULL;

  /* find tail of list */
  ptr = head;
  while ( ptr->next )
    ptr = ptr->next;

  /* construct new list-element */
  dim = entry->length;
  if ( (newElement = LALCalloc (1, sizeof (*newElement))) == NULL)
    return NULL;
  if ( (newElement->entry.data = LALCalloc (dim, sizeof(entry->data[0]))) == NULL ) {
    LALFree (newElement);
    return NULL;
  }
  newElement->entry.length = dim;
  memcpy (newElement->entry.data, entry->data, dim * sizeof(entry->data[0]) );

  /* link this to the tail of list */
  ptr->next = newElement;
  newElement->prev = ptr;

  return newElement;

} /* XLALREAL8VectorListAddEntry() */

static int frEvent2snglInspiral(SnglInspiralTable **snglInspiralEvent, 
    FrEvent *frameEvent )
{
  FrEvent              *frEvt   = NULL;
  SnglInspiralTable    *snglEvt = NULL;
  int                   numEvt  = 0;
  double                timeAfter = 0;

  /* If we already have events in snglInspiralEvent, 
   * wind on to the end of the list */
  for( snglEvt = *snglInspiralEvent; snglEvt; snglEvt=snglEvt->next);

  /* store the frameEvents in the snglInspiral linked list */
  for(frEvt=frameEvent; frEvt; frEvt=frEvt->next, ++numEvt) 
  {
    if ( !(*snglInspiralEvent) )
    {
      *snglInspiralEvent = snglEvt = (SnglInspiralTable * ) 
        LALCalloc( 1, sizeof(SnglInspiralTable) );
    }
    else
    {
      snglEvt = snglEvt->next = (SnglInspiralTable * ) 
        LALCalloc( 1, sizeof(SnglInspiralTable) );
    }

    /* read data from the frEvt */
    snprintf(snglEvt->search, LIGOMETA_SEARCH_MAX, "%s", frEvt->name);
    snglEvt->snr = frEvt->amplitude;
    snglEvt->end_time.gpsSeconds = frEvt->GTimeS;
    snglEvt->end_time.gpsNanoSeconds = frEvt->GTimeN;
    timeAfter = frEvt->timeAfter;
    XLALGPSAdd(&snglEvt->end_time,timeAfter);
    snglEvt->eff_distance = FrEventGetParam ( frEvt, ignore_const("distance (Mpc)") );
    snglEvt->mass1 = FrEventGetParam ( frEvt, ignore_const("mass1") );
    snglEvt->mass2 = FrEventGetParam ( frEvt, ignore_const("mass2") );
    snglEvt->tau0 =FrEventGetParam ( frEvt, ignore_const("tau0") );
    snglEvt->tau3 = FrEventGetParam ( frEvt, ignore_const("tau1p5") );
    snglEvt->coa_phase = FrEventGetParam ( frEvt, ignore_const("phase") );
    snglEvt->chisq = FrEventGetParam ( frEvt, ignore_const("chi2") );

    /* populate additional colums */
    snglEvt->mtotal = snglEvt->mass1 + snglEvt->mass2;
    snglEvt->eta = (snglEvt->mass1 * snglEvt->mass2) /
      (snglEvt->mtotal * snglEvt->mtotal);
    snglEvt->mchirp = pow( snglEvt->eta, 0.6) * snglEvt->mtotal;
    snprintf(snglEvt->ifo, LIGOMETA_IFO_MAX, "%s", ifo);
  }
  return( numEvt );
}

LALFrameUFrDetector *XLALFrameUFrDetectorAlloc_FrameL_(const char *name,
    const char *prefix, double latitude, double longitude, double elevation,
    double azimuthX, double azimuthY, double altitudeX, double altitudeY, double midpointX, double midpointY, int localTime)
{
    LALFrameUFrDetector *detector;
    detector = LALCalloc(1, sizeof(*detector));
    if (!detector)
        XLAL_ERROR_NULL(XLAL_ENOMEM);
    detector->handle = calloc(1, sizeof(*detector->handle));
    if (!detector->handle)
        XLAL_ERROR_NULL(XLAL_ENOMEM);
    detector->handle->classe = FrDetectorDef();
    detector->handle->name = strdup(name);
    if (!detector->handle->name) {
        XLALFrameUFrDetectorFree(detector);
        XLAL_ERROR_NULL(XLAL_ENOMEM);
    }
    if (prefix) {
        memcpy(detector->prefix, prefix, 2);
        memcpy(detector->handle->prefix, prefix, 2);
    }

    detector->handle->longitude = longitude;    /* longitude (east of greenwich) in radians      */
    detector->handle->latitude = latitude;      /* latitude (north of equator) in radians        */
    detector->handle->elevation = elevation;    /* detector altitude (meter)                     */
    detector->handle->armXazimuth = azimuthX;   /* orientation of X arm in radians CW from North */
    detector->handle->armYazimuth = azimuthY;   /* orientation of Y arm in radians CW from North */
    /* Azimuth values should be in the range 0 to 2pi */
    detector->handle->armXaltitude = altitudeX; /* altitude (pitch) of the X arm                 */
    detector->handle->armYaltitude = altitudeY; /* altitude (pitch) of the Y arm                 */
    detector->handle->armXmidpoint = midpointX; /* vertex to  middle of the X arm distance       */
    detector->handle->armYmidpoint = midpointY; /* vertex to  middle of the Y arm distance       */
    detector->handle->localTime = localTime;    /* local time - UTC time (second)                */
    return detector;
}

/* Does no error checking, so check immediately after calling. */
static void
allocate(
    REAL4               x,
    REAL4               y,
    REAL4               z,
    REAL4               f0,
    SnglInspiralTable **tmplt,
    INT4               *ntiles,
    BOOLEAN             havePsi )
{
  REAL4 mass, eta, m1, m2;

  *tmplt = (*tmplt)->next = (SnglInspiralTable *) LALCalloc( 1,
           sizeof(SnglInspiralTable) );

  mass = -y/x / (16.0*LAL_PI*LAL_PI*f0);
  eta = 16.0457 * pow( -x*x/y/y/y/y/y, 0.3333333 );
  m1 = 0.5*mass* (1 + sqrt(1 - 4*eta));
  m2 = 0.5*mass* (1 - sqrt(1 - 4*eta));

  if ( ! havePsi )
  {
    (*tmplt)->mass1 = m1;
    (*tmplt)->mass2 = m2;
    (*tmplt)->eta = eta;
    (*tmplt)->mchirp = pow(m1*m2,0.6)/pow(m1+m2,0.2);
  }
  (*tmplt)->psi0 = x*pow(f0,5./3);
  (*tmplt)->psi3 = y*pow(f0,2./3);
  (*tmplt)->beta = z*pow(f0,2./3);
  ++(*ntiles);
} /* allocate() */

/**
 * Function for interpolating PSD to a given sample rate
 */
REAL8FrequencySeries *
XLALInterpolatePSD( REAL8FrequencySeries *in,      /**< input strain time series */
                    REAL8                deltaFout /**< sample rate of time series */)
{
  REAL8FrequencySeries *ret=NULL;
  REAL8 deltaFin, r, y_1, y_2;
  UINT4 k, lo, numPoints;

  deltaFin = in->deltaF;

  /* length of output vector */
  numPoints = (UINT4) (in->data->length * deltaFin / deltaFout);

  /* allocate memory */
  ret = LALCalloc(1, sizeof(*ret));
  if (!ret)
  {
    XLAL_ERROR_NULL( XLAL_ENOMEM );
  }

  ret->data = XLALCreateREAL8Vector( numPoints );
  if (! ret->data)
  {
    XLAL_ERROR_NULL( XLAL_ENOMEM );
  }

  ret->deltaF = deltaFout;

  /* copy values from in which should be the same */
  ret->epoch       = in->epoch;
  ret->f0          = in->f0;
  ret->sampleUnits = in->sampleUnits;
  strcpy(ret->name, in->name);

  /* go over points of output vector and interpolate linearly
     using closest points of input */
  for (k = 0; k < numPoints; k++) {
    lo = (UINT4)( k*deltaFout / deltaFin);

    /* y_1 and y_2 are the input values at x1 and x2 */
    /* here we need to make sure that we don't exceed
       bounds of input vector */
    if ( lo < in->data->length - 1) {
      y_1 = in->data->data[lo];
      y_2 = in->data->data[lo+1];

      /* we want to calculate y_2*r + y_1*(1-r) where
         r = (x-x1)/(x2-x1) */
      r = k*deltaFout / deltaFin - lo;

      ret->data->data[k] = y_2 * r + y_1 * (1 - r);
    }
    else {
      ret->data->data[k] = 0.0;
    }
  }

  return ret;
}

REAL4Vector *
XLALComputeFreq(
		REAL4TimeSeries *hp,
		REAL4TimeSeries *hc)
{
  REAL4Vector *Freq = NULL;
  REAL4Vector *hpDot = NULL, *hcDot = NULL;
  UINT4 k, len;
  REAL8 dt;

  len = hp->data->length;
  dt = hp->deltaT;
  Freq = LALCalloc(1, sizeof(*Freq));
  Freq= XLALCreateREAL4Vector(len);

  hpDot = LALCalloc(1, sizeof(*hpDot));
  hpDot= XLALCreateREAL4Vector(len);
  hcDot = LALCalloc(1, sizeof(*hcDot));
  hcDot= XLALCreateREAL4Vector(len);

  /* Construct the dot vectors (2nd order differencing) */
  hpDot->data[0] = 0.0;
  hpDot->data[len] = 0.0;
  hcDot->data[0] = 0.0;
  hcDot->data[len] = 0.0;
  for( k = 1; k < len-1; k++)
    {
      hpDot->data[k] = 1./(2.*dt) *
	(hp->data->data[k+1]-hp->data->data[k-1]);
      hcDot->data[k] = 1./(2.*dt) *
	(hc->data->data[k+1]-hc->data->data[k-1]);
    }

  /* Compute frequency using the fact that  */
  /*h(t) = A(t) e^(i Phi) = Re(h) + i Im(h) */
  for( k = 0; k < len; k++)
    {
      Freq->data[k] = hcDot->data[k] * hp->data->data[k] -
	hpDot->data[k] * hc->data->data[k];
      Freq->data[k] /= LAL_TWOPI;
      Freq->data[k] /= (pow(hp->data->data[k],2.) + pow(hc->data->data[k], 2.));
    }

  return Freq;

}

int main(int argc, char *argv[])
{
	char tstr[32]; // string to hold GPS time -- 31 characters is enough
	const double H0 = 0.72 * LAL_H0FAC_SI; // Hubble's constant in seconds
	const size_t length = 65536; // number of points in a segment
	const size_t stride = length / 2; // number of points in a stride
	size_t i, n;
	REAL8FrequencySeries *OmegaGW = NULL;
	REAL8TimeSeries **seg = NULL;
	LIGOTimeGPS epoch;
	gsl_rng *rng;

	XLALSetErrorHandler(XLALAbortErrorHandler);

	parseargs(argc, argv);

	XLALGPSSetREAL8(&epoch, tstart);
	gsl_rng_env_setup();
	rng = gsl_rng_alloc(gsl_rng_default);
	OmegaGW = XLALSimSGWBOmegaGWFlatSpectrum(Omega0, flow, srate/length, length/2 + 1);

	n = duration * srate;
	seg = LALCalloc(numDetectors, sizeof(*seg));
	printf("# time (s)");
	for (i = 0; i < numDetectors; ++i) {
		char name[LALNameLength];
		snprintf(name, sizeof(name), "%s:STRAIN", detectors[i].frDetector.prefix);
		seg[i] = XLALCreateREAL8TimeSeries(name, &epoch, 0.0, 1.0/srate, &lalStrainUnit, length);
		printf("\t%s (strain)", name);
	}
	printf("\n");

	XLALSimSGWB(seg, detectors, numDetectors, 0, OmegaGW, H0, rng); // first time to initilize

	while (1) { // infinite loop
		size_t j;
		for (j = 0; j < stride; ++j, --n) { // output first stride points
			LIGOTimeGPS t = seg[0]->epoch;
			if (n == 0) // check if we're done
				goto end;
			printf("%s", XLALGPSToStr(tstr, XLALGPSAdd(&t, j * seg[0]->deltaT)));
			for (i = 0; i < numDetectors; ++i)
				printf("\t%e", seg[i]->data->data[j]);
			printf("\n");
		}
		XLALSimSGWB(seg, detectors, numDetectors, stride, OmegaGW, H0, rng); // make more data
	}

end:
	for (i = 0; i < numDetectors; ++i)
		XLALDestroyREAL8TimeSeries(seg[i]);
	XLALFree(seg);
	XLALDestroyREAL8FrequencySeries(OmegaGW);
	LALCheckMemoryLeaks();

	return 0;
}

static int frSimEvent2simInspiral (SimInspiralTable **simInspiralEvent,
    FrSimEvent       *frSimEvent )
{
  FrSimEvent           *frSimEvt   = NULL;
  SimInspiralTable     *simEvt     = NULL;
  int                   numSim     = 0;

  /* If we already have events in snglInspiralEvent, 
   * wind on to the end of the list */
  for( simEvt = *simInspiralEvent; simEvt; simEvt = simEvt->next);

  /* store the frameEvents in the snglInspiral linked list */
  for( frSimEvt = frSimEvent; frSimEvt; frSimEvt = frSimEvt->next, ++numSim) 
  {
    if ( !(*simInspiralEvent) )
    {
      *simInspiralEvent = simEvt = (SimInspiralTable * ) 
        LALCalloc( 1, sizeof(SimInspiralTable) );
    }
    else
    {
      simEvt = simEvt->next = (SimInspiralTable * ) 
        LALCalloc( 1, sizeof(SimInspiralTable) );
    }

    /* read data from the frSimEvt */
    snprintf(simEvt->waveform, LIGOMETA_SEARCH_MAX, "%s", frSimEvt->name);
    simEvt->geocent_end_time.gpsSeconds = frSimEvt->GTimeS;
    simEvt->geocent_end_time.gpsNanoSeconds = frSimEvt->GTimeN;
    simEvt->v_end_time = simEvt->geocent_end_time;


    simEvt->distance = FrSimEventGetParam ( frSimEvt, ignore_const("distance") );
    simEvt->eff_dist_v = simEvt->distance;

    simEvt->mass1 = FrSimEventGetParam ( frSimEvt, ignore_const("m1") );
    simEvt->mass2 = FrSimEventGetParam ( frSimEvt, ignore_const("m2") );
  }

  return ( numSim );
}

/** Create a DetectorStateSeries with length entries */
DetectorStateSeries*
XLALCreateDetectorStateSeries ( UINT4 length )		/**< number of entries */
{
  DetectorStateSeries *ret = NULL;

  if ( (ret = LALCalloc(1, sizeof(DetectorStateSeries) )) == NULL ) {
    XLALPrintError ("%s: failed to LALCalloc(1, %zu)\n", __func__, sizeof(DetectorStateSeries) );
    XLAL_ERROR_NULL ( XLAL_ENOMEM );
  }

  if ( (ret->data = LALCalloc (length, sizeof(DetectorState) )) == NULL ) {
    XLALFree (ret);
    XLALPrintError ("%s: failed to LALCalloc(%d, %zu)\n", __func__, length, sizeof(DetectorState) );
    XLAL_ERROR_NULL ( XLAL_ENOMEM );
  }

  ret->length = length;

  return ret;

} /* XLALCreateDetectorStateSeries() */

int main(int argc, char *argv[])
{
	const double H0 = 0.72 * LAL_H0FAC_SI; // Hubble's constant in seconds
	const double srate = 16384.0; // sampling rate in Hertz
	const size_t length = 65536; // number of points in a segment
	const size_t stride = length / 2; // number of points in a stride
	size_t i, n;
	REAL8FrequencySeries *OmegaGW = NULL;
	REAL8TimeSeries **seg = NULL;
	LIGOTimeGPS epoch;
	gsl_rng *rng;

	XLALSetErrorHandler(XLALAbortErrorHandler);

	parseargs(argc, argv);

	XLALGPSSetREAL8(&epoch, tstart);
	gsl_rng_env_setup();
	rng = gsl_rng_alloc(gsl_rng_default);
	OmegaGW = XLALSimSGWBOmegaGWFlatSpectrum(Omega0, flow, srate/length, length/2 + 1);

	n = duration * srate;
	seg = LALCalloc(numDetectors, sizeof(*seg));
	for (i = 0; i < numDetectors; ++i)
		seg[i] = XLALCreateREAL8TimeSeries("STRAIN", &epoch, 0.0, 1.0/srate, &lalStrainUnit, length);

	XLALSimSGWB(seg, detectors, numDetectors, 0, OmegaGW, H0, rng); // first time to initilize
	while (1) { // infinite loop
		double t0 = XLALGPSGetREAL8(&seg[0]->epoch);
		size_t j;
		for (j = 0; j < stride; ++j, --n) { // output first stride points
			if (n == 0) // check if we're done
				goto end;
			printf("%.9f", t0 + j * seg[0]->deltaT);
			for (i = 0; i < numDetectors; ++i)
				printf("\t%e", seg[i]->data->data[j]);
			printf("\n");
		}
		XLALSimSGWB(seg, detectors, numDetectors, stride, OmegaGW, H0, rng); // make more data
	}

end:
	for (i = 0; i < numDetectors; ++i)
		XLALDestroyREAL8TimeSeries(seg[i]);
	XLALFree(seg);
	XLALDestroyREAL8FrequencySeries(OmegaGW);
	LALCheckMemoryLeaks();

	return 0;
}

void
XLALComputeComplexVector(
			 COMPLEX8Vector **outPlus,
			 COMPLEX8Vector **outCross,
			 REAL4FrequencySeries *Ampl,
			 REAL4FrequencySeries *Phase)
{
  COMPLEX8Vector *uPlus = NULL, *uCross = NULL;
  COMPLEX8 num;
  REAL4 Re, Im;
  UINT4 k, n;

  n = Ampl->data->length;

  uPlus = LALCalloc(1, sizeof(*uPlus));
  uPlus = XLALCreateCOMPLEX8Vector(n);
  uCross = LALCalloc(1, sizeof(*uCross));
  uCross = XLALCreateCOMPLEX8Vector(n);

  for( k = 0 ; k < n ; k++ ) {
    Re = Ampl->data->data[k] * cos(Phase->data->data[k]);
    Im = Ampl->data->data[k] * sin(Phase->data->data[k]);
    num.re = Re;
    num.im = Im;
    uPlus->data[k] = num;

    Re = Ampl->data->data[k] * cos(Phase->data->data[k] - LAL_PI/2.);
    Im = Ampl->data->data[k] * sin(Phase->data->data[k] - LAL_PI/2.);
    num.re = Re;
    num.im = Im;
    uCross->data[k] = num;
  }

  *outPlus = uPlus;
  *outCross = uCross;

}

/**
 * @brief Opens a LALFrStream associated with a LALCache
 * @details
 * This routine creates a #LALFrStream that is a stream associated with
 * the frame files contained in a LALCache.
 * @param cache Pointer to a LALCache structure describing the frame files to stream.
 * @returns Pointer to a newly created #LALFrStream structure.
 * @retval NULL Failure.
 */
LALFrStream *XLALFrStreamCacheOpen(LALCache * cache)
{
    LALFrStream *stream;
    size_t i;

    if (!cache)
        XLAL_ERROR_NULL(XLAL_EFAULT);

    stream = LALCalloc(1, sizeof(*stream));
    if (!stream)
        XLAL_ERROR_NULL(XLAL_ENOMEM);
    stream->cache = XLALCacheDuplicate(cache);

    /* check cache entries for t0 and dt; if these are not set then read
     * the framefile to try to get them */
    for (i = 0; i < stream->cache->length; ++i) {
        if (stream->cache->list[i].t0 == 0 || stream->cache->list[i].dt == 0) {
            LIGOTimeGPS end;
            size_t nFrame;
            if (XLALFrStreamFileOpen(stream, i) < 0) {
                XLALFrStreamClose(stream);
                XLAL_ERROR_NULL(XLAL_EIO);
            }
            nFrame = XLALFrFileQueryNFrame(stream->file);
            stream->cache->list[i].t0 = stream->epoch.gpsSeconds;
            XLALFrFileQueryGTime(&end, stream->file, nFrame - 1);
            XLALGPSAdd(&end, XLALFrFileQueryDt(stream->file, nFrame - 1));
            stream->cache->list[i].dt =
                ceil(XLALGPSGetREAL8(&end)) - stream->cache->list[i].t0;
            XLALFrStreamFileClose(stream);
        }
    }

    /* sort and uniqify the cache */
    if (XLALCacheSort(stream->cache) || XLALCacheUniq(stream->cache)) {
        XLALFrStreamClose(stream);
        XLAL_ERROR_NULL(XLAL_EFUNC);
    }

    stream->mode = LAL_FR_STREAM_DEFAULT_MODE;

    /* open up the first file */
    if (XLALFrStreamFileOpen(stream, 0) < 0) {
        XLALFrStreamClose(stream);
        XLAL_ERROR_NULL(XLAL_EFUNC);
    }
    return stream;
}

REAL8 start_freq_from_frame_url(CHAR  *url)
{

  FrameH *frame=NULL;
  FrFile *frFile=NULL;
  FrHistory *frHist=NULL;
  FrHistory *thisHist;
  CHAR *comment=NULL;
  CHAR *token=NULL;
  REAL8 ret=0;
  CHAR *path;

  /* convert url to path by skipping protocol part of protocol:path */
  path = strchr(url, ':');
  if (path == NULL)
    path = url;
  else
    path+=strlen("://localhost"); /* skip the '://localhost' -- now on the path */

  frFile =  FrFileINew( path );
  frame = FrameRead (frFile);
  frHist = frame->history;
  thisHist = frHist;
  while (thisHist) {

    /* get history comment string and parse it */
    comment = LALCalloc(1, (strlen(thisHist->comment)+1)*sizeof(CHAR));
    strcpy(comment, thisHist->comment);

    token = strtok(comment,":");

    if (strstr(token,"freqStart22") || strstr(token,"freq_start_22")) {
      token = strtok(NULL,":");
      ret = atof(token);
    }

    LALFree(comment);
    comment = NULL;

    thisHist = thisHist->next;
  }

  FrFileIEnd( frFile );
  return ret;
}

void
XLALSortSimRingdown(
    SimRingdownTable **head,
    int (*comparefunc)(const SimRingdownTable * const *,
      const SimRingdownTable * const *)
    )

{
  INT4 i;
  INT4 length;
  SimRingdownTable *event;
  SimRingdownTable **array;

  /* empty list --> no-op */
  if(!head || !*head)
    return;

  /* count the number of events in the list */
  for(length = 0, event = *head; event; event = event->next)
    length++;

  /* construct an array of pointers into the list */
  array = LALCalloc(length, sizeof(*array));
  for(i = 0, event = *head; event; event = event->next)
    array[i++] = event;

  /* sort the array using the specified function */
  qsort(array, length, sizeof(*array),
      (int(*)(const void *, const void *)) comparefunc);

  /* re-link the list according to the sorted array */
  for(i = 0; i < length; i++, head = &(*head)->next)
    *head = array[i];
  *head = NULL;

  /* free the array */
  LALFree(array);
}

/**
 * Function for reading a numerical relativity metadata file.
 * It returns a list of numrel wave parameters.  It uses
 * XLALParseDataFile() for reading the data file.  This automatically
 * takes care of removing comment lines starting with # and other details.
 */
void
LALNRDataFind( LALStatus *status,   /**< pointer to LALStatus structure */
	       NRWaveCatalog *out,  /**< [out] list of numrel metadata */
	       const CHAR *dir,     /**< [in] directory with data files */
	       const CHAR *filename /**< [in] File with metadata information */)
{
  LALParsedDataFile *cfgdata=NULL;
  UINT4 k, numWaves;

  INITSTATUS(status);
  ATTATCHSTATUSPTR (status);

  ASSERT (filename != NULL, status, NRWAVEIO_ENULL, NRWAVEIO_MSGENULL );
  ASSERT ( out != NULL, status, NRWAVEIO_ENULL, NRWAVEIO_MSGENULL );
  ASSERT ( dir != NULL, status, NRWAVEIO_ENULL, NRWAVEIO_MSGENULL );

  if ( XLALParseDataFile ( &cfgdata, filename ) != XLAL_SUCCESS ) {
    ABORT( status, NRWAVEIO_EFILE, NRWAVEIO_MSGEFILE );
  }
  numWaves = cfgdata->lines->nTokens; /*number of waves */

  /* allocate memory for output catalog */
  out->length = numWaves;
  out->data = LALCalloc(1, out->length * sizeof(NRWaveMetaData));
  if ( out->data == NULL) {
    ABORT( status, NRWAVEIO_ENOMEM, NRWAVEIO_MSGENOMEM );
  }

  /* now get wave parameters from each line of data */
  for (k = 0; k < numWaves; k++) {
    TRY(LALGetSingleNRMetaData( status->statusPtr, out->data + k, dir, cfgdata->lines->tokens[k]), status);
  }

  XLALDestroyParsedDataFile (cfgdata);

  DETATCHSTATUSPTR(status);
  RETURN(status);
}

int
XLALCoincSimRingdownTest (
    SimRingdownTable   **simHead,
    CoincRingdownTable **coincHead,
    SimRingdownTable   **missedSimHead,
    CoincRingdownTable **missedCoincHead
    )

{
  CoincRingdownTable    *thisCoinc       = *coincHead;
  CoincRingdownTable    *prevCoinc       = NULL;
  CoincRingdownTable    *thisMissedCoinc = NULL;
  SimRingdownTable      *thisSim         = NULL;
  SimRingdownTable      *prevSim         = NULL;
  SimRingdownTable      *thisMissedSim   = NULL;
  SnglRingdownTable     *thisSngl        = NULL;
  EventIDColumn         *thisId          = NULL;

  InterferometerNumber   ifoInCoinc = LAL_UNKNOWN_IFO;
  int                    numSimFound = 0;

  if ( !*coincHead )
  {
    XLALPrintInfo(
        "XLALCoincSimRingdown: Empty coincInspiral passed as input" );
    *missedSimHead = *simHead;
    *simHead = NULL;
    return( 0 );
  }

  *coincHead = NULL;

   while( thisCoinc )
   {
     thisSim = NULL;
     /* loop over the interferometers to get the event_id*/

     for ( ifoInCoinc = (InterferometerNumber) 0; ifoInCoinc < LAL_NUM_IFO; ifoInCoinc++)
     {
       if ( (thisSngl = thisCoinc->snglRingdown[ifoInCoinc]) )
       {
         thisSim = thisSngl->event_id->simRingdownTable;
         break;
       }
     }

     for ( ; ifoInCoinc < LAL_NUM_IFO; ifoInCoinc++)
     {
       if ( (thisSngl = thisCoinc->snglRingdown[ifoInCoinc]) &&
           (thisSim != thisSngl->event_id->simRingdownTable) )
       {
         thisSim = NULL;
         break;
       }
     }

     if ( thisSim )

     {
       /* thisCoinc is coincident with a thisSim */
       thisCoinc->simRingdown = thisSim;

       /* set the event_id's */
       if ( !thisSim->event_id )
       {
         thisId = thisSim->event_id = LALCalloc( 1, sizeof(EventIDColumn) );
       }
       else
       {
         for ( thisId = thisSim->event_id; thisId->next; thisId = thisId->next);
         thisId = thisId->next = LALCalloc( 1, sizeof(EventIDColumn) );
       }
       thisId->simRingdownTable = thisSim;
       thisId->coincRingdownTable = thisCoinc;

       if ( ! *coincHead )
       {
         *coincHead = thisCoinc;
       }

       XLALPrintInfo( "+" );
       /* move on to
        * * next coinc */
       prevCoinc = thisCoinc;
       thisCoinc = thisCoinc->next;
     }
     else
     {
       /* discard this event and move on to the next one */
       if ( ! *missedCoincHead )
       {
         *missedCoincHead = thisMissedCoinc = thisCoinc;
       }
       else
       {
         thisMissedCoinc = thisMissedCoinc->next = thisCoinc;
       }

       if ( prevCoinc ) prevCoinc->next = thisCoinc->next;
       thisCoinc = thisCoinc->next;
//.........这里部分代码省略.........

void
XLALMCMCMetro (
  LALMCMCParameter  **paraPtr,
  LALMCMCInput       *inputMCMC
) {

  static LALStatus status;

  LALMCMCParameter *parameter   = NULL;
  LALMCMCParam     *paraHead    = NULL;
  gsl_matrix       *startCovMat = NULL;
  gsl_matrix       *covMat      = NULL;
  REAL4Vector      *vector      = NULL;
  StateMCMC state;
  REAL4     logPrior, logLikelihood, currentLogPosterior;
  REAL8     range=0.0;
  int       foundBurnin, loopMCMC, move, dim, i,j;

  /* initialize the parameter structure */
  parameter = *paraPtr;
  inputMCMC->funcInit( parameter, inputMCMC->inspiralTable );

  dim            = parameter->dimension;
  inputMCMC->dim = dim;

  /* initialize the internal vectors */
  for (paraHead=parameter->param; paraHead;
       paraHead=paraHead->next)
  {
    paraHead->core->chain=NULL;
    LALSCreateVector( &status, &paraHead->core->chain, inputMCMC->numberDraw+1);
  }

  /* initialize vectors in the MCMC structure */
  LALSCreateVector( &status, &vector, dim);
  inputMCMC->mean = (REAL8*)LALCalloc( sizeof(REAL8), dim );
  inputMCMC->xdiff= (REAL8*)LALCalloc( sizeof(REAL8), dim );
  inputMCMC->ndiff= (REAL8*)LALCalloc( sizeof(REAL8), dim );

  /* prepare the gsl matrices */
  startCovMat = gsl_matrix_alloc( dim, dim);
  covMat      = gsl_matrix_alloc( dim, dim);

  /* set everything to zeros in the matrix */
  for (i=0;i<dim;i++)
    for (j=0;j<dim;j++)
      gsl_matrix_set( startCovMat, i, j, 0.0 );

  /* populate the startCovMat with initial values */
  printf("populating the covariance matrix:\n");
  for (paraHead=parameter->param,i=0; paraHead; paraHead=paraHead->next,i++)
  {
    range=paraHead->core->maxVal - paraHead->core->minVal;
    gsl_matrix_set( startCovMat, i, i, 0.001*range*range );
    printf("element %d: %f\n", i, 0.001*range*range);
  }
  gsl_matrix_memcpy( covMat, startCovMat );

  /* initialize the state */
  inputMCMC->counter=inputMCMC->counterState=inputMCMC->counterAccept=inputMCMC->counterAcceptDraw=0;
  inputMCMC->burninCounter=0;
  loopMCMC=1;
  if (inputMCMC->useScaling)
  {
    state = doScaling;
    inputMCMC->scaling = inputMCMC->scalePeak;
  }
  else if (inputMCMC->useAnnealing)
  {
    state = doAnnealing;
    inputMCMC->scaling = inputMCMC->scaleNormal;
  }
  else
  {
    state = unknownState; /* TODO: error*/
    printf("ERROR: unknown state\n");
    exit(1);
  }


  /* get the first values */
  inputMCMC->funcPrior( inputMCMC, parameter );
  logLikelihood = inputMCMC->funcLikelihood( inputMCMC, parameter);
  logPrior = parameter->logPrior;
  currentLogPosterior=logLikelihood+logPrior;

  do
  {
    /* increase the counters */
    inputMCMC->counter++;      /*k*/
    inputMCMC->counterState++; /*c*/

    /* print generel information */
    printf("####################################################\n");
    printf("\nMCMCSTATUS k: %d c: %d  state: ",
           inputMCMC->counter, inputMCMC->counterState);
    printState( state);
    printf("\n");


//.........这里部分代码省略.........

/** The main function of binary2sft.c
 *
 */
int main( int argc, char *argv[] )  {

  UserInput_t uvar = empty_UserInput;           /* user input variables */
  INT4 i,j;                                     /* counter */
  SFTVector *SFTvect = NULL;
  char *noisestr = XLALCalloc(1,sizeof(char));

  /**********************************************************************************/
  /* register and read all user-variables */
  if (XLALReadUserVars(argc,argv,&uvar)) {
    LogPrintf(LOG_CRITICAL,"%s : XLALReadUserVars() failed with error = %d\n",__func__,xlalErrno);
    return 1;
  }
  LogPrintf(LOG_DEBUG,"%s : read in uservars\n",__func__);

  /**********************************************************************************/
  /* read in the cache file */
  FILE *cachefp = NULL;
  if ((cachefp = fopen(uvar.cachefile,"r")) == NULL) {
    LogPrintf(LOG_CRITICAL,"%s : failed to open binary input file %s\n",__func__,uvar.cachefile);
    return 1;
  }
  i = 0;
  while (fscanf(cachefp,"%*s %*d %*d")!=EOF) i++;
  INT4 Nfiles = i;
  fclose(cachefp);
  LogPrintf(LOG_DEBUG,"%s : counted %d files listed in the cache file.\n",__func__,Nfiles);

  /* allocate memory */
  char **filenames = LALCalloc(Nfiles,sizeof(char*));
  LIGOTimeGPSVector fileStart;
  fileStart.data = LALCalloc(Nfiles,sizeof(LIGOTimeGPS));
  for (i=0;i<Nfiles;i++) filenames[i] = LALCalloc(512,sizeof(char));

  if ((cachefp = fopen(uvar.cachefile,"r")) == NULL) {
    LogPrintf(LOG_CRITICAL,"%s : failed to open binary input file %s\n",__func__,uvar.cachefile);
    return 1;
  }

  for (i=0;i<Nfiles;i++) {
    fscanf(cachefp,"%s %d %d %*d",filenames[i],&(fileStart.data[i].gpsSeconds),&(fileStart.data[i].gpsNanoSeconds));
  }
  fclose(cachefp);

  /* initialise the random number generator */
  gsl_rng * r;
  if (XLALInitgslrand(&r,uvar.seed)) {
    LogPrintf(LOG_CRITICAL,"%s: XLALinitgslrand() failed with error = %d\n",__func__,xlalErrno);
    XLAL_ERROR(XLAL_EFAULT);
  }

  /* setup the binaryToSFT parameters */
  BinaryToSFTparams par;
  par.tsft = uvar.tsft;
  par.freq = uvar.freq;
  par.freqband = uvar.freqband;
  par.tsamp = uvar.tsamp;
  par.highpassf = uvar.highpassf;
  par.amp_inj = uvar.amp_inj;
  par.f_inj = uvar.f_inj;
  par.asini_inj = uvar.asini_inj;
  XLALGPSSetREAL8(&(par.tasc_inj),uvar.tasc_inj);
  par.tref = fileStart.data[0];
  par.P_inj = uvar.P_inj;
  par.phi_inj = uvar.phi_inj;
  par.r = r;

  /**********************************************************************************/
  /* loop over the input files */
  long int ntot = 0;
  for (j=0;j<Nfiles;j++) {

    UINT4 k = 0;
    INT8Vector *np = NULL;
    REAL8Vector *R = NULL;
    par.tstart = fileStart.data[j];
    REAL8 norm1 = par.tsamp/par.tsft;
    REAL8 norm2 = 1.0/(par.tsamp*par.tsft);
    UINT4 oldlen;
    if (SFTvect==NULL) oldlen = 0;
    else oldlen = SFTvect->length;
    LogPrintf(LOG_DEBUG,"%s : working on file %s\n",__func__,filenames[j]);

    if (XLALBinaryToSFTVector(&SFTvect,filenames[j],&(fileStart.data[j]),&par,&np,&R)) {
      LogPrintf(LOG_CRITICAL,"%s : failed to convert binary input file %s to sfts\n",__func__,filenames[j]);
      return 1;
    }
    if ((np!=NULL) && (R!=NULL)) {
      for (k=0;k<np->length;k++) {
        ntot += np->data[k];
        char temp[64];
        sprintf(temp,"%d %e %e\n",SFTvect->data[oldlen+k].epoch.gpsSeconds,(REAL8)np->data[k]*norm1,R->data[k]*norm2);
        noisestr = (char *)XLALRealloc(noisestr,sizeof(char)*(1+strlen(noisestr)+strlen(temp)));
        strcat(noisestr,temp);
      }
      XLALDestroyINT8Vector(np);
      XLALDestroyREAL8Vector(R);
//.........这里部分代码省略.........

/**
 * \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];
//.........这里部分代码省略.........

/**
 * Main function for injecting numetrical relativity waveforms.
 * Takes as input a list of injections, and adds h(t) to a given
 * timeseries for a specified ifo and a dynamic range factor.
 */
void InjectNumRelWaveformsREAL8 (LALStatus      *status,       /**< pointer to LALStatus structure */
                            REAL8TimeSeries     *chan,         /**< [out] the output time series */
                            SimInspiralTable    *injections,   /**< [in] list of injections */
                            CHAR                ifo[3],        /**< [in] 2 char code for interferometer */
                            REAL8               freqLowCutoff, /**< [in] Lower cutoff frequency */
                            REAL8               snrLow,        /**< [in] lower cutoff value of snr */
                            REAL8               snrHigh,       /**< TO BE DOCUMENTED */
                            CHAR                *fname)       /**< [in] higher cutoff value of snr */
{
  SimInspiralTable *thisInj = NULL;
  REAL8 startFreq, startFreqHz, massTotal;
  REAL8 thisSNR;
  SimInspiralTable *simTableOut=NULL;
  SimInspiralTable *thisInjOut=NULL;

  INITSTATUS(status);
  ATTATCHSTATUSPTR (status);
  ASSERT( chan, status, INSPIRALH_ENULL, INSPIRALH_MSGENULL );
  ASSERT( ifo, status, INSPIRALH_ENULL, INSPIRALH_MSGENULL );


  /* loop over injections */
  for ( thisInj = injections; thisInj; thisInj = thisInj->next )
    {

      startFreq = start_freq_from_frame_url(thisInj->numrel_data);
      massTotal = (thisInj->mass1 + thisInj->mass2) * LAL_MTSUN_SI;
      startFreqHz = startFreq / ( LAL_TWOPI * massTotal);

      if (startFreqHz < freqLowCutoff)
        {
          REAL8TimeSeries *strain = NULL;
          strain  = XLALNRInjectionStrain(ifo, thisInj);
          thisSNR = calculate_ligo_snr_from_strain_real8(strain, ifo);

           /* set channel name */
           snprintf( chan->name, LALNameLength * sizeof( CHAR ),
                    "%s:STRAIN", ifo );

          if ((thisSNR < snrHigh) && (thisSNR > snrLow))
            {
              /* simTableOut will be null only the