23 static const double pi = 3.14159265358979312;
27 return (1+73.0/24*e2+37.0/96*e2*e2)/pow(1.0+e,3.5);
31 {
return (1.0+7.0/8*e*e)/(1.0+
e)/(1+e);
34 double bisect(
double (*
f)(
double),
double a,
double b)
37 if (fa == 0.)
return a;
38 if (fb == 0.)
return b;
62 double F = ( p*p*p - 2*Mbh*(3.0+e*
e)*p*p + Mbh*Mbh*pow(3.0+e*e, 2)*p -
63 4*Mbh*a*a*pow(1.0-e*e, 2) )/p/p/p;
64 double N = (2.0/p)*(-Mbh*p*p+(Mbh*Mbh*(3.0+e*e)-a*
a)*p-Mbh*a*a*(1+3*e*e) );
67 double xsq, des =
fabs(N*N-4.0*F*C);
69 if (a>0)xsq = (-N - sqrt(des))/2.0/F;
70 else xsq = (-N + sqrt(des))/2.0/F;
72 return p*p-xsq*(1.0+
e)*(3.0-e);
78 {
return zw_func(p, rw_e, rw_a);
94 else return (6.0+2.0*e)/(1.0+
e);
104 double w, R0 =sqrt(r*r + a*a*(1+2*m/r));
105 if (a>0)w = m/(m*a + r*sqrt(m*r));
106 else w = m/(m*a- r*sqrt(m*r));
107 double Delta = r*r + a*a - 2*m*
r;
108 return r*r/(2*
pi)/sqrt(
fabs( 3.0*r*r*Delta + 4.0*m/w/w*(r*R0*R0*w*w - 4*m*a*w - r + 2*m) ) );
112 {
double gamma_parabolic = 0.19;
121 double E = a*(e-1.0)/2.0;
122 double const b = (64.0*3.14159265358979312)/5.0;
123 double L = mu*sqrt(rmin*(1.0+e));
124 double deltaE = b*mu*mu*pow(1.0/rmin, 3.5)*
E_ecc(e);
125 double deltaL = b*a*a*
J_ecc(e);
129 e_new = sqrt(1.0+2*E*L*L/mu/mu/mu);
130 rmin_new = (e_new-1.0)/2.0/(E/mu);
133 void gwave_data(
double e,
double rmin,
double mu,
double& e_new,
double& rmin_new)
135 double E = mu*(e-1.0)/2.0/rmin;
136 double L = mu*sqrt(rmin*(1.0+e));
137 double rp[7] = {6.95, 7.22, 7.5, 8.75, 10, 12.5, 15};
138 double E_data[7] = { 0.006753, 0.003573, 0.002420,
139 0.000730, 0.000334, 0.000105, 3.081720e-05};
140 double J_data[7] = { 0.071858, 0.044782, 0.035159,
141 0.015815, 0.009692, 0.004601, 2.144661e-03};
148 e_new = rmin_new = 0;
149 if ((1.0+2*E*L*L/mu/mu/mu)<0)
return;
150 e_new = sqrt(1.0+2*E*L*L/mu/mu/mu);
151 rmin_new = (e_new-1.0)/2.0/(E/mu);
158 double getE0(
double e,
double rmin,
double mu,
double rc)
159 {
return mu*(e-1.0)/2.0/rmin + mu/2.0/rc;
165 double rparabolic = 6.89;
173 double E0 =
getE0(e,rmin,mu,rc);
176 return E0*( 1 - pow((rmin-rc)/Delta,gamma) );
181 double L0 = mu*(sqrt(rmin*(1.0+e)) - sqrt(rc));
184 return L0*( 1 - pow((rmin-rc)/Delta, gamma) );
189 double E = mu*(e-1.0)/2.0/rmin;
190 double L = mu*sqrt(rmin*(1.0+e));
197 if ((rmin-rc)>=Delta or (64.0*
pi)/5.0*mu*mu*pow(1.0/rmin,3.5)*
E_ecc(e)>deltaE)
202 double f = 1.0+2*E*L*L/mu/mu/mu;
204 printf(
"E=%e deltaE=%e rmin=%e e=%e\n", E, deltaE, rmin, e);
208 rmin_new = (e_new-1.0)/2.0/(E/mu);
216 double E0 =
getE0(e,rmin,mu,rc);
218 if (rmin>rc)deltaE = E0*(1- pow((rmin-rc)/Delta, gamma) );
220 double deltaEpm = (64.0*
pi)/5.0*mu*mu*pow(1.0/rmin, 3.5)*
E_ecc(e);
221 double ratio = deltaE/deltaEpm;
222 if (
fabs(rmin-rc)>Delta || ratio < 1)ratio = 1;
231 double F = ( p*p*p - 2*Mbh*(3+e*
e)*p*p + p*Mbh*Mbh*pow(3.0+e*e, 2)-4*Mbh*a*a*pow(1.0-e*e,2) )/p/p/p;
232 double N = (2/p)*( -Mbh*p*p + p*( Mbh*Mbh*(3.0+e*e)-a*
a) - Mbh*a*a*(1+3*e*e) );
233 double C = a*a-Mbh*p;
236 double xsq, des =
fabs(N*N-4*F*C);
238 if (a>=0)xsq = (-N - sqrt(des))/2/F;
239 else xsq = (-N + sqrt(des))/2/F;
240 return sqrt(
fabs( 1 - (Mbh/p)*(1.0-e*e)*(1.0-xsq/p/p*(1.0-e*e)) ) );
247 double F = ( p*p*p - 2*Mbh*(3+e*
e)*p*p + p*Mbh*Mbh*pow(3.0+e*e, 2)-4*Mbh*a*a*pow(1.0-e*e,2) )/p/p/p;
248 double N = (2/p)*( -Mbh*p*p + p*( Mbh*Mbh*(3.0+e*e)-a*
a) - Mbh*a*a*(1+3*e*e) );
249 double C = a*a-Mbh*p;
252 double xsq, des =
fabs(N*N-4*F*C);
254 if (a>=0)xsq = (-N - sqrt(des))/2/F;
255 else xsq = (-N + sqrt(des))/2/F;
256 double E = sqrt( 1 - (Mbh/p)*(1.0-e*e)*(1.0-xsq/p/p*(1.0-e*e)) );
257 return sqrt(
fabs(xsq)) + a*E;
261 {
double a = Jspin+mu*sqrt((1+e)*rp);
263 while(
fabs(a-a1)>1.0e-5){
278 double des = pow(p/x,4) + 4*p*p*p/x/x*(E*E-1);
280 double y = 0.5*(pow(p/x, 2)-sqrt(des));
281 return sqrt(
fabs(1-y));
298 if (E>1)
printf(
"Error in orbital_param_geo, E>0\n");
312 double E0 =
getE0(e,rmin,mu,rc);
315 return E0*(1.0 - pow ((rmin-rc)/Delta, gamma) );
321 double L0 = mu*(sqrt(rmin*(1.0+e)) - sqrt(rc));
324 return L0*(1.0 - pow( (rmin-rc)/Delta, gamma ) );
336 if ((rmin-rc)>=Delta || (64.0*
pi)/5.0*mu*mu*pow(1.0/rmin, 3.5)*
E_ecc(e)>deltaE){
337 deltaE = (64.0*
pi)/5.0*mu*mu*pow(1.0/rmin,3.5)*
E_ecc(e);
338 deltaL = (64.0*
pi)/5.0*mu*mu/rmin/rmin*
J_ecc(e);
344 printf(
"%lf %lf %lf %lf\n", E/mu, L/mu,
351 double rparabolic = 6.89;
364 double E0 =
getE0(e,rmin,mu,rc);
366 if (rmin>rc)deltaE = E0*(1.0-pow((rmin-rc)/Delta, gamma) );
368 double deltaEpm = (64.0*
pi)/5.0*mu*mu*pow(1.0/rmin,3.5)*
E_ecc(e);
369 double ratio = deltaE/deltaEpm;
370 if (
fabs(rmin-rc)>Delta || ratio < 1)ratio = 1.0;
376 {
double E = 1+(e-1.0)/2.0/rmin;
377 double L = sqrt(rmin*(1.0+e));
384 double f = 1.0+2*E*L*
L;
386 if (E>0)rmin_new = (e_new-1.0)/2.0/E;
387 else if (E==0.) rmin_new = L*L/2.0;
391 double a_eff(
double e,
double rp,
double mu,
double Jbh)
392 {
return 0.5*sqrt( (1+e)*rp/(2*6.89) )+Jbh;
double geo_crit_radius(double e, double a)
double newt_crit_radius(double e, double a)
double bisect(double(*f)(double), double a, double b)
wavearray< double > a(hp.size())
double rzoom_whirl(double e, double a)
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
double geo_amp_enhance(double e, double rmin, double mu, double a)
double ang_mom_geo(double rp, double e, double a)
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
void orbital_param_geo_to_newt(double rmin, double e, double a, double &rmin_new, double &e_new)
double newt_getDeltaE(double e, double rmin, double mu, double a)
double scaling_exp(double e, double a)
double geo_getDeltaE(double e, double rmin, double mu, double a)
double newt_amp_enhance(double e, double rmin, double mu, double a)
double crit_radius(double e, double a)
double getE0(double e, double rmin, double mu, double rc)
double amp_enhance(double e, double rmin, double mu, double a)
double geo_getDeltaL(double e, double rmin, double mu, double a)
printf("total live time: non-zero lags = %10.1f \, liveTot)
void gwave_data(double e, double rmin, double mu, double &e_new, double &rmin_new)
double scaling_exp_zw(double a, double r0)
void newt_close_encounter(double e, double rmin, double mu, double a, double &e_new, double &rmin_new)
void close_encounter_pm(double e, double rmin, double mu, double &e_new, double &rmin_new)
void orbital_param_geo(double E, double L, double a, double &rp, double &e)
double newt_getDeltaL(double e, double rmin, double mu, double a)
double a_eff(double e, double rp, double mu, double Jbh)
double interp(double v, double *x, double *y, int n)
double fabs(const Complex &x)
void geo_close_encounter(double e, double rmin, double mu, double a, double &e_new, double &rmin_new)
double energy_geo(double rp, double e, double a)
double ang_mom_eff_geo(double rp, double e, double Jspin, double mu)
void orbital_param_newt_to_geo(double rmin, double e, double a, double &rp, double &ee)
double zw_func(double p, double e, double a)