C++ DETATCHSTATUSPTR函数代码示例

void
LALStringToZ(LALStatus * stat, COMPLEX16 * value, const CHAR * string,
             CHAR ** endptr)
{
    REAL8 re, im;       /* real and imaginary parts */
    CHAR *end;  /* substring following parsed numbers */

    INITSTATUS(stat);
    ATTATCHSTATUSPTR(stat);

    /* Check for valid input arguments. */
    ASSERT(value, stat, STRINGINPUTH_ENUL, STRINGINPUTH_MSGENUL);
    ASSERT(string, stat, STRINGINPUTH_ENUL, STRINGINPUTH_MSGENUL);

    /* Parse string.  Return if nothing was parsed. */
    TRY(LALStringToD(stat->statusPtr, &re, string, &end), stat);
    TRY(LALStringToD(stat->statusPtr, &im, end, &end), stat);
    if (string == end) {
        if (endptr)
            *endptr = end;
        DETATCHSTATUSPTR(stat);
        RETURN(stat);
    }

    /* Set values and return. */
    *value = re;
    *value += im * I;
    if (endptr)
        *endptr = end;
    DETATCHSTATUSPTR(stat);
    RETURN(stat);
}

void LALForwardRealDFT(
    LALStatus      *status,
    COMPLEX8Vector *output,
    REAL4Vector    *input
    )
{
  COMPLEX8Vector *a = NULL;
  COMPLEX8Vector *b = NULL;
  UINT4 n;
  UINT4 j;
  UINT4 k;

  INITSTATUS(status);
  ATTATCHSTATUSPTR( status );

  n = input->length;

  TRY( LALCCreateVector( status->statusPtr, &a, n ), status );
  TRY( LALCCreateVector( status->statusPtr, &b, n ), status );

  for ( j = 0; j < n; ++j )
    a->data[j] = input->data[j];

  TRY( LALDFT( status->statusPtr, b, a, -1 ), status );

  for ( k = 0; k <= n / 2; ++k )
    output->data[k] = b->data[k];

  TRY( LALCDestroyVector( status->statusPtr, &a ), status );
  TRY( LALCDestroyVector( status->statusPtr, &b ), status );

  DETATCHSTATUSPTR( status );
  RETURN( status );
}

void ComputeFoft(LALStatus   *status,
		 REAL8Vector          *foft,
                 HoughTemplate        *pulsarTemplate,
		 REAL8Vector          *timeDiffV,
		 REAL8Cart3CoorVector *velV,
                 REAL8                 timeBase){
  
  INT4   mObsCoh;
  REAL8   f0new, vcProdn, timeDiffN;
  INT4    f0newBin;
  REAL8   sourceDelta, sourceAlpha, cosDelta;
  INT4    j,i, nspin, factorialN; 
  REAL8Cart3Coor  sourceLocation;
  
  /* --------------------------------------------- */
  INITSTATUS(status);
  ATTATCHSTATUSPTR (status);
  
  /*   Make sure the arguments are not NULL: */
  ASSERT (foft,  status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL);
  ASSERT (pulsarTemplate,  status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL);
  ASSERT (timeDiffV,  status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL);
  ASSERT (velV,  status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL);
  
  ASSERT (foft->data,  status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL);
  ASSERT (timeDiffV->data,  status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL);
  ASSERT (velV->data,  status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL);
  
  sourceDelta = pulsarTemplate->latitude;
  sourceAlpha = pulsarTemplate->longitude;
  cosDelta = cos(sourceDelta);
  
  sourceLocation.x = cosDelta* cos(sourceAlpha);
  sourceLocation.y = cosDelta* sin(sourceAlpha);
  sourceLocation.z = sin(sourceDelta);
    
  mObsCoh = foft->length;    
  nspin = pulsarTemplate->spindown.length;
  
  for (j=0; j<mObsCoh; ++j){  /* loop for all different time stamps */
    vcProdn = velV->data[j].x * sourceLocation.x
      + velV->data[j].y * sourceLocation.y
      + velV->data[j].z * sourceLocation.z;
    f0new = pulsarTemplate->f0;
    factorialN = 1;
    timeDiffN = timeDiffV->data[j];
    
    for (i=0; i<nspin;++i){ /* loop for spin-down values */
      factorialN *=(i+1);
      f0new += pulsarTemplate->spindown.data[i]* timeDiffN / factorialN;
      timeDiffN *= timeDiffN;
    }
    f0newBin = floor( f0new * timeBase + 0.5);
    foft->data[j] = f0newBin * (1.0 +vcProdn) / timeBase;
  }    
    
  DETATCHSTATUSPTR (status);
  /* normal exit */
  RETURN (status);
}

void
LALInspiralEccentricityTemplates(
   LALStatus        *status,
   REAL4Vector      *signalvec1,
   REAL4Vector      *signalvec2,
   InspiralTemplate *params
   )
 {

   INT4 count;

   INITSTATUS(status);
   ATTATCHSTATUSPTR(status);

   ASSERT(signalvec1, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signalvec2, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signalvec1->data, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);
   ASSERT(signalvec2->data, status, LALINSPIRALH_ENULL, LALINSPIRALH_MSGENULL);

   /* Initially the waveforms are empty */
   memset(signalvec1->data, 0, signalvec1->length * sizeof(REAL4));
   memset(signalvec2->data, 0, signalvec2->length * sizeof(REAL4));

   /* Call the engine function */
   LALInspiralEccentricityEngine(status->statusPtr, signalvec1, signalvec2, NULL, NULL, NULL, &count, params);
   CHECKSTATUSPTR(status);

   DETATCHSTATUSPTR(status);
   RETURN (status);

}

/** \see See \ref Eigen_c for documentation */
void
LALDSymmetricEigenValues( LALStatus *stat, REAL8Vector *values, REAL8Array *matrix )
{
  REAL8Vector *offDiag = NULL; /* off-diagonal line of
                                  tri-diagonalized matrix */

  INITSTATUS(stat);
  ATTATCHSTATUSPTR( stat );

  /* Check dimension length.  All other argument testing is done by
     the subroutines. */
  ASSERT( values, stat, MATRIXUTILSH_ENUL, MATRIXUTILSH_MSGENUL );
  ASSERT( values->length, stat, MATRIXUTILSH_ENUL, MATRIXUTILSH_MSGENUL );

  /* Allocate an off-diagonal vector for the tri-diagonal matrix. */
  TRY( LALDCreateVector( stat->statusPtr, &offDiag, values->length ),
       stat );

  /* Call the subroutines. */
  LALDSymmetricToTriDiagonal2( stat->statusPtr, values, matrix,
			       offDiag );
  BEGINFAIL( stat ) {
    TRY( LALDDestroyVector( stat->statusPtr, &offDiag ), stat );
  } ENDFAIL( stat );
  LALDTriDiagonalToDiagonal2( stat->statusPtr, values, matrix,
			      offDiag );
  BEGINFAIL( stat ) {
    TRY( LALDDestroyVector( stat->statusPtr, &offDiag ), stat );
  } ENDFAIL( stat );

  /* Clean up. */
  TRY( LALDDestroyVector( stat->statusPtr, &offDiag ), stat );
  DETATCHSTATUSPTR( stat );
  RETURN( stat );
}

void LALSelectPeakWhiteNoise(LALStatus  *status,
			     UCHARPeakGram        *pg,
			     REAL8                *thr,
			     REAL8Periodogram1    *peri)
{ 
  

  INT4  length;
  INT4  nPeaks;
  
  /* --------------------------------------------- */
  INITSTATUS(status);
  ATTATCHSTATUSPTR (status);
  
  /*   Make sure the arguments are not NULL: */
  ASSERT (peri, status, PEAKSELECTH_ENULL, PEAKSELECTH_MSGENULL);
  ASSERT (pg,   status, PEAKSELECTH_ENULL, PEAKSELECTH_MSGENULL);
  ASSERT (thr,  status, PEAKSELECTH_ENULL, PEAKSELECTH_MSGENULL);

  pg->epoch.gpsSeconds     = peri->epoch.gpsSeconds;
  pg->epoch.gpsNanoSeconds = peri->epoch.gpsNanoSeconds;
  pg->timeBase     = peri->timeBase;
  pg->fminBinIndex = peri->fminBinIndex;
  
  length = peri->length;
  nPeaks = 0;
  
  ASSERT (length==pg->length,  status, PEAKSELECTH_EVAL, PEAKSELECTH_MSGEVAL);
  
  if (length > 0){
    UCHAR  *out;
    REAL8  *in1;
    REAL8  threshold;
    INT4      n;

    ASSERT (peri->data, status, PEAKSELECTH_ENULL, PEAKSELECTH_MSGENULL);
    ASSERT (pg->data,   status, PEAKSELECTH_ENULL, PEAKSELECTH_MSGENULL);
    out = pg->data;
    in1 = peri->data;
    threshold = *thr;
    n = length;
    
    while (n-- >0){
      if ( *in1 > threshold){
        *out = 1 ;
	nPeaks++ ;
      } else {
        *out = 0;
      }     
      ++out;
      ++in1;    
    }         
  }
  
  pg->nPeaks = nPeaks;
  
  DETATCHSTATUSPTR (status);
  /* normal exit */
  RETURN (status);
}

/**
 * \ingroup LALInspiralBank_h
 * \brief Module to check whether a point with coordinates (x,y) is inside
 * a polygon defined by the vectors (vx, vy), which size (n) must be
 * provided. The functions returns 1 if the point is inside or 0 otherwise.
 *
 * \author Cokelaer. T
 *
 * ### Notes ###
 *
 * Tested in matlab codes and some BCV tests within lal/lalapps.
 *
 */
void LALInsidePolygon(  LALStatus          *status,	/**< LAL status pointer */
                        REAL4              *inputx,	/**< [in] two arrays of floats defining the polygon */
                        REAL4              *inputy,	/**< [in] two arrays of floats defining the polygon */
                        INT4               n,		/**< [in] the size of the vectors */
                        REAL4              x,		/**< [in] the coordinate of the point */
                        REAL4              y,		/**< [in] the coordinate of the point */
                        INT4               *valid	/**< [out] 0 if outside and 1 if inside */
                        )
{


   INITSTATUS(status);
   ATTATCHSTATUSPTR(status);
   ASSERT (n>=3,  status, LALINSPIRALBANKH_ENULL, LALINSPIRALBANKH_MSGENULL);

   {
     int i, j, c = 0;
     for (i = 0, j = n-1; i < n; j = i++) {
       if ((((inputy[i] <= y) && (y < inputy[j])) ||
	    ((inputy[j] <= y) && (y < inputy[i]))) &&
	   (x < (inputx[j] - inputx[i]) * (y - inputy[i]) / (inputy[j] - inputy[i]) + inputx[i]))
	 c = !c;
     }
     *valid = c;
   }

   DETATCHSTATUSPTR(status);
   RETURN(status);

}

/**
 * register all our "user-variables"
 */
void
InitUserVars (LALStatus *status, struct CommandLineArgsTag *CLA)
{

  INITSTATUS(status);
  ATTATCHSTATUSPTR (status);

  /* Initialize default values */
  CLA->Tsft=1800;
  CLA->nTsft=0;            
  CLA->timestamps=NULL;
  CLA->gpsStart=-1;
  CLA->sqrtSh=1.0;
  
  /** Default year-span of ephemeris-files to be used */
  CLA->ephemEarth = XLALStringDuplicate("earth00-19-DE405.dat.gz");
  CLA->ephemSun = XLALStringDuplicate("sun00-19-DE405.dat.gz");
  
  /* ---------- register all our user-variable ---------- */
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->Alpha),      "Alpha",          REAL8,  'a', OPTIONAL,  "Sky position Alpha (equatorial coordinates) in radians") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->Alpha),      "longitude",      REAL8,  0,   DEVELOPER, "[DEPRECATED] Use --Alpha instead!") == XLAL_SUCCESS, XLAL_EFUNC);

  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->Delta),      "Delta",          REAL8,  'd', OPTIONAL,  "Sky position Delta (equatorial coordinates) in radians") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->Delta),      "latitude",       REAL8,  0,   DEVELOPER, "[DEPRECATED] Use --Delta instead!") == XLAL_SUCCESS, XLAL_EFUNC);
 
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->phi0),       "phi0",           REAL8,  'Q', OPTIONAL,  "Phi_0: Initial phase in radians") == XLAL_SUCCESS, XLAL_EFUNC);

  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->psi),        "psi",            REAL8,  'Y', OPTIONAL,  "Polarisation in radians") == XLAL_SUCCESS, XLAL_EFUNC);

  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->cosi),       "cosi",           REAL8,  'i', OPTIONAL,  "Cos(iota)") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->cosiota),    "cosiota",        REAL8,  0,   DEVELOPER, "[DEPRECATED] Use --cosi instead") == XLAL_SUCCESS, XLAL_EFUNC);

  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->h0),         "h0",             REAL8,  's', OPTIONAL,  "Strain amplitude h_0") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->sqrtSh),     "sqrtSh",         REAL8,  'N', OPTIONAL,  "Noise floor: one-sided sqrt(Sh) in 1/sqrt(Hz)") == XLAL_SUCCESS, XLAL_EFUNC);
  
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->timestamps), "timestampsFile", STRING, 'T', OPTIONAL,  "Name of timestamps file") == XLAL_SUCCESS, XLAL_EFUNC);
  
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->gpsStart),   "startTime",      INT4,   'S', OPTIONAL,  "GPS start time of continuous observation") == XLAL_SUCCESS, XLAL_EFUNC);
  
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->Tsft),       "Tsft",           REAL8,  't', OPTIONAL,  "Length of an SFT in seconds") == XLAL_SUCCESS, XLAL_EFUNC);
  
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->nTsft),      "nTsft",          INT4,   'n', OPTIONAL,  "Number of SFTs") == XLAL_SUCCESS, XLAL_EFUNC);
  
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->IFO),        "IFO",            STRING, 'D', OPTIONAL,  "Detector: H1, H2, L1, G1, ... ") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->detector),   "detector",       STRING, 0,   DEVELOPER, "[DEPRECATED] Use --IFO instead!") == XLAL_SUCCESS, XLAL_EFUNC);
  
  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->ephemEarth), "ephemEarth",     STRING, 0,   OPTIONAL,  "Earth ephemeris file to use") == XLAL_SUCCESS, XLAL_EFUNC);

  XLAL_CHECK_LAL( status, XLALRegisterNamedUvar(&(CLA->ephemSun),   "ephemSun",       STRING, 0,   OPTIONAL,  "Sun ephemeris file to use") == XLAL_SUCCESS, XLAL_EFUNC);

  /* ----- added for mfd_v4 compatibility ---------- */
  XLAL_CHECK_LAL ( status, XLALRegisterNamedUvar(&(CLA->duration),  "duration",       REAL8, 0,    OPTIONAL,  "Duration of requested signal in seconds") == XLAL_SUCCESS, XLAL_EFUNC);
  
  XLAL_CHECK_LAL ( status, XLALRegisterNamedUvar(&(CLA->aPlus),     "aPlus",          REAL8, 0,    OPTIONAL,  "Plus polarization amplitude aPlus") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_LAL ( status, XLALRegisterNamedUvar(&(CLA->aCross),    "aCross",         REAL8, 0,    OPTIONAL,  "Cross polarization amplitude aCross") == XLAL_SUCCESS, XLAL_EFUNC);


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

/**
 * \deprecated Use XLALComputeDetAMResponse() instead.
 */
void LALComputeDetAMResponse(LALStatus * status, LALDetAMResponse * pResponse, const LALDetAndSource * pDetAndSrc, const LIGOTimeGPS * gps)
{
	double fplus, fcross;

	INITSTATUS(status);
	ATTATCHSTATUSPTR(status);

	ASSERT(pResponse != (LALDetAMResponse *) NULL, status, DETRESPONSEH_ENULLOUTPUT, DETRESPONSEH_MSGENULLOUTPUT);

	ASSERT(pDetAndSrc != NULL, status, DETRESPONSEH_ENULLINPUT, DETRESPONSEH_MSGENULLINPUT);

	ASSERT(gps != (LIGOTimeGPS *) NULL, status, DETRESPONSEH_ENULLINPUT, DETRESPONSEH_MSGENULLINPUT);

	/* source coordinates must be in equatorial system */
	ASSERT(pDetAndSrc->pSource->equatorialCoords.system == COORDINATESYSTEM_EQUATORIAL, status, DETRESPONSEH_ESRCNOTEQUATORIAL, DETRESPONSEH_MSGESRCNOTEQUATORIAL);

	XLALComputeDetAMResponse(&fplus, &fcross, pDetAndSrc->pDetector->response, pDetAndSrc->pSource->equatorialCoords.longitude, pDetAndSrc->pSource->equatorialCoords.latitude, pDetAndSrc->pSource->orientation, XLALGreenwichMeanSiderealTime(gps));

	pResponse->plus = fplus;
	pResponse->cross = fcross;
	pResponse->scalar = 0.0;	/* not implemented */

	DETATCHSTATUSPTR(status);
	RETURN(status);
}

/**
 * Calculate the binary system time delay using the pulsar parameters in
 * \c params
 */
void
LALBinaryPulsarDeltaT( LALStatus            *status,
                       BinaryPulsarOutput   *output,
                       BinaryPulsarInput    *input,
                       BinaryPulsarParams   *params ){
  INITSTATUS(status);
  ATTATCHSTATUSPTR(status);

  /* Check input arguments */
  ASSERT(input != (BinaryPulsarInput *)NULL, status,
  BINARYPULSARTIMINGH_ENULLINPUT, BINARYPULSARTIMINGH_MSGENULLINPUT);

  ASSERT(output != (BinaryPulsarOutput *)NULL, status,
  BINARYPULSARTIMINGH_ENULLOUTPUT, BINARYPULSARTIMINGH_MSGENULLOUTPUT);

  ASSERT(params != (BinaryPulsarParams *)NULL, status,
  BINARYPULSARTIMINGH_ENULLPARAMS, BINARYPULSARTIMINGH_MSGENULLPARAMS);

  ASSERT((!strcmp(params->model, "BT")) ||
         (!strcmp(params->model, "BT1P")) ||
         (!strcmp(params->model, "BT2P")) ||
         (!strcmp(params->model, "BTX")) ||
         (!strcmp(params->model, "ELL1")) ||
         (!strcmp(params->model, "DD")) ||
         (!strcmp(params->model, "DDS")) ||
         (!strcmp(params->model, "MSS")) ||
         (!strcmp(params->model, "T2")), status,
         BINARYPULSARTIMINGH_ENULLBINARYMODEL,
         BINARYPULSARTIMINGH_MSGNULLBINARYMODEL);

  XLALBinaryPulsarDeltaT( output, input, params );

  DETATCHSTATUSPTR(status);
  RETURN(status);
}

/**
 * Function to clean a sft vector -- calls LALCleanCOMPLEX8SFT repeatedly for each
 * sft in vector
 */
void LALCleanMultiSFTVect (LALStatus       *status,   /**< pointer to LALStatus structure */
			   MultiSFTVector  *multVect, /**< SFTVector to be cleaned */
			   INT4            width,     /**< maximum width to be cleaned -- set sufficiently large if all bins in each line are to be cleaned*/
			   INT4            window,    /**< window size for noise floor estimation in vicinity of a line */
			   LineNoiseInfo   *lineInfo, /**< list of lines */
			   RandomParams    *randPar   /**< parameters for generating random noise */)
{

  UINT4 k;

  INITSTATUS(status);
  ATTATCHSTATUSPTR (status);

  ASSERT (multVect, status, SFTCLEANH_ENULL, SFTCLEANH_MSGENULL);
  ASSERT (multVect->data, status, SFTCLEANH_ENULL, SFTCLEANH_MSGENULL);
  ASSERT (multVect->length > 0, status, SFTCLEANH_EVAL, SFTCLEANH_MSGEVAL);
  ASSERT (lineInfo, status, SFTCLEANH_ENULL, SFTCLEANH_MSGENULL);
  ASSERT (lineInfo->nLines, status, SFTCLEANH_EVAL, SFTCLEANH_MSGEVAL);
  ASSERT (lineInfo->lineFreq, status, SFTCLEANH_ENULL, SFTCLEANH_MSGENULL);
  ASSERT (lineInfo->leftWing, status, SFTCLEANH_ENULL, SFTCLEANH_MSGENULL);
  ASSERT (lineInfo->rightWing, status, SFTCLEANH_ENULL, SFTCLEANH_MSGENULL);
  ASSERT (window > 0, status, SFTCLEANH_EVAL, SFTCLEANH_MSGEVAL);
  ASSERT (width >= 0, status, SFTCLEANH_EVAL, SFTCLEANH_MSGEVAL);


  for ( k = 0; k < multVect->length; k++) {
    TRY (LALCleanSFTVector (status->statusPtr, multVect->data[k], width, window, lineInfo, randPar), status);
  }

  DETATCHSTATUSPTR (status);
  /* normal exit */
  RETURN (status);
}

void
LALFreqTimeRealFFT(
    LALStatus               *status,
    REAL4TimeSeries         *time,
    COMPLEX8FrequencySeries *freq,
    RealFFTPlan             *plan
    )
{
  INITSTATUS(status);
  XLAL_PRINT_DEPRECATION_WARNING("XLALREAL4FreqTimeFFT");
  ATTATCHSTATUSPTR( status );

  ASSERT( plan, status, TIMEFREQFFTH_ENULL, TIMEFREQFFTH_MSGENULL );
  ASSERT( freq, status, TIMEFREQFFTH_ENULL, TIMEFREQFFTH_MSGENULL );
  ASSERT( time, status, TIMEFREQFFTH_ENULL, TIMEFREQFFTH_MSGENULL );
  ASSERT( time->data, status, TIMEFREQFFTH_ENULL, TIMEFREQFFTH_MSGENULL );
  ASSERT( time->data->length, status,
      TIMEFREQFFTH_ESIZE, TIMEFREQFFTH_MSGESIZE );
  ASSERT( freq->deltaF > 0, status, TIMEFREQFFTH_ERATE, TIMEFREQFFTH_MSGERATE );

  /* call the XLAL function */
  if ( XLALREAL4FreqTimeFFT( time, freq, plan ) == XLAL_FAILURE )
  {
    XLALClearErrno();
    ABORTXLAL( status );
  }

  DETATCHSTATUSPTR( status );
  RETURN( status );
}

/** Main driver funtion for doing Numerical Relativity Injections */
void LALDriveNRInject( LALStatus *status,	/**< pointer to LALStatus structure */
		       REAL4TimeSeries *injData, /**< The time series to inject into */
		       SimInspiralTable *injections, /**< The list of injections to perform */
		       NumRelInjectParams *params /**< Parameters */
		       )
{

  REAL4TimeVectorSeries *sumStrain = NULL;
  SimInspiralTable *thisInj    = NULL;   /* current injection */

  INITSTATUS(status);
  ATTATCHSTATUSPTR (status);

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

      TRY( LALAddStrainModes(status->statusPtr, &sumStrain, params->nrCatalog,
				  params->modeLlo, params->modeLhi, thisInj), status);

      TRY( LALInjectStrainGW( status->statusPtr, injData, sumStrain, thisInj, params->ifo, params->dynRange), status);

      XLALDestroyREAL4VectorSequence ( sumStrain->data );
      LALFree( sumStrain );
      sumStrain = NULL;

    } /* end loop over injections */


  DETATCHSTATUSPTR(status);
  RETURN(status);

}

static void dAstroOmega (LALStatus *s, REAL8 *domegaz, REAL8 z, void *p)
 {
  AstroOmegaParams params;
  AstroOmegaSourceParams sourcep;
  /*AstroOmegaCosmoParams cosmop;*/
  REAL8LALSDensity *SDensitySource;
  REAL8 Rc, dEgw, nu, nuz;

  INITSTATUS(s);
  ATTATCHSTATUSPTR (s);
  params = *((AstroOmegaParams *)p);
  sourcep = params.sourceparams;
  SDensitySource = sourcep.SDensitySource;
  nu = *((REAL8 *)params.extraparams);

  /*frequency in the source frame*/
  nuz = (1. + z) * nu;
  /*single spectral energy density in the source frame*/
  SDensitySource(&dEgw, nuz);
  /*cosmic formation rate*/
  SFR(&Rc, z);
  *domegaz = dEgw * Rc / pow((1.+z),3.5);
  CHECKSTATUSPTR (s);
  DETATCHSTATUSPTR (s);
  RETURN (s);
 }

void LALAstroOmega (LALStatus *s, REAL8 *omeganu, REAL8 nu, void *p)
 {
   DIntegrateIn  zint;
   AstroOmegaParams params;
   AstroOmegaSourceParams sourcep;
   AstroOmegaCosmoParams cosmop;
   REAL8 omegaz, zmax, numax, lambda;
   INITSTATUS(s);
   ATTATCHSTATUSPTR (s);

   params = *((AstroOmegaParams *)p);
   cosmop=params.cosmoparams;
   sourcep = params.sourceparams;
   numax= sourcep.numax;
   lambda = sourcep.lambda;

   if((nu >= numax)||(nu <= 0.)){*omeganu = 0.;}
   else
    {
     if (nu < (numax / 6.)) {zmax = 5.;}
     else {zmax = (numax / nu) - 1.;}

     zint.function = dAstroOmega;
     zint.xmin     = 0;
     zint.xmax     = zmax;
     zint.type     = ClosedInterval;

     LALDRombergIntegrate (s->statusPtr, &omegaz, &zint, &params);
     *omeganu = 4.66e-56 * lambda / (cosmop.ho * cosmop.ho) * nu * omegaz;
     }
  CHECKSTATUSPTR (s);
  DETATCHSTATUSPTR (s);
  RETURN (s);
 }