31 inline double min(
double x,
double y)
35 {
const double C = 299792458;
36 const double G = 6.67428e-11;
38 const double L = G/C/C*SolarMass*totalMass;
41 double* tt =
new double[
N];
43 for(
int i=0;
i<
N; ++
i){
49 double tMax = tt[N-1];
51 if(tMax*Rate<2.1e9)nSamples = tMax*
Rate;
58 for(
int i=0;
i<nSamples; ++
i){
62 double frac = (ts - tt[j-1])/(tt[j] - tt[j-1]);
63 hps[
i] = H[j-1].
Re() + (H[
j].
Re() - H[j-1].
Re())*frac;
64 hxs[
i] = H[j-1].
Im() + (H[
j].
Im() - H[j-1].
Im())*frac;
73 printf(
"Error! masses must be >=0\n");
78 printf(
"Error! rmin must be >=0\n");
82 printf(
"Error! Orbit must be initially parabolic or bound.\n");
99 double t_post_merger = 400;
101 double Jspin = abh*Mbh*Mbh;
103 double a0 = Jspin+mu*L0;
106 if (e0<1 && r0 > rmin0*(1+e0)/(1-e0)) r0 = rmin0*(1+
e0)/(1-e0);
109 int N_plot = t_end/
dt;
112 geodesic geo(r0, phi0, E0, L0, dir, Jspin, mu);
114 bool success = geo.
integrate(t_end, N_plot, 0,
true);
122 double* r_LR =
new double[N_plot];
123 for(
int i=0;
i<N_plot; ++
i)
124 r_LR[
i] = 2*(1+cos(2*(acos(-
min(mu*geo.
Lt[
i] + Jspin, 1.))/3))) - geo.
r[
i];
125 tmerge =
interp(0, r_LR, geo.
t, N_plot, 1., 1.);
129 int Ntsig = (geo.
t[N_plot-1]+t_post_merger)/dt;
130 double* tsig =
new double[Ntsig];
131 for(
int i=0;
i<Ntsig; ++
i) tsig[
i] = dt*
i;
133 double* hr =
interp(tsig, Ntsig, geo.
t, geo.
hplus, N_plot, 0, 0);
134 double* hi =
interp(tsig, Ntsig, geo.
t, geo.
hcross, N_plot, 0, 0);
137 for(
int i=0;
i<Ntsig; ++
i)hsig[
i] =
Complex(hr[
i], hi[i]);
139 a_final = mu*geo.
Lt[N_plot-1] + abh*Mbh*Mbh;
141 double tstart_merger = tmerge;
142 if (tmerge>0)
irs_merger(dt, tmerge, tsig, hsig, Ntsig, a_final, tstart_merger);
144 int retVal =
Sample(hsig, tsig, Ntsig, m1+m2, 16384., Hp, Hx)==0;
wavearray< double > t(hp.size())
double min(double x, double y)
virtual void rate(double r)
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\ layers : "<< nLAYERS<< "\ dF(hz) : "<< dF<< "\ dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1) *itime+ifreq;double time=itime *dT;double freq=(ifreq >0) ? ifreq *dF :dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
bool integrate(double D_t, int &N, double rmin=0, bool stop_alr=false)
int Sample(Complex *H, double *t, int N, double totalMass, double Rate, wavearray< double > &hps, wavearray< double > &hxs)
printf("total live time: non-zero lags = %10.1f \, liveTot)
wavearray< double > ts(N)
double interp(double v, double *x, double *y, int n)
int getEBBH(double m1, double m2, double rmin0, double e0, wavearray< double > &Hp, wavearray< double > &Hx, double t_end)
double energy_geo(double rp, double e, double a)
double ang_mom_eff_geo(double rp, double e, double Jspin, double mu)
virtual void resize(unsigned int)
void irs_merger(double dt, double tmerge, double *tsig, Complex *hsig, int lentsig, double a, double tstart=-1)