28char blackhole_killing_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_killing.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $" ;
59#include "utilitaires.h"
67 const double& theta_i,
const int& nrk_phi,
68 const int& nrk_theta)
const {
86 cout <<
"Not yet prepared!!!" << endl ;
95 double dp = 2. * M_PI / double(np) ;
106 xi_t.
set(0) = xi_i(0) ;
107 xi_p.
set(0) = xi_i(1) ;
108 xi_l.
set(0) = xi_i(2) ;
115 xi_ini.
set(0) = xi_i(0) ;
116 xi_ini.
set(1) = xi_i(1) ;
117 xi_ini.
set(2) = xi_i(2) ;
119 double pp_0 = phi_i ;
121 for (
int i=1; i<np+1; i++) {
125 xi_t.
set(i) = xi(0) ;
126 xi_p.
set(i) = xi(1) ;
127 xi_l.
set(i) = xi(2) ;
155 for (
int k=0; k<np; k++) {
187 source_phi =
pow(
confo, 2.) * rr * st / xi_phi ;
193 double hh = 2. * M_PI / double(nn) ;
197 double t1, t2, t3, t4, t5 ;
205 for (
int i=0; i<mm; i++) {
207 t1 = hh * double(4*i) ;
208 t2 = hh * double(4*i+1) ;
209 t3 = hh * double(4*i+2) ;
210 t4 = hh * double(4*i+3) ;
211 t5 = hh * double(4*i+4) ;
213 integ += (hh/45.) * (14.*source_phi.
val_point(rah,M_PI/2.,t1)
214 + 64.*source_phi.
val_point(rah,M_PI/2.,t2)
215 + 24.*source_phi.
val_point(rah,M_PI/2.,t3)
216 + 64.*source_phi.
val_point(rah,M_PI/2.,t4)
217 + 14.*source_phi.
val_point(rah,M_PI/2.,t5)
222 cout <<
"Black_hole:: t_f = " << integ << endl ;
223 double ratio = 0.5 * integ / M_PI ;
225 cout <<
"Black_hole:: t_f / 2M_PI = " << ratio << endl ;
227 for (
int k=0; k<np; k++) {
240 double dt = 0.5 * M_PI / double(nt-1) ;
251 for (
int i=0; i<np; i++) {
253 xi_th.
set(nt-1, i) = xi_t(i) ;
254 xi_ph.
set(nt-1, i) = xi_p(i) ;
255 xi_ll.
set(nt-1, i) = xi_l(i) ;
259 for (
int i=0; i<np; i++) {
262 xi_ini.
set(0) = xi_t(i) ;
263 xi_ini.
set(1) = xi_p(i) ;
264 xi_ini.
set(2) = xi_l(i) ;
266 double pp = double(i) * dp ;
267 double tt_0 = theta_i ;
269 for (
int j=1; j<nt; j++) {
273 xi_th.
set(nt-1-j, i) = xi(0) ;
274 xi_ph.
set(nt-1-j, i) = xi(1) ;
275 xi_ll.
set(nt-1-j, i) = xi(2) ;
292 killing.
set(1) = 0. ;
293 killing.
set(2) = 0.5 ;
294 killing.
set(3) = 0.5 ;
296 for (
int l=0; l<2; l++) {
297 for (
int i=0; i<nr; i++) {
298 for (
int j=0; j<nt; j++) {
299 for (
int k=0; k<np; k++) {
300 (killing.
set(2)).set_grid_point(l, k, j, i) = xi_th(j, k) ;
301 (killing.
set(3)).set_grid_point(l, k, j, i) = xi_ph(j, k) ;
311 double check_norm = 0. ;
312 source_phi =
pow(
confo, 2.) * rr * st / killing(3) ;
315 for (
int i=0; i<mm; i++) {
317 t1 = hh * double(4*i) ;
318 t2 = hh * double(4*i+1) ;
319 t3 = hh * double(4*i+2) ;
320 t4 = hh * double(4*i+3) ;
321 t5 = hh * double(4*i+4) ;
323 check_norm += (hh/45.) *
324 ( 14.*source_phi.
val_point(rah,M_PI/4.,t1)
325 + 64.*source_phi.
val_point(rah,M_PI/4.,t2)
326 + 24.*source_phi.
val_point(rah,M_PI/4.,t3)
327 + 64.*source_phi.
val_point(rah,M_PI/4.,t4)
328 + 14.*source_phi.
val_point(rah,M_PI/4.,t5) ) ;
332 cout <<
"Black_hole:: t_f for M_PI/4 = " << check_norm / M_PI
Tbl runge_kutta_phi_bh(const Tbl &xi_i, const double &phi_i, const int &nrk) const
Compute a forth-order Runge-Kutta integration to the phi direction for the solution of the Killing ve...
Map & mp
Mapping associated with the black hole.
virtual double rad_ah() const
Radius of the apparent horizon.
Tbl runge_kutta_theta_bh(const Tbl &xi_i, const double &theta_i, const double &phi, const int &nrk) const
Compute a forth-order Runge-Kutta integration to the theta direction for the solution of the Killing ...
Vector killing_vect_bh(const Tbl &xi_i, const double &phi_i, const double &theta_i, const int &nrk_phi, const int &nrk_theta) const
Compute the Killing vector of a black hole normalized so that its affine length is 2 M_PI.
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Scalar confo
Conformal factor generated by the black hole.
Coord r
r coordinate centered on the grid
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Tensor field of valence 0 (or component of a tensorial field).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Tensor field of valence 1.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Scalar & set(int)
Read/write access to a component.
Cmp pow(const Cmp &, int)
Power .
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Standard units of space, time and mass.