LORENE
binary_anashift_xcts.C
1/*
2 * Method of class Binary_xcts to set some analytical form
3 * to the shift vector (see file binary_xcts.h for documentation).
4 */
5
6/*
7 * Copyright (c) 2010 Michal Bejger
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License version 2
13 * as published by the Free Software Foundation.
14 *
15 * LORENE is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with LORENE; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 *
24 */
25
26char binary_anashift_xcts_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_anashift_xcts.C,v 1.3 2014/10/13 08:52:45 j_novak Exp $" ;
27
28/*
29 * $Id: binary_anashift_xcts.C,v 1.3 2014/10/13 08:52:45 j_novak Exp $
30 * $Log: binary_anashift_xcts.C,v $
31 * Revision 1.3 2014/10/13 08:52:45 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.2 2010/06/15 07:58:32 m_bejger
35 * Minor corrections
36 *
37 * Revision 1.1 2010/05/04 07:35:54 m_bejger
38 * Initial version
39 *
40 * $Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_anashift_xcts.C,v 1.3 2014/10/13 08:52:45 j_novak Exp $
41 *
42 */
43
44// Headers C
45#include "math.h"
46
47// Headers Lorene
48#include "binary_xcts.h"
49#include "tenseur.h"
50#include "unites.h"
51
52namespace Lorene {
54
55 using namespace Unites ;
56
57 for (int i=0; i<2; i++) {
58
59 // Radius of the star:
60 double a0 = et[i]->ray_eq() ;
61
62 // Mass ratio
63 double p_mass = et[i]->mass_g() / et[1-i]->mass_g() ;
64
65 // G M Omega R / (1 + mass_ratio)
66 double www = ggrav * et[i]->mass_g() * omega
67 * separation() / (1. + p_mass) ;
68
69 const Map& mp = et[i]->get_mp() ;
70 Scalar tmp(mp) ;
71 Scalar tmp_ext(mp) ;
72 int nzet = et[i]->get_nzet() ;
73 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
74
75 Vector w_beta (mp, CON, mp.get_bvect_cart()) ;
76 Scalar khi_beta (mp) ;
77
78 // Computation of w_beta
79 // ----------------------
80 // X component
81 // -----------
82
83 w_beta.set(1) = 0 ;
84
85 // Y component
86 // -----------
87
88 // For the incompressible case :
89 tmp = - 6 * www / a0 * ( 1 - (mp.r)*(mp.r) / (3*a0*a0) ) ;
90
91 tmp.annule(nzet, nzm1) ;
92 tmp_ext = - 4 * www / mp.r ;
93 tmp_ext.annule(0, nzet-1) ;
94
95 w_beta.set(2) = tmp + tmp_ext ;
96
97 // Z component
98 // -----------
99 w_beta.set(3) = 0 ;
100
101 w_beta.std_spectral_base() ;
102
103 // Computation of khi_beta
104 // ------------------------
105
106 tmp = 2 * www / a0 * (mp.y) * ( 1 - 3 * (mp.r)*(mp.r) / (5*a0*a0) ) ;
107 tmp.annule(nzet, nzm1) ;
108 tmp_ext = 0.8 * www * a0*a0 * (mp.sint) * (mp.sinp)
109 / (mp.r * mp.r) ;
110 tmp_ext.annule(0, nzet-1) ;
111
112 khi_beta = tmp + tmp_ext ;
113
114 // Sets the standard spectral bases for a scalar field
115 khi_beta.std_spectral_base() ;
116
117 // Computation of beta auto.
118 // --------------------------
119
120 Tensor xdw_temp (w_beta.derive_con(et[i]->get_flat())) ;
121
122 Tenseur x_d_w_temp (et[i]->get_mp(),2,CON,et[i]->get_mp().get_bvect_cart()) ;
123 x_d_w_temp.set_etat_qcq() ;
124
125 for (int j=0; j<3; j++)
126 for (int k=0; k<3; k++)
127 x_d_w_temp.set(j,k) = xdw_temp(k+1, j+1) ;
128
129 Tenseur x_d_w = skxk (x_d_w_temp) ;
130 x_d_w.dec_dzpuis() ;
131
132 Vector xdw (et[i]->get_mp(), CON, et[i]->get_mp().get_bvect_cart()) ;
133 for (int j=0; j<3; j++)
134 xdw.set(j+1) = x_d_w(j) ;
135
136 // See Eq (92) from Gourgoulhon et al.(2001) and with the new
137 // convention for shift = - N^i
138
139 Vector d_khi = khi_beta.derive_con(et[i]->get_flat()) ;
140 d_khi.dec_dzpuis(2) ;
141
142 et[i]->set_beta_auto() = - 7./8. * w_beta + 1./8. *
143 (d_khi + xdw) ;
144
146
147 }
148
149}
150}
void analytical_shift()
Sets some analytical template for the shift vector (via the members w_shift and khi_shift of the two ...
double separation() const
Returns the coordinate separation of the two stellar centers [r_unit].
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
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...
Vector & set_beta_auto()
Read/write 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:1400
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.