28 char scalarBH_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Compobj/scalarBH.C,v 1.6 2016/05/10 12:52:32 f_vincent Exp $" ;
82 ifstream file(file_name) ;
84 cerr <<
"Problem in opening the file " << file_name << endl ;
90 int nrfile, nthetafile;
91 file >> nrfile >> nthetafile ;
95 if (rHor<0. || nrfile<0. || nthetafile<0.){
96 cerr <<
"In scalarBH::scalarBH(Map,char*): " 97 <<
"Bad parameters rHor, nrfile, nthetafile" << endl;
104 cout <<
"nr, ntheta from file = " << nrfile <<
" " << nthetafile << endl;
105 cout <<
"rHor from file = " << rHor << endl;
107 cout <<
"Scalar field values provided" << endl;
109 cout <<
"Scalar field values not provided, put to zero" << endl;
111 cerr <<
"In scalarBH::scalarBH(Map,char*): " 112 <<
"Bad parameter isphi" << endl;
116 double* Xfile =
new double[nrfile * nthetafile] ;
117 double* thetafile =
new double[nrfile * nthetafile] ;
118 double* f0file =
new double[nrfile * nthetafile] ;
119 double* f1file =
new double[nrfile * nthetafile] ;
120 double* f2file =
new double[nrfile * nthetafile] ;
121 double* wwfile =
new double[nrfile * nthetafile] ;
122 double* sfieldfile =
new double[nrfile * nthetafile] ;
124 cout <<
"Reading metric data... ";
125 for (
int ii=0;ii<nrfile*nthetafile;ii++){
128 file >> thetafile[ii] ;
133 file >> sfieldfile[ii] ;
140 cout <<
"done" << endl;
143 double Xbefmax = Xfile[nrfile-2];
153 double delta_theta = 0.1*fabs(thetafile[nrfile] - thetafile[0]);
155 cout <<
"Starting interpolating Lorene grid... " ;
164 for (
int l=l_min; l<nz; l++) {
169 for (
int k=0; k<np; k++){
170 for (
int j=0; j<nt; j++){
171 double th0 = theta(l, k, j, 0);
172 for (
int i=0; i<nr; i++){
173 double r0 = rr(l, k, j, i);
175 if (r0 < __infinity){
176 x0 =
sqrt(r0*r0 - rH2);
187 double thc = thetafile[ith];
188 while (fabs(th0 - thc) > delta_theta){
191 thc = thetafile[ith];
194 double xc = Xfile[ir];
196 cerr <<
"In scalarBH::ScalarBH(): " 197 "r should be zero here" << endl;
207 }
else if(xx0 > Xbefmax && xx0< 1.){
220 double f0interp, f1interp, f2interp, winterp, sfieldinterp;
222 double xcinf = Xfile[ir-1], xcsup = Xfile[ir+1];
224 if (ir-3>0) {irext1=ir-2; irext2=ir-3;}
225 else if (ir+3<nrfile) {irext1=ir+2; irext2=ir+3;}
227 cerr <<
"scalarBH::scalarBH(): bad radial indice" << endl;
230 double xcext1 = Xfile[irext1], xcext2 = Xfile[irext2];
238 if (fabs(thc-th0)>delta_theta){
239 cerr <<
"scalarBH::ScalarBH(): theta problem in grid" << endl;
240 cerr <<
"Theta info: " << thc <<
" " << th0 << endl;
243 if (xx0 <= Xbefmax &&
244 (xx0 < xcinf || xx0 > xc || xx0 > xcsup)){
245 cerr <<
"scalarBH::ScalarBH(): rad problem in grid" << endl;
246 cerr <<
"Radial info: " << xcinf <<
" " << xx0 <<
" " 247 << xc <<
" " << xcsup << endl;
249 }
else if (xx0 > Xbefmax &&
250 (xx0 < xcinf || xx0 < xc || xx0 > xcsup)){
251 cerr <<
"scalarBH::ScalarBH(): special rad " 252 "problem in grid" << endl;
253 cerr <<
"Radial info: " << xcinf <<
" " << xx0 <<
" " 254 << xc <<
" " << xcsup << endl;
258 double polyrinf = (xx0-xc)*(xx0-xcsup)*(xx0-xcext1)*(xx0-xcext2)/((xcinf-xc)*(xcinf-xcsup)*(xcinf-xcext1)*(xcinf-xcext2));
259 double polyrmid = (xx0-xcinf)*(xx0-xcsup)*(xx0-xcext1)*(xx0-xcext2)/((xc-xcinf)*(xc-xcsup)*(xc-xcext1)*(xc-xcext2));
260 double polyrsup = (xx0-xcinf)*(xx0-xc)*(xx0-xcext1)*(xx0-xcext2)/((xcsup-xcinf)*(xcsup-xc)*(xcsup-xcext1)*(xcsup-xcext2));
261 double polyrext1 = (xx0-xcinf)*(xx0-xc)*(xx0-xcsup)*(xx0-xcext2)/((xcext1-xcinf)*(xcext1-xc)*(xcext1-xcsup)*(xcext1-xcext2));
262 double polyrext2 = (xx0-xcinf)*(xx0-xc)*(xx0-xcsup)*(xx0-xcext1)/((xcext2-xcinf)*(xcext2-xc)*(xcext2-xcsup)*(xcext2-xcext1));
265 double f0ext1 = f0file[irext1], f0ext2 = f0file[irext2],
266 f0inf = f0file[ir-1],
268 f0sup = f0file[ir+1], f1ext1=f1file[irext1],
269 f1ext2=f1file[irext2],
270 f1inf = f1file[ir-1], f1mid = f1file[ir],
271 f1sup = f1file[ir+1], f2ext1=f2file[irext1],
272 f2ext2=f2file[irext2],
273 f2inf = f2file[ir-1], f2mid = f2file[ir],
274 f2sup = f2file[ir+1], wext1=wwfile[irext1],
275 wext2=wwfile[irext2],
276 winf = wwfile[ir-1], wmid = wwfile[ir],
277 wsup = wwfile[ir+1], sfext1=sfieldfile[irext1],
278 sfext2=sfieldfile[irext2],
279 sfinf = sfieldfile[ir-1],
280 sfmid = sfieldfile[ir], sfsup = sfieldfile[ir+1];
287 f0interp = f0ext1*polyrext1 + f0ext2*polyrext2 + f0inf*polyrinf
288 + f0mid*polyrmid + f0sup*polyrsup;
289 f1interp = f1ext1*polyrext1 + f1ext2*polyrext2 + f1inf*polyrinf
290 + f1mid*polyrmid + f1sup*polyrsup;
291 f2interp = f2ext1*polyrext1 + f2ext2*polyrext2 + f2inf*polyrinf
292 + f2mid*polyrmid + f2sup*polyrsup;
293 winterp = wext1*polyrext1 + wext2*polyrext2 + winf*polyrinf
294 + wmid*polyrmid + wsup*polyrsup;
295 sfieldinterp = sfext1*polyrext1 + sfext2*polyrext2
297 + sfmid*polyrmid + sfsup*polyrsup;
304 f0interp = f0file[ir] ;
305 f1interp = f1file[ir] ;
306 f2interp = f2file[ir] ;
307 winterp = wwfile[ir] ;
308 sfieldinterp = sfieldfile[ir] ;
337 cout <<
"done." << endl;
342 cout <<
"Starting updating metric... " ;
344 cout <<
"done." << endl;
447 ost << endl <<
"Black hole with scalar hair (class ScalarBH) " << endl ;
501 void ScalarBH::update_metric() {
518 gam.set(2,2) =
exp(2*
ff1);
519 gam.set(2,2).std_spectral_base() ;
521 gam.set(3,3) =
exp(2*
ff2) ;
522 gam.set(3,3).std_spectral_base() ;
virtual ~ScalarBH()
Destructor.
Scalar ff2
Metric field F_2 of Herdeiro & Radu (2015)
Cmp exp(const Cmp &)
Exponential.
virtual void sauve(FILE *) const
Save in a file.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Cmp sqrt(const Cmp &)
Square root.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tensor field of valence 0 (or component of a tensorial field).
Base class for coordinate mappings.
Scalar ww
Metric field W of Herdeiro & Radu (2015)
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
virtual void del_deriv() const
Deletes all the derived quantities.
Tbl & set_domain(int l)
Read/write of the value in a given domain.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Coord tet
coordinate centered on the grid
void operator=(const ScalarBH &)
Assignment to another ScalarBH.
Scalar ff1
Metric field F_1 of Herdeiro & Radu (2015)
Base class for stationary compact objects (***under development***).
void operator=(const Compobj &)
Assignment to another Compobj.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Vector beta
Shift vector .
ScalarBH(Map &mp_i, const char *file_name)
Standard constructor.
int get_nzone() const
Returns the number of domains.
Scalar ff0
Metric field F_0 of Herdeiro & Radu (2015)
double rHor
Event horizon coordinate radius.
virtual void del_deriv() const
Deletes all the derived quantities.
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Black hole with scalar hair spacetime (***under development***).
Scalar nn
Lapse function N .
Scalar sfield
Scalar field (modulus of Phi)
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Scalar & set(int)
Read/write access to a component.
Class intended to describe valence-2 symmetric tensors.
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Coord r
r coordinate centered on the grid