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 
29 
30 
31 /*
32  * $Id: binary_anashift.C,v 1.10 2016/12/05 16:17:47 j_novak Exp $
33  * $Log: binary_anashift.C,v $
34  * Revision 1.10 2016/12/05 16:17:47 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.9 2014/10/13 08:52:44 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.8 2005/09/13 19:38:31 f_limousin
41  * Reintroduction of the resolution of the equations in cartesian coordinates.
42  *
43  * Revision 1.7 2005/02/17 17:34:50 f_limousin
44  * Change the name of some quantities to be consistent with other classes
45  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
46  *
47  * Revision 1.6 2004/03/25 10:29:01 j_novak
48  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
49  *
50  * Revision 1.5 2004/02/27 10:01:32 f_limousin
51  * Correct sign of shift_auto to agree with the new convention
52  * for shift.
53  *
54  * Revision 1.4 2004/01/22 10:09:41 f_limousin
55  * First executable version
56  *
57  * Revision 1.3 2004/01/20 15:21:23 f_limousin
58  * First version
59  *
60  *
61  * $Header: /cvsroot/Lorene/C++/Source/Binary/binary_anashift.C,v 1.10 2016/12/05 16:17:47 j_novak Exp $
62  *
63  */
64 
65 // Headers C
66 #include "math.h"
67 
68 // Headers Lorene
69 #include "binary.h"
70 #include "tenseur.h"
71 #include "unites.h"
72 
73 namespace Lorene {
75 
76  using namespace Unites ;
77 
78  for (int i=0; i<2; i++) {
79 
80  // Radius of the star:
81  double a0 = et[i]->ray_eq() ;
82 
83  // Mass ratio
84  double p_mass = et[i]->mass_g() / et[1-i]->mass_g() ;
85 
86  // G M Omega R / (1+p)
87  double www = ggrav * et[i]->mass_g() * omega
88  * separation() / (1. + p_mass) ;
89 
90  const Map& mp = et[i]->get_mp() ;
91  Scalar tmp(mp) ;
92  Scalar tmp_ext(mp) ;
93  int nzet = et[i]->get_nzet() ;
94  int nzm1 = mp.get_mg()->get_nzone() - 1 ;
95 
96  Vector w_beta (mp, CON, mp.get_bvect_cart()) ;
97  Scalar khi_beta (mp) ;
98 
99  // Computation of w_beta
100  // ----------------------
101  // X component
102  // -----------
103 
104  w_beta.set(1) = 0 ;
105 
106  // Y component
107  // -----------
108 
109  // For the incompressible case :
110  tmp = - 6 * www / a0 * ( 1 - (mp.r)*(mp.r) / (3*a0*a0) ) ;
111 
112  tmp.annule(nzet, nzm1) ;
113  tmp_ext = - 4 * www / mp.r ;
114  tmp_ext.annule(0, nzet-1) ;
115 
116  w_beta.set(2) = tmp + tmp_ext ;
117 
118  // Z component
119  // -----------
120  w_beta.set(3) = 0 ;
121 
122  w_beta.std_spectral_base() ;
123 
124  // Computation of khi_beta
125  // ------------------------
126 
127  tmp = 2 * www / a0 * (mp.y) * ( 1 - 3 * (mp.r)*(mp.r) / (5*a0*a0) ) ;
128  tmp.annule(nzet, nzm1) ;
129  tmp_ext = 0.8 * www * a0*a0 * (mp.sint) * (mp.sinp)
130  / (mp.r * mp.r) ;
131  tmp_ext.annule(0, nzet-1) ;
132 
133  khi_beta = tmp + tmp_ext ;
134 
135  // Sets the standard spectral bases for a scalar field
136  khi_beta.std_spectral_base() ;
137 
138 
139  // Computation of beta auto.
140  // --------------------------
141 
142  Tensor xdw_temp (w_beta.derive_con(et[i]->get_flat())) ;
143 
144  Tenseur x_d_w_temp (et[i]->get_mp(),2,CON,et[i]->get_mp().get_bvect_cart()) ;
145  x_d_w_temp.set_etat_qcq() ;
146  for (int j=0; j<3; j++)
147  for (int k=0; k<3; k++)
148  x_d_w_temp.set(j,k) = xdw_temp(k+1, j+1) ;
149 
150  Tenseur x_d_w = skxk (x_d_w_temp) ;
151  x_d_w.dec_dzpuis() ;
152 
153  Vector xdw (et[i]->get_mp(), CON, et[i]->get_mp().get_bvect_cart()) ;
154  for (int j=0; j<3; j++)
155  xdw.set(j+1) = x_d_w(j) ;
156 
157  // See Eq (92) from Gourgoulhon et al.(2001) and with the new
158  // convention for shift = - N^i
159 
160  Vector d_khi = khi_beta.derive_con(et[i]->get_flat()) ;
161  d_khi.dec_dzpuis(2) ;
162 
163  et[i]->set_beta_auto() = - 7./8. * w_beta + 1./8. *
164  (d_khi + xdw) ;
165 
167 
168  }
169 
170 }
171 }
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Lorene prototypes.
Definition: app_hor.h:67
void analytical_shift()
Sets some analytical template for the shift vector (via the members w_shift and khi_shift of the two ...
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: binary.h:94
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
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:322
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:817
Coord sint
Definition: map.h:739
double separation() const
Returns the coordinate separation of the two stellar centers [r_unit].
Definition: binary.C:610
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
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:859
Coord sinp
Definition: map.h:741
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Tensor handling.
Definition: tensor.h:294
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
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:1023
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:809
void dec_dzpuis()
dzpuis -= 1 ;
Definition: tenseur.C:1120
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
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:745
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:111
virtual double mass_g() const
Gravitational mass.
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
Coord r
r coordinate centered on the grid
Definition: map.h:736
Vector & set_beta_auto()
Read/write of .
Definition: star_bin.C:471