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
30char 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
68namespace 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].
const Map & get_mp() const
Returns the mapping.
Definition etoile.h:659
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Definition etoile.h:667
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
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
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.