LORENE
binaire_orbite.C
1 /*
2 * Method of class Binaire to compute the orbital angular velocity {\tt omega}
3 * and the position of the rotation axis {\tt x_axe}.
4 *
5 * (See file binaire.h for documentation)
6 *
7 */
8
9/*
10 * Copyright (c) 2000-2003 Eric Gourgoulhon
11 * Copyright (c) 2000-2001 Keisuke Taniguchi
12 *
13 * This file is part of LORENE.
14 *
15 * LORENE is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License as published by
17 * the Free Software Foundation; either version 2 of the License, or
18 * (at your option) any later version.
19 *
20 * LORENE is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with LORENE; if not, write to the Free Software
27 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28 *
29 */
30
31
32
33
34/*
35 * $Id: binaire_orbite.C,v 1.9 2016/12/05 16:17:47 j_novak Exp $
36 * $Log: binaire_orbite.C,v $
37 * Revision 1.9 2016/12/05 16:17:47 j_novak
38 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
39 *
40 * Revision 1.8 2014/10/24 11:27:49 j_novak
41 * Minor change in the setting of parameters for zero_list, to avoid problems with g++-4.8. To be further explored...
42 *
43 * Revision 1.7 2014/10/13 08:52:44 j_novak
44 * Lorene classes and functions now belong to the namespace Lorene.
45 *
46 * Revision 1.6 2014/10/06 15:12:59 j_novak
47 * Modified #include directives to use c++ syntax.
48 *
49 * Revision 1.5 2011/03/27 16:42:21 e_gourgoulhon
50 * Added output via new function save_profile for graphics.
51 *
52 * Revision 1.4 2009/06/18 18:40:57 k_taniguchi
53 * Added a slightly modified code to determine
54 * the orbital angular velocity.
55 *
56 * Revision 1.3 2004/03/25 10:28:59 j_novak
57 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
58 *
59 * Revision 1.2 2003/09/16 13:41:27 e_gourgoulhon
60 * Search of sub-intervals containing the zero(s) via zero_list
61 * Selection of the sub-interval as the closest one to previous value of omega.
62 * Replaced zerosec by zerosec_b.
63 *
64 * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
65 * LORENE
66 *
67 * Revision 2.9 2001/08/14 12:36:45 keisuke
68 * Change of the minimum and maximum values for searching a new position
69 * of the rotation axis.
70 *
71 * Revision 2.8 2001/03/20 16:06:55 keisuke
72 * Addition of the method directly to calculate
73 * the orbital angular velocity (we do not use this method!).
74 *
75 * Revision 2.7 2001/03/19 17:10:28 keisuke
76 * Change of the definition of the canter of mass.
77 *
78 * Revision 2.6 2001/03/19 13:37:04 keisuke
79 * Set x_axe to be zero for the case of identical stars.
80 *
81 * Revision 2.5 2001/03/17 16:17:20 keisuke
82 * Input a subroutine to determine the position of the rotation axis
83 * in the case of binary systems composed of different stars.
84 *
85 * Revision 2.4 2000/10/02 08:29:34 keisuke
86 * Change the method to construct d_logn_auto in the calculation of dnulg[i].
87 *
88 * Revision 2.3 2000/09/22 15:56:32 keisuke
89 * *** empty log message ***
90 *
91 * Revision 2.2 2000/09/22 15:52:12 keisuke
92 * Changement calcul de dlogn (prise en compte d'une partie divergente
93 * dans logn_auto par l'appel a d_logn_auto).
94 *
95 * Revision 2.1 2000/02/12 18:36:09 eric
96 * *** empty log message ***
97 *
98 * Revision 2.0 2000/02/12 17:09:21 eric
99 * *** empty log message ***
100 *
101 *
102 * $Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_orbite.C,v 1.9 2016/12/05 16:17:47 j_novak Exp $
103 *
104 */
105
106// Headers C
107#include <cmath>
108
109// Headers Lorene
110#include "binaire.h"
111#include "eos.h"
112#include "param.h"
113#include "utilitaires.h"
114#include "unites.h"
115
116#include "scalar.h"
117#include "graphique.h"
118
119namespace Lorene {
120double fonc_binaire_axe(double , const Param& ) ;
121double fonc_binaire_orbit(double , const Param& ) ;
122
123//******************************************************************************
124
125void Binaire::orbit(double fact_omeg_min, double fact_omeg_max, double& xgg1,
126 double& xgg2) {
127
128 using namespace Unites ;
129
130 //-------------------------------------------------------------
131 // Evaluation of various quantities at the center of each star
132 //-------------------------------------------------------------
133
134 double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
135 double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
136
137 for (int i=0; i<2; i++) {
138
139 const Map& mp = et[i]->get_mp() ;
140
141 const Cmp& logn_auto_regu = et[i]->get_logn_auto_regu()() ;
142 const Cmp& logn_comp = et[i]->get_logn_comp()() ;
143 const Cmp& loggam = et[i]->get_loggam()() ;
144 const Cmp& nnn = et[i]->get_nnn()() ;
145 const Cmp& a_car = et[i]->get_a_car()() ;
146 const Tenseur& shift = et[i]->get_shift() ;
147 const Tenseur& d_logn_auto_div = et[i]->get_d_logn_auto_div() ;
148
149 Tenseur dln_auto_div = d_logn_auto_div ;
150
151 if ( *(dln_auto_div.get_triad()) != ref_triad ) {
152
153 // Change the basis from spherical coordinate to Cartesian one
154 dln_auto_div.change_triad( mp.get_bvect_cart() ) ;
155
156 // Change the basis from mapping coordinate to absolute one
157 dln_auto_div.change_triad( ref_triad ) ;
158
159 }
160
161 //----------------------------------
162 // Calcul de d/dX( nu + ln(Gamma) ) au centre de l'etoile ---> dnulg[i]
163 //----------------------------------
164
165 // Facteur de passage x --> X :
166 double factx ;
167 if (fabs(mp.get_rot_phi()) < 1.e-14) {
168 factx = 1. ;
169 }
170 else {
171 if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
172 factx = - 1. ;
173 }
174 else {
175 cout << "Binaire::orbit : unknown value of rot_phi !" << endl ;
176 abort() ;
177 }
178 }
179
180 Cmp tmp = logn_auto_regu + logn_comp + loggam ;
181
182 // ... gradient suivant X :
183 dnulg[i] = dln_auto_div(0)(0, 0, 0, 0)
184 + factx * tmp.dsdx()(0, 0, 0, 0) ;
185
186 tmp = logn_auto_regu + logn_comp ;
187 cout << "dlnndx_div_c : " << dln_auto_div(0)(0, 0, 0, 0) << endl ;
188 cout << "dlnndx_c : " << dln_auto_div(0)(0, 0, 0, 0) + factx*tmp.dsdx()(0, 0, 0, 0) << endl ;
189
190 cout << "dloggamdx_c : " << factx*loggam.dsdx()(0, 0, 0, 0) << endl ;
191 Scalar stmp(logn_comp) ;
192 save_profile(stmp, 0., 10., 0.5*M_PI, 0., "prof_logn.d") ;
193 stmp = loggam ;
194 save_profile(stmp, 0., 1.8, 0.5*M_PI, 0., "prof_loggam.d") ;
195
196 //----------------------------------
197 // Calcul de A^2/N^2 au centre de l'etoile ---> asn2[i]
198 //----------------------------------
199
200 double nc = nnn(0, 0, 0, 0) ;
201 double a2c = a_car(0, 0, 0, 0) ;
202 asn2[i] = a2c / (nc * nc) ;
203
204 if ( et[i]->is_relativistic() ) {
205
206 //----------------------------------
207 // Calcul de d/dX(A^2/N^2) au centre de l'etoile ---> dasn2[i]
208 //----------------------------------
209
210 double da2c = factx * a_car.dsdx()(0, 0, 0, 0) ;
211 double dnc = factx * nnn.dsdx()(0, 0, 0, 0) ;
212
213 dasn2[i] = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
214
215 //----------------------------------
216 // Calcul de N^Y au centre de l'etoile ---> ny[i]
217 //----------------------------------
218
219 ny[i] = shift(1)(0, 0, 0, 0) ;
220 nyso[i] = ny[i] / omega ;
221
222 //----------------------------------
223 // Calcul de dN^Y/dX au centre de l'etoile ---> dny[i]
224 //----------------------------------
225
226 dny[i] = factx * shift(1).dsdx()(0, 0, 0, 0) ;
227 dnyso[i] = dny[i] / omega ;
228
229 //----------------------------------
230 // Calcul de (N^X)^2 + (N^Y)^2 + (N^Z)^2
231 // au centre de l'etoile ---> npn[i]
232 //----------------------------------
233
234 tmp = flat_scalar_prod(shift, shift)() ;
235
236 npn[i] = tmp(0, 0, 0, 0) ;
237 npnso2[i] = npn[i] / omega / omega ;
238
239 //----------------------------------
240 // Calcul de d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
241 // au centre de l'etoile ---> dnpn[i]
242 //----------------------------------
243
244 dnpn[i] = factx * tmp.dsdx()(0, 0, 0, 0) ;
245 dnpnso2[i] = dnpn[i] / omega / omega ;
246
247 } // fin du cas relativiste
248 else {
249 dasn2[i] = 0 ;
250 ny[i] = 0 ;
251 nyso[i] = 0 ;
252 dny[i] = 0 ;
253 dnyso[i] = 0 ;
254 npn[i] = 0 ;
255 npnso2[i] = 0 ;
256 dnpn[i] = 0 ;
257 dnpnso2[i] = 0 ;
258 }
259
260 cout << "Binaire::orbit: central d(nu+log(Gam))/dX : "
261 << dnulg[i] << endl ;
262 cout << "Binaire::orbit: central A^2/N^2 : " << asn2[i] << endl ;
263 cout << "Binaire::orbit: central d(A^2/N^2)/dX : " << dasn2[i] << endl ;
264 cout << "Binaire::orbit: central N^Y : " << ny[i] << endl ;
265 cout << "Binaire::orbit: central dN^Y/dX : " << dny[i] << endl ;
266 cout << "Binaire::orbit: central N.N : " << npn[i] << endl ;
267 cout << "Binaire::orbit: central d(N.N)/dX : " << dnpn[i] << endl ;
268
269 //----------------------
270 // Pour information seulement : 1/ calcul des positions des "centres de
271 // de masse"
272 // 2/ calcul de dH/dX en r=0
273 //-----------------------
274
275 ori_x[i] = (et[i]->get_mp()).get_ori_x() ;
276
277 xgg[i] = factx * (et[i]->xa_barycenter() - ori_x[i]) ;
278
279 } // fin de la boucle sur les etoiles
280
281 xgg1 = xgg[0] ;
282 xgg2 = xgg[1] ;
283
284//---------------------------------
285// Position de l'axe de rotation
286//---------------------------------
287
288 int relat = ( et[0]->is_relativistic() ) ? 1 : 0 ;
289 double ori_x1 = ori_x[0] ;
290 double ori_x2 = ori_x[1] ;
291
292 if ( et[0]->get_eos() == et[1]->get_eos() &&
293 et[0]->get_ent()()(0,0,0,0) == et[1]->get_ent()()(0,0,0,0) ) {
294
295 x_axe = 0. ;
296
297 }
298 else {
299
300 Param paraxe ;
301 paraxe.add_int(relat) ;
302 paraxe.add_double( ori_x1, 0) ;
303 paraxe.add_double( ori_x2, 1) ;
304 paraxe.add_double( dnulg[0], 2) ;
305 paraxe.add_double( dnulg[1], 3) ;
306 paraxe.add_double( asn2[0], 4) ;
307 paraxe.add_double( asn2[1], 5) ;
308 paraxe.add_double( dasn2[0], 6) ;
309 paraxe.add_double( dasn2[1], 7) ;
310 paraxe.add_double( nyso[0], 8) ;
311 paraxe.add_double( nyso[1], 9) ;
312 paraxe.add_double( dnyso[0], 10) ;
313 paraxe.add_double( dnyso[1], 11) ;
314 paraxe.add_double( npnso2[0], 12) ;
315 paraxe.add_double( npnso2[1], 13) ;
316 paraxe.add_double( dnpnso2[0], 14) ;
317 paraxe.add_double( dnpnso2[1], 15) ;
318
319 int nitmax_axe = 200 ;
320 int nit_axe ;
321 double precis_axe = 1.e-13 ;
322
323 x_axe = zerosec(fonc_binaire_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
324 precis_axe, nitmax_axe, nit_axe) ;
325
326 cout << "Binaire::orbit : Number of iterations in zerosec for x_axe : "
327 << nit_axe << endl ;
328 }
329
330 cout << "Binaire::orbit : x_axe [km] : " << x_axe / km << endl ;
331
332//-------------------------------------
333// Calcul de la vitesse orbitale
334//-------------------------------------
335
336 Param parf ;
337 parf.add_int(relat) ;
338 parf.add_double( ori_x1, 0) ;
339 parf.add_double( dnulg[0], 1) ;
340 parf.add_double( asn2[0], 2) ;
341 parf.add_double( dasn2[0], 3) ;
342 parf.add_double( ny[0], 4) ;
343 parf.add_double( dny[0], 5) ;
344 parf.add_double( npn[0], 6) ;
345 parf.add_double( dnpn[0], 7) ;
346 parf.add_double( x_axe, 8) ;
347
348 double omega1 = fact_omeg_min * omega ;
349 double omega2 = fact_omeg_max * omega ;
350 cout << "Binaire::orbit: omega1, omega2 [rad/s] : "
351 << omega1 * f_unit << " " << omega2 * f_unit << endl ;
352
353 // Search for the various zeros in the interval [omega1,omega2]
354 // ------------------------------------------------------------
355 int nsub = 50 ; // total number of subdivisions of the interval
356 Tbl* azer = 0x0 ;
357 Tbl* bzer = 0x0 ;
358 zero_list(fonc_binaire_orbit, parf, omega1, omega2, nsub,
359 azer, bzer) ;
360
361 // Search for the zero closest to the previous value of omega
362 // ----------------------------------------------------------
363 double omeg_min, omeg_max ;
364 int nzer = azer->get_taille() ; // number of zeros found by zero_list
365 cout << "Binaire:orbit : " << nzer <<
366 " zero(s) found in the interval [omega1, omega2]." << endl ;
367 cout << "omega, omega1, omega2 : " << omega << " " << omega1
368 << " " << omega2 << endl ;
369 cout << "azer : " << *azer << endl ;
370 cout << "bzer : " << *bzer << endl ;
371
372 if (nzer == 0) {
373 cout <<
374 "Binaire::orbit: WARNING : no zero detected in the interval"
375 << endl << " [" << omega1 * f_unit << ", "
376 << omega2 * f_unit << "] rad/s !" << endl ;
377 omeg_min = omega1 ;
378 omeg_max = omega2 ;
379 }
380 else {
381 double dist_min = fabs(omega2 - omega1) ;
382 int i_dist_min = -1 ;
383 for (int i=0; i<nzer; i++) {
384 // Distance of previous value of omega from the center of the
385 // interval [azer(i), bzer(i)]
386 double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
387 if (dist < dist_min) {
388 dist_min = dist ;
389 i_dist_min = i ;
390 }
391 }
392 omeg_min = (*azer)(i_dist_min) ;
393 omeg_max = (*bzer)(i_dist_min) ;
394 }
395
396 delete azer ; // Tbl allocated by zero_list
397 delete bzer ; //
398
399 cout << "Binaire:orbit : interval selected for the search of the zero : "
400 << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
401 << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
402
403 // Computation of the zero in the selected interval by the secant method
404 // ---------------------------------------------------------------------
405 int nitermax = 200 ;
406 int niter ;
407 double precis = 1.e-13 ;
408 omega = zerosec_b(fonc_binaire_orbit, parf, omeg_min, omeg_max,
409 precis, nitermax, niter) ;
410
411 cout << "Binaire::orbit : Number of iterations in zerosec for omega : "
412 << niter << endl ;
413
414 cout << "Binaire::orbit : omega [rad/s] : "
415 << omega * f_unit << endl ;
416
417
418 /* We do not use the method below.
419 // Direct calculation
420 //--------------------
421
422 double om2_et1 ;
423 double om2_et2 ;
424
425 double x_et1 = ori_x1 - x_axe ;
426 double x_et2 = ori_x2 - x_axe ;
427
428 if (relat == 1) {
429
430 double andan_et1 = 0.5 * dasn2[0] + asn2[0] * dnulg[0] ;
431 double andan_et2 = 0.5 * dasn2[1] + asn2[1] * dnulg[1] ;
432
433 double bpb_et1 = x_et1 * x_et1 - 2. * nyso[0] * x_et1 + npnso2[0] ;
434 double bpb_et2 = x_et2 * x_et2 - 2. * nyso[1] * x_et2 + npnso2[1] ;
435
436 double cpc_et1 = 0.5 * dnpnso2[0] + x_et1 * (1. - dnyso[0]) - nyso[0] ;
437 double cpc_et2 = 0.5 * dnpnso2[1] + x_et2 * (1. - dnyso[1]) - nyso[1] ;
438
439 om2_et1 = dnulg[0] / (andan_et1 * bpb_et1 + asn2[0] * cpc_et1) ;
440 om2_et2 = dnulg[1] / (andan_et2 * bpb_et2 + asn2[1] * cpc_et2) ;
441
442 }
443 else {
444
445 om2_et1 = dnulg[0] / x_et1 ;
446 om2_et2 = dnulg[1] / x_et2 ;
447
448 }
449
450 double ome_et1 = sqrt( om2_et1 ) ;
451 double ome_et2 = sqrt( om2_et2 ) ;
452
453 double diff_om1 = 1. - ome_et1 / ome_et2 ;
454 double diff_om2 = 1. - ome_et1 / omega ;
455
456 cout << "Binaire::orbit : omega (direct) [rad/s] : "
457 << ome_et1 * f_unit << endl ;
458
459 cout << "Binaire::orbit : relative difference between " << endl
460 << " omega_star1 and omega_star2 (direct) : " << diff_om1 << endl
461 << " omega (direct) and omega (iretation) : " << diff_om2 << endl ;
462 */
463
464}
465
466
467void Binaire::orbit_eqmass(double fact_omeg_min, double fact_omeg_max,
468 double mass1, double mass2,
469 double& xgg1, double& xgg2) {
470
471 using namespace Unites ;
472
473 //-------------------------------------------------------------
474 // Evaluation of various quantities at the center of each star
475 //-------------------------------------------------------------
476
477 double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
478 double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
479
480 for (int i=0; i<2; i++) {
481
482 const Map& mp = et[i]->get_mp() ;
483
484 const Cmp& logn_auto_regu = et[i]->get_logn_auto_regu()() ;
485 const Cmp& logn_comp = et[i]->get_logn_comp()() ;
486 const Cmp& loggam = et[i]->get_loggam()() ;
487 const Cmp& nnn = et[i]->get_nnn()() ;
488 const Cmp& a_car = et[i]->get_a_car()() ;
489 const Tenseur& shift = et[i]->get_shift() ;
490 const Tenseur& d_logn_auto_div = et[i]->get_d_logn_auto_div() ;
491
492 Tenseur dln_auto_div = d_logn_auto_div ;
493
494 if ( *(dln_auto_div.get_triad()) != ref_triad ) {
495
496 // Change the basis from spherical coordinate to Cartesian one
497 dln_auto_div.change_triad( mp.get_bvect_cart() ) ;
498
499 // Change the basis from mapping coordinate to absolute one
500 dln_auto_div.change_triad( ref_triad ) ;
501
502 }
503
504 //----------------------------------
505 // Calcul de d/dX( nu + ln(Gamma) ) au centre de l'etoile ---> dnulg[i]
506 //----------------------------------
507
508 // Facteur de passage x --> X :
509 double factx ;
510 if (fabs(mp.get_rot_phi()) < 1.e-14) {
511 factx = 1. ;
512 }
513 else {
514 if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
515 factx = - 1. ;
516 }
517 else {
518 cout << "Binaire::orbit : unknown value of rot_phi !" << endl ;
519 abort() ;
520 }
521 }
522
523 Cmp tmp = logn_auto_regu + logn_comp + loggam ;
524
525 // ... gradient suivant X :
526 dnulg[i] = dln_auto_div(0)(0, 0, 0, 0)
527 + factx * tmp.dsdx()(0, 0, 0, 0) ;
528
529 //----------------------------------
530 // Calcul de A^2/N^2 au centre de l'etoile ---> asn2[i]
531 //----------------------------------
532
533 double nc = nnn(0, 0, 0, 0) ;
534 double a2c = a_car(0, 0, 0, 0) ;
535 asn2[i] = a2c / (nc * nc) ;
536
537 if ( et[i]->is_relativistic() ) {
538
539 //----------------------------------
540 // Calcul de d/dX(A^2/N^2) au centre de l'etoile ---> dasn2[i]
541 //----------------------------------
542
543 double da2c = factx * a_car.dsdx()(0, 0, 0, 0) ;
544 double dnc = factx * nnn.dsdx()(0, 0, 0, 0) ;
545
546 dasn2[i] = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
547
548 //----------------------------------
549 // Calcul de N^Y au centre de l'etoile ---> ny[i]
550 //----------------------------------
551
552 ny[i] = shift(1)(0, 0, 0, 0) ;
553 nyso[i] = ny[i] / omega ;
554
555 //----------------------------------
556 // Calcul de dN^Y/dX au centre de l'etoile ---> dny[i]
557 //----------------------------------
558
559 dny[i] = factx * shift(1).dsdx()(0, 0, 0, 0) ;
560 dnyso[i] = dny[i] / omega ;
561
562 //----------------------------------
563 // Calcul de (N^X)^2 + (N^Y)^2 + (N^Z)^2
564 // au centre de l'etoile ---> npn[i]
565 //----------------------------------
566
567 tmp = flat_scalar_prod(shift, shift)() ;
568
569 npn[i] = tmp(0, 0, 0, 0) ;
570 npnso2[i] = npn[i] / omega / omega ;
571
572 //----------------------------------
573 // Calcul de d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
574 // au centre de l'etoile ---> dnpn[i]
575 //----------------------------------
576
577 dnpn[i] = factx * tmp.dsdx()(0, 0, 0, 0) ;
578 dnpnso2[i] = dnpn[i] / omega / omega ;
579
580 } // fin du cas relativiste
581 else {
582 dasn2[i] = 0 ;
583 ny[i] = 0 ;
584 nyso[i] = 0 ;
585 dny[i] = 0 ;
586 dnyso[i] = 0 ;
587 npn[i] = 0 ;
588 npnso2[i] = 0 ;
589 dnpn[i] = 0 ;
590 dnpnso2[i] = 0 ;
591 }
592
593 cout << "Binaire::orbit: central d(nu+log(Gam))/dX : "
594 << dnulg[i] << endl ;
595 cout << "Binaire::orbit: central A^2/N^2 : " << asn2[i] << endl ;
596 cout << "Binaire::orbit: central d(A^2/N^2)/dX : " << dasn2[i] << endl ;
597 cout << "Binaire::orbit: central N^Y : " << ny[i] << endl ;
598 cout << "Binaire::orbit: central dN^Y/dX : " << dny[i] << endl ;
599 cout << "Binaire::orbit: central N.N : " << npn[i] << endl ;
600 cout << "Binaire::orbit: central d(N.N)/dX : " << dnpn[i] << endl ;
601
602 //----------------------
603 // Pour information seulement : 1/ calcul des positions des "centres de
604 // de masse"
605 // 2/ calcul de dH/dX en r=0
606 //-----------------------
607
608 ori_x[i] = (et[i]->get_mp()).get_ori_x() ;
609
610 xgg[i] = factx * (et[i]->xa_barycenter() - ori_x[i]) ;
611
612 } // fin de la boucle sur les etoiles
613
614 xgg1 = xgg[0] ;
615 xgg2 = xgg[1] ;
616
617//---------------------------------
618// Position de l'axe de rotation
619//---------------------------------
620
621 int relat = ( et[0]->is_relativistic() ) ? 1 : 0 ;
622 double ori_x1 = ori_x[0] ;
623 double ori_x2 = ori_x[1] ;
624
625 if ( et[0]->get_eos() == et[1]->get_eos() && mass1 == mass2 ) {
626
627 x_axe = 0. ;
628
629 }
630 else {
631
632 Param paraxe ;
633 paraxe.add_int(relat) ;
634 paraxe.add_double( ori_x1, 0) ;
635 paraxe.add_double( ori_x2, 1) ;
636 paraxe.add_double( dnulg[0], 2) ;
637 paraxe.add_double( dnulg[1], 3) ;
638 paraxe.add_double( asn2[0], 4) ;
639 paraxe.add_double( asn2[1], 5) ;
640 paraxe.add_double( dasn2[0], 6) ;
641 paraxe.add_double( dasn2[1], 7) ;
642 paraxe.add_double( nyso[0], 8) ;
643 paraxe.add_double( nyso[1], 9) ;
644 paraxe.add_double( dnyso[0], 10) ;
645 paraxe.add_double( dnyso[1], 11) ;
646 paraxe.add_double( npnso2[0], 12) ;
647 paraxe.add_double( npnso2[1], 13) ;
648 paraxe.add_double( dnpnso2[0], 14) ;
649 paraxe.add_double( dnpnso2[1], 15) ;
650
651 int nitmax_axe = 200 ;
652 int nit_axe ;
653 double precis_axe = 1.e-13 ;
654
655 x_axe = zerosec(fonc_binaire_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
656 precis_axe, nitmax_axe, nit_axe) ;
657
658 cout << "Binaire::orbit : Number of iterations in zerosec for x_axe : "
659 << nit_axe << endl ;
660 }
661
662 cout << "Binaire::orbit : x_axe [km] : " << x_axe / km << endl ;
663
664//-------------------------------------
665// Calcul de la vitesse orbitale
666//-------------------------------------
667
668 Param parf ;
669 parf.add_int(relat) ;
670 parf.add_double( ori_x1, 0) ;
671 parf.add_double( dnulg[0], 1) ;
672 parf.add_double( asn2[0], 2) ;
673 parf.add_double( dasn2[0], 3) ;
674 parf.add_double( ny[0], 4) ;
675 parf.add_double( dny[0], 5) ;
676 parf.add_double( npn[0], 6) ;
677 parf.add_double( dnpn[0], 7) ;
678 parf.add_double( x_axe, 8) ;
679
680 double omega1 = fact_omeg_min * omega ;
681 double omega2 = fact_omeg_max * omega ;
682 cout << "Binaire::orbit: omega1, omega2 [rad/s] : "
683 << omega1 * f_unit << " " << omega2 * f_unit << endl ;
684
685 // Search for the various zeros in the interval [omega1,omega2]
686 // ------------------------------------------------------------
687 int nsub = 50 ; // total number of subdivisions of the interval
688 Tbl* azer = 0x0 ;
689 Tbl* bzer = 0x0 ;
690 zero_list(fonc_binaire_orbit, parf, omega1, omega2, nsub,
691 azer, bzer) ;
692
693 // Search for the zero closest to the previous value of omega
694 // ----------------------------------------------------------
695 double omeg_min, omeg_max ;
696 int nzer = azer->get_taille() ; // number of zeros found by zero_list
697 cout << "Binaire:orbit : " << nzer <<
698 " zero(s) found in the interval [omega1, omega2]." << endl ;
699 cout << "omega, omega1, omega2 : " << omega << " " << omega1
700 << " " << omega2 << endl ;
701 cout << "azer : " << *azer << endl ;
702 cout << "bzer : " << *bzer << endl ;
703
704 if (nzer == 0) {
705 cout <<
706 "Binaire::orbit: WARNING : no zero detected in the interval"
707 << endl << " [" << omega1 * f_unit << ", "
708 << omega2 * f_unit << "] rad/s !" << endl ;
709 omeg_min = omega1 ;
710 omeg_max = omega2 ;
711 }
712 else {
713 double dist_min = fabs(omega2 - omega1) ;
714 int i_dist_min = -1 ;
715 for (int i=0; i<nzer; i++) {
716 // Distance of previous value of omega from the center of the
717 // interval [azer(i), bzer(i)]
718 double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
719 if (dist < dist_min) {
720 dist_min = dist ;
721 i_dist_min = i ;
722 }
723 }
724 omeg_min = (*azer)(i_dist_min) ;
725 omeg_max = (*bzer)(i_dist_min) ;
726 }
727
728 delete azer ; // Tbl allocated by zero_list
729 delete bzer ; //
730
731 cout << "Binaire:orbit : interval selected for the search of the zero : "
732 << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
733 << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
734
735 // Computation of the zero in the selected interval by the secant method
736 // ---------------------------------------------------------------------
737 int nitermax = 200 ;
738 int niter ;
739 double precis = 1.e-13 ;
740 omega = zerosec_b(fonc_binaire_orbit, parf, omeg_min, omeg_max,
741 precis, nitermax, niter) ;
742
743 cout << "Binaire::orbit : Number of iterations in zerosec for omega : "
744 << niter << endl ;
745
746 cout << "Binaire::orbit : omega [rad/s] : "
747 << omega * f_unit << endl ;
748
749
750 /* We do not use the method below.
751 // Direct calculation
752 //--------------------
753
754 double om2_et1 ;
755 double om2_et2 ;
756
757 double x_et1 = ori_x1 - x_axe ;
758 double x_et2 = ori_x2 - x_axe ;
759
760 if (relat == 1) {
761
762 double andan_et1 = 0.5 * dasn2[0] + asn2[0] * dnulg[0] ;
763 double andan_et2 = 0.5 * dasn2[1] + asn2[1] * dnulg[1] ;
764
765 double bpb_et1 = x_et1 * x_et1 - 2. * nyso[0] * x_et1 + npnso2[0] ;
766 double bpb_et2 = x_et2 * x_et2 - 2. * nyso[1] * x_et2 + npnso2[1] ;
767
768 double cpc_et1 = 0.5 * dnpnso2[0] + x_et1 * (1. - dnyso[0]) - nyso[0] ;
769 double cpc_et2 = 0.5 * dnpnso2[1] + x_et2 * (1. - dnyso[1]) - nyso[1] ;
770
771 om2_et1 = dnulg[0] / (andan_et1 * bpb_et1 + asn2[0] * cpc_et1) ;
772 om2_et2 = dnulg[1] / (andan_et2 * bpb_et2 + asn2[1] * cpc_et2) ;
773
774 }
775 else {
776
777 om2_et1 = dnulg[0] / x_et1 ;
778 om2_et2 = dnulg[1] / x_et2 ;
779
780 }
781
782 double ome_et1 = sqrt( om2_et1 ) ;
783 double ome_et2 = sqrt( om2_et2 ) ;
784
785 double diff_om1 = 1. - ome_et1 / ome_et2 ;
786 double diff_om2 = 1. - ome_et1 / omega ;
787
788 cout << "Binaire::orbit : omega (direct) [rad/s] : "
789 << ome_et1 * f_unit << endl ;
790
791 cout << "Binaire::orbit : relative difference between " << endl
792 << " omega_star1 and omega_star2 (direct) : " << diff_om1 << endl
793 << " omega (direct) and omega (iretation) : " << diff_om2 << endl ;
794 */
795
796}
797
798
799//*************************************************
800// Function used for search of the rotation axis
801//*************************************************
802
803double fonc_binaire_axe(double x_rot, const Param& paraxe) {
804
805 int relat = paraxe.get_int() ;
806
807 double ori_x1 = paraxe.get_double(0) ;
808 double ori_x2 = paraxe.get_double(1) ;
809 double dnulg_1 = paraxe.get_double(2) ;
810 double dnulg_2 = paraxe.get_double(3) ;
811 double asn2_1 = paraxe.get_double(4) ;
812 double asn2_2 = paraxe.get_double(5) ;
813 double dasn2_1 = paraxe.get_double(6) ;
814 double dasn2_2 = paraxe.get_double(7) ;
815 double nyso_1 = paraxe.get_double(8) ;
816 double nyso_2 = paraxe.get_double(9) ;
817 double dnyso_1 = paraxe.get_double(10) ;
818 double dnyso_2 = paraxe.get_double(11) ;
819 double npnso2_1 = paraxe.get_double(12) ;
820 double npnso2_2 = paraxe.get_double(13) ;
821 double dnpnso2_1 = paraxe.get_double(14) ;
822 double dnpnso2_2 = paraxe.get_double(15) ;
823
824 double om2_star1 ;
825 double om2_star2 ;
826
827 double x1 = ori_x1 - x_rot ;
828 double x2 = ori_x2 - x_rot ;
829
830 if (relat == 1) {
831
832 double andan_1 = 0.5 * dasn2_1 + asn2_1 * dnulg_1 ;
833 double andan_2 = 0.5 * dasn2_2 + asn2_2 * dnulg_2 ;
834
835 double bpb_1 = x1 * x1 - 2. * nyso_1 * x1 + npnso2_1 ;
836 double bpb_2 = x2 * x2 - 2. * nyso_2 * x2 + npnso2_2 ;
837
838 double cpc_1 = 0.5 * dnpnso2_1 + x1 * (1. - dnyso_1) - nyso_1 ;
839 double cpc_2 = 0.5 * dnpnso2_2 + x2 * (1. - dnyso_2) - nyso_2 ;
840
841 om2_star1 = dnulg_1 / (andan_1 * bpb_1 + asn2_1 * cpc_1) ;
842 om2_star2 = dnulg_2 / (andan_2 * bpb_2 + asn2_2 * cpc_2) ;
843
844 }
845 else {
846
847 om2_star1 = dnulg_1 / x1 ;
848 om2_star2 = dnulg_2 / x2 ;
849
850 }
851
852 return om2_star1 - om2_star2 ;
853
854}
855
856//*****************************************************************************
857// Fonction utilisee pour la recherche de omega par la methode de la secante
858//*****************************************************************************
859
860double fonc_binaire_orbit(double om, const Param& parf) {
861
862 int relat = parf.get_int() ;
863
864 double xc = parf.get_double(0) ;
865 double dnulg = parf.get_double(1) ;
866 double asn2 = parf.get_double(2) ;
867 double dasn2 = parf.get_double(3) ;
868 double ny = parf.get_double(4) ;
869 double dny = parf.get_double(5) ;
870 double npn = parf.get_double(6) ;
871 double dnpn = parf.get_double(7) ;
872 double x_axe = parf.get_double(8) ;
873
874 double xx = xc - x_axe ;
875 double om2 = om*om ;
876
877 double dphi_cent ;
878
879 if (relat == 1) {
880 double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
881
882 dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
883 - 0.5*bpb* dasn2 )
884 / ( 1 - asn2 * bpb ) ;
885 }
886 else {
887 dphi_cent = - om2*xx ;
888 }
889
890 return dnulg + dphi_cent ;
891
892}
893
894
895}
Lorene::Binaire::et
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): {\tt et[0]} contains the address of {\tt star...
Definition binaire.h:127
Lorene::Binaire::orbit_eqmass
void orbit_eqmass(double fact_omeg_min, double fact_omeg_max, double mass1, double mass2, double &xgg1, double &xgg2)
Computes the orbital angular velocity {\tt omega} and the position of the rotation axis {\tt x_axe}...
Definition binaire_orbite.C:467
Lorene::Binaire::ref_triad
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Definition binaire.h:115
Lorene::Binaire::omega
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition binaire.h:132
Lorene::Binaire::orbit
void orbit(double fact_omeg_min, double fact_omeg_max, double &xgg1, double &xgg2)
Computes the orbital angular velocity {\tt omega} and the position of the rotation axis {\tt x_axe}...
Definition binaire_orbite.C:125
Lorene::Binaire::x_axe
double x_axe
Absolute X coordinate of the rotation axis.
Definition binaire.h:136
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::dsdx
const Cmp & dsdx() const
Returns of *this , where .
Definition cmp_deriv.C:151
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_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:295
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::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:249
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
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::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:304
Lorene::Tenseur::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:707
Lorene::Tenseur::change_triad
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tenseur.C:674
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::zerosec
double zerosec(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter, bool abort=true)
Finding the zero a function.
Definition zerosec.C:92
Lorene::flat_scalar_prod
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Definition tenseur_operateur.C:656
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.