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