LORENE
binary_orbit_xcts.C
1/*
2 * Method of class Binary_xcts to compute the orbital angular velocity
3 * {\tt omega} and the position of the rotation axis {\tt x_axe}.
4 * (See file binary_xcts.h for documentation)
5 */
6
7/*
8 * Copyright (c) 2010 Michal Bejger
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License version 2
14 * as published by the Free Software Foundation.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27char binary_orbit_xcts_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_orbit_xcts.C,v 1.14 2014/10/24 14:10:24 j_novak Exp $" ;
28
29/*
30 * $Id: binary_orbit_xcts.C,v 1.14 2014/10/24 14:10:24 j_novak Exp $
31 * $Log: binary_orbit_xcts.C,v $
32 * Revision 1.14 2014/10/24 14:10:24 j_novak
33 * Minor change to prevent weird error from g++-4.8...
34 *
35 * Revision 1.13 2014/10/13 08:52:45 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.12 2014/10/06 15:12:59 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.11 2011/03/30 13:14:27 m_bejger
42 * Psi and chi rewritten using auto and comp parts to improve the convergence (in all the remaining fields, not only logn)
43 *
44 * Revision 1.10 2011/03/28 17:13:37 m_bejger
45 * logn in dnulg stated using Psi1,2 and chi1,2
46 *
47 * Revision 1.9 2011/03/27 16:41:19 e_gourgoulhon
48 * -- Corrected sign of ny and dny for star no. 2
49 * -- Added output via new function save_profile for graphics
50 *
51 * Revision 1.8 2011/03/27 14:58:48 m_bejger
52 * dnulg by means of dsdx(); rearrangements to use primary variables
53 *
54 * Revision 1.7 2011/03/25 16:28:36 e_gourgoulhon
55 * Still in progress
56 *
57 * Revision 1.6 2010/12/09 10:41:20 m_bejger
58 * For testing; not sure if working properly
59 *
60 * Revision 1.5 2010/10/26 19:45:45 m_bejger
61 * Cleanup
62 *
63 * Revision 1.4 2010/07/16 16:27:19 m_bejger
64 * This version is basically a copy of the one used by Binaire (binaire_orbite.C)
65 *
66 * Revision 1.3 2010/06/17 14:15:41 m_bejger
67 * Using method get_Psi()
68 *
69 * Revision 1.2 2010/06/15 07:57:30 m_bejger
70 * Minor corrections
71 *
72 * Revision 1.1 2010/05/04 07:35:54 m_bejger
73 * Initial version
74 *
75 * $Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_orbit_xcts.C,v 1.14 2014/10/24 14:10:24 j_novak Exp $
76 *
77 */
78
79// Headers C
80#include <cmath>
81
82// Headers Lorene
83#include "binary_xcts.h"
84#include "eos.h"
85#include "param.h"
86#include "utilitaires.h"
87#include "unites.h"
88
89#include "graphique.h"
90
91namespace Lorene {
92double fonc_binary_xcts_axe(double , const Param& ) ;
93double fonc_binary_xcts_orbit(double , const Param& ) ;
94
95//******************************************************************************
96
97void Binary_xcts::orbit(double fact_omeg_min,
98 double fact_omeg_max,
99 double& xgg1,
100 double& xgg2) {
101
102 using namespace Unites ;
103
104 //-------------------------------------------------------------
105 // Evaluation of various quantities at the center of each star
106 //-------------------------------------------------------------
107
108 double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
109 double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
110
111 for (int i=0; i<2; i++) {
112
113 const Map& mp = et[i]->get_mp() ;
114 const Metric& flat = et[i]->get_flat() ;
115
116 //------------------------------------------------------------------
117 // Recasting Phi and chi to manifestly equal auto and comp part
118 // - more fortunate from the point of view of Omega computation
119 //------------------------------------------------------------------
120 Scalar chi = et[i]->get_chi_auto() + et[i]->get_chi_comp() + 1 ;
121 Scalar Psi = et[i]->get_Psi_auto() + et[i]->get_Psi_comp() + 1 ;
122
123 Scalar logn = log(chi) - log(Psi) ;
124 logn.std_spectral_base() ;
125
126 // Sign convention for shift (beta^i = - N^i)
127 Vector shift = - ( et[i]->get_beta() ) ;
128 shift.change_triad(et[i]->mp.get_bvect_cart()) ;
129
130 //------------------------------------------------------------------
131 // d/dX( log(N) + log(Gamma) )
132 // in the center of the star ---> dnulg[i]
133 //------------------------------------------------------------------
134
135 // Facteur de passage x --> X :
136 double factx ;
137
138 if (fabs(mp.get_rot_phi()) < 1.e-14) factx = 1. ;
139 else {
140
141 if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
142 factx = - 1. ;
143
144 } else {
145
146 cout << "Binary_xcts::orbit : unknown value of rot_phi !" << endl ;
147 abort() ;
148 }
149 }
150
151 Scalar tmp = logn + et[i]->get_loggam() ;
152 dnulg[i] = factx*tmp.dsdx().val_grid_point(0, 0, 0, 0) ;
153
154 // For graphical outputs:
155 Scalar tgraph = logn - log( (1. + et[i]->get_chi_auto()) / (1. + et[i]->get_Psi_auto()) ) ;
156 // tmp = log( (1. + et[i]->get_chi_comp()) / (1. + et[i]->get_Psi_comp()) ) ;
157 tgraph.std_spectral_base() ;
158 save_profile(tgraph, 0., 10., 0.5*M_PI, 0., "prof_logn.d") ;
159 save_profile(et[i]->get_loggam(), 0., 1.8, 0.5*M_PI, 0., "prof_loggam.d") ;
160
161 //------------------------------------------------------------------
162 // Psi^4/N^2 = in the center of the star ---> asn2[i]
163 //------------------------------------------------------------------
164
165 Scalar Psi6schi2 = pow(Psi, 6)/(chi % chi) ;
166 Psi6schi2.std_spectral_base() ;
167 asn2[i] = Psi6schi2.val_grid_point(0, 0, 0, 0) ;
168
169 //------------------------------------------------------------------
170 // d/dX(A^2/N^2) in the center of the star ---> dasn2[i]
171 //------------------------------------------------------------------
172
173 dasn2[i] = Psi6schi2.dsdx().val_grid_point(0, 0, 0, 0) ;
174
175 //------------------------------------------------------------------
176 // N^Y in the center of the star ---> ny[i]
177 //------------------------------------------------------------------
178
179 ny[i] = factx*shift(2).val_grid_point(0, 0, 0, 0) ;
180
181 nyso[i] = ny[i] / omega ;
182
183 //------------------------------------------------------------------
184 // dN^Y/dX in the center of the star ---> dny[i]
185 //------------------------------------------------------------------
186
187 dny[i] = shift(2).dsdx().val_grid_point(0, 0, 0, 0) ;
188
189 dnyso[i] = dny[i] / omega ;
190
191 //------------------------------------------------------------------
192 // (N^X)^2 + (N^Y)^2 + (N^Z)^2
193 // in the center of the star ---> npn[i]
194 //------------------------------------------------------------------
195
196 tmp = contract(shift, 0, shift.up_down(flat), 0) ;
197
198 npn[i] = tmp.val_grid_point(0, 0, 0, 0) ;
199 npnso2[i] = npn[i] / omega / omega ;
200
201 //------------------------------------------------------------------
202 // d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
203 // in the center of the star ---> dnpn[i]
204 //------------------------------------------------------------------
205
206 dnpn[i] = factx * tmp.dsdx().val_grid_point(0, 0, 0, 0) ;
207 dnpnso2[i] = dnpn[i] / omega / omega ;
208
209 cout << "Binary_xcts::orbit: central d(nu+log(Gam))/dX : "
210 << dnulg[i] << endl ;
211 cout << "Binary_xcts::orbit: central A^2/N^2 : "
212 << asn2[i] << endl ;
213 cout << "Binary_xcts::orbit: central d(A^2/N^2)/dX : "
214 << dasn2[i] << endl ;
215 cout << "Binary_xcts::orbit: central N^Y : "
216 << ny[i] << endl ;
217 cout << "Binary_xcts::orbit: central dN^Y/dX : "
218 << dny[i] << endl ;
219 cout << "Binary_xcts::orbit: central N.N : "
220 << npn[i] << endl ;
221 cout << "Binary_xcts::orbit: central d(N.N)/dX : "
222 << dnpn[i] << endl ;
223
224
225 ori_x[i] = (et[i]->get_mp()).get_ori_x() ;
226 xgg[i] = factx * (et[i]->xa_barycenter() - ori_x[i]) ;
227
228 }
229
230 xgg1 = xgg[0] ;
231 xgg2 = xgg[1] ;
232
233//---------------------------------
234// axis of rotation
235//---------------------------------
236
237 int relat = 1 ;
238
239 double ori_x1 = ori_x[0] ;
240 double ori_x2 = ori_x[1] ;
241
242 if ( et[0]->get_eos() == et[1]->get_eos() &&
243 fabs( et[0]->get_ent().val_grid_point(0,0,0,0) -
244 et[1]->get_ent().val_grid_point(0,0,0,0) ) < 1.e-14 ) {
245
246 x_axe = 0. ;
247
248 } else {
249
250 Param paraxe ;
251 paraxe.add_int(relat) ;
252 paraxe.add_double( ori_x1, 0) ;
253 paraxe.add_double( ori_x2, 1) ;
254 paraxe.add_double( dnulg[0], 2) ;
255 paraxe.add_double( dnulg[1], 3) ;
256 paraxe.add_double( asn2[0], 4) ;
257 paraxe.add_double( asn2[1], 5) ;
258 paraxe.add_double( dasn2[0], 6) ;
259 paraxe.add_double( dasn2[1], 7) ;
260 paraxe.add_double( nyso[0], 8) ;
261 paraxe.add_double( nyso[1], 9) ;
262 paraxe.add_double( dnyso[0], 10) ;
263 paraxe.add_double( dnyso[1], 11) ;
264 paraxe.add_double( npnso2[0], 12) ;
265 paraxe.add_double( npnso2[1], 13) ;
266 paraxe.add_double( dnpnso2[0], 14) ;
267 paraxe.add_double( dnpnso2[1], 15) ;
268
269 int nitmax_axe = 200 ;
270 int nit_axe ;
271 double precis_axe = 1.e-13 ;
272
273 x_axe = zerosec(fonc_binary_xcts_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
274 precis_axe, nitmax_axe, nit_axe) ;
275
276 cout << "Binary_xcts::orbit : Number of iterations in zerosec for x_axe : "
277 << nit_axe << endl ;
278 }
279
280 cout << "Binary_xcts::orbit: x_axe [km] : " << x_axe / km << endl ;
281
282//-------------------------------------
283// Calcul de la vitesse orbitale
284//-------------------------------------
285
286 Param parf ;
287 parf.add_int(relat) ;
288 parf.add_double( ori_x1, 0) ;
289 parf.add_double( dnulg[0], 1) ;
290 parf.add_double( asn2[0], 2) ;
291 parf.add_double( dasn2[0], 3) ;
292 parf.add_double( ny[0], 4) ;
293 parf.add_double( dny[0], 5) ;
294 parf.add_double( npn[0], 6) ;
295 parf.add_double( dnpn[0], 7) ;
296 parf.add_double( x_axe, 8) ;
297
298 double omega1 = fact_omeg_min * omega ;
299 double omega2 = fact_omeg_max * omega ;
300 cout << "Binary_xcts::orbit: omega1, omega2 [rad/s] : "
301 << omega1 * f_unit << " " << omega2 * f_unit << endl ;
302
303 // Search for the various zeros in the interval [omega1, omega2]
304 // ------------------------------------------------------------
305 int nsub = 50 ; // total number of subdivisions of the interval
306 Tbl* azer = 0x0 ;
307 Tbl* bzer = 0x0 ;
308 zero_list(fonc_binary_xcts_orbit, parf, omega1, omega2, nsub,
309 azer, bzer) ;
310
311 // Search for the zero closest to the previous value of omega
312 // ----------------------------------------------------------
313 double omeg_min, omeg_max ;
314 int nzer = azer->get_taille() ; // number of zeros found by zero_list
315 cout << "Binary_xcts:orbit : " << nzer <<
316 " zero(s) found in the interval [omega1, omega2]." << endl ;
317 cout << "omega, omega1, omega2 : " << omega << " " << omega1
318 << " " << omega2 << endl ;
319 cout << "azer : " << *azer << endl ;
320 cout << "bzer : " << *bzer << endl ;
321
322 if (nzer == 0) {
323 cout <<
324 "Binary_xcts::orbit: WARNING : no zero detected in the interval"
325 << endl << " [" << omega1 * f_unit << ", "
326 << omega2 * f_unit << "] rad/s !" << endl ;
327 omeg_min = omega1 ;
328 omeg_max = omega2 ;
329 }
330 else {
331 double dist_min = fabs(omega2 - omega1) ;
332 int i_dist_min = -1 ;
333 for (int i=0; i<nzer; i++) {
334 // Distance of previous value of omega from the center of the
335 // interval [azer(i), bzer(i)]
336 double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
337 if (dist < dist_min) {
338 dist_min = dist ;
339 i_dist_min = i ;
340 }
341 }
342 omeg_min = (*azer)(i_dist_min) ;
343 omeg_max = (*bzer)(i_dist_min) ;
344 }
345
346 delete azer ; // Tbl allocated by zero_list
347 delete bzer ; //
348
349 cout << "Binary_xcts::orbit : interval selected for the search of the zero : "
350 << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
351 << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
352
353 // Computation of the zero in the selected interval by the secant method
354 // ---------------------------------------------------------------------
355 int nitermax = 200 ;
356 int niter ;
357 double precis = 1.e-13 ;
358 omega = zerosec_b(fonc_binary_xcts_orbit, parf, omeg_min, omeg_max,
359 precis, nitermax, niter) ;
360
361 cout << "Binary_xcts::orbit : Number of iterations in zerosec for omega : "
362 << niter << endl ;
363
364 cout << "Binary_xcts::orbit : omega [rad/s] : "
365 << omega * f_unit << endl ;
366
367}
368
369//*************************************************
370// Function used for search of the rotation axis
371//*************************************************
372
373double fonc_binary_xcts_axe(double x_rot, const Param& paraxe) {
374
375 int relat = paraxe.get_int() ;
376
377 double ori_x1 = paraxe.get_double(0) ;
378 double ori_x2 = paraxe.get_double(1) ;
379 double dnulg_1 = paraxe.get_double(2) ;
380 double dnulg_2 = paraxe.get_double(3) ;
381 double asn2_1 = paraxe.get_double(4) ;
382 double asn2_2 = paraxe.get_double(5) ;
383 double dasn2_1 = paraxe.get_double(6) ;
384 double dasn2_2 = paraxe.get_double(7) ;
385 double nyso_1 = paraxe.get_double(8) ;
386 double nyso_2 = paraxe.get_double(9) ;
387 double dnyso_1 = paraxe.get_double(10) ;
388 double dnyso_2 = paraxe.get_double(11) ;
389 double npnso2_1 = paraxe.get_double(12) ;
390 double npnso2_2 = paraxe.get_double(13) ;
391 double dnpnso2_1 = paraxe.get_double(14) ;
392 double dnpnso2_2 = paraxe.get_double(15) ;
393
394 double om2_star1 ;
395 double om2_star2 ;
396
397 double x1 = ori_x1 - x_rot ;
398 double x2 = ori_x2 - x_rot ;
399
400 if (relat == 1) {
401
402 double andan_1 = 0.5 * dasn2_1 + asn2_1 * dnulg_1 ;
403 double andan_2 = 0.5 * dasn2_2 + asn2_2 * dnulg_2 ;
404
405 double bpb_1 = x1 * x1 - 2. * nyso_1 * x1 + npnso2_1 ;
406 double bpb_2 = x2 * x2 - 2. * nyso_2 * x2 + npnso2_2 ;
407
408 double cpc_1 = 0.5 * dnpnso2_1 + x1 * (1. - dnyso_1) - nyso_1 ;
409 double cpc_2 = 0.5 * dnpnso2_2 + x2 * (1. - dnyso_2) - nyso_2 ;
410
411 om2_star1 = dnulg_1 / (andan_1 * bpb_1 + asn2_1 * cpc_1) ;
412 om2_star2 = dnulg_2 / (andan_2 * bpb_2 + asn2_2 * cpc_2) ;
413
414 }
415 else {
416
417 om2_star1 = dnulg_1 / x1 ;
418 om2_star2 = dnulg_2 / x2 ;
419
420 }
421
422 return om2_star1 - om2_star2 ;
423
424}
425
426//*****************************************************************************
427// Fonction utilisee pour la recherche de omega par la methode de la secante
428//*****************************************************************************
429
430double fonc_binary_xcts_orbit(double om, const Param& parf) {
431
432 int relat = parf.get_int() ;
433
434 double xc = parf.get_double(0) ;
435 double dnulg = parf.get_double(1) ;
436 double asn2 = parf.get_double(2) ;
437 double dasn2 = parf.get_double(3) ;
438 double ny = parf.get_double(4) ;
439 double dny = parf.get_double(5) ;
440 double npn = parf.get_double(6) ;
441 double dnpn = parf.get_double(7) ;
442 double x_axe = parf.get_double(8) ;
443
444 double xx = xc - x_axe ;
445 double om2 = om*om ;
446
447 double dphi_cent ;
448
449 if (relat == 1) {
450 double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
451
452 dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
453 - 0.5*bpb* dasn2 )
454 / ( 1 - asn2 * bpb ) ;
455 }
456 else {
457 dphi_cent = - om2*xx ;
458 }
459
460 return dnulg + dphi_cent ;
461
462}
463
464
465}
double x_axe
Absolute X coordinate of the rotation axis.
Definition binary_xcts.h:84
void orbit(double fact_omeg_min, double fact_omeg_max, double &xgg1, double &xgg2)
Computes the orbital angular velocity omega and the position of the rotation axis x_axe.
Star_bin_xcts * et[2]
Array of the two stars (to perform loops on the stars): et[0] contains the address of star1 and et[1]...
Definition binary_xcts.h:75
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition binary_xcts.h:80
Base class for coordinate mappings.
Definition map.h:670
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Metric for tensor calculation.
Definition metric.h:90
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 int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:292
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
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 .
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density,.
const Scalar & get_chi_auto() const
Returns the scalar field generated principally by the star.
Definition star.h:1372
const Scalar & get_chi_comp() const
Returns the scalar field generated principally by the companion star.
Definition star.h:1377
const Scalar & get_Psi_comp() const
Returns the scalar field generated principally by the companion star.
Definition star.h:1367
const Metric & get_flat() const
Return the flat metric defined on the mapping (Spherical components with respect to the mapping of th...
Definition star.h:1400
const Scalar & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition star.h:1336
const Scalar & get_Psi_auto() const
Returns the scalar field generated principally by the star.
Definition star.h:1362
const Map & get_mp() const
Returns the mapping.
Definition star.h:355
const Eos & get_eos() const
Returns the equation of state.
Definition star.h:361
const Vector & get_beta() const
Returns the shift vector .
Definition star.h:402
Basic array class.
Definition tbl.h:161
int get_taille() const
Gives the total size (ie dim.taille)
Definition tbl.h:397
Tensor field of valence 1.
Definition vector.h:188
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
void save_profile(const Scalar &uu, double r_min, double r_max, double theta, double phi, const char *filename)
Saves in a file the profile of a Scalar along some radial axis determined by a fixed value of .
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.
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:89
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.