Logo coherent WaveBurst  
Library Reference Guide
Logo
egw_model.cc
Go to the documentation of this file.
1 /*
2 # Copyright (C) 2019 Sergey Klimenko, Valentin Necula
3 #
4 # This program is free software: you can redistribute it and/or modify
5 # it under the terms of the GNU General Public License as published by
6 # the Free Software Foundation, either version 3 of the License, or
7 # (at your option) any later version.
8 #
9 # This program is distributed in the hope that it will be useful,
10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 # GNU General Public License for more details.
13 #
14 # You should have received a copy of the GNU General Public License
15 # along with this program. If not, see <https://www.gnu.org/licenses/>.
16 */
17 
18 
19 #include "numpy.hh"
20 
21 #include "stdio.h"
22 
23 static const double pi = 3.14159265358979312;
24 
25 double E_ecc(double e)
26 { double e2 = e*e;
27  return (1+73.0/24*e2+37.0/96*e2*e2)/pow(1.0+e,3.5);
28 }
29 
30 double J_ecc(double e)
31 { return (1.0+7.0/8*e*e)/(1.0+e)/(1+e);
32 }
33 
34 double bisect(double (*f)(double), double a, double b)
35 { double fa = f(a);
36  double fb = f(b);
37  if (fa == 0.) return a;
38  if (fb == 0.) return b;
39 
40  double tol = 1.0e-14;
41 
42  while(fabs(b-a)>tol){
43  double c = 0.5*(a+b);
44  double fc = f(c);
45  if (fc==0.)return c;
46  if (fb*fc < 0.0){
47  a = c;
48  fa = fc;
49  }
50  else{
51  b = c;
52  fb = fc;
53  }
54  }
55  return 0.5*(a+b);
56 }
57 
58 double zw_func(double p, double e, double a)
59 { //"Equals zero for zoom-whirl parameters"
60  // Formulae from Glampedakis and Kennefick
61  double Mbh = 1.0;
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) );
65  double C = a*a-Mbh*p;
66  C *= C;
67  double xsq, des = fabs(N*N-4.0*F*C);
68  if (a==0.)des = 0;
69  if (a>0)xsq = (-N - sqrt(des))/2.0/F; // Prograde motion
70  else xsq = (-N + sqrt(des))/2.0/F; // Retrograde motion
71 
72  return p*p-xsq*(1.0+e)*(3.0-e);
73 }
74 
75 double rw_e, rw_a;
76 
77 double rw_f(double p)
78 { return zw_func(p, rw_e, rw_a);
79 }
80 
81 double rzoom_whirl(double e, double a)
82 { //"Value of rp for zoom-whirl orbit in Kerr spacetime"
83  rw_e = e;
84  rw_a = a;
85  double plow, phigh;
86  if (a>0){
87  plow = 1.0+e;
88  phigh = 6.0 + 2.0*e;
89  }
90  else if(a<0){
91  plow = 6.0+2.0*e;
92  phigh = 15;
93  }
94  else return (6.0+2.0*e)/(1.0+e);
95 
96  double p = bisect(rw_f, plow, phigh);
97  return p/(1.0+e);
98 }
99 
100 double scaling_exp_zw(double a, double r0)
101 { // "Instability exponent as a function of BH spin and radius of unstable circular orbit"
102  double r = r0;
103  double m = 1;
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) ) );
109 }
110 
111 double scaling_exp(double e, double a)
112 { double gamma_parabolic = 0.19; //# Exponent for e=1 a=a_para
113  double a_para = 0.5;
114  return scaling_exp_zw(a, rzoom_whirl(e,a))/scaling_exp_zw(a_para, rzoom_whirl(1,a_para))*gamma_parabolic;
115 }
116 
117 
118 void close_encounter_pm(double e, double rmin, double mu, double& e_new, double& rmin_new)
119 { //"Change in orbital parameters due to a close encounter according to Peters & Mathews"
120  double a = mu/rmin;
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);
126 
127  E -= deltaE;
128  L -= deltaL;
129  e_new = sqrt(1.0+2*E*L*L/mu/mu/mu);
130  rmin_new = (e_new-1.0)/2.0/(E/mu);
131 }
132 
133 void gwave_data(double e, double rmin, double mu, double& e_new, double& rmin_new)
134 { // "Change in orbital parameters due to g-wave results from simulations"
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};
142 
143 
144  double deltaE = interp(rmin, rp, E_data, 7)*E_ecc(e)/E_ecc(1);
145  double deltaL = interp(rmin, rp, J_data, 7)*J_ecc(e)/J_ecc(1);
146  E -= deltaE;
147  L -= deltaL;
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);
152 }
153 
154 double getDelta()
155 { return 3;
156 }
157 
158 double getE0(double e, double rmin, double mu, double rc)
159 { return mu*(e-1.0)/2.0/rmin + mu/2.0/rc;
160 }
161 
162 
163 double newt_crit_radius(double e, double a)
164 { // "Effective periastron distance at large separation for the BH-NS system to merge function of eccentricity and BH spin"
165  double rparabolic = 6.89; // # Capture for e=1, a=0.5 orbit
166  double a_para = 0.5;
167  double r = (rparabolic/rzoom_whirl(1.0,a_para))*rzoom_whirl(e,a);
168  return r;
169 }
170 
171 double newt_getDeltaE(double e, double rmin, double mu, double a)
172 { double rc = newt_crit_radius(e, a);
173  double E0 = getE0(e,rmin,mu,rc);
174  double Delta = getDelta();
175  double gamma = scaling_exp(e,a);
176  return E0*( 1 - pow((rmin-rc)/Delta,gamma) );
177 }
178 
179 double newt_getDeltaL(double e, double rmin, double mu, double a)
180 { double rc = newt_crit_radius(e, a);
181  double L0 = mu*(sqrt(rmin*(1.0+e)) - sqrt(rc));
182  double Delta = getDelta();
183  double gamma = scaling_exp(e,a);
184  return L0*( 1 - pow((rmin-rc)/Delta, gamma) );
185 }
186 
187 void newt_close_encounter(double e, double rmin, double mu, double a, double& e_new, double& rmin_new)
188 { //"Change in orbital parameters due to a close encounter according to ZW model"
189  double E = mu*(e-1.0)/2.0/rmin;
190  double L = mu*sqrt(rmin*(1.0+e));
191  double deltaE = newt_getDeltaE(e, rmin, mu, a);
192  double deltaL = newt_getDeltaL(e, rmin, mu, a);
193  double Delta = getDelta();
194  double rc = newt_crit_radius(e, a);
195  // Stitch on P&M for large rmin
196 
197  if ((rmin-rc)>=Delta or (64.0*pi)/5.0*mu*mu*pow(1.0/rmin,3.5)*E_ecc(e)>deltaE)
198  return close_encounter_pm(e, rmin, mu, e_new, rmin_new);
199 
200  E -= deltaE;
201  L -= deltaL;
202  double f = 1.0+2*E*L*L/mu/mu/mu;
203  if (f<0){
204  printf("E=%e deltaE=%e rmin=%e e=%e\n", E, deltaE, rmin, e);
205  f = 0;
206  }
207  e_new = sqrt(f);
208  rmin_new = (e_new-1.0)/2.0/(E/mu);
209 }
210 
211 
212 double newt_amp_enhance(double e, double rmin, double mu, double a)
213 { double gamma = scaling_exp(e,a);
214  double Delta = getDelta();
215  double rc = newt_crit_radius(e, a);
216  double E0 = getE0(e,rmin,mu,rc);
217  double deltaE = E0;
218  if (rmin>rc)deltaE = E0*(1- pow((rmin-rc)/Delta, gamma) );
219 
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;
223  return sqrt(ratio);
224 }
225 
226 
227 double energy_geo(double rp, double e, double a)
228 {
229  double p = rp*(1+e);
230  double Mbh = 1.0;
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;
234  C *=C;
235 
236  double xsq, des = fabs(N*N-4*F*C);
237  if(a==0)des = 0;
238  if (a>=0)xsq = (-N - sqrt(des))/2/F; // Prograde motion
239  else xsq = (-N + sqrt(des))/2/F; // Retrograde motion
240  return sqrt( fabs( 1 - (Mbh/p)*(1.0-e*e)*(1.0-xsq/p/p*(1.0-e*e)) ) );
241 }
242 
243 double ang_mom_geo(double rp, double e, double a)
244 {
245  double p = rp*(1+e);
246  double Mbh = 1.0;
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;
250  C *=C;
251 
252  double xsq, des = fabs(N*N-4*F*C);
253  if(a==0)des = 0;
254  if (a>=0)xsq = (-N - sqrt(des))/2/F; // Prograde motion
255  else xsq = (-N + sqrt(des))/2/F; // Retrograde motion
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;
258 }
259 
260 double ang_mom_eff_geo(double rp, double e, double Jspin, double mu)
261 { double a = Jspin+mu*sqrt((1+e)*rp);
262  double a1 = Jspin+mu*ang_mom_geo(rp, e, a);
263  while(fabs(a-a1)>1.0e-5){
264  a = a1;
265  a1 = Jspin+mu*ang_mom_geo(rp,e,a);
266  }
267  return ang_mom_geo(rp,e,a);
268 }
269 
270 double opg_x;
271 double opg_E;
272 double opg_L;
273 double opg_a;
274 
275 double opg_ecc(double p)
276 { double x = opg_x;
277  double E = opg_E;
278  double des = pow(p/x,4) + 4*p*p*p/x/x*(E*E-1);
279  if (des<0)des = 0;
280  double y = 0.5*(pow(p/x, 2)-sqrt(des));
281  return sqrt(fabs(1-y));
282 }
283 
284 double opg_f(double p)
285 { double e = opg_ecc(p);
286  return opg_E - energy_geo(p/(1+e),e,opg_a);
287 }
288 
289 double opg_g(double p)
290 { return ang_mom_geo(p/2.0,1.0, opg_a)- opg_L;
291 }
292 
293 void orbital_param_geo(double E, double L, double a, double& rp, double& e)
294 { opg_x = L-a*E;
295  opg_E = E;
296  opg_a = a;
297  opg_L = L;
298  if (E>1)printf("Error in orbital_param_geo, E>0\n");
299  double p;
300  if (E==1.)p = bisect(opg_g,1.0,100);
301  else p = bisect(opg_f, 1.0, 100);
302  e = opg_ecc(p);
303  rp = p/(1.0+e);
304 }
305 
306 double geo_crit_radius(double e, double a)
307 { return rzoom_whirl(e, a);
308 }
309 
310 double geo_getDeltaE(double e, double rmin, double mu, double a)
311 { double rc = geo_crit_radius(e, a);
312  double E0 = getE0(e,rmin,mu,rc);
313  double Delta = getDelta();
314  double gamma = scaling_exp(e,a);
315  return E0*(1.0 - pow ((rmin-rc)/Delta, gamma) );
316 }
317 
318 
319 double geo_getDeltaL(double e, double rmin, double mu, double a)
320 { double rc = geo_crit_radius(e, a);
321  double L0 = mu*(sqrt(rmin*(1.0+e)) - sqrt(rc));
322  double Delta = getDelta();
323  double gamma = scaling_exp(e,a);
324  return L0*(1.0 - pow( (rmin-rc)/Delta, gamma ) );
325 }
326 
327 void geo_close_encounter(double e, double rmin, double mu, double a, double& e_new, double& rmin_new)
328 { double E = mu*energy_geo(rmin,e,a);
329  double L = mu*ang_mom_geo(rmin,e,a);
330 
331  double deltaE = geo_getDeltaE(e, rmin, mu, a);
332  double deltaL = geo_getDeltaL(e, rmin, mu, a);
333  double Delta = getDelta();
334  double rc = geo_crit_radius(e, a);
335 
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);
339  }
340  E = E - deltaE;
341  L = L - deltaL;
342  orbital_param_geo(E/mu,L/mu,a,rmin_new, e_new);
343  if (e_new==1)
344  printf("%lf %lf %lf %lf\n", E/mu, L/mu,
345  energy_geo(rmin_new,e_new,a), ang_mom_geo(rmin_new,e_new,a) );
346 }
347 
348 
349 double crit_radius(double e, double a)
350 { // "Effective periastron distance at large separation for the BH-NS system to merge function of eccentricity and BH spin"
351  double rparabolic = 6.89; // Capture for e=1, a=0.5 orbit
352  double a_para = 0.5;
353  return ( rparabolic/rzoom_whirl(1, a_para) )*rzoom_whirl(e,a);
354 }
355 
356 
357 double geo_amp_enhance(double e, double rmin, double mu, double a)
358 { return 1;}
359 
360 double amp_enhance(double e, double rmin, double mu, double a)
361 { double gamma = scaling_exp(e,a);
362  double Delta = getDelta();
363  double rc = crit_radius(e, a);
364  double E0 = getE0(e,rmin,mu,rc);
365  double deltaE = E0;
366  if (rmin>rc)deltaE = E0*(1.0-pow((rmin-rc)/Delta, gamma) );
367 
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;
371  return sqrt(ratio);
372 }
373 
374 
375 void orbital_param_newt_to_geo(double rmin,double e, double a, double& rp, double& ee)
376 { double E = 1+(e-1.0)/2.0/rmin;
377  double L = sqrt(rmin*(1.0+e));
378  orbital_param_geo(E,L,a, rp, ee);
379 }
380 
381 void orbital_param_geo_to_newt(double rmin, double e, double a, double& rmin_new, double& e_new)
382 { double E = energy_geo(rmin,e,a) - 1.0;
383  double L = ang_mom_geo(rmin,e,a);
384  double f = 1.0+2*E*L*L;
385  e_new = sqrt(f);
386  if (E>0)rmin_new = (e_new-1.0)/2.0/E;
387  else if (E==0.) rmin_new = L*L/2.0;
388 }
389 
390 
391 double a_eff(double e, double rp, double mu, double Jbh)
392 { return 0.5*sqrt( (1+e)*rp/(2*6.89) )+Jbh;
393 }
double E_ecc(double e)
Definition: egw_model.cc:25
static const double C
Definition: GNGen.cc:28
double geo_crit_radius(double e, double a)
Definition: egw_model.cc:306
double newt_crit_radius(double e, double a)
Definition: egw_model.cc:163
double bisect(double(*f)(double), double a, double b)
Definition: egw_model.cc:34
double rw_a
Definition: egw_model.cc:75
wavearray< double > a(hp.size())
double rzoom_whirl(double e, double a)
Definition: egw_model.cc:81
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)
Definition: egw_model.cc:357
double J_ecc(double e)
Definition: egw_model.cc:30
double opg_x
Definition: egw_model.cc:270
double ang_mom_geo(double rp, double e, double a)
Definition: egw_model.cc:243
int m
Definition: cwb_net.C:28
#define N
cout<< "SNR "<< xsnr<< endl;wavearray< double > f
Definition: ComputeSNR.C:75
void orbital_param_geo_to_newt(double rmin, double e, double a, double &rmin_new, double &e_new)
Definition: egw_model.cc:381
double newt_getDeltaE(double e, double rmin, double mu, double a)
Definition: egw_model.cc:171
wavearray< double > w
Definition: Test1.C:27
double scaling_exp(double e, double a)
Definition: egw_model.cc:111
double opg_L
Definition: egw_model.cc:272
double geo_getDeltaE(double e, double rmin, double mu, double a)
Definition: egw_model.cc:310
double newt_amp_enhance(double e, double rmin, double mu, double a)
Definition: egw_model.cc:212
double opg_E
Definition: egw_model.cc:271
double rw_e
Definition: egw_model.cc:75
double crit_radius(double e, double a)
Definition: egw_model.cc:349
double getE0(double e, double rmin, double mu, double rc)
Definition: egw_model.cc:158
double amp_enhance(double e, double rmin, double mu, double a)
Definition: egw_model.cc:360
double geo_getDeltaL(double e, double rmin, double mu, double a)
Definition: egw_model.cc:319
printf("total live time: non-zero lags = %10.1f \, liveTot)
double rp
void gwave_data(double e, double rmin, double mu, double &e_new, double &rmin_new)
Definition: egw_model.cc:133
double scaling_exp_zw(double a, double r0)
Definition: egw_model.cc:100
void newt_close_encounter(double e, double rmin, double mu, double a, double &e_new, double &rmin_new)
Definition: egw_model.cc:187
double rw_f(double p)
Definition: egw_model.cc:77
double F
double e
double opg_g(double p)
Definition: egw_model.cc:289
regression r
Definition: Regression_H1.C:44
void close_encounter_pm(double e, double rmin, double mu, double &e_new, double &rmin_new)
Definition: egw_model.cc:118
double getDelta()
Definition: egw_model.cc:154
void orbital_param_geo(double E, double L, double a, double &rp, double &e)
Definition: egw_model.cc:293
double opg_a
Definition: egw_model.cc:273
double newt_getDeltaL(double e, double rmin, double mu, double a)
Definition: egw_model.cc:179
double opg_ecc(double p)
Definition: egw_model.cc:275
double a_eff(double e, double rp, double mu, double Jbh)
Definition: egw_model.cc:391
double interp(double v, double *x, double *y, int n)
Definition: numpy.cc:21
double fabs(const Complex &x)
Definition: numpy.cc:55
void geo_close_encounter(double e, double rmin, double mu, double a, double &e_new, double &rmin_new)
Definition: egw_model.cc:327
double energy_geo(double rp, double e, double a)
Definition: egw_model.cc:227
static const double pi
Definition: egw_model.cc:23
double ang_mom_eff_geo(double rp, double e, double Jspin, double mu)
Definition: egw_model.cc:260
void orbital_param_newt_to_geo(double rmin, double e, double a, double &rp, double &ee)
Definition: egw_model.cc:375
double zw_func(double p, double e, double a)
Definition: egw_model.cc:58
double opg_f(double p)
Definition: egw_model.cc:284
wavearray< double > y
Definition: Test10.C:31
static double Delta
Definition: geodesics.cc:26