LORENE
binaire_ana_shift.C
1 /*
2  * Method of class Binaire to set some analytical form to the shift vector.
3  *
4  * (see file binaire.h for documentation).
5  */
6 
7 /*
8  * Copyright (c) 2000-2001 Eric Gourgoulhon
9  * Copyright (c) 2000-2001 Keisuke Taniguchi
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 char binaire_ana_shift_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_ana_shift.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $" ;
31 
32 /*
33  * $Id: binaire_ana_shift.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $
34  * $Log: binaire_ana_shift.C,v $
35  * Revision 1.3 2014/10/13 08:52:44 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.2 2004/03/25 10:28:59 j_novak
39  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
40  *
41  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
42  * LORENE
43  *
44  * Revision 2.3 2000/03/17 15:36:26 eric
45  * Suppression de l'appel a analytical_omega().
46  *
47  * Revision 2.2 2000/03/17 15:27:11 eric
48  * Appel de la fonction analytical_omega() pour fixer la valeur de omega.
49  *
50  * Revision 2.1 2000/03/16 09:37:28 eric
51  * Utilisation du cas incompressible plutot que n=1.
52  *
53  * Revision 2.0 2000/03/15 16:43:35 eric
54  * *** empty log message ***
55  *
56  *
57  * $Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_ana_shift.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $
58  *
59  */
60 
61 // Headers C
62 #include "math.h"
63 
64 // Headers Lorene
65 #include "binaire.h"
66 #include "unites.h"
67 
68 namespace Lorene {
70 
71  // Does nothing for a Newtonian star
72  // ---------------------------------
73  if ( !star1.is_relativistic() ){
74  assert( !star2.is_relativistic() ) ;
75  return ;
76  }
77 
78 
79  using namespace Unites ;
80 
81  for (int i=0; i<2; i++) {
82 
83  // Radius of the star:
84  double a0 = et[i]->ray_eq() ;
85 
86  // Mass ratio
87  double p_mass = et[i]->mass_g() / et[1-i]->mass_g() ;
88 
89  // G M Omega R / (1+p)
90  double www = ggrav * et[i]->mass_g() * omega
91  * separation() / (1. + p_mass) ;
92 
93  const Map& mp = et[i]->get_mp() ;
94  Cmp tmp(mp) ;
95  Cmp tmp_ext(mp) ;
96  int nzet = et[i]->get_nzet() ;
97  int nzm1 = mp.get_mg()->get_nzone() - 1 ;
98 
99  // Computation of w_shift
100  // ----------------------
101  et[i]->set_w_shift().set_etat_qcq() ;
102 
103  // X component
104  // -----------
105  et[i]->set_w_shift().set(0) = 0 ;
106 
107  // Y component
108  // -----------
109 
110 // For the incompressible case :
111  tmp = - 6 * www / a0 * ( 1 - (mp.r)*(mp.r) / (3*a0*a0) ) ;
112 
113 // For the compressible (n=1) case :
114 // Mtbl xi = M_PI * mp.r / a0 ;
115 // Mtbl sinc = sin(xi) / xi ;
116 // The value of sinc is set to 1 at the origin
117 // for (int k=0; k<mp.get_mg()->get_np(0); k++) {
118 // for (int j=0; j<mp.get_mg()->get_nt(0); j++) {
119 // sinc.set(0, k, j, 0) = 1 ;
120 // }
121 // }
122 // tmp = - 4 * www / a0 * ( 1 + sinc ) ;
123 
124  tmp.annule(nzet, nzm1) ;
125  tmp_ext = - 4 * www / mp.r ;
126  tmp_ext.annule(0, nzet-1) ;
127 
128  et[i]->set_w_shift().set(1) = tmp + tmp_ext ;
129 
130  // Z component
131  // -----------
132  et[i]->set_w_shift().set(2) = 0 ;
133 
134  // Sets the standard spectral bases for Cartesian components
135  et[i]->set_w_shift().set_std_base() ;
136 
137  // Computation of khi_shift
138  // ------------------------
139 
140  tmp = 2 * www / a0 * (mp.y) * ( 1 - 3 * (mp.r)*(mp.r) / (5*a0*a0) ) ;
141  tmp.annule(nzet, nzm1) ;
142  tmp_ext = 0.8 * www * a0*a0 * (mp.sint) * (mp.sinp)
143  / (mp.r * mp.r) ;
144  tmp_ext.annule(0, nzet-1) ;
145 
146  et[i]->set_khi_shift() = tmp + tmp_ext ;
147 
148  // Sets the standard spectral bases for a scalar field
149  et[i]->set_khi_shift().set_std_base() ;
150 
151  }
152 
153 }
154 }
void analytical_shift()
Sets some analytical template for the shift vector (via the members {\tt w_shift} and {\tt khi_shift}...
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): {\tt et[0]} contains the address of {\tt star...
Definition: binaire.h:124
double separation() const
Returns the coordinate separation of the two stellar centers [{\tt r_unit}].
Definition: binaire.C:457
Etoile_bin star2
Second star of the system.
Definition: binaire.h:118
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: binaire.h:129
Etoile_bin star1
First star of the system.
Definition: binaire.h:115
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
Tenseur & set_w_shift()
Read/write of w_shift.
Definition: etoile_bin.C:538
virtual double mass_g() const
Gravitational mass.
Tenseur & set_khi_shift()
Read/write of khi_shift.
Definition: etoile_bin.C:545
int get_nzet() const
Returns the number of domains occupied by the star.
Definition: etoile.h:662
double ray_eq() const
Coordinate radius at , [r_unit].
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Definition: etoile.h:667
const Map & get_mp() const
Returns the mapping.
Definition: etoile.h:659
Base class for coordinate mappings.
Definition: map.h:670
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
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Coord sinp
Definition: map.h:723
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
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 set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.