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 char star_bhns_kinema_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_kinema.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $" ;
29 
30 /*
31  * $Id: star_bhns_kinema.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $
32  * $Log: star_bhns_kinema.C,v $
33  * Revision 1.4 2014/10/13 08:53:41 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.3 2014/10/06 15:13:16 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.2 2008/05/15 19:16:06 k_taniguchi
40  * Change of a parameter.
41  *
42  * Revision 1.1 2007/06/22 01:32:00 k_taniguchi
43  * *** empty log message ***
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_kinema.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $
47  *
48  */
49 
50 // C++ headers
51 //#include <>
52 
53 // C headers
54 #include <cmath>
55 
56 // Lorene headers
57 #include "star_bhns.h"
58 #include "unites.h"
59 
60 namespace Lorene {
61 void Star_bhns::kinema_bhns(bool kerrschild, const double& mass_bh,
62  const double& sepa, double omega,
63  double x_rot, double y_rot) {
64 
65  // Fundamental constants and units
66  // -------------------------------
67  using namespace Unites ;
68 
69  int nz = mp.get_mg()->get_nzone() ;
70  int nzm1 = nz - 1 ;
71 
72  //----------------------
73  // Computation of B^i/N
74  //----------------------
75 
76  // 1/ Computation of omega m^i
77 
78  const Coord& xa = mp.xa ;
79  const Coord& ya = mp.ya ;
80 
81  // bsn.change_triad(mp.get_bvect_cart()) ;
82 
83  if (fabs(mp.get_rot_phi()) < 1.e-10) {
84 
85  bsn.set(1) = - omega * (ya - y_rot) ;
86  bsn.set(2) = omega * (xa - x_rot) ;
87  bsn.set(3) = 0. ;
88 
89  }
90  else {
91 
92  bsn.set(1) = omega * (ya - y_rot) ;
93  bsn.set(2) = - omega * (xa - x_rot) ;
94  bsn.set(3) = 0. ;
95 
96  }
97 
99  bsn.annule_domain(nzm1) ;
100 
101  // 2/ Addition of shift_tot and division by lapse
102 
103  for (int i=1; i<=3; i++) {
104  bsn.set(i) = confo_tot * ( bsn(i) + shift_tot(i) ) / lapconf_tot ;
105  }
107  bsn.annule_domain(nzm1) ;
108 
109  //-----------------------------------------------------
110  // Lorentz factor between the co-orbiting ---> gam0
111  // observer and the Eulerian one
112  // See Eq (23) and (24) from Gourgoulhon et al. (2001)
113  //-----------------------------------------------------
114 
115  Scalar bsn2(mp) ;
116  bsn2 = bsn(1)%bsn(1) + bsn(2)%bsn(2) + bsn(3)%bsn(3) ;
117  bsn2.std_spectral_base() ;
118 
119  if (kerrschild) {
120 
121  double mass = ggrav * mass_bh ;
122 
123  Scalar xx(mp) ;
124  xx = mp.x ;
125  xx.std_spectral_base() ;
126  Scalar yy(mp) ;
127  yy = mp.y ;
128  yy.std_spectral_base() ;
129  Scalar zz(mp) ;
130  zz = mp.z ;
131  zz.std_spectral_base() ;
132 
133  double yns = mp.get_ori_y() ;
134 
135  Scalar rbh(mp) ;
136  rbh = sqrt( (xx+sepa)*(xx+sepa) + (yy+yns)*(yy+yns) + zz*zz ) ;
137  rbh.std_spectral_base() ;
138 
139  Vector ll(mp, CON, mp.get_bvect_cart()) ;
140  ll.set_etat_qcq() ;
141  ll.set(1) = (xx+sepa) / rbh ;
142  ll.set(2) = (yy+yns) / rbh ;
143  ll.set(3) = zz / rbh ;
144  ll.std_spectral_base() ;
145 
146  Scalar msr(mp) ;
147  msr = mass / rbh ;
148  msr.std_spectral_base() ;
149 
150  Scalar llbsn(mp) ;
151  llbsn = ll(1)%bsn(1) + ll(2)%bsn(2) + ll(3)%bsn(3) ;
152  llbsn.std_spectral_base() ;
153 
154  Scalar tmp1(mp) ;
155  tmp1 = 2. * msr % llbsn % llbsn ;
156  tmp1.std_spectral_base() ;
157 
158  gam0 = 1. / sqrt(1. - psi4*(bsn2+tmp1)) ;
160 
161  }
162  else { // Isotropic coordinates with the maximal slicing
163 
164  gam0 = 1. / sqrt(1. - psi4%bsn2) ;
166 
167  }
168 
169  //-----------------------
170  // Centrifugal potential
171  //-----------------------
172 
173  pot_centri = - log( gam0 ) ;
174  pot_centri.annule_domain(nzm1) ;
176 
177  // The derived quantities are obsolete
178  // -----------------------------------
179 
180  del_deriv() ;
181 
182 }
183 }
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
Coord y
y coordinate centered on the grid
Definition: map.h:727
Coord ya
Absolute y coordinate.
Definition: map.h:731
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:770
Coord z
z coordinate centered on the grid
Definition: map.h:728
Coord x
x coordinate centered on the grid
Definition: map.h:726
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
Coord xa
Absolute x coordinate.
Definition: map.h:730
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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
Vector shift_tot
Total shift vector.
Definition: star_bhns.h:144
Vector bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition: star_bhns.h:99
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: star_bhns.C:341
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 .
Scalar confo_tot
Total conformal factor.
Definition: star_bhns.h:163
Scalar psi4
Fourth power of the total conformal factor.
Definition: star_bhns.h:176
Scalar gam0
Lorentz factor between the co-orbiting observer and the Eulerian one.
Definition: star_bhns.h:107
Scalar lapconf_tot
Total lapconf function.
Definition: star_bhns.h:119
Scalar pot_centri
Centrifugal potential.
Definition: star_bhns.h:110
Map & mp
Mapping associated with the star.
Definition: star.h:180
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
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.