LORENE
phys_param.C
1 /*
2  * Method of class Isol_hor to compute physical parameters of the horizon
3  *
4  * (see file isol_hor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Jose Luis Jaramillo
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: phys_param.C,v 1.14 2016/12/05 16:17:56 j_novak Exp $
32  * $Log: phys_param.C,v $
33  * Revision 1.14 2016/12/05 16:17:56 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.13 2014/10/13 08:53:01 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.12 2014/10/06 15:13:11 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.11 2005/11/02 16:09:44 jl_jaramillo
43  * changes in boundary_nn_Dir_lapl
44  *
45  * Revision 1.10 2005/04/15 11:54:21 jl_jaramillo
46  * function to compute the expansion of spherical surfaces
47  *
48  * Revision 1.9 2005/03/22 13:25:36 f_limousin
49  * Small changes. The angular velocity and A^{ij} are computed
50  * with a differnet sign.
51  *
52  * Revision 1.8 2005/03/03 10:10:14 f_limousin
53  * Add the function area_hor().
54  *
55  * Revision 1.7 2005/02/07 10:35:42 f_limousin
56  * Minor changes.
57  *
58  * Revision 1.6 2004/12/22 18:16:16 f_limousin
59  * Mny different changes.
60  *
61  * Revision 1.5 2004/11/18 12:30:01 jl_jaramillo
62  * Definition of b_tilde
63  *
64  * Revision 1.4 2004/10/29 15:44:13 jl_jaramillo
65  * ADM angular momentum added.
66  *
67  * Revision 1.3 2004/09/17 13:37:21 f_limousin
68  * Correction of an error in calculation of the radius
69  *
70  * Revision 1.2 2004/09/09 16:54:53 f_limousin
71  * Add the 2 lines $Id: phys_param.C,v 1.14 2016/12/05 16:17:56 j_novak Exp $Log: for CVS
72  *
73  *
74  *
75  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/phys_param.C,v 1.14 2016/12/05 16:17:56 j_novak Exp $
76  *
77  */
78 
79 // C++ headers
80 #include "headcpp.h"
81 
82 // C headers
83 #include <cstdlib>
84 #include <cassert>
85 
86 // Lorene headers
87 #include "isol_hor.h"
88 #include "metric.h"
89 #include "evolution.h"
90 #include "unites.h"
91 #include "scalar.h"
92 #include "vector.h"
93 #include "graphique.h"
94 #include "utilitaires.h"
95 
96 
97 namespace Lorene {
99 
100  Vector get_radial_vect (ff.get_mp(), CON, *(ff.get_triad()) ) ;
101 
102  get_radial_vect.set(1) = gam_uu()(1,1) ;
103 
104  get_radial_vect.set(2) = gam_uu()(1,2) ;
105 
106  get_radial_vect.set(3) = gam_uu()(1,3) ;
107 
108  get_radial_vect = get_radial_vect / sqrt(gam_uu()(1,1)) ;
109 
110  get_radial_vect.std_spectral_base() ;
111 
112 
113  return get_radial_vect ;
114 
115 }
116 
117 
118 // Think of defining this as a pointer
120 
121  Vector get_radial_vect (ff.get_mp(), CON, *(ff.get_triad()) ) ;
122 
123  get_radial_vect.set(1) = (met_gamt.con())(1,1) ;
124 
125  get_radial_vect.set(2) = (met_gamt.con())(1,2) ;
126 
127  get_radial_vect.set(3) = (met_gamt.con())(1,3) ;
128 
129  get_radial_vect = get_radial_vect / sqrt((met_gamt.con())(1,1)) ;
130 
131  get_radial_vect.std_spectral_base() ;
132 
133 
134  return get_radial_vect ;
135 
136 }
137 
138 
139 const Scalar Isol_hor::b_tilde()const {
140 
141  Scalar tmp = contract( beta(), 0, met_gamt.radial_vect()
142  .down(0, met_gamt), 0) ;
143 
144  return tmp ;
145 
146 }
147 
148 
150 
151  Scalar tmp = sqrt( gam_dd()(2,2) * gam_dd()(3,3) - gam_dd()(2,3)
152  * gam_dd()(2,3)) ;
153 
154  tmp.std_spectral_base() ;
155 
156  return tmp ;
157 
158 }
159 
160 double Isol_hor::area_hor() const {
161 
162  Scalar integrand (darea_hor()) ;
163  integrand.raccord(1) ;
164 
165  return mp.integrale_surface(integrand, radius + 1e-15) ;
166 
167 }
168 
169 
170 double Isol_hor::radius_hor() const {
171 
172  double resu = area_hor() / (4. * M_PI);
173 
174  resu = pow(resu, 1./2.) ;
175 
176  return resu ;
177 
178 }
179 
180 
181 double Isol_hor::ang_mom_hor()const {
182 
183  // Vector \partial_phi
184  Vector phi (ff.get_mp(), CON, *(ff.get_triad()) ) ;
185 
186  Scalar tmp (ff.get_mp() ) ;
187  tmp = 1 ;
188  tmp.std_spectral_base() ;
189  tmp.mult_rsint() ;
190 
191  phi.set(1) = 0. ;
192  phi.set(2) = 0. ;
193  phi.set(3) = tmp ;
194 
195  Scalar k_rphi = contract(contract( radial_vect_hor(), 0, k_dd(), 0), 0,
196  phi, 0) / (8. * M_PI) ;
197 
198  Scalar integrand = k_rphi * darea_hor() ; // we correct with the curved
199  // element of area
200 
201  double ang_mom = mp.integrale_surface(integrand, radius + 1e-15) ;
202 
203  return ang_mom ;
204 
205 }
206 
207 
208 // Mass (fundamental constants made 1)
209 double Isol_hor::mass_hor()const {
210 
211  double rr = radius_hor() ;
212 
213  double tmp = sqrt( pow( rr, 4) + 4 * pow( ang_mom_hor(), 2) ) / ( 2 * rr ) ;
214 
215  return tmp ;
216 
217 }
218 
219 
220 // Surface gravity
221 double Isol_hor::kappa_hor() const{
222 
223  double rr = radius_hor() ;
224 
225  double jj = ang_mom_hor() ;
226 
227  double tmp = (pow( rr, 4) - 4 * pow( jj, 2)) / ( 2 * pow( rr, 3)
228  * sqrt( pow( rr, 4) + 4 * pow( jj, 2) ) ) ;
229 
230  return tmp ;
231 
232 }
233 
234 
235 // Orbital velocity
236 double Isol_hor::omega_hor()const {
237 
238  double rr = radius_hor() ;
239 
240  double jj = ang_mom_hor() ;
241 
242  double tmp = 2 * jj / ( rr * sqrt( pow( rr, 4) + 4 * pow( jj, 2) ) ) ;
243 
244  return tmp ;
245 
246 }
247 
248 
249 // ADM angular momentum
250 
251 double Isol_hor::ang_mom_adm()const {
252 
253  Scalar integrand = (k_dd()(1,3) - gam_dd()(1,3) * trk()) / (8. * M_PI) ;
254 
255  integrand.mult_rsint() ; // in order to pass from the triad
256  // component to the coordinate basis
257 
258  double tmp = mp.integrale_surface_infini(integrand) ;
259 
260  return tmp ;
261 
262 }
263 
264 // Expansion
265 
267 
268  Scalar expa = contract(gam().radial_vect().derive_cov(gam()), 0,1)
269  + contract(contract(k_dd(), 0, gam().radial_vect(), 0),
270  0, gam().radial_vect(), 0) - trk() ;
271 
272  return expa ;
273 }
274 }
double area_hor() const
Area of the horizon.
Definition: phys_param.C:160
virtual const Vector & beta() const
shift vector at the current time step (jtime )
const Vector radial_vect_hor() const
Vector radial normal.
Definition: phys_param.C:98
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
Lorene prototypes.
Definition: app_hor.h:67
double ang_mom_hor() const
Angular momentum (modulo)
Definition: phys_param.C:181
double radius
Radius of the horizon in LORENE&#39;s units.
Definition: isol_hor.h:266
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
Definition: metric.h:309
double kappa_hor() const
Surface gravity.
Definition: phys_param.C:221
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
Definition: metric.C:365
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
Tensor field of valence 1.
Definition: vector.h:188
const Map & get_mp() const
Returns the mapping.
Definition: metric.h:202
double omega_hor() const
Orbital velocity.
Definition: phys_param.C:236
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
virtual const Sym_tensor & gam_dd() const
Induced metric (covariant components ) at the current time step (jtime )
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
const Scalar darea_hor() const
Element of area of the horizon.
Definition: phys_param.C:149
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
Scalar expansion() const
Expansion of the outgoing null normal ( )
Definition: phys_param.C:266
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
Definition: phys_param.C:139
Metric met_gamt
3 metric tilde
Definition: isol_hor.h:326
double ang_mom_adm() const
ADM angular Momentum.
Definition: phys_param.C:251
const Metric & gam() const
Induced metric at the current time step (jtime )
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
double radius_hor() const
Radius of the horizon.
Definition: phys_param.C:170
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime ) ...
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:510
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
const Vector tradial_vect_hor() const
Vector radial normal tilde.
Definition: phys_param.C:119
double mass_hor() const
Mass computed at the horizon.
Definition: phys_param.C:209
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.