78 double dnulg, p6sl2, dp6sl2 ;
79 double shiftx, shifty, dshiftx, dshifty, shift2, dshift2 ;
80 double x_orb, y_orb, y_separ, xbh_orb, mhsr ;
97 double mass = ggrav * massbh ;
114 cout <<
"Bin_bhns::orbit_omega : unknown value of rot_phi !"
125 dnulg = factx * ( ((lapconf_auto.
dsdx()).val_grid_point(0,0,0,0)
126 + dlapconf_comp(1).val_grid_point(0,0,0,0))
128 - ((confo_auto.
dsdx()).val_grid_point(0,0,0,0)
129 + dconfo_comp(1).val_grid_point(0,0,0,0))
131 + (tmp1.
dsdx()).val_grid_point(0,0,0,0) ) ;
141 p6sl2 =
pow(confo_c,6.) / lapconf_c / lapconf_c ;
148 double dlapconf_c = factx *
149 ( (lapconf_auto.
dsdx()).val_grid_point(0,0,0,0)
150 + dlapconf_comp(1).val_grid_point(0,0,0,0) ) ;
152 double dpsi6_c = 6. * factx *
pow(confo_c,5.)
153 * ((confo_auto.
dsdx()).val_grid_point(0,0,0,0)
154 + dconfo_comp(1).val_grid_point(0,0,0,0)) ;
156 dp6sl2 = (dpsi6_c - 2.*
pow(confo_c,6.)*dlapconf_c/lapconf_c)
157 / lapconf_c / lapconf_c ;
164 shiftx = shift(1).val_grid_point(0,0,0,0) ;
165 shifty = shift(2).val_grid_point(0,0,0,0) ;
172 dshiftx = factx * ((shift_auto(1).dsdx()).val_grid_point(0,0,0,0)
173 + dshift_comp(1,1).val_grid_point(0,0,0,0)) ;
175 dshifty = factx * ((shift_auto(2).dsdx()).val_grid_point(0,0,0,0)
176 + dshift_comp(1,2).val_grid_point(0,0,0,0)) ;
185 tmp2 = shift(1)%shift(1) + shift(2)%shift(2) + shift(3)%shift(3) ;
186 tmp2.std_spectral_base() ;
188 shift2 = tmp2.val_grid_point(0,0,0,0) ;
196 dshift2 = 2.*factx*((shift(1).val_grid_point(0,0,0,0))
197 * ((shift_auto(1).dsdx()).val_grid_point(0,0,0,0)
198 + dshift_comp(1,1).val_grid_point(0,0,0,0))
199 +(shift(2).val_grid_point(0,0,0,0))
200 * ((shift_auto(2).dsdx()).val_grid_point(0,0,0,0)
201 + dshift_comp(1,2).val_grid_point(0,0,0,0))
202 +(shift(3).val_grid_point(0,0,0,0))
203 * ((shift_auto(3).dsdx()).val_grid_point(0,0,0,0)
204 + dshift_comp(1,3).val_grid_point(0,0,0,0)) ) ;
226 cout <<
"Bin_bhns::orbit_omega: central d(log(lap)+log(Gam))/dX : "
228 cout <<
"Bin_bhns::orbit_omega: central psi^6/lapconf^2 : "
230 cout <<
"Bin_bhns::orbit_omega: central d(psi^6/lapconf^2)/dX : "
232 cout <<
"Bin_bhns::orbit_omega: central shift^X : "
234 cout <<
"Bin_bhns::orbit_omega: central shift^Y : "
236 cout <<
"Bin_bhns::orbit_omega: central d(shift^X)/dX : "
238 cout <<
"Bin_bhns::orbit_omega: central d(shift^Y)/dX : "
240 cout <<
"Bin_bhns::orbit_omega: central shift^i shift_i : "
242 cout <<
"Bin_bhns::orbit_omega: central d(shift^i shift_i)/dX : "
267 double omega1 = fact_omeg_min *
omega ;
268 double omega2 = fact_omeg_max *
omega ;
270 cout <<
"Bin_bhns::orbit_omega: omega1, omega2 [rad/s] : "
271 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
279 zero_list(func_binbhns_orbit_ks, parorb, omega1, omega2, nsub,
285 double omeg_min, omeg_max ;
288 cout <<
"Bin_bhns::orbit_omega: " << nzer
289 <<
" zero(s) found in the interval [omega1, omega2]." << endl ;
290 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
291 <<
" " << omega2 << endl ;
292 cout <<
"azer : " << *azer << endl ;
293 cout <<
"bzer : " << *bzer << endl ;
297 "Bin_bhns::orbit_omega: WARNING : no zero detected in the interval"
298 << endl <<
" [" << omega1 * f_unit <<
", "
299 << omega2 * f_unit <<
"] rad/s !" << endl ;
304 double dist_min = fabs(omega2 - omega1) ;
305 int i_dist_min = -1 ;
306 for (
int i=0; i<nzer; i++) {
310 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
312 if (dist < dist_min) {
317 omeg_min = (*azer)(i_dist_min) ;
318 omeg_max = (*bzer)(i_dist_min) ;
324 cout <<
"Bin_bhns::orbit_omega: interval selected for the search of the zero : "
325 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
326 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s " << endl ;
333 double precis = 1.e-13 ;
334 omega =
zerosec_b(func_binbhns_orbit_ks, parorb, omeg_min, omeg_max,
335 precis, nitermax, niter) ;
337 cout <<
"Bin_bhns::orbit_omega: Number of iterations in zerosec for omega : "
340 cout <<
"Bin_bhns::orbit_omega: omega [rad/s] : "
341 <<
omega * f_unit << endl ;
351 double dnulg, p6sl2, dp6sl2 ;
352 double shiftx, shifty, dshiftx, dshifty, shift2, dshift2 ;
353 double x_orb, y_orb ;
384 cout <<
"Bin_bhns::orbit_omega: unknown value of rot_phi !"
395 dnulg = factx * ( ((lapconf_auto.
dsdx()).val_grid_point(0,0,0,0)
396 + dlapconf_comp(1).val_grid_point(0,0,0,0))
398 - ((confo_auto.
dsdx()).val_grid_point(0,0,0,0)
399 + dconfo_comp(1).val_grid_point(0,0,0,0))
401 + (tmp1.
dsdx()).val_grid_point(0,0,0,0) ) ;
411 p6sl2 =
pow(confo_c,6.) / lapconf_c / lapconf_c ;
418 double dlapconf_c = factx *
419 ( (lapconf_auto.
dsdx()).val_grid_point(0,0,0,0)
420 + dlapconf_comp(1).val_grid_point(0,0,0,0) ) ;
422 double dpsi6_c = 6. * factx *
pow(confo_c,5.)
423 * ((confo_auto.
dsdx()).val_grid_point(0,0,0,0)
424 + dconfo_comp(1).val_grid_point(0,0,0,0)) ;
426 dp6sl2 = (dpsi6_c - 2.*
pow(confo_c,6.)*dlapconf_c/lapconf_c)
427 / lapconf_c / lapconf_c ;
434 shiftx = shift(1).val_grid_point(0,0,0,0) ;
435 shifty = shift(2).val_grid_point(0,0,0,0) ;
442 dshiftx = factx * ( (shift_auto(1).dsdx()).val_grid_point(0,0,0,0)
443 + dshift_comp(1,1).val_grid_point(0,0,0,0) ) ;
445 dshifty = factx * ( (shift_auto(2).dsdx()).val_grid_point(0,0,0,0)
446 + dshift_comp(1,2).val_grid_point(0,0,0,0) ) ;
455 tmp2 = shift(1)%shift(1) + shift(2)%shift(2) + shift(3)%shift(3) ;
456 tmp2.std_spectral_base() ;
458 shift2 = tmp2.val_grid_point(0,0,0,0) ;
466 dshift2 = 2.*factx*( (shift(1).val_grid_point(0,0,0,0))
467 * ((shift_auto(1).dsdx()).val_grid_point(0,0,0,0)
468 + dshift_comp(1,1).val_grid_point(0,0,0,0))
469 +(shift(2).val_grid_point(0,0,0,0))
470 * ((shift_auto(2).dsdx()).val_grid_point(0,0,0,0)
471 + dshift_comp(1,2).val_grid_point(0,0,0,0))
472 +(shift(3).val_grid_point(0,0,0,0))
473 * ((shift_auto(3).dsdx()).val_grid_point(0,0,0,0)
474 + dshift_comp(1,3).val_grid_point(0,0,0,0)) ) ;
488 cout <<
"Bin_bhns::orbit_omega: central d(log(lap)+log(Gam))/dX : "
490 cout <<
"Bin_bhns::orbit_omega: central psi^6/lapconf^2 : "
492 cout <<
"Bin_bhns::orbit_omega: central d(psi^6/lapconf^2)/dX : "
494 cout <<
"Bin_bhns::orbit_omega: central shift^X : "
496 cout <<
"Bin_bhns::orbit_omega: central shift^Y : "
498 cout <<
"Bin_bhns::orbit_omega: central d(shift^X)/dX : "
500 cout <<
"Bin_bhns::orbit_omega: central d(shift^Y)/dX : "
502 cout <<
"Bin_bhns::orbit_omega: central shift^i shift_i : "
504 cout <<
"Bin_bhns::orbit_omega: central d(shift^i shift_i)/dX : "
525 double omega1 = fact_omeg_min *
omega ;
526 double omega2 = fact_omeg_max *
omega ;
528 cout <<
"Bin_bhns::orbit_omega: omega1, omega2 [rad/s] : "
529 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
537 zero_list(func_binbhns_orbit_is, parorb, omega1, omega2, nsub,
543 double omeg_min, omeg_max ;
546 cout <<
"Bin_bhns::orbit_omega: " << nzer
547 <<
" zero(s) found in the interval [omega1, omega2]." << endl ;
548 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
549 <<
" " << omega2 << endl ;
550 cout <<
"azer : " << *azer << endl ;
551 cout <<
"bzer : " << *bzer << endl ;
555 "Bin_bhns::orbit_omega: WARNING : no zero detected in the interval"
556 << endl <<
" [" << omega1 * f_unit <<
", "
557 << omega2 * f_unit <<
"] rad/s !" << endl ;
562 double dist_min = fabs(omega2 - omega1) ;
563 int i_dist_min = -1 ;
564 for (
int i=0; i<nzer; i++) {
568 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
570 if (dist < dist_min) {
575 omeg_min = (*azer)(i_dist_min) ;
576 omeg_max = (*bzer)(i_dist_min) ;
582 cout <<
"Bin_bhns::orbit_omega: interval selected for the search of the zero : "
583 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
584 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s " << endl ;
591 double precis = 1.e-13 ;
592 omega =
zerosec_b(func_binbhns_orbit_is, parorb, omeg_min, omeg_max,
593 precis, nitermax, niter) ;
595 cout <<
"Bin_bhns::orbit_omega: Number of iterations in zerosec for omega : "
598 cout <<
"Bin_bhns::orbit_omega: omega [rad/s] : "
599 <<
omega * f_unit << endl ;