LORENE
bin_bhns_orbit.C
1/*
2 * Methods of class Bin_bhns to compute the orbital angular velocity
3 *
4 * (see file bin_bhns.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005-2007 Keisuke Taniguchi
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char bin_bhns_orbit_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_orbit.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $" ;
29
30/*
31 * $Id: bin_bhns_orbit.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
32 * $Log: bin_bhns_orbit.C,v $
33 * Revision 1.4 2014/10/13 08:52:41 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.3 2014/10/06 15:13:00 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.2 2008/05/15 19:01:28 k_taniguchi
40 * Change of some parameters.
41 *
42 * Revision 1.1 2007/06/22 01:10:20 k_taniguchi
43 * *** empty log message ***
44 *
45 *
46 * $Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_orbit.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
47 *
48 */
49
50// C++ headers
51//#include <>
52
53// C headers
54#include <cmath>
55
56// Lorene headers
57#include "bin_bhns.h"
58#include "param.h"
59#include "utilitaires.h"
60#include "unites.h"
61
62namespace Lorene {
63double func_binbhns_orbit_ks(double , const Param& ) ;
64double func_binbhns_orbit_is(double , const Param& ) ;
65
66//**********************************************************************
67
68void Bin_bhns::orbit_omega(double fact_omeg_min, double fact_omeg_max) {
69
70 using namespace Unites ;
71
72 if (hole.is_kerrschild()) {
73
74 //--------------------------------------------------------------------
75 // Evaluation of various quantities at the center of the neutron star
76 //--------------------------------------------------------------------
77
78 double dnulg, p6sl2, dp6sl2 ;
79 double shiftx, shifty, dshiftx, dshifty, shift2, dshift2 ;
80 double x_orb, y_orb, y_separ, xbh_orb, mhsr ;
81
82 const Map& mp = star.get_mp() ;
83
84 const Scalar& lapconf = star.get_lapconf_tot() ;
85 const Scalar& lapconf_auto = star.get_lapconf_auto() ;
86 const Scalar& confo = star.get_confo_tot() ;
87 const Scalar& confo_auto = star.get_confo_auto() ;
88 const Scalar& loggam = star.get_loggam() ;
89 const Vector& shift = star.get_shift_tot() ;
90 const Vector& shift_auto = star.get_shift_auto() ;
91
92 const Vector& dlapconf_comp = star.get_d_lapconf_comp() ;
93 const Vector& dconfo_comp = star.get_d_confo_comp() ;
94 const Tensor& dshift_comp = star.get_d_shift_comp() ;
95
96 const double& massbh = hole.get_mass_bh() ;
97 double mass = ggrav * massbh ;
98
99 //----------------------------------------------------------
100 // Calculation of d/dX( ln(lapconf) - ln(psi) + ln(Gamma) )
101 // at the center of NS
102 //----------------------------------------------------------
103
104 // Factor to translate x --> X
105 double factx ;
106 if (fabs(mp.get_rot_phi()) < 1.e-14) {
107 factx = 1. ;
108 }
109 else {
110 if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
111 factx = - 1. ;
112 }
113 else {
114 cout << "Bin_bhns::orbit_omega : unknown value of rot_phi !"
115 << endl ;
116 abort() ;
117 }
118 }
119
120 Scalar tmp1(mp) ;
121 tmp1 = loggam ;
122 tmp1.std_spectral_base() ;
123
124 // d/dX tmp1
125 dnulg = factx * ( ((lapconf_auto.dsdx()).val_grid_point(0,0,0,0)
126 + dlapconf_comp(1).val_grid_point(0,0,0,0))
127 / lapconf.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))
130 / confo.val_grid_point(0,0,0,0)
131 + (tmp1.dsdx()).val_grid_point(0,0,0,0) ) ;
132
133
134 //----------------------------------------------------
135 // Calculation of psi^6/lapconf^2 at the center of NS
136 //----------------------------------------------------
137
138 double lapconf_c = lapconf.val_grid_point(0,0,0,0) ;
139 double confo_c = confo.val_grid_point(0,0,0,0) ;
140
141 p6sl2 = pow(confo_c,6.) / lapconf_c / lapconf_c ;
142
143
144 //----------------------------------------------------------
145 // Calculation of d/dX(psi^6/lapconf^2) at the center of NS
146 //----------------------------------------------------------
147
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) ) ;
151
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)) ;
155
156 dp6sl2 = (dpsi6_c - 2.*pow(confo_c,6.)*dlapconf_c/lapconf_c)
157 / lapconf_c / lapconf_c ;
158
159
160 //--------------------------------------------------------
161 // Calculation of shift^X and shift^Y at the center of NS
162 //--------------------------------------------------------
163
164 shiftx = shift(1).val_grid_point(0,0,0,0) ;
165 shifty = shift(2).val_grid_point(0,0,0,0) ;
166
167
168 //------------------------------------------------------------------
169 // Calculation of d shift^X/dX and d shift^Y/dX at the center of NS
170 //------------------------------------------------------------------
171
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)) ;
174
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)) ;
177
178
179 //--------------------------------------------------------
180 // Calculation of (shift^X)^2 + (shift^Y)^2 + (shift^Z)^2
181 // at the center of NS
182 //--------------------------------------------------------
183
184 Scalar tmp2(mp) ;
185 tmp2 = shift(1)%shift(1) + shift(2)%shift(2) + shift(3)%shift(3) ;
186 tmp2.std_spectral_base() ;
187
188 shift2 = tmp2.val_grid_point(0,0,0,0) ;
189
190
191 //----------------------------------------------------------------
192 // Calculation of d/dX( (shift^X)^2 + (shift^Y)^2 + (shift^Z)^2 )
193 // at the center of NS
194 //----------------------------------------------------------------
195
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)) ) ;
205
206
207 //-------------------------
208 // Information of position
209 //-------------------------
210
211 x_orb = (star.get_mp()).get_ori_x() - x_rot ;
212 y_orb = (star.get_mp()).get_ori_y() - y_rot ;
213 y_separ = (star.get_mp()).get_ori_y() - (hole.get_mp()).get_ori_y() ;
214
215 xbh_orb = (hole.get_mp()).get_ori_x() - x_rot ;
216
217 //------------------------------
218 // Calculation of H_BH / r_bh^4
219 //------------------------------
220
221 mhsr = mass / pow( separ*separ+y_separ*y_separ, 2.5 ) ;
222
223 // Output
224 // ------
225
226 cout << "Bin_bhns::orbit_omega: central d(log(lap)+log(Gam))/dX : "
227 << dnulg << endl ;
228 cout << "Bin_bhns::orbit_omega: central psi^6/lapconf^2 : "
229 << p6sl2 << endl ;
230 cout << "Bin_bhns::orbit_omega: central d(psi^6/lapconf^2)/dX : "
231 << dp6sl2 << endl ;
232 cout << "Bin_bhns::orbit_omega: central shift^X : "
233 << shiftx << endl ;
234 cout << "Bin_bhns::orbit_omega: central shift^Y : "
235 << shifty << endl ;
236 cout << "Bin_bhns::orbit_omega: central d(shift^X)/dX : "
237 << dshiftx << endl ;
238 cout << "Bin_bhns::orbit_omega: central d(shift^Y)/dX : "
239 << dshifty << endl ;
240 cout << "Bin_bhns::orbit_omega: central shift^i shift_i : "
241 << shift2 << endl ;
242 cout << "Bin_bhns::orbit_omega: central d(shift^i shift_i)/dX : "
243 << dshift2 << endl ;
244
245
246 //---------------------------------------------------------------//
247 // Calculation of the orbital angular velocity //
248 //---------------------------------------------------------------//
249
250 Param parorb ;
251 parorb.add_double(dnulg, 0) ;
252 parorb.add_double(p6sl2, 1) ;
253 parorb.add_double(dp6sl2, 2) ;
254 parorb.add_double(shiftx, 3) ;
255 parorb.add_double(shifty, 4) ;
256 parorb.add_double(dshiftx, 5) ;
257 parorb.add_double(dshifty, 6) ;
258 parorb.add_double(shift2, 7) ;
259 parorb.add_double(dshift2, 8) ;
260 parorb.add_double(x_orb, 9) ;
261 parorb.add_double(y_orb, 10) ;
262 parorb.add_double(separ, 11) ;
263 parorb.add_double(y_separ, 12) ;
264 parorb.add_double(xbh_orb, 13) ;
265 parorb.add_double(mhsr, 14) ;
266
267 double omega1 = fact_omeg_min * omega ;
268 double omega2 = fact_omeg_max * omega ;
269
270 cout << "Bin_bhns::orbit_omega: omega1, omega2 [rad/s] : "
271 << omega1 * f_unit << " " << omega2 * f_unit << endl ;
272
273 // Search for the various zeros in the interval [omega1,omega2]
274 // ------------------------------------------------------------
275
276 int nsub = 50 ; // total number of subdivisions of the interval
277 Tbl* azer = 0x0 ;
278 Tbl* bzer = 0x0 ;
279 zero_list(func_binbhns_orbit_ks, parorb, omega1, omega2, nsub,
280 azer, bzer) ;
281
282 // Search for the zero closest to the previous value of omega
283 // ----------------------------------------------------------
284
285 double omeg_min, omeg_max ;
286 int nzer = azer->get_taille() ; // number of zeros found by zero_list
287
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 ;
294
295 if (nzer == 0) {
296 cout <<
297 "Bin_bhns::orbit_omega: WARNING : no zero detected in the interval"
298 << endl << " [" << omega1 * f_unit << ", "
299 << omega2 * f_unit << "] rad/s !" << endl ;
300 omeg_min = omega1 ;
301 omeg_max = omega2 ;
302 }
303 else {
304 double dist_min = fabs(omega2 - omega1) ;
305 int i_dist_min = -1 ;
306 for (int i=0; i<nzer; i++) {
307 // Distance of previous value of omega from the center of the
308 // interval [azer(i), bzer(i)]
309
310 double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
311
312 if (dist < dist_min) {
313 dist_min = dist ;
314 i_dist_min = i ;
315 }
316 }
317 omeg_min = (*azer)(i_dist_min) ;
318 omeg_max = (*bzer)(i_dist_min) ;
319 }
320
321 delete azer ; // Tbl allocated by zero_list
322 delete bzer ;
323
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 ;
327
328 // Computation of the zero in the selected interval by the secant method
329 // ---------------------------------------------------------------------
330
331 int nitermax = 200 ;
332 int niter ;
333 double precis = 1.e-13 ;
334 omega = zerosec_b(func_binbhns_orbit_ks, parorb, omeg_min, omeg_max,
335 precis, nitermax, niter) ;
336
337 cout << "Bin_bhns::orbit_omega: Number of iterations in zerosec for omega : "
338 << niter << endl ;
339
340 cout << "Bin_bhns::orbit_omega: omega [rad/s] : "
341 << omega * f_unit << endl ;
342
343
344 }
345 else { // Isotropic coordinates with the maximal slicing
346
347 //--------------------------------------------------------------------
348 // Evaluation of various quantities at the center of the neutron star
349 //--------------------------------------------------------------------
350
351 double dnulg, p6sl2, dp6sl2 ;
352 double shiftx, shifty, dshiftx, dshifty, shift2, dshift2 ;
353 double x_orb, y_orb ;
354
355 const Map& mp = star.get_mp() ;
356
357 const Scalar& lapconf = star.get_lapconf_tot() ;
358 const Scalar& lapconf_auto = star.get_lapconf_auto() ;
359 const Scalar& confo = star.get_confo_tot() ;
360 const Scalar& confo_auto = star.get_confo_auto() ;
361 const Scalar& loggam = star.get_loggam() ;
362 const Vector& shift = star.get_shift_tot() ;
363 const Vector& shift_auto = star.get_shift_auto() ;
364
365 const Vector& dlapconf_comp = star.get_d_lapconf_comp() ;
366 const Vector& dconfo_comp = star.get_d_confo_comp() ;
367 const Tensor& dshift_comp = star.get_d_shift_comp() ;
368
369 //----------------------------------------------------------
370 // Calculation of d/dX( ln(lapconf) - ln(psi) + ln(Gamma) )
371 // at the center of NS
372 //----------------------------------------------------------
373
374 // Factor to translate x --> X
375 double factx ;
376 if (fabs(mp.get_rot_phi()) < 1.e-14) {
377 factx = 1. ;
378 }
379 else {
380 if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
381 factx = - 1. ;
382 }
383 else {
384 cout << "Bin_bhns::orbit_omega: unknown value of rot_phi !"
385 << endl ;
386 abort() ;
387 }
388 }
389
390 Scalar tmp1(mp) ;
391 tmp1 = loggam ;
392 tmp1.std_spectral_base() ;
393
394 // d/dX tmp1
395 dnulg = factx * ( ((lapconf_auto.dsdx()).val_grid_point(0,0,0,0)
396 + dlapconf_comp(1).val_grid_point(0,0,0,0))
397 / lapconf.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))
400 / confo.val_grid_point(0,0,0,0)
401 + (tmp1.dsdx()).val_grid_point(0,0,0,0) ) ;
402
403
404 //----------------------------------------------------
405 // Calculation of psi^6/lapconf^2 at the center of NS
406 //----------------------------------------------------
407
408 double lapconf_c = lapconf.val_grid_point(0,0,0,0) ;
409 double confo_c = confo.val_grid_point(0,0,0,0) ;
410
411 p6sl2 = pow(confo_c,6.) / lapconf_c / lapconf_c ;
412
413
414 //----------------------------------------------------------
415 // Calculation of d/dX(psi^6/lapconf^2) at the center of NS
416 //----------------------------------------------------------
417
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) ) ;
421
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)) ;
425
426 dp6sl2 = (dpsi6_c - 2.*pow(confo_c,6.)*dlapconf_c/lapconf_c)
427 / lapconf_c / lapconf_c ;
428
429
430 //--------------------------------------------------------
431 // Calculation of shift^X and shift^Y at the center of NS
432 //--------------------------------------------------------
433
434 shiftx = shift(1).val_grid_point(0,0,0,0) ;
435 shifty = shift(2).val_grid_point(0,0,0,0) ;
436
437
438 //------------------------------------------------------------------
439 // Calculation of d shift^X/dX and d shift^Y/dX at the center of NS
440 //------------------------------------------------------------------
441
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) ) ;
444
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) ) ;
447
448
449 //--------------------------------------------------------
450 // Calculation of (shift^X)^2 + (shift^Y)^2 + (shift^Z)^2
451 // at the center of NS
452 //--------------------------------------------------------
453
454 Scalar tmp2(mp) ;
455 tmp2 = shift(1)%shift(1) + shift(2)%shift(2) + shift(3)%shift(3) ;
456 tmp2.std_spectral_base() ;
457
458 shift2 = tmp2.val_grid_point(0,0,0,0) ;
459
460
461 //----------------------------------------------------------------
462 // Calculation of d/dX( (shift^X)^2 + (shift^Y)^2 + (shift^Z)^2 )
463 // at the center of NS
464 //----------------------------------------------------------------
465
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)) ) ;
475
476
477 //-------------------------
478 // Information of position
479 //-------------------------
480
481 x_orb = (star.get_mp()).get_ori_x() - x_rot ;
482 y_orb = (star.get_mp()).get_ori_y() - y_rot ;
483
484
485 // Output
486 // ------
487
488 cout << "Bin_bhns::orbit_omega: central d(log(lap)+log(Gam))/dX : "
489 << dnulg << endl ;
490 cout << "Bin_bhns::orbit_omega: central psi^6/lapconf^2 : "
491 << p6sl2 << endl ;
492 cout << "Bin_bhns::orbit_omega: central d(psi^6/lapconf^2)/dX : "
493 << dp6sl2 << endl ;
494 cout << "Bin_bhns::orbit_omega: central shift^X : "
495 << shiftx << endl ;
496 cout << "Bin_bhns::orbit_omega: central shift^Y : "
497 << shifty << endl ;
498 cout << "Bin_bhns::orbit_omega: central d(shift^X)/dX : "
499 << dshiftx << endl ;
500 cout << "Bin_bhns::orbit_omega: central d(shift^Y)/dX : "
501 << dshifty << endl ;
502 cout << "Bin_bhns::orbit_omega: central shift^i shift_i : "
503 << shift2 << endl ;
504 cout << "Bin_bhns::orbit_omega: central d(shift^i shift_i)/dX : "
505 << dshift2 << endl ;
506
507
508 //---------------------------------------------------------------//
509 // Calculation of the orbital angular velocity //
510 //---------------------------------------------------------------//
511
512 Param parorb ;
513 parorb.add_double(dnulg, 0) ;
514 parorb.add_double(p6sl2, 1) ;
515 parorb.add_double(dp6sl2, 2) ;
516 parorb.add_double(shiftx, 3) ;
517 parorb.add_double(shifty, 4) ;
518 parorb.add_double(dshiftx, 5) ;
519 parorb.add_double(dshifty, 6) ;
520 parorb.add_double(shift2, 7) ;
521 parorb.add_double(dshift2, 8) ;
522 parorb.add_double(x_orb, 9) ;
523 parorb.add_double(y_orb, 10) ;
524
525 double omega1 = fact_omeg_min * omega ;
526 double omega2 = fact_omeg_max * omega ;
527
528 cout << "Bin_bhns::orbit_omega: omega1, omega2 [rad/s] : "
529 << omega1 * f_unit << " " << omega2 * f_unit << endl ;
530
531 // Search for the various zeros in the interval [omega1,omega2]
532 // ------------------------------------------------------------
533
534 int nsub = 50 ; // total number of subdivisions of the interval
535 Tbl* azer = 0x0 ;
536 Tbl* bzer = 0x0 ;
537 zero_list(func_binbhns_orbit_is, parorb, omega1, omega2, nsub,
538 azer, bzer) ;
539
540 // Search for the zero closest to the previous value of omega
541 // ----------------------------------------------------------
542
543 double omeg_min, omeg_max ;
544 int nzer = azer->get_taille() ; // number of zeros found by zero_list
545
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 ;
552
553 if (nzer == 0) {
554 cout <<
555 "Bin_bhns::orbit_omega: WARNING : no zero detected in the interval"
556 << endl << " [" << omega1 * f_unit << ", "
557 << omega2 * f_unit << "] rad/s !" << endl ;
558 omeg_min = omega1 ;
559 omeg_max = omega2 ;
560 }
561 else {
562 double dist_min = fabs(omega2 - omega1) ;
563 int i_dist_min = -1 ;
564 for (int i=0; i<nzer; i++) {
565 // Distance of previous value of omega from the center of the
566 // interval [azer(i), bzer(i)]
567
568 double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
569
570 if (dist < dist_min) {
571 dist_min = dist ;
572 i_dist_min = i ;
573 }
574 }
575 omeg_min = (*azer)(i_dist_min) ;
576 omeg_max = (*bzer)(i_dist_min) ;
577 }
578
579 delete azer ; // Tbl allocated by zero_list
580 delete bzer ;
581
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 ;
585
586 // Computation of the zero in the selected interval by the secant method
587 // ---------------------------------------------------------------------
588
589 int nitermax = 200 ;
590 int niter ;
591 double precis = 1.e-13 ;
592 omega = zerosec_b(func_binbhns_orbit_is, parorb, omeg_min, omeg_max,
593 precis, nitermax, niter) ;
594
595 cout << "Bin_bhns::orbit_omega: Number of iterations in zerosec for omega : "
596 << niter << endl ;
597
598 cout << "Bin_bhns::orbit_omega: omega [rad/s] : "
599 << omega * f_unit << endl ;
600
601 }
602}
603
604//******************************
605// Function for searching omega
606//******************************
607
608double func_binbhns_orbit_ks(double om, const Param& parorb) {
609
610 double dnulg = parorb.get_double(0) ;
611 double p6sl2 = parorb.get_double(1) ;
612 double dp6sl2 = parorb.get_double(2) ;
613 double shiftx = parorb.get_double(3) ;
614 double shifty = parorb.get_double(4) ;
615 double dshiftx = parorb.get_double(5) ;
616 double dshifty = parorb.get_double(6) ;
617 double shift2 = parorb.get_double(7) ;
618 double dshift2 = parorb.get_double(8) ;
619 double x_orb = parorb.get_double(9) ;
620 double y_orb = parorb.get_double(10) ;
621 double x_separ = parorb.get_double(11) ;
622 double y_separ = parorb.get_double(12) ;
623 double xbh_orb = parorb.get_double(13) ;
624 double mhsr = parorb.get_double(14) ;
625
626 double om2 = om * om ;
627 /*
628 double bpb = om2 * (x_orb * x_orb + y_orb * y_orb)
629 + 2. * om * (shifty * x_orb - shiftx * y_orb) + shift2
630 + 2. * mhsr * (x_separ*x_separ+y_separ*y_separ)
631 * (x_separ*shiftx + y_separ*shifty + om * y_separ * xbh_orb)
632 * (x_separ*shiftx + y_separ*shifty + om * y_separ * xbh_orb) ;
633 */
634
635 double bpb = om2 * (x_orb * x_orb + y_orb * y_orb
636 + 2.*mhsr*(x_separ*x_separ+y_separ*y_separ)
637 * y_separ * y_separ * xbh_orb * xbh_orb)
638 + 2.*om*(shifty * x_orb - shiftx * y_orb
639 + 2.*mhsr*(x_separ*x_separ+y_separ*y_separ)
640 * (shiftx*x_separ + shifty*y_separ) * y_separ * xbh_orb)
641 + shift2 + 2.*mhsr*(x_separ*x_separ+y_separ*y_separ)
642 * (shiftx*x_separ + shifty*y_separ)
643 * (shiftx*x_separ + shifty*y_separ) ;
644
645 double dlngam0 =
646 ( 0.5 * dp6sl2 * bpb
647 + p6sl2 * (0.5*dshift2
648 + om * (shifty - dshiftx*y_orb + dshifty*x_orb)
649 + om2 * x_orb
650 - mhsr * x_separ * (x_separ*shiftx+y_separ*shifty)
651 * (x_separ*shiftx+y_separ*shifty)
652 + 2. * mhsr * (x_separ*shiftx+y_separ*shifty)
653 * (y_separ*y_separ*shiftx - x_separ*y_separ*shifty
654 +(x_separ*x_separ+y_separ*y_separ)*(x_separ*dshiftx
655 +y_separ*dshifty))
656 + 2. * mhsr * om * y_separ * xbh_orb
657 * ( (y_separ*y_separ-2.*x_separ*x_separ)*shiftx
658 - 3. * x_separ * y_separ * shifty
659 +(x_separ*x_separ+y_separ*y_separ)*(x_separ*dshiftx
660 +y_separ*dshifty) )
661 - 3. * mhsr * om2 * x_separ * y_separ * y_separ * xbh_orb
662 * xbh_orb)
663 ) / (1 - p6sl2 * bpb) ;
664
665 return dnulg - dlngam0 ;
666
667}
668
669double func_binbhns_orbit_is(double om, const Param& parorb) {
670
671 double dnulg = parorb.get_double(0) ;
672 double p6sl2 = parorb.get_double(1) ;
673 double dp6sl2 = parorb.get_double(2) ;
674 double shiftx = parorb.get_double(3) ;
675 double shifty = parorb.get_double(4) ;
676 double dshiftx = parorb.get_double(5) ;
677 double dshifty = parorb.get_double(6) ;
678 double shift2 = parorb.get_double(7) ;
679 double dshift2 = parorb.get_double(8) ;
680 double x_orb = parorb.get_double(9) ;
681 double y_orb = parorb.get_double(10) ;
682
683 double om2 = om * om ;
684
685 double bpb = om2 * (x_orb * x_orb + y_orb * y_orb)
686 + 2. * om * (shifty * x_orb - shiftx * y_orb) + shift2 ;
687
688 double dlngam0 = ( 0.5 * dp6sl2 * bpb
689 + p6sl2 * (0.5*dshift2
690 + om *
691 (shifty - dshiftx*y_orb + dshifty*x_orb)
692 + om2 * x_orb)
693 ) / (1 - p6sl2 * bpb) ;
694
695 return dnulg - dlngam0 ;
696
697}
698}
double y_rot
Absolute Y coordinate of the rotation axis.
Definition bin_bhns.h:89
Hole_bhns hole
Black hole.
Definition bin_bhns.h:72
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition bin_bhns.h:80
double separ
Absolute orbital separation between two centers of BH and NS.
Definition bin_bhns.h:83
void orbit_omega(double fact_omeg_min, double fact_omeg_max)
Estimates the relative error on the Hamiltonian constraint.
double x_rot
Absolute X coordinate of the rotation axis.
Definition bin_bhns.h:86
Star_bhns star
Neutron star.
Definition bin_bhns.h:75
bool is_kerrschild() const
Returns true for a Kerr-Schild background, false for a Conformally flat one.
Definition blackhole.h:218
const Map & get_mp() const
Returns the mapping.
Definition blackhole.h:213
double get_mass_bh() const
Returns the gravitational mass of BH [{\tt m_unit}].
Definition blackhole.h:221
Base class for coordinate mappings.
Definition map.h:670
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Parameter storage.
Definition param.h:125
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
const Scalar & dsdx() const
Returns of *this , where .
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
Definition star_bhns.h:347
const Tensor & get_d_shift_comp() const
Returns the derivative of the shift vector generated by the companion black hole.
Definition star_bhns.h:380
const Scalar & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition star_bhns.h:321
const Scalar & get_confo_tot() const
Returns the total conformal factor.
Definition star_bhns.h:391
const Vector & get_shift_tot() const
Returns the part total shift vector.
Definition star_bhns.h:372
const Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the star.
Definition star_bhns.h:383
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion black hole.
Definition star_bhns.h:401
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the star.
Definition star_bhns.h:364
const Scalar & get_lapconf_auto() const
Returns the part of the lapconf function generated by the star.
Definition star_bhns.h:339
const Vector & get_d_lapconf_comp() const
Returns the derivative of the lapse function generated by the companion black hole.
Definition star_bhns.h:361
const Map & get_mp() const
Returns the mapping.
Definition star.h:355
Basic array class.
Definition tbl.h:161
int get_taille() const
Gives the total size (ie dim.taille)
Definition tbl.h:397
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
void zero_list(double(*f)(double, const Param &), const Param &par, double xmin, double xmax, int nsub, Tbl *&az, Tbl *&bz)
Locates approximatively all the zeros of a function in a given interval.
Definition zero_list.C:57
double zerosec_b(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter)
Finding the zero a function on a bounded domain.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.