LORENE
binary_anashift.C
1/*
2 * Method of class Binary to set some analytical form to the shift vector.
3 *
4 * (see file binary.h for documentation).
5 */
6
7/*
8 * Copyright (c) 2004 Francois Limousin
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 as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
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
29char binary_anashift_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binary/binary_anashift.C,v 1.9 2014/10/13 08:52:44 j_novak Exp $" ;
30
31/*
32 * $Id: binary_anashift.C,v 1.9 2014/10/13 08:52:44 j_novak Exp $
33 * $Log: binary_anashift.C,v $
34 * Revision 1.9 2014/10/13 08:52:44 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.8 2005/09/13 19:38:31 f_limousin
38 * Reintroduction of the resolution of the equations in cartesian coordinates.
39 *
40 * Revision 1.7 2005/02/17 17:34:50 f_limousin
41 * Change the name of some quantities to be consistent with other classes
42 * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
43 *
44 * Revision 1.6 2004/03/25 10:29:01 j_novak
45 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
46 *
47 * Revision 1.5 2004/02/27 10:01:32 f_limousin
48 * Correct sign of shift_auto to agree with the new convention
49 * for shift.
50 *
51 * Revision 1.4 2004/01/22 10:09:41 f_limousin
52 * First executable version
53 *
54 * Revision 1.3 2004/01/20 15:21:23 f_limousin
55 * First version
56 *
57 *
58 * $Header: /cvsroot/Lorene/C++/Source/Binary/binary_anashift.C,v 1.9 2014/10/13 08:52:44 j_novak Exp $
59 *
60 */
61
62// Headers C
63#include "math.h"
64
65// Headers Lorene
66#include "binary.h"
67#include "tenseur.h"
68#include "unites.h"
69
70namespace Lorene {
72
73 using namespace Unites ;
74
75 for (int i=0; i<2; i++) {
76
77 // Radius of the star:
78 double a0 = et[i]->ray_eq() ;
79
80 // Mass ratio
81 double p_mass = et[i]->mass_g() / et[1-i]->mass_g() ;
82
83 // G M Omega R / (1+p)
84 double www = ggrav * et[i]->mass_g() * omega
85 * separation() / (1. + p_mass) ;
86
87 const Map& mp = et[i]->get_mp() ;
88 Scalar tmp(mp) ;
89 Scalar tmp_ext(mp) ;
90 int nzet = et[i]->get_nzet() ;
91 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
92
93 Vector w_beta (mp, CON, mp.get_bvect_cart()) ;
94 Scalar khi_beta (mp) ;
95
96 // Computation of w_beta
97 // ----------------------
98 // X component
99 // -----------
100
101 w_beta.set(1) = 0 ;
102
103 // Y component
104 // -----------
105
106 // For the incompressible case :
107 tmp = - 6 * www / a0 * ( 1 - (mp.r)*(mp.r) / (3*a0*a0) ) ;
108
109 tmp.annule(nzet, nzm1) ;
110 tmp_ext = - 4 * www / mp.r ;
111 tmp_ext.annule(0, nzet-1) ;
112
113 w_beta.set(2) = tmp + tmp_ext ;
114
115 // Z component
116 // -----------
117 w_beta.set(3) = 0 ;
118
119 w_beta.std_spectral_base() ;
120
121 // Computation of khi_beta
122 // ------------------------
123
124 tmp = 2 * www / a0 * (mp.y) * ( 1 - 3 * (mp.r)*(mp.r) / (5*a0*a0) ) ;
125 tmp.annule(nzet, nzm1) ;
126 tmp_ext = 0.8 * www * a0*a0 * (mp.sint) * (mp.sinp)
127 / (mp.r * mp.r) ;
128 tmp_ext.annule(0, nzet-1) ;
129
130 khi_beta = tmp + tmp_ext ;
131
132 // Sets the standard spectral bases for a scalar field
133 khi_beta.std_spectral_base() ;
134
135
136 // Computation of beta auto.
137 // --------------------------
138
139 Tensor xdw_temp (w_beta.derive_con(et[i]->get_flat())) ;
140
141 Tenseur x_d_w_temp (et[i]->get_mp(),2,CON,et[i]->get_mp().get_bvect_cart()) ;
142 x_d_w_temp.set_etat_qcq() ;
143 for (int j=0; j<3; j++)
144 for (int k=0; k<3; k++)
145 x_d_w_temp.set(j,k) = xdw_temp(k+1, j+1) ;
146
147 Tenseur x_d_w = skxk (x_d_w_temp) ;
148 x_d_w.dec_dzpuis() ;
149
150 Vector xdw (et[i]->get_mp(), CON, et[i]->get_mp().get_bvect_cart()) ;
151 for (int j=0; j<3; j++)
152 xdw.set(j+1) = x_d_w(j) ;
153
154 // See Eq (92) from Gourgoulhon et al.(2001) and with the new
155 // convention for shift = - N^i
156
157 Vector d_khi = khi_beta.derive_con(et[i]->get_flat()) ;
158 d_khi.dec_dzpuis(2) ;
159
160 et[i]->set_beta_auto() = - 7./8. * w_beta + 1./8. *
161 (d_khi + xdw) ;
162
164
165 }
166
167}
168}
Star_bin * 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.h:89
void analytical_shift()
Sets some analytical template for the shift vector (via the members w_shift and khi_shift of the two ...
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition binary.h:94
double separation() const
Returns the coordinate separation of the two stellar centers [r_unit].
Definition binary.C:604
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
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord sint
Definition map.h:721
Coord r
r coordinate centered on the grid
Definition map.h:718
Coord sinp
Definition map.h:723
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
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
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
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:859
Vector & set_beta_auto()
Read/write of .
Definition star_bin.C:468
virtual double mass_g() const
Gravitational mass.
const Map & get_mp() const
Returns the mapping.
Definition star.h:355
double ray_eq() const
Coordinate radius at , [r_unit].
int get_nzet() const
Returns the number of domains occupied by the star.
Definition star.h:358
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
void dec_dzpuis()
dzpuis -= 1 ;
Definition tenseur.C:1104
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.