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;