29 char bin_bhns_extr_orbit_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bin_bhns_extr/bin_bhns_extr_orbit.C,v 1.5 2014/10/24 14:10:23 j_novak Exp $" ;
58 #include "bin_bhns_extr.h"
61 #include "utilitaires.h"
65 double fonc_bhns_orbit_ks(
double,
const Param& ) ;
66 double fonc_bhns_orbit_cf(
double,
const Param& ) ;
71 double fact_omeg_max) {
79 double dnulg, asn2, dasn2, nx, dnx, ny, dny, npn, dnpn ;
91 Tenseur dln_auto_div = d_logn_auto_div ;
105 Tenseur dln_comp = d_logn_comp ;
138 cout <<
"Bin_bhns_extr::orbit_omega : unknown value of rot_phi !"
144 Cmp tmp = logn_auto_regu + loggam ;
147 dnulg = dln_auto_div(0)(0, 0, 0, 0) + dln_comp(0)(0, 0, 0, 0)
148 + factx * tmp.
dsdx()(0, 0, 0, 0) ;
153 double nc = nnn(0, 0, 0, 0) ;
154 double a2c = a_car(0, 0, 0, 0) ;
155 asn2 = a2c / (nc * nc) ;
163 double da2c = factx * a_car.
dsdx()(0, 0, 0, 0) ;
164 double dnc = factx * nnn.
dsdx()(0, 0, 0, 0) ;
166 dasn2 = ( da2c - 2. * a2c / nc * dnc ) / (nc*nc) ;
171 nx = shift(0)(0, 0, 0, 0) ;
176 dnx = factx * shift(0).dsdx()(0, 0, 0, 0) ;
181 ny = shift(1)(0, 0, 0, 0) ;
186 dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
193 npn = tmp(0, 0, 0, 0) ;
199 dnpn = factx * tmp.
dsdx()(0, 0, 0, 0) ;
203 cout <<
"Bin_bhns_extr::orbit_omega : "
204 <<
"It should be the relativistic calculation !" << endl ;
208 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. d(nu+log(Gam))/dX : "
210 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. A^2/N^2 : "
212 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. d(A^2/N^2)/dX : "
214 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. N^X : " << nx << endl ;
215 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. dN^X/dX : "
217 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. N^Y : " << ny << endl ;
218 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. dN^Y/dX : "
220 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. N.N : " << npn << endl ;
221 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. d(N.N)/dX : "
244 double omega1 = fact_omeg_min *
omega ;
245 double omega2 = fact_omeg_max *
omega ;
247 cout <<
"Bin_bhns_extr::orbit_omega: omega1, omega2 [rad/s] : "
248 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
255 zero_list(fonc_bhns_orbit_ks, parf, omega1, omega2, nsub,
260 double omeg_min, omeg_max ;
262 cout <<
"Bin_bhns_extr::orbit_omega : " << nzer <<
263 "zero(s) found in the interval [omega1, omega2]." << endl ;
264 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
265 <<
" " << omega2 << endl ;
266 cout <<
"azer : " << *azer << endl ;
267 cout <<
"bzer : " << *bzer << endl ;
270 cout <<
"Bin_bhns_extr::orbit_omega: WARNING : "
271 <<
"no zero detected in the interval" << endl
272 <<
" [" << omega1 * f_unit <<
", "
273 << omega2 * f_unit <<
"] rad/s !" << endl ;
278 double dist_min = fabs(omega2 - omega1) ;
279 int i_dist_min = -1 ;
280 for (
int i=0; i<nzer; i++) {
283 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
285 if (dist < dist_min) {
290 omeg_min = (*azer)(i_dist_min) ;
291 omeg_max = (*bzer)(i_dist_min) ;
297 cout <<
"Bin_bhns_extr:orbit_omega : "
298 <<
"interval selected for the search of the zero : "
299 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
300 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s "
308 double precis = 1.e-13 ;
310 precis, nitermax, niter) ;
312 cout <<
"Bin_bhns_extr::orbit_omega : "
313 <<
"Number of iterations in zerosec for omega : "
316 cout <<
"Bin_bhns_extr::orbit_omega : omega [rad/s] : "
317 <<
omega * f_unit << endl ;
323 double fact_omeg_max) {
331 double dnulg, asn2, dasn2, ny, dny, npn, dnpn ;
343 Tenseur dln_auto_div = d_logn_auto_div ;
371 "Bin_bhns_extr::orbit_omega_cfc : unknown value of rot_phi !"
377 Cmp tmp = logn_auto_regu + logn_comp + loggam ;
380 dnulg = dln_auto_div(0)(0, 0, 0, 0)
381 + factx * tmp.
dsdx()(0, 0, 0, 0) ;
386 double nc = nnn(0, 0, 0, 0) ;
387 double a2c = a_car(0, 0, 0, 0) ;
388 asn2 = a2c / (nc * nc) ;
396 double da2c = factx * a_car.
dsdx()(0, 0, 0, 0) ;
397 double dnc = factx * nnn.
dsdx()(0, 0, 0, 0) ;
399 dasn2 = ( da2c - 2. * a2c / nc * dnc ) / (nc*nc) ;
404 ny = shift(1)(0, 0, 0, 0) ;
409 dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
416 npn = tmp(0, 0, 0, 0) ;
422 dnpn = factx * tmp.
dsdx()(0, 0, 0, 0) ;
426 cout <<
"Bin_bhns_extr::orbit_omega_cfc : "
427 <<
"It should be the relativistic calculation !" << endl ;
431 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(nu+log(Gam))/dX : "
433 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. A^2/N^2 : "
435 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(A^2/N^2)/dX : "
437 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. N^Y : "
439 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. dN^Y/dX : "
441 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. N.N : "
443 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(N.N)/dX : "
463 double omega1 = fact_omeg_min *
omega ;
464 double omega2 = fact_omeg_max *
omega ;
466 cout <<
"Bin_bhns_extr::orbit_omega_cfc: omega1, omega2 [rad/s] : "
467 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
474 zero_list(fonc_bhns_orbit_cf, parf, omega1, omega2, nsub,
479 double omeg_min, omeg_max ;
481 cout <<
"Bin_bhns_extr::orbit_omega_cfc : " << nzer <<
482 "zero(s) found in the interval [omega1, omega2]." << endl ;
483 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
484 <<
" " << omega2 << endl ;
485 cout <<
"azer : " << *azer << endl ;
486 cout <<
"bzer : " << *bzer << endl ;
489 cout <<
"Bin_bhns_extr::orbit_omega_cfc: WARNING : "
490 <<
"no zero detected in the interval" << endl
491 <<
" [" << omega1 * f_unit <<
", "
492 << omega2 * f_unit <<
"] rad/s !" << endl ;
497 double dist_min = fabs(omega2 - omega1) ;
498 int i_dist_min = -1 ;
499 for (
int i=0; i<nzer; i++) {
502 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
504 if (dist < dist_min) {
509 omeg_min = (*azer)(i_dist_min) ;
510 omeg_max = (*bzer)(i_dist_min) ;
516 cout <<
"Bin_bhns_extr:orbit_omega_cfc : "
517 <<
"interval selected for the search of the zero : "
518 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
519 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s "
527 double precis = 1.e-13 ;
529 precis, nitermax, niter) ;
531 cout <<
"Bin_bhns_extr::orbit_omega_cfc : "
532 <<
"Number of iterations in zerosec for omega : "
535 cout <<
"Bin_bhns_extr::orbit_omega_cfc : omega [rad/s] : "
536 <<
omega * f_unit << endl ;
545 double fonc_bhns_orbit_ks(
double om,
const Param& parf) {
566 double bpb = om2 * xc*xc - 2.*om * ny * xc + npn + 2.*msr * nx*nx;
568 dphi_cent = ( asn2* ( om* (ny + xc*dny) - om2*xc - 0.5*dnpn
569 - 2.*msr * nx * dnx + msr * nx*nx / xc )
570 - 0.5*bpb* dasn2 ) / ( 1 - asn2 * bpb ) ;
573 cout <<
"Bin_bhns_extr::orbit_omega_ks : "
574 <<
"It should be the relativistic calculation !" << endl ;
578 return dnulg + dphi_cent ;
582 double fonc_bhns_orbit_cf(
double om,
const Param& parf) {
584 int relat = parf.get_int() ;
586 double xc = parf.get_double(0) ;
587 double dnulg = parf.get_double(1) ;
588 double asn2 = parf.get_double(2) ;
589 double dasn2 = parf.get_double(3) ;
590 double ny = parf.get_double(4) ;
591 double dny = parf.get_double(5) ;
592 double npn = parf.get_double(6) ;
593 double dnpn = parf.get_double(7) ;
600 double bpb = om2 * xc*xc - 2.*om * ny * xc + npn ;
602 dphi_cent = ( asn2* ( om* (ny + xc*dny) - om2*xc - 0.5*dnpn )
603 - 0.5*bpb* dasn2 ) / ( 1 - asn2 * bpb ) ;
606 cout <<
"Bin_bhns_extr::orbit_omega_cf : "
607 <<
"It should be the relativistic calculation !" << endl ;
611 return dnulg + dphi_cent ;