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 
26 char 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 
52 namespace 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 }
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:391
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
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Base class for coordinate mappings.
Definition: map.h:670
double separation() const
Returns the coordinate separation of the two stellar centers [r_unit].
Definition: binary_xcts.C:434
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
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
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
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
Vector & set_beta_auto()
Read/write of .
Coord sint
Definition: map.h:721
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
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
void analytical_shift()
Sets some analytical template for the shift vector (via the members w_shift and khi_shift of the two ...
virtual double mass_g() const
Gravitational mass.
Coord sinp
Definition: map.h:723
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: binary_xcts.h:80
Tensor handling.
Definition: tensor.h:288
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
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
void dec_dzpuis()
dzpuis -= 1 ;
Definition: tenseur.C:1104
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
int get_nzet() const
Returns the number of domains occupied by the star.
Definition: star.h:358
Coord y
y coordinate centered on the grid
Definition: map.h:727
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Coord r
r coordinate centered on the grid
Definition: map.h:718