Merge branch 'glitch-subtraction-fix' into 'master'

tweaks to glitch subtraction steps for O3a catalog

See merge request lscsoft/bayeswave!167

src/BayesWaveMCMC.c

@@ -96,11 +96,12 @@ void kickstart_glitch_model(struct Data *data, struct Prior *prior, struct Model
   FILE *waveletStartFile;
   double *params = malloc(data->NW*(sizeof(double)));
   int size;
-  int d=0;
+  int d;
   int d_counter;
   
   for(ifo=0; ifo<data->NI; ifo++)
   {
+    d = 0;
     sprintf(filename,"glitch_%s.dat",data->ifos[ifo]);
     if( (waveletStartFile = fopen(filename,"r")) == NULL )
     {
@@ -147,7 +148,16 @@ void kickstart_glitch_model(struct Data *data, struct Prior *prior, struct Model
       fclose(waveletStartFile);
       
       //fill residual
+      FILE *tempfile = fopen("tempfile.dat","w");
       for(i=0; i<data->N; i++) data->r[ifo][i] = data->s[ifo][i] - model->glitch[ifo]->templates[i];
+      for(i=0; i<data->N/2; i++)
+      {
+        data->r[ifo][i] = data->s[ifo][i] - model->glitch[ifo]->templates[i];
+        fprintf(tempfile,"%lg %lg %lg\n",data->r[ifo][2*i]*data->r[ifo][2*i]+data->r[ifo][2*i+1]*data->r[ifo][2*i+1] ,
+                data->s[ifo][2*i]*data->s[ifo][2*i]+data->s[ifo][2*i+1]*data->s[ifo][2*i+1] ,
+                model->glitch[ifo]->templates[2*i]*model->glitch[ifo]->templates[2*i]+model->glitch[ifo]->templates[2*i+1]*model->glitch[ifo]->templates[2*i+1]);
+      }
+      fclose(tempfile);
     }
   }
   free(params);

src/BayesWaveWavelet.c

@@ -97,14 +97,14 @@ void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs)
   double t = imin*dt;
 
   //brute force
-  /*
+  
   t = 0.0;
   for(n=0; n<N; n++)
   {
     arg = (t-t0)/tau;
     hs[n] += Amp * exp(-arg*arg) * cos(LAL_TWOPI*f0*(t-t0) + phi);
     t += dt;
-  }*/
+  }
   
   //recursive phase
   double phase    =  LAL_TWOPI*f0*(t-t0) + phi;
@@ -120,7 +120,7 @@ void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs)
   double dexpdt  = exp( -(2.*dt2)                / tau2 );
   double env     = exp( -((t0-t)*(t0-t))         / tau2 );
   double dexp2it = exp( -((double)(imin)*2.*dt2) / tau2 );
-  
+  /*
   for(n=imin; n<imax; n++)
   {
     //wavelet
@@ -132,7 +132,7 @@ void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs)
     //iterate wavelet amplitude
     env     *= dexpt*dexp2it;
     dexp2it *= dexpdt;
-  }
+  }*/
 }
 
 double SineGaussianDistance(double *params0, double *params1)