Commit 1d703941 authored by Tyson Littenberg's avatar Tyson Littenberg

Merge branch 'fix-td-wavelet' into 'master'

fixed td wavelet

See merge request lscsoft/bayeswave!152
parents 5be47e55 82388338
......@@ -221,9 +221,12 @@ int main(int argc, char *argv[])
for(int n=0; n<N_burnin; n++)
{
fscanf(paramsfile, "%i", &glitch_model_size);
for(int d=1; d<=glitch_model_size; d++)
if(glitch_model_size>0)
{
for(i=0; i<5; i++) fscanf(paramsfile,"%lg", &glitch_model_params[i]);
for(int d=1; d<=glitch_model_size; d++)
{
for(i=0; i<5; i++) fscanf(paramsfile,"%lg", &glitch_model_params[i]);
}
}
}
......@@ -245,17 +248,21 @@ int main(int argc, char *argv[])
fscanf(paramsfile, "%i", &glitch_model_size);
for(int d=1; d<=glitch_model_size; d++)
if(glitch_model_size>0)
{
for(i=0; i<NP; i++) fscanf(paramsfile,"%lg", &glitch_model_params[i]);
for(int d=1; d<=glitch_model_size; d++)
{
for(i=0; i<NP; i++) fscanf(paramsfile,"%lg", &glitch_model_params[i]);
SineGaussianTime(glitch_model, glitch_model_params, NBW, 1, data->bw_seglen);
}
SineGaussianTime(glitch_model, glitch_model_params, NBW, 1, data->bw_seglen);
//copy glitch_model to master glitch array for safe keeping
for(i=0; i<NBW; i++) glitch_model_posterior[i][n] = glitch_model[i];
}
//copy glitch_model to master glitch array for safe keeping
for(i=0; i<NBW; i++) glitch_model_posterior[i][n] = glitch_model[i];
printProgress((double)n/(double)N_burnin);
printProgress((double)(n+1)/(double)N_burnin);
}
printf("\n");
......
......@@ -55,9 +55,22 @@ void SineGaussianBandwidth(double *sigpar, double *fmin, double *fmax)
}
void SineGaussianDuration(double *sigpar, double *tmin, double *tmax)
{
double t0 = sigpar[0];
double f0 = sigpar[1];
double Q = sigpar[2];
double tau = Q/(LAL_TWOPI*f0);
*tmin = t0 - 3.*tau;
*tmax = t0 + 3.*tau;
}
void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs)
{
double f0, t0, Q, Amp, phi;
double tmin, tmax;
t0 = sigpar[0];
f0 = sigpar[1];
......@@ -71,9 +84,18 @@ void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs)
double dt2 = dt*dt;
double arg;
double t = 0.0;
int n;
SineGaussianBandwidth(sigpar,&tmin,&tmax);
int imin = (int)floor(tmin/dt);
int imax = (int)floor(tmax/dt);
if(imin < 1 ) imin = 1;
if(imax > N-1) imax = N;
//recompute start time (might be a bit off from tmin because of sampling)
double t = imin*dt;
//brute force
/*
for(n=0; n<N; n++)
......@@ -84,30 +106,29 @@ void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs)
}*/
//recursive phase
double phase = -LAL_TWOPI*f0*(t0) + phi;
double dphase = LAL_TWOPI*f0*dt;
double dre = cos(dphase)-1;
double dim = sin(dphase);
double cosphase = cos(phase);
double sinphase = sin(phase);
double phase = -LAL_TWOPI*f0*(t-t0) + phi;
double dphase = LAL_TWOPI*f0*dt;
double dre = cos(dphase)-1;
double dim = sin(dphase);
double cosphase = cos(phase);
double sinphase = sin(phase);
//recursive amplitude
double dexpt = exp( (2.*t0*dt - dt2)/(tau2) );
double dexpdt = exp( -2.*dt2/(tau2) );
double env = exp(t0*t0/(tau2));
double dexp2it = 1.0;
double dexpt = exp( (2.*t0*dt - dt2) / tau2 );
double dexpdt = exp( -(2.*dt2) / tau2 );
double env = exp( -((t0-t)*(t0-t)) / tau2 );
double dexp2it = exp( -((double)(imin)*2.*dt2) / tau2 );
for(n=0; n<N; n++)
for(n=imin; n<imax; n++)
{
//arg = (t-t0)/tau;
//hs[n] += Amp * exp(-arg*arg) * cosphase;
//wavelet
hs[n] += Amp * env * cosphase;
//t += dt;
//iterate wavelet phase
recursive_phase_evolution(dre, dim, &cosphase, &sinphase);
//recursive gaussian envelope
//iterate wavelet amplitude
env *= dexpt*dexp2it;
dexp2it *= dexpdt;
}
......
......@@ -26,6 +26,7 @@
double SineGaussianSNR(double *params, double *Snf, double Tobs);
void SineGaussianBandwidth(double *sigpar, double *fmin, double *fmax);
void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs);
void SineGaussianFourier(double *hs, double *sigpar, int N, int flag, double Tobs);
void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs);
void SineGaussianFisher(double *params, struct FisherMatrix *fisher);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment