LAL: StochasticMC.c Source File

00001 /*
00002 *  Copyright (C) 2007 Jolien Creighton, Kaice T. Reilly, Robert Adam Mercer, Tania Regimbau
00003 *
00004 *  This program is free software; you can redistribute it and/or modify
00005 *  it under the terms of the GNU General Public License as published by
00006 *  the Free Software Foundation; either version 2 of the License, or
00007 *  (at your option) any later version.
00008 *
00009 *  This program is distributed in the hope that it will be useful,
00010 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 *  GNU General Public License for more details.
00013 *
00014 *  You should have received a copy of the GNU General Public License
00015 *  along with with program; see the file COPYING. If not, write to the
00016 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
00017 *  MA  02111-1307  USA
00018 */
00019 
00020 
00021 /*
00022 <lalVerbatim file="StochasticMCCV">
00023 Author: Tania Regimbau, Sukanta Bose, Jeff Noel
00024 $Id: StochasticMC.c,v 3.24 2007/06/08 14:41:47 bema Exp $
00025 </lalVerbatim> 
00026 <lalLaTeX>
00027 \subsection{Module \texttt{StochasticMC.c}}
00028 \label{ss:StochasticMC.c}
00029 
00030 Routine used by the stochastic DSO to do software injection.
00031 
00032 \subsubsection*{Prototypes}
00033 \input{StochasticMCCP}
00034 \idx{LALStochasticMCDso()}
00035 \subsubsection*{Description}
00036 This routine simulates time-domain signal in a pair 
00037 of detectors using Sukanta Bose's code SimulateSB.c, whitened with the adequate response function that can be used in LALwrapper.  
00038 \idx{LALStochasticMCDsoSplice()}
00039 \subsubsection*{Description}
00040 In this version, long time-series are constructed by concatenating short segments of simulated data whitened with the adequate response function. Segments are sinusoidally spliced into consecutive segments using Jeff Noel's function SinusoidalSplice, in order to avoid discontinuities in the final time serie.  
00041 \subsubsection*{Algorithm}
00042 
00043 The following program shows how to use the routines LALStochasticMCDso and LALStochasticMCDsoSplice
00044 
00045 \begin verbatim
00046 #include <math.h>
00047 #include <string.h>
00048 #include <stdio.h>
00049 #ifdef HAVE_UNISTD_H
00050 #include <unistd.h>
00051 #endif
00052 #ifdef HAVE_GETOPT_H
00053 #include <getopt.h>
00054 #endif
00055 #include <FrameL.h>
00056 #include <lal/LALStdio.h>
00057 #include <lal/LALStdlib.h>
00058 #include <lal/AVFactories.h>
00059 #include <lal/PrintFTSeries.h>
00060 #include <lal/FrameStream.h>
00061 #include <lal/FrameCalibration.h>
00062 #include <lal/Calibration.h>
00063 #include <lal/LALConstants.h>
00064 #include <lal/LALStatusMacros.h>
00065 #include <lal/StochasticCrossCorrelation.h>
00066 #include <lal/RealFFT.h>
00067 #include <lal/ComplexFFT.h>
00068 #include <lal/Units.h>
00069 #include <lal/StreamInput.h>
00070 #include <lal/PrintVector.h>
00071 #include <lal/Random.h>
00072 #include <lal/SimulateSB.h>
00073 #include <lal/StochasticMC.h>
00074 
00075 NRCSID(STOCHASTICTESTMCC, "$Id: StochasticMC.c,v 3.24 2007/06/08 14:41:47 bema Exp $");
00076 int main( ){
00077   
00078   static LALStatus status;
00079   SSSimStochBGOutput MCoutput;
00080   StochasticMCParams MCparams;
00081   StochasticMCSInput MCinput;
00082 
00083   //output structure
00084   REAL4TimeSeries SimStochBG1;
00085   REAL4TimeSeries SimStochBG2;
00086 
00087   //input structure
00088   CalibrationFunctions calfuncs1,calfuncs2;
00089   CalibrationUpdateParams  calfacts1,calfacts2;
00090   COMPLEX8FrequencySeries responseFunction1,responseFunction2,sensingFunction1,sensingFunction2; 
00091   LIGOTimeGPS epoch;
00092   COMPLEX8TimeSeries oloopfactor1,sfactor1,oloopfactor2,sfactor2;
00093   COMPLEX8Vector *response[2]={NULL,NULL};
00094   COMPLEX8Vector *sensing[2]={NULL,NULL};
00095   COMPLEX8Vector *oloopfactor[2]={NULL,NULL};
00096   COMPLEX8Vector *sfactor[2]={NULL,NULL};
00097   
00098   //parameters  
00099   UINT4 timeref, starttime, caltime;
00100   UINT4 freqlen,length,lengthseg,caloffset;
00101   REAL8 sRate,deltaT, deltaF;
00102 
00103   INT4 i, j;
00104   LALUnit countPerStrain = {0,{0,0,0,0,0,-1,1},{0,0,0,0,0,0,0}};
00105   
00106   CHAR Outfile1[LALNameLength];
00107   CHAR Outfile2[LALNameLength]; 
00108   FILE *pfone, *pftwo;
00109   
00110   status.statusPtr = NULL;
00111   lengthseg = 61440;
00112   sRate = 1024.;
00113   starttime = 729331170;
00114   
00115   //write parameters 
00116   MCparams.lengthseg = lengthseg; MCparams.numseg = 3;
00117   MCparams.sRate = sRate;
00118   MCparams.starttime = starttime;
00119   MCparams.seed = 173;
00120   MCparams.fRef = 100.; MCparams.f0 = 0.;
00121   MCparams.omegaRef = 1.; MCparams.alpha = 0;
00122   MCparams.site1 = 0; MCparams.site2 = 0; 
00123 
00124   //set other parameters
00125   length = MCparams.numseg * lengthseg;
00126   deltaT = 1./sRate;
00127   freqlen = lengthseg / 2 + 1;
00128   caloffset =  lengthseg / (2 * sRate);
00129   deltaF = 1.0/(deltaT*lengthseg);
00130   timeref = 0; 
00131   caltime = starttime + caloffset;
00132   for (i=0;i<2;i++)
00133   {
00134     LALCCreateVector(&status,&response[i],freqlen);
00135     LALCCreateVector(&status,&sensing[i],freqlen);
00136   }
00137  
00138  for (i=0;i<2;i++)
00139   {
00140     LALCCreateVector(&status,&oloopfactor[i],length);
00141     LALCCreateVector(&status,&sfactor[i],length);
00142   }
00143 
00144  for (i=0;i<2;i++)
00145   {
00146     for (j=0; j < freqlen; j++) 
00147      {
00148       response[i]->data[j].re = 1.;
00149       response[i]->data[j].im = 1.;
00150       sensing[i]->data[j].re = 1.;
00151       sensing[i]->data[j].im = 1.;
00152      }
00153   }
00154 
00155   for (i=0;i<2;i++)
00156   {
00157     for (j=0; j < length; j++) 
00158      {
00159       oloopfactor[i]->data[j].re = 1.;
00160       oloopfactor[i]->data[j].im = 1.;
00161       sfactor[i]->data[j].re = 1.;
00162       sfactor[i]->data[j].im = 1.;
00163      }
00164   }
00165     
00166   memset(&responseFunction1,0,sizeof(COMPLEX8FrequencySeries));
00167   memset(&responseFunction2,0,sizeof(COMPLEX8FrequencySeries));
00168   memset(&sensingFunction1,0,sizeof(COMPLEX8FrequencySeries));
00169   memset(&sensingFunction2,0,sizeof(COMPLEX8FrequencySeries));
00170   
00171   responseFunction1.epoch.gpsSeconds = timeref;
00172   responseFunction1.epoch.gpsNanoSeconds = 0;
00173   responseFunction1.deltaF = deltaF;
00174   responseFunction1.f0 = 0;
00175   responseFunction1.sampleUnits = countPerStrain;
00176   responseFunction1.data = response[0];
00177 
00178   responseFunction2.epoch.gpsSeconds = timeref;
00179   responseFunction2.epoch.gpsNanoSeconds = 0;
00180   responseFunction2.deltaF = deltaF;
00181   responseFunction2.f0 = 0;
00182   responseFunction2.sampleUnits = countPerStrain;
00183   responseFunction2.data = response[1];
00184   
00185   sensingFunction1.epoch.gpsSeconds = timeref;
00186   sensingFunction1.epoch.gpsNanoSeconds = 0;
00187   sensingFunction1.deltaF = deltaF;
00188   sensingFunction1.f0 = 0;
00189   sensingFunction1.sampleUnits = countPerStrain;
00190   sensingFunction1.data = sensing[0];
00191 
00192   sensingFunction2.epoch.gpsSeconds = timeref;
00193   sensingFunction2.epoch.gpsNanoSeconds = 0;
00194   sensingFunction2.deltaF = deltaF;
00195   sensingFunction2.f0 = 0;
00196   sensingFunction2.sampleUnits = countPerStrain;
00197   sensingFunction2.data = sensing[1];
00198    
00199   memset(&oloopfactor1,0,sizeof(COMPLEX8TimeSeries));
00200   memset(&oloopfactor2,0,sizeof(COMPLEX8TimeSeries));
00201   memset(&sfactor1,0,sizeof(COMPLEX8TimeSeries));
00202   memset(&sfactor2,0,sizeof(COMPLEX8TimeSeries));
00203 
00204   oloopfactor1.epoch.gpsSeconds = timeref;
00205   oloopfactor1.epoch.gpsNanoSeconds = 0;
00206   oloopfactor1.deltaT = 1. / sRate;
00207   oloopfactor1.f0 = 0;
00208   oloopfactor1.sampleUnits = lalADCCountUnit;
00209   oloopfactor1.data = oloopfactor[0];
00210 
00211   oloopfactor2.epoch.gpsSeconds = timeref;
00212   oloopfactor2.epoch.gpsNanoSeconds = 0;
00213   oloopfactor2.deltaT = 1. / sRate;
00214   oloopfactor2.f0 = 0;
00215   oloopfactor2.sampleUnits = lalADCCountUnit;
00216   oloopfactor2.data = oloopfactor[1];
00217 
00218   sfactor1.epoch.gpsSeconds = timeref;
00219   sfactor1.epoch.gpsNanoSeconds = 0;
00220   sfactor1.deltaT = 1. / sRate;
00221   sfactor1.f0 = 0;  
00222   sfactor1.sampleUnits = lalADCCountUnit;
00223   sfactor1.data = sfactor[0];
00224 
00225   sfactor2.epoch.gpsSeconds = timeref;
00226   sfactor2.epoch.gpsNanoSeconds = 0;
00227   sfactor2.deltaT = 1. / sRate;
00228   sfactor2.f0 = 0;
00229   sfactor2.sampleUnits = lalADCCountUnit;
00230   sfactor2.data = sfactor[1];
00231   
00232   calfuncs1.responseFunction = &responseFunction1; 
00233   calfuncs1.sensingFunction = &sensingFunction1; 
00234   calfuncs2.responseFunction = &responseFunction2;
00235   calfuncs2.sensingFunction = &sensingFunction2;
00236   
00237   calfacts1.openLoopFactor = &oloopfactor1; 
00238   calfacts1.sensingFactor = &sfactor1;
00239   calfacts2.openLoopFactor = &oloopfactor2;
00240   calfacts2.sensingFactor = &sfactor2;
00241   
00242   epoch.gpsSeconds = caltime;
00243   epoch.gpsNanoSeconds = 0;
00244   calfacts1.epoch = epoch;calfacts2.epoch = epoch;
00245 
00246   //write input parameters    
00247   MCinput.calfuncs1 = calfuncs1; 
00248   MCinput.calfuncs2 = calfuncs2;
00249   MCinput.calfacts1 = calfacts1; 
00250   MCinput.calfacts2 = calfacts2;
00251    
00252   SimStochBG1.data = NULL;
00253   LALSCreateVector(&status, &(SimStochBG1.data), length);
00254   SimStochBG2.data = NULL;
00255   LALSCreateVector(&status, &(SimStochBG2.data), length);
00256   
00257   MCoutput.SSimStochBG1 = &SimStochBG1;
00258   MCoutput.SSimStochBG2 = &SimStochBG2;
00259 
00260   LALStochasticMCDso (&status,&MCoutput,&MCinput,&MCparams);
00261 
00262   //create output files 
00263 
00264   if (print){
00265   LALSnprintf( Outfile1, LALNameLength * sizeof(CHAR),
00266             "%s_1_%d.ilwd",_IFO1,starttime);
00267   LALSnprintf( Outfile2, LALNameLength * sizeof(CHAR),
00268             "%s_2_%d.ilwd",_IFO2,starttime);
00269 
00270   //print output to ilwds
00271      
00272    if (print)
00273     {
00274      pfone=LALFopen(Outfile1,"w");pftwo=LALFopen(Outfile2,"w");
00275      fprintf(pfone,"<?ilwd?>\n");fprintf(pftwo,"<?ilwd?>\n");
00276      fprintf(pfone,"<ilwd comment='%s'",Outfile1);
00277      fprintf(pfone," name='%s' size='1'>\n",Outfile1);
00278      fprintf(pfone," <real_8 dims='%d' name='%s'>",length,Outfile1);
00279      fprintf(pftwo,"<ilwd comment='%s'",Outfile2);
00280      fprintf(pftwo," name='%s' size='1'>\n",Outfile2);
00281      fprintf(pftwo," <real_8 dims='%d' name='%s'>",length,Outfile2);
00282   
00283      for(i=0;(UINT4)i<length;i++)
00284        {
00285         fprintf(pfone,"% e",SimStochBG1.data->data[i]);
00286         fprintf(pftwo,"% e",SimStochBG2.data->data[i]);
00287        } 
00288      fprintf(pfone,"</real_8>");fprintf(pftwo,"</real_8>");
00289      fprintf(pfone," </ilwd>");fprintf(pftwo," </ilwd>");       
00290      LALFclose(pfone);LALFclose(pftwo);
00291     }
00292   }
00293 }
00294 \end{verbatim}
00295 
00296 \subsubsection*{Uses}
00297 \begin{verbatim}
00298 LALSSSimStochBGTimeSeries()
00299 LALUpdateCalibration()
00300 LALResponseConvert()
00301 \end{verbatim}
00302 </lalLaTeX> */
00303 
00304 
00305 #include <math.h>
00306 #include <string.h>
00307 #include <stdio.h>
00308 #ifdef HAVE_UNISTD_H
00309 #include <unistd.h>
00310 #endif
00311 #ifdef HAVE_GETOPT_H
00312 #include <getopt.h>
00313 #endif
00314 #include <lal/LALStdio.h>
00315 #include <lal/LALStdlib.h>
00316 #include <lal/AVFactories.h>
00317 #include <lal/Calibration.h>
00318 #include <lal/LALConstants.h>
00319 #include <lal/LALStatusMacros.h>
00320 #include <lal/StochasticCrossCorrelation.h>
00321 #include <lal/RealFFT.h>
00322 #include <lal/ComplexFFT.h>
00323 #include <lal/Units.h>
00324 #include <lal/Random.h>
00325 #include <lal/SimulateSB.h>
00326 #include <lal/StochasticMC.h>
00327 
00328 
00329 NRCSID(STOCHASTICMCC, "$Id: StochasticMC.c,v 3.24 2007/06/08 14:41:47 bema Exp $");
00330 
00331 static int verbose = 0;
00332 
00333 static void SinusoidalSplice(REAL4Vector **longData, REAL4Vector **shortData, REAL4Vector *output, UINT4 nSpliceSegs, UINT4 offset);
00334  
00335 void LALStochasticMCDso (LALStatus *status,
00336                          SSSimStochBGOutput *MCoutput,
00337                          StochasticMCInput  *MCinput,
00338                          StochasticMCParams *MCparams)
00339 {
00340   
00341   /* This stores the parameters of functions used */
00342   StochasticOmegaGWParameters        parametersOmega;
00343   
00344   /* This stores the structures of SimulateSB */
00345   SSSimStochBGParams                 SBParams;
00346   SSSimStochBGInput                  SBInput; 
00347   SSSimStochBGOutput                 SBOutput;
00348   REAL4TimeSeries                    whitenedSSimStochBG1;
00349   REAL4TimeSeries                    whitenedSSimStochBG2;
00350   
00351   REAL4FrequencySeries               omegaGW;
00352   COMPLEX8FrequencySeries            wFilter1;
00353   COMPLEX8FrequencySeries            wFilter2;
00354 
00355   
00356   /* vector to store response functions of a pair of detectors  */
00357   COMPLEX8Vector                    *response[2]={NULL,NULL};
00358   
00359   /* Counters & output */
00360   INT4                               i,loop;
00361   INT4                               seed;
00362   REAL8                              totnorm;
00363   REAL8                              totnorm2;
00364 
00365   /* various times, frequencies, sample rates */
00366         INT4                               caltime;
00367   UINT4                              starttime;
00368   UINT4                              length;
00369         INT4                               lengthseg;
00370   INT4                               numseg;
00371   UINT4                              caliboffset;
00372   UINT4                              freqlen;
00373   REAL8                              deltaF, deltaT, sRate;
00374   INT4 site1, site2;
00375   REAL8 omegaRef, f0, fRef, alpha;
00376   
00377   /*calibration information*/
00378   CalibrationFunctions          calfuncs1,calfuncs2 ;
00379   CalibrationUpdateParams       calfacts1,calfacts2 ; 
00380 
00381   const LALUnit countPerStrain = {0,{0,0,0,0,0,-1,1},{0,0,0,0,0,0,0}};
00382  
00383   /* initialize status pointer */
00384    INITSTATUS(status, "LALStochasticMC", STOCHASTICMCC);
00385    ATTATCHSTATUSPTR(status);
00386 
00387     /* ERROR CHECKING */
00388   
00389  /***** check input/output structures exist *****/
00390 
00391  /* output structure */
00392  
00393  ASSERT(MCoutput !=NULL,status, 
00394         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00395  ASSERT(MCoutput->SSimStochBG1->data !=NULL,status, 
00396         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00397  ASSERT(MCoutput->SSimStochBG2->data !=NULL,status, 
00398         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00399  ASSERT(MCoutput->SSimStochBG1->data->data !=NULL,status, 
00400         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00401  ASSERT(MCoutput->SSimStochBG2->data->data !=NULL,status, 
00402  STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00403 
00404  /* input structure */
00405  ASSERT(MCinput != NULL, status, 
00406         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00407   
00408  ASSERT(MCinput->calfuncs1.responseFunction->data !=NULL,status, 
00409          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00410  ASSERT(MCinput->calfuncs2.responseFunction->data !=NULL,status, 
00411          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00412  ASSERT(MCinput->calfuncs1.responseFunction->data->data !=NULL,
00413          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00414  ASSERT(MCinput->calfuncs2.responseFunction->data->data !=NULL,
00415         status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00416 
00417  ASSERT(MCinput->calfuncs1.sensingFunction->data !=NULL,status, 
00418          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00419  ASSERT(MCinput->calfuncs2.sensingFunction->data !=NULL,status, 
00420          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00421  ASSERT(MCinput->calfuncs1.sensingFunction->data->data !=NULL,
00422         status, STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00423  ASSERT(MCinput->calfuncs2.sensingFunction->data->data !=NULL,
00424         status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00425 
00426   ASSERT(MCinput->calfacts1.openLoopFactor->data !=NULL,status, 
00427          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00428   ASSERT(MCinput->calfacts2.openLoopFactor->data !=NULL,status, 
00429          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00430   ASSERT(MCinput->calfacts1.openLoopFactor->data->data !=NULL, 
00431          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00432   ASSERT(MCinput->calfacts2.openLoopFactor->data->data !=NULL,
00433          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00434 
00435   ASSERT(MCinput->calfacts1.sensingFactor->data !=NULL,status, 
00436          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00437   ASSERT(MCinput->calfacts2.sensingFactor->data !=NULL,status, 
00438          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00439   ASSERT(MCinput->calfacts1.sensingFactor->data->data !=NULL, 
00440          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00441   ASSERT(MCinput->calfacts2.sensingFactor->data->data !=NULL, 
00442          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00443 
00444 
00445   /************* check parameter structures ***********/
00446 
00447   /* lengthseg is non-zero */
00448   ASSERT(MCparams->lengthseg > 0, status, 
00449          STOCHASTICMCH_ENULLLEN,STOCHASTICMCH_MSGENULLLEN);
00450 
00451   /* numseg is non-zero */
00452   ASSERT(MCparams->numseg > 0, status, 
00453          STOCHASTICMCH_ENULLSEG,STOCHASTICMCH_MSGENULLSEG);
00454 
00455   /* sampling rate is non zero */
00456   ASSERT(MCparams->sRate > 0, status, 
00457          STOCHASTICMCH_ENULLSRATE,STOCHASTICMCH_MSGENULLSRATE);
00458  
00459  
00460   /************* done with null pointers *****************/
00461   
00462   
00463   /**** check for legality ****/
00464   
00465   /* start frequency must not be negative */
00466   if (MCparams->f0 < 0)
00467     {
00468       ABORT(status,STOCHASTICMCH_ENEGFMIN,STOCHASTICMCH_MSGENEGFMIN );
00469     }
00470 
00471   
00472   /* read input stracture for calibration info*/
00473    
00474   calfuncs1 = MCinput->calfuncs1; 
00475   calfuncs2 = MCinput->calfuncs2;
00476   calfacts1 = MCinput->calfacts1; 
00477   calfacts2 = MCinput->calfacts2;
00478 
00479   /*read parameters*/
00480   
00481   fRef = MCparams->fRef;f0 = MCparams->f0;
00482   omegaRef = MCparams->omegaRef; alpha = MCparams->alpha;
00483   site1 = MCparams->site1; site2 = MCparams->site2;
00484   starttime = MCparams->starttime;
00485   lengthseg = MCparams->lengthseg; numseg = MCparams->numseg;
00486   sRate =  MCparams->sRate;
00487   seed = MCparams->seed;
00488 
00489   /*derive other parameters*/
00490   
00491   length = lengthseg * numseg;
00492   freqlen = lengthseg / 2 + 1;
00493   caliboffset =  lengthseg / (2 * sRate);
00494   deltaT = 1.0 / sRate; 
00495   deltaF = 1.0/(deltaT*lengthseg);
00496   caltime = starttime + caliboffset;
00497   
00498    
00499 /** check for mismatches **/
00500   /* epoch */
00501     
00502   if (calfacts1.epoch.gpsSeconds != caltime) 
00503     {
00504       ABORT(status,STOCHASTICMCH_EMMEPOCH,STOCHASTICMCH_MSGEMMEPOCH);
00505     }
00506   if (calfacts2.epoch.gpsSeconds != caltime) 
00507     {
00508       ABORT(status,STOCHASTICMCH_EMMEPOCH,STOCHASTICMCH_MSGEMMEPOCH);
00509     }
00510    if (calfacts1.epoch.gpsNanoSeconds != 0) 
00511     {
00512       ABORT(status,STOCHASTICMCH_EMMEPOCH,STOCHASTICMCH_MSGEMMEPOCH);
00513     }
00514   if (calfacts2.epoch.gpsNanoSeconds != 0) 
00515     {
00516       ABORT(status,STOCHASTICMCH_EMMEPOCH,STOCHASTICMCH_MSGEMMEPOCH);
00517     }
00518   
00519 
00520   /* Create vectors */
00521   
00522   omegaGW.data = NULL;
00523   LALSCreateVector(status->statusPtr, &(omegaGW.data), freqlen);
00524   
00525   whitenedSSimStochBG1.data = NULL;
00526   LALSCreateVector(status->statusPtr, &(whitenedSSimStochBG1.data), lengthseg);
00527   whitenedSSimStochBG2.data = NULL;
00528   LALSCreateVector(status->statusPtr, &(whitenedSSimStochBG2.data), lengthseg);
00529 
00530   
00531   /* define SimulateSBParams */
00532   SBParams.length         = lengthseg; 
00533   SBParams.deltaT         = 1.0/sRate;
00534   SBParams.detectorOne    = lalCachedDetectors[site1];
00535   SBParams.detectorTwo    = lalCachedDetectors[site2]; 
00536   SBParams.SSimStochBGTimeSeries1Unit = lalADCCountUnit;
00537   SBParams.SSimStochBGTimeSeries2Unit = lalADCCountUnit;
00538 
00539   
00540   /* find omegaGW */
00541   parametersOmega.length   = freqlen;
00542   parametersOmega.f0       = f0;
00543   parametersOmega.deltaF   = deltaF;
00544   parametersOmega.alpha    = alpha;
00545   parametersOmega.fRef     = fRef;
00546   parametersOmega.omegaRef = omegaRef;
00547   LALStochasticOmegaGW(status->statusPtr, &omegaGW, &parametersOmega);
00548     
00549   memset( &wFilter1, 0, sizeof(COMPLEX8FrequencySeries) );
00550   memset( &wFilter2, 0, sizeof(COMPLEX8FrequencySeries) );
00551 
00552   for (i=0;i<2;i++)
00553     {
00554       LALCCreateVector(status->statusPtr, &response[i],freqlen);
00555     }
00556   /*wFilter1.epoch.gpsSeconds = caltime;*/
00557   wFilter1.epoch.gpsNanoSeconds = 0;
00558   wFilter1.deltaF = deltaF ;
00559   wFilter1.f0 = f0;
00560   wFilter1.sampleUnits = countPerStrain;
00561   wFilter1.data = response[0];
00562 
00563   /*wFilter2.epoch.gpsSeconds = caltime;*/
00564   wFilter2.epoch.gpsNanoSeconds = 0;
00565   wFilter2.deltaF = deltaF  ;
00566   wFilter2.f0 = f0;
00567   wFilter2.sampleUnits = countPerStrain;
00568   wFilter2.data = response[1];
00569 
00570 
00571   for (loop = 0; loop < numseg; loop ++)
00572    
00573    {
00574     seed = seed + loop*2;
00575     SBParams.seed = seed;
00576          
00577     wFilter1.epoch.gpsSeconds = caltime;
00578     wFilter2.epoch.gpsSeconds = caltime;
00579 
00580     calfacts1.epoch.gpsSeconds = caltime; 
00581     calfacts2.epoch.gpsSeconds = caltime;
00582     
00583     /* compute response function from calibration parameters */
00584     LALUpdateCalibration(status->statusPtr,&calfuncs1,&calfuncs1,&calfacts1);
00585     LALUpdateCalibration(status->statusPtr,&calfuncs2,&calfuncs2,&calfacts2);
00586     LALResponseConvert(status->statusPtr,&wFilter1,calfuncs1.responseFunction);
00587     LALResponseConvert(status->statusPtr,&wFilter2,calfuncs2.responseFunction);
00588      
00589    
00590     /* force DC to be 0 and nyquist to be real */
00591     response[0]->data[0].re = 0.0;
00592     response[0]->data[0].im = 0.0;
00593     response[1]->data[0].re = 0.0;
00594     response[1]->data[0].im = 0.0;
00595     response[0]->data[freqlen-1].im = 0;
00596     response[1]->data[freqlen-1].im = 0;
00597 
00598 
00599     /* define SSSimStochBGInput */
00600     SBInput.omegaGW                  = &omegaGW;
00601     SBInput.whiteningFilter1         = &wFilter1;
00602     SBInput.whiteningFilter2         = &wFilter2;
00603     
00604     SBOutput.SSimStochBG1 = &whitenedSSimStochBG1;
00605     SBOutput.SSimStochBG2 = &whitenedSSimStochBG2;
00606 
00607     /* generate whitened simulated SB data */     
00608     LALSSSimStochBGTimeSeries(status->statusPtr, &SBOutput, &SBInput, &SBParams);
00609      if ( verbose ) 
00610        {
00611 
00612         /* Mean square */
00613         totnorm2=0.0;
00614         for (i=0;i<lengthseg;i++)
00615         totnorm2+=((whitenedSSimStochBG1.data->data[i])*(whitenedSSimStochBG1.data->data[i]));
00616     
00617         totnorm2/=lengthseg;
00618         printf("Mean square of whitened output is: %e\n",totnorm2);
00619   
00620         /* check normalizations */
00621         totnorm=0.0;
00622         for (i=1;(UINT4)i<freqlen;i++){
00623          REAL8 freq=i*sRate/lengthseg;
00624          REAL8 resp = response[0]->data[i].re;
00625          totnorm+=resp*resp*(omegaGW.data->data[i])/(freq*freq*freq);
00626          }
00627         totnorm*=0.3*LAL_H0FAC_SI*LAL_H0FAC_SI*sRate/(LAL_PI*LAL_PI*lengthseg);
00628         printf("Mean square of whitened output should be: %e.  Ratio is %e\n",totnorm,totnorm/totnorm2);
00629         }
00630 
00631       for (i = 0; i < lengthseg; i ++)
00632         {
00633           MCoutput->SSimStochBG1->data->data[i + loop * lengthseg] = whitenedSSimStochBG1.data->data[i];
00634           MCoutput->SSimStochBG2->data->data[i + loop * lengthseg] = whitenedSSimStochBG2.data->data[i];
00635         }
00636      caltime = caltime + caliboffset;            
00637    }
00638    
00639      /* assign parameters and data to output */
00640       
00641    MCoutput->SSimStochBG1->f0 = f0;
00642    MCoutput->SSimStochBG1->deltaT = deltaT;
00643    MCoutput->SSimStochBG1->epoch.gpsSeconds = starttime;
00644    MCoutput->SSimStochBG1->epoch.gpsNanoSeconds = 0;
00645    MCoutput->SSimStochBG1->sampleUnits = lalADCCountUnit;
00646    strncpy( MCoutput->SSimStochBG1->name, 
00647                "Whitened-SimulatedSBOne", LALNameLength );
00648 
00649    MCoutput->SSimStochBG2->f0 = f0;
00650    MCoutput->SSimStochBG2->deltaT = deltaT;
00651    MCoutput->SSimStochBG2->epoch.gpsSeconds = starttime;
00652    MCoutput->SSimStochBG2->epoch.gpsNanoSeconds = 0;
00653    MCoutput->SSimStochBG2->sampleUnits =lalADCCountUnit;
00654    strncpy( MCoutput->SSimStochBG2->name, 
00655                "Whitened-SimulatedSBTwo", LALNameLength );
00656     
00657   
00658   /* clean up, and exit */  
00659 
00660   /* deallocate memory for output vector */ 
00661 
00662   LALSDestroyVector(status->statusPtr, &(omegaGW.data));
00663   for (i=0;i<2;i++){
00664     LALCDestroyVector(status->statusPtr, &response[i]);
00665   }
00666   /* clean up valid data */
00667   LALSDestroyVector(status->statusPtr, &(whitenedSSimStochBG1.data));
00668   LALSDestroyVector(status->statusPtr, &(whitenedSSimStochBG2.data));
00669 
00670   
00671   DETATCHSTATUSPTR(status);
00672   RETURN(status);
00673 }
00674     
00675    
00676 
00677 void LALStochasticMCDsoSplice (LALStatus *status,
00678                          SSSimStochBGOutput *MCoutput,
00679                          StochasticMCInput  *MCinput,
00680                          StochasticMCParams *MCparams)
00681 {
00682   
00683   /* This stores the parameters of functions used */
00684   StochasticOmegaGWParameters        parametersOmega;
00685   
00686   /* This stores the structures of SimulateSB */
00687   SSSimStochBGParams                 SBParams;
00688   SSSimStochBGInput                  SBInput; 
00689   SSSimStochBGOutput                 SBOutput;
00690   REAL4TimeSeries                    whitenedSSimStochBG1;
00691   REAL4TimeSeries                    whitenedSSimStochBG2;
00692   
00693   REAL4FrequencySeries               omegaGW;
00694   COMPLEX8FrequencySeries            wFilter1;
00695   COMPLEX8FrequencySeries            wFilter2;
00696 
00697   REAL4Vector **longTrain1 = NULL, **shortTrain1 = NULL, **longTrain2 = NULL, **shortTrain2 = NULL;
00698   /* vector to store response functions of a pair of detectors  */
00699   COMPLEX8Vector                    *response[2]={NULL,NULL};
00700   
00701   /* Counters & output */
00702   INT4                               i,m,k,loop;
00703   INT4                               seed;
00704 
00705   /* various times, frequencies, sample rates */
00706         INT4                               caltime;
00707   UINT4                              starttime;
00708         INT4                               lengthseg;
00709   UINT4                              length;
00710         INT4                               numseg;
00711         INT4                               numsegsplice;
00712         INT4                               numsegtot;
00713   UINT4                              spliceoffset;
00714   UINT4                              freqlen;
00715   REAL8                              deltaF, deltaT, sRate;
00716   INT4 site1, site2;
00717   REAL8 omegaRef, f0, fRef, alpha;
00718   
00719   /*calibration information*/
00720   CalibrationFunctions          calfuncs1,calfuncs2 ;
00721   CalibrationUpdateParams       calfacts1,calfacts2 ; 
00722 
00723   const LALUnit countPerStrain = {0,{0,0,0,0,0,-1,1},{0,0,0,0,0,0,0}};
00724  
00725   /* initialize status pointer */
00726    INITSTATUS(status, "LALStochasticMC", STOCHASTICMCC);
00727    ATTATCHSTATUSPTR(status);
00728 
00729     /* ERROR CHECKING */
00730   
00731  /***** check input/output structures exist *****/
00732 
00733  /* output structure */
00734  
00735  ASSERT(MCoutput !=NULL,status, 
00736         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00737  ASSERT(MCoutput->SSimStochBG1->data !=NULL,status, 
00738         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00739  ASSERT(MCoutput->SSimStochBG2->data !=NULL,status, 
00740         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00741  ASSERT(MCoutput->SSimStochBG1->data->data !=NULL,status, 
00742         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00743  ASSERT(MCoutput->SSimStochBG2->data->data !=NULL,status, 
00744  STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00745 
00746  /* input structure */
00747  ASSERT(MCinput != NULL, status, 
00748         STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00749   
00750  ASSERT(MCinput->calfuncs1.responseFunction->data !=NULL,status, 
00751          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00752  ASSERT(MCinput->calfuncs2.responseFunction->data !=NULL,status, 
00753          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00754  ASSERT(MCinput->calfuncs1.responseFunction->data->data !=NULL,
00755          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00756  ASSERT(MCinput->calfuncs2.responseFunction->data->data !=NULL,
00757         status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00758 
00759  ASSERT(MCinput->calfuncs1.sensingFunction->data !=NULL,status, 
00760          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00761  ASSERT(MCinput->calfuncs2.sensingFunction->data !=NULL,status, 
00762          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00763  ASSERT(MCinput->calfuncs1.sensingFunction->data->data !=NULL,
00764         status, STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00765  ASSERT(MCinput->calfuncs2.sensingFunction->data->data !=NULL,
00766         status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00767 
00768   ASSERT(MCinput->calfacts1.openLoopFactor->data !=NULL,status, 
00769          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00770   ASSERT(MCinput->calfacts2.openLoopFactor->data !=NULL,status, 
00771          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00772   ASSERT(MCinput->calfacts1.openLoopFactor->data->data !=NULL, 
00773          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00774   ASSERT(MCinput->calfacts2.openLoopFactor->data->data !=NULL,
00775          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00776 
00777   ASSERT(MCinput->calfacts1.sensingFactor->data !=NULL,status, 
00778          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00779   ASSERT(MCinput->calfacts2.sensingFactor->data !=NULL,status, 
00780          STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00781   ASSERT(MCinput->calfacts1.sensingFactor->data->data !=NULL, 
00782          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00783   ASSERT(MCinput->calfacts2.sensingFactor->data->data !=NULL, 
00784          status,STOCHASTICMCH_ENULLP,STOCHASTICMCH_MSGENULLP);
00785 
00786 
00787   /************* check parameter structures ***********/
00788 
00789   /* lengthseg is non-zero */
00790   ASSERT(MCparams->lengthseg > 0, status, 
00791          STOCHASTICMCH_ENULLLEN,STOCHASTICMCH_MSGENULLLEN);
00792 
00793   /* numseg is non-zero */
00794   ASSERT(MCparams->numseg > 0, status, 
00795          STOCHASTICMCH_ENULLSEG,STOCHASTICMCH_MSGENULLSEG);
00796 
00797   /* sampling rate is non zero */
00798   ASSERT(MCparams->sRate > 0, status, 
00799          STOCHASTICMCH_ENULLSRATE,STOCHASTICMCH_MSGENULLSRATE);
00800  
00801  
00802   /************* done with null pointers *****************/
00803   
00804   
00805   /**** check for legality ****/
00806   
00807   /* start frequency must not be negative */
00808   if (MCparams->f0 < 0)
00809     {
00810       ABORT(status,STOCHASTICMCH_ENEGFMIN,STOCHASTICMCH_MSGENEGFMIN );
00811     }
00812 
00813   
00814   /* read input stracture for calibration info*/
00815    
00816   calfuncs1 = MCinput->calfuncs1; 
00817   calfuncs2 = MCinput->calfuncs2;
00818   calfacts1 = MCinput->calfacts1; 
00819   calfacts2 = MCinput->calfacts2;
00820 
00821   /*read parameters*/
00822   
00823   fRef = MCparams->fRef;f0 = MCparams->f0;
00824   omegaRef = MCparams->omegaRef; alpha = MCparams->alpha;
00825   site1 = MCparams->site1; site2 = MCparams->site2;
00826   starttime = MCparams->starttime;
00827   lengthseg = MCparams->