LORENE
star_bhns_kinema.C
1 /*
2  * Method of class Star_bhns to compute kinematic quantities
3  *
4  * (see file star_bhns.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005-2007 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 version 2
15  * as published by the Free Software Foundation.
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  * $Id: star_bhns_kinema.C,v 1.5 2016/12/05 16:18:16 j_novak Exp $
32  * $Log: star_bhns_kinema.C,v $
33  * Revision 1.5 2016/12/05 16:18:16 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.4 2014/10/13 08:53:41 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.3 2014/10/06 15:13:16 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.2 2008/05/15 19:16:06 k_taniguchi
43  * Change of a parameter.
44  *
45  * Revision 1.1 2007/06/22 01:32:00 k_taniguchi
46  * *** empty log message ***
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_kinema.C,v 1.5 2016/12/05 16:18:16 j_novak Exp $
50  *
51  */
52 
53 // C++ headers
54 //#include <>
55 
56 // C headers
57 #include <cmath>
58 
59 // Lorene headers
60 #include "star_bhns.h"
61 #include "unites.h"
62 
63 namespace Lorene {
64 void Star_bhns::kinema_bhns(bool kerrschild, const double& mass_bh,
65  const double& sepa, double omega,
66  double x_rot, double y_rot) {
67 
68  // Fundamental constants and units
69  // -------------------------------
70  using namespace Unites ;
71 
72  int nz = mp.get_mg()->get_nzone() ;
73  int nzm1 = nz - 1 ;
74 
75  //----------------------
76  // Computation of B^i/N
77  //----------------------
78 
79  // 1/ Computation of omega m^i
80 
81  const Coord& xa = mp.xa ;
82  const Coord& ya = mp.ya ;
83 
84  // bsn.change_triad(mp.get_bvect_cart()) ;
85 
86  if (fabs(mp.get_rot_phi()) < 1.e-10) {
87 
88  bsn.set(1) = - omega * (ya - y_rot) ;
89  bsn.set(2) = omega * (xa - x_rot) ;
90  bsn.set(3) = 0. ;
91 
92  }
93  else {
94 
95  bsn.set(1) = omega * (ya - y_rot) ;
96  bsn.set(2) = - omega * (xa - x_rot) ;
97  bsn.set(3) = 0. ;
98 
99  }
100 
102  bsn.annule_domain(nzm1) ;
103 
104  // 2/ Addition of shift_tot and division by lapse
105 
106  for (int i=1; i<=3; i++) {
107  bsn.set(i) = confo_tot * ( bsn(i) + shift_tot(i) ) / lapconf_tot ;
108  }
110  bsn.annule_domain(nzm1) ;
111 
112  //-----------------------------------------------------
113  // Lorentz factor between the co-orbiting ---> gam0
114  // observer and the Eulerian one
115  // See Eq (23) and (24) from Gourgoulhon et al. (2001)
116  //-----------------------------------------------------
117 
118  Scalar bsn2(mp) ;
119  bsn2 = bsn(1)%bsn(1) + bsn(2)%bsn(2) + bsn(3)%bsn(3) ;
120  bsn2.std_spectral_base() ;
121 
122  if (kerrschild) {
123 
124  double mass = ggrav * mass_bh ;
125 
126  Scalar xx(mp) ;
127  xx = mp.x ;
128  xx.std_spectral_base() ;
129  Scalar yy(mp) ;
130  yy = mp.y ;
131  yy.std_spectral_base() ;
132  Scalar zz(mp) ;
133  zz = mp.z ;
134  zz.std_spectral_base() ;
135 
136  double yns = mp.get_ori_y() ;
137 
138  Scalar rbh(mp) ;
139  rbh = sqrt( (xx+sepa)*(xx+sepa) + (yy+yns)*(yy+yns) + zz*zz ) ;
140  rbh.std_spectral_base() ;
141 
142  Vector ll(mp, CON, mp.get_bvect_cart()) ;
143  ll.set_etat_qcq() ;
144  ll.set(1) = (xx+sepa) / rbh ;
145  ll.set(2) = (yy+yns) / rbh ;
146  ll.set(3) = zz / rbh ;
147  ll.std_spectral_base() ;
148 
149  Scalar msr(mp) ;
150  msr = mass / rbh ;
151  msr.std_spectral_base() ;
152 
153  Scalar llbsn(mp) ;
154  llbsn = ll(1)%bsn(1) + ll(2)%bsn(2) + ll(3)%bsn(3) ;
155  llbsn.std_spectral_base() ;
156 
157  Scalar tmp1(mp) ;
158  tmp1 = 2. * msr % llbsn % llbsn ;
159  tmp1.std_spectral_base() ;
160 
161  gam0 = 1. / sqrt(1. - psi4*(bsn2+tmp1)) ;
163 
164  }
165  else { // Isotropic coordinates with the maximal slicing
166 
167  gam0 = 1. / sqrt(1. - psi4%bsn2) ;
169 
170  }
171 
172  //-----------------------
173  // Centrifugal potential
174  //-----------------------
175 
176  pot_centri = - log( gam0 ) ;
177  pot_centri.annule_domain(nzm1) ;
179 
180  // The derived quantities are obsolete
181  // -----------------------------------
182 
183  del_deriv() ;
184 
185 }
186 }
Scalar psi4
Fourth power of the total conformal factor.
Definition: star_bhns.h:176
Coord xa
Absolute x coordinate.
Definition: map.h:748
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Vector shift_tot
Total shift vector.
Definition: star_bhns.h:144
Scalar lapconf_tot
Total lapconf function.
Definition: star_bhns.h:119
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
Map & mp
Mapping associated with the star.
Definition: star.h:180
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: star_bhns.C:344
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Vector bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition: star_bhns.h:99
double get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:788
Lorene prototypes.
Definition: app_hor.h:67
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
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
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
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
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
Scalar confo_tot
Total conformal factor.
Definition: star_bhns.h:163
void kinema_bhns(bool kerrschild, const double &mass_bh, const double &sepa, double omega, double x_rot, double y_rot)
Computes the quantities bsn and pot_centri .
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Scalar gam0
Lorentz factor between the co-orbiting observer and the Eulerian one.
Definition: star_bhns.h:107
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
Scalar pot_centri
Centrifugal potential.
Definition: star_bhns.h:110
Coord ya
Absolute y coordinate.
Definition: map.h:749
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
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Coord y
y coordinate centered on the grid
Definition: map.h:745
Coord x
x coordinate centered on the grid
Definition: map.h:744
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
Coord z
z coordinate centered on the grid
Definition: map.h:746