Skip to content

GitLab

Commit 85850767 authored by Michal Dovčiak's avatar Michal Dovčiak
Browse files

Add new reverberation model based on XILLVER tables by Javier Garcia. 

Add kynxilrev model.
Rename the package of kynrefrev and kynxilrev models to kynreverb.
parent b33445e0
Expand all Hide whitespace changes
Inline Side-by-side
README.md
View file @ 85850767
This diff is collapsed.
fft_reverberation.c
View file @ 85850767
......@@ -6,9 +6,16 @@
#define PI 3.14159265358979
#define PI2 6.28318530717959
extern int xs_write(char* wrtstr, int idest);
extern double t0, fwrap;
extern int exclude_energy;
//extern int xs_write(char* wrtstr, int idest);
//extern double t0, fwrap;
//extern int exclude_energy;
// Let's declare variables that can be seen from ide and FFT routine
double *del=NULL, t0, fwrap=0.;
float *radius=NULL;
long int nradius;
int exclude_energy=0;
void fft_reverberation(const double *ear, int ne, double *photar,
double frequency1, double frequency2, int photar_sw,
......@@ -16,7 +23,7 @@ void fft_reverberation(const double *ear, int ne, double *photar,
double *far, double *far_prim,
double *flux_bands, double *flux_bands_prim,
int nbands, double tot_time_sec,
double flare_duration_sec, char *cparam){
double flare_duration_sec, char *cparam, char *cname){
void four(double *data_r, double *data_im, long n1, int isign);
......@@ -70,12 +77,12 @@ FILE *fw1, *fw2, *fw3, *fw4, *fw5, *fw6;
//heap instead of on an available size-limited stack!
if ( ( data_r = (double *) malloc( nn * sizeof(double) ) ) == NULL ||
( data_im = (double *) malloc( nn * sizeof(double) ) ) == NULL ) {
xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5);
xs_write("fft_reverberation: Failed to allocate memory for tmp arrays.", 5);
for (ie = 0; ie < ne; ie++) photar[ie] = 0.;
goto error;
}
if ( ( freq = (double *) malloc( (nn/2+1) * sizeof(double) ) ) == NULL ) {
xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5);
xs_write("fft_reverberation: Failed to allocate memory for tmp arrays.", 5);
for (ie = 0; ie < ne; ie++) photar[ie] = 0.;
goto error;
}
......@@ -96,7 +103,7 @@ if( ( r_part = (double *) malloc( ne * (nn/2+1) * sizeof(double) ) ) == NULL |
( phase_tot_u2 = (double *) malloc( ne * (nn/2+1) * sizeof(double) ) ) == NULL ||
( ampl_tot_etot2 = (double *) malloc( ne * (nn/2+1) * sizeof(double) ) ) == NULL ||
( phase_tot_etot2 = (double *) malloc( ne * (nn/2+1) * sizeof(double) ) ) == NULL ) {
xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5);
xs_write("fft_reverberation: Failed to allocate memory for tmp arrays.", 5);
for (ie = 0; ie < ne; ie++) photar[ie] = 0.;
goto error;
}
......@@ -112,7 +119,7 @@ if( nbands > 0 && (
( ampl_bands_tot = (double *) malloc( nbands * (nn/2+1) * sizeof(double) ) ) == NULL ||
( phase_bands_tot = (double *) malloc( nbands * (nn/2+1) * sizeof(double) ) ) == NULL ||
( phase_bands_tot_u1 = (double *) malloc( nbands * (nn/2+1) * sizeof(double) ) ) == NULL ) ){
xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5);
xs_write("fft_reverberation: Failed to allocate memory for tmp arrays.", 5);
for (ie = 0; ie < ne; ie++) photar[ie] = 0.;
goto error;
}
......@@ -1319,17 +1326,17 @@ if( abs(photar_sw) > 10 && nbands > 0){
// Write the tables
// energy dependent Fourier transform
if(abs(photar_sw)<=10){
sprintf(filename, "kynrefrev_%s_real.dat",cparam);
sprintf(filename, "%s_%s_real.dat",cname,cparam);
fw1 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_imag.dat",cparam);
sprintf(filename, "%s_%s_imag.dat",cname,cparam);
fw2 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_ampl.dat",cparam);
sprintf(filename, "%s_%s_ampl.dat",cname,cparam);
fw3 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_phase.dat",cparam);
sprintf(filename, "%s_%s_phase.dat",cname,cparam);
fw4 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_phase_u1.dat",cparam);
sprintf(filename, "%s_%s_phase_u1.dat",cname,cparam);
fw5 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_phase_u2.dat",cparam);
sprintf(filename, "%s_%s_phase_u2.dat",cname,cparam);
fw6 = fopen(filename, "w");
for(ifr=0;ifr<=nn/2;ifr++){
for(ie=0;ie<ne;ie++){
......@@ -1355,15 +1362,15 @@ if(abs(photar_sw)<=10){
fclose(fw6);
}
// Fourier transform for chosen energy bands
sprintf(filename, "kynrefrev_%s_bands_real.dat",cparam);
sprintf(filename, "%s_%s_bands_real.dat",cname,cparam);
fw1 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_bands_imag.dat",cparam);
sprintf(filename, "%s_%s_bands_imag.dat",cname,cparam);
fw2 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_bands_ampl.dat",cparam);
sprintf(filename, "%s_%s_bands_ampl.dat",cname,cparam);
fw3 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_bands_phase.dat",cparam);
sprintf(filename, "%s_%s_bands_phase.dat",cname,cparam);
fw4 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_bands_phase_u1.dat",cparam);
sprintf(filename, "%s_%s_bands_phase_u1.dat",cname,cparam);
fw5 = fopen(filename, "w");
//sprintf(filename, "kynrefrev_%s_bands_phase_u2.dat",cparam);
//fw4 = fopen(filename, "w");
......@@ -1393,17 +1400,17 @@ fclose(fw4);
fclose(fw5);
// energy dependent Fourier transform when primary is included
if(abs(photar_sw)<=10){
sprintf(filename, "kynrefrev_%s_real_tot.dat",cparam);
sprintf(filename, "%s_%s_real_tot.dat",cname,cparam);
fw1 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_imag_tot.dat",cparam);
sprintf(filename, "%s_%s_imag_tot.dat",cname,cparam);
fw2 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_ampl_tot.dat",cparam);
sprintf(filename, "%s_%s_ampl_tot.dat",cname,cparam);
fw3 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_phase_tot.dat",cparam);
sprintf(filename, "%s_%s_phase_tot.dat",cname,cparam);
fw4 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_phase_tot_u1.dat",cparam);
sprintf(filename, "%s_%s_phase_tot_u1.dat",cname,cparam);
fw5 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_phase_tot_u2.dat",cparam);
sprintf(filename, "%s_%s_phase_tot_u2.dat",cname,cparam);
fw6 = fopen(filename, "w");
for(ifr=0;ifr<=nn/2;ifr++){
for(ie=0;ie<ne;ie++){
......@@ -1429,15 +1436,15 @@ if(abs(photar_sw)<=10){
fclose(fw6);
}
// Fourier transform for chosen energy bands when primary is included
sprintf(filename, "kynrefrev_%s_bands_real_tot.dat",cparam);
sprintf(filename, "%s_%s_bands_real_tot.dat",cname,cparam);
fw1 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_bands_imag_tot.dat",cparam);
sprintf(filename, "%s_%s_bands_imag_tot.dat",cname,cparam);
fw2 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_bands_ampl_tot.dat",cparam);
sprintf(filename, "%s_%s_bands_ampl_tot.dat",cname,cparam);
fw3 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_bands_phase_tot.dat",cparam);
sprintf(filename, "%s_%s_bands_phase_tot.dat",cname,cparam);
fw4 = fopen(filename, "w");
sprintf(filename, "kynrefrev_%s_bands_phase_tot_u1.dat",cparam);
sprintf(filename, "%s_%s_bands_phase_tot_u1.dat",cname,cparam);
fw5 = fopen(filename, "w");
for(ifr=0;ifr<=nn/2;ifr++){
fprintf(fw1, "%E\t", freq[ifr]);
......@@ -1465,7 +1472,7 @@ fclose(fw4);
fclose(fw5);
// frequency integrated energy dependent functions
if(abs(photar_sw)<=10){
sprintf(filename, "kynrefrev_%s_fft_tot_int.dat",cparam);
sprintf(filename, "%s_%s_fft_tot_int.dat",cname,cparam);
fw1 = fopen(filename, "w");
for(ie=0;ie<ne;ie++)
fprintf(fw1, "%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\n",
......@@ -1485,13 +1492,13 @@ if(abs(photar_sw)<=10){
im_part_tot_int_etot-im_part_tot_int[ie],
r_part_tot_int_etot- r_part_tot_int[ie]))/PI/freq_wrap0);
fclose(fw1);
sprintf(filename, "kynrefrev_%s_freq_wrap.dat",cparam);
sprintf(filename, "%s_%s_freq_wrap.dat",cname,cparam);
fw1 = fopen(filename, "w");
fprintf(fw1, "%E", freq_wrap0);
fclose(fw1);
// frequency band integrated Fourier transform
if( frequency1 > 0. && frequency1 < frequency2){
sprintf(filename, "kynrefrev_%s_fft_tot_fband.dat",cparam);
sprintf(filename, "%s_%s_fft_tot_fband.dat",cname,cparam);
fw1 = fopen(filename, "w");
for(ie=0;ie<ne;ie++)
fprintf(fw1, "%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\n",
......@@ -1871,3 +1878,17 @@ for (i=0;i<n1;i++) {
free((char*) (data));
}
#undef SWAP
/*******************************************************************************
*******************************************************************************/
double incgamma(double a, double x){
double gfunc;
long i;
gfunc=0.;
for(i=100000;i>=1;i--) gfunc=(a-i)*i/(2.*i+1.+x-a+gfunc);
gfunc = pow(x,a) * exp(-x) / (1.+x-a+gfunc);
return(gfunc);
}
/*******************************************************************************
*******************************************************************************/
lmodel-kynrefrev.dat → lmodel-kynreverb.dat
View file @ 85850767
......@@ -38,3 +38,44 @@ dt/Af s/- 0. -1e8 -1e8 1e8 1e8 -1.
k/qf " " 1. -1e8 -1e8 1e8 1e8 -1.
$xsw 16.
nthreads " " 4. 1. 1. 100. 100. -1.
kynxilrev 39 0. 1.0e20 c_KYNxilrev add 0
a/M GM/c 1. 0. 0. 1. 1. 0.2
theta_o deg 30. 0. 0. 89. 89. 5.
rin GM/c^2 1. 1. 1. 1000. 1000. -0.5
$ms 1.
rout GM/c^2 1000. 1. 1. 1000. 1000 -1.0
phi deg 0. -180. -180. 180. 180. -1.
dphi deg 360. 0. 0. 360. 360. -1.
M/M8 " " 0.1 1e-8 1e-8 1e+3 1e+3 5.
height GM/c^2 3. -20. -20. 100. 100. -1.
PhoIndex " " 2.0 1.2 1.2 3.4 3.4 -0.1
L/Ledd " " 0.001 0. 0. 1e+10 1e+10 10.
Np:Nr " " 1.0 0. 0. 10. 10. -0.1
density " " 1. 1. 1. 1e+4 1e+4 5.
den_prof " " 0. -5. -5. 0. 0. 0.1
abun " " 1. 0.5 0.5 20. 20. 0.1
Ecut keV 300. 5. 5. 1000. 1000. -10.
therm " " 0. -2. -2. 2. 2. -0.1
arate Ledd 0.1 1.e-5 1.e-5 100. 100. -0.1
f_col " " 2.4 1. 1. 3. 3. -0.1
alpha GM/c^2 -6. -100. -100. 100. 100. 1.
beta GM/c^2 0. -100. -100. 100. 100. 1.
rcloud GM/c^2 0. 0. 0. 100. 100. 1.
zshift " " 0. -0.999 -0.999 10. 10. -0.1
limb " " 0. 0. 0. 2. 2. -1.
$tab 11.
ntable " " 80. 0. 0. 99. 99. -1.
nrad " " -4488. -10000. -10000. 10000. 10000. -100.
division " " -1. 0. 0. 1. 1. -1.
nphi " " 180. 1. 1. 20000. 20000. -100.
deltaT GM/c^3 1. 1e-3 1e-3 10. 10. -0.05
nt_ratio " " 1. 1. 1. 10. 10. -1.
t1/f1/E1 s/Hz/keV 0.3 -1e8 -1e8 1e16 1e16 -1.
t2/f2/E2 s/Hz/keV 0.8 -1e8 -1e8 1e16 1e16 -1.
Eref1 keV 1. -1. -1. 100. 100. -1.
Eref2 keV 3. 0. 0. 100. 100. -1.
dt/Af s/- 0. -1e8 -1e8 1e8 1e8 -1.
k/qf " " 1. -1e8 -1e8 1e8 1e8 -1.
$xsw 16.
nthreads " " 4. 1. 1. 100. 100. -1.
other/Makefile
View file @ 85850767
SRC_KYNREFREV=xside_threads.c xside.c libxspec.c fft_reverberation.c xskynrefrev.c
SRC_KYNXILREV=xside_threads.c xside.c libxspec.c fft_reverberation.c xskynxilrev.c
CC=gcc
INCLUDE=-I/usr/local/share/heasoft/x86_64-unknown-linux-gnu-libc2.19/include
......@@ -8,3 +9,6 @@ CFLAGS=-fPIC -O3 -Wall -DOUTSIDE_XSPEC -lm
kynrefrev: $(SRC_KYNREFREV)
$(CC) $(CFLAGS) $(INCLUDE) -o $@ $(SRC_KYNREFREV) $(LIBRARY_PATH) $(LIBRARY)
kynxilrev: $(SRC_KYNXILREV)
$(CC) $(CFLAGS) $(INCLUDE) -o $@ $(SRC_KYNXILREV) $(LIBRARY_PATH) $(LIBRARY)
xside.c
View file @ 85850767
......@@ -2069,13 +2069,12 @@ while(iiloop<(iloop+nloops) && iiloop<=(nrad*nphi)){
}//next iiloop
error:
if(far_loc != NULL); free((void *) far_loc); far_loc = NULL;
if(far_loc != NULL){free((void *) far_loc); far_loc = NULL;}
if(polar){
if(qar_loc != NULL); free((void *) qar_loc); qar_loc = NULL;
if(uar_loc != NULL); free((void *) uar_loc); uar_loc = NULL;
if(var_loc != NULL); free((void *) var_loc); var_loc = NULL;
if(qar_loc != NULL){free((void *) qar_loc); qar_loc = NULL;}
if(uar_loc != NULL){free((void *) uar_loc); uar_loc = NULL;}
if(var_loc != NULL){free((void *) var_loc); var_loc = NULL;}
}
return;
}
/*******************************************************************************
......
xskynrefrev.c
View file @ 85850767
......@@ -567,10 +567,10 @@ return(0);
#define EC 300.
// Let's declare variables that can be seen from ide and FFT routine
double *del=NULL, t0, fwrap=0.;
float *radius=NULL;
long int nradius;
int exclude_energy=0;
extern double *del, t0, fwrap;
extern float *radius;
extern long int nradius;
extern int exclude_energy;
/* Let's declare variables that are common for the main and emissivity
subroutines */
......@@ -588,7 +588,7 @@ extern char* FGMODF(void);
extern char* FGMSTR(char* dname);
extern void FPMSTR(const char* value1, const char* value2);
extern int xs_write(char* wrtstr, int idest);
//extern double incgamma(double a, double x);
extern double incgamma(double a, double x);
//extern void cutoffpl(double *ear, int ne, double *param, double *photar);
void KYNrefrev(double *ear_xspec, int ne_xspec, const double *param,
......@@ -605,7 +605,8 @@ extern void fft_reverberation(const double *ear, int ne, double *photar,
long nn, long nt, double *far, double *far_prim,
double *flux_bands, double *flux_bands_prim,
int nbands, double tot_time_sec,
double flare_duration_sec, char *cparam);
double flare_duration_sec, char *cparam,
char *cname);
void emis_KYNrefrev(double** ear_loc, const int ne_loc, const long nt,
......@@ -615,7 +616,7 @@ void emis_KYNrefrev(double** ear_loc, const int ne_loc, const long nt,
const double alpha_o, const double beta_o,
const double delay, const double g);
double incgamma(double a, double x);
//double incgamma(double a, double x);
/* Let's declare static variables whose values should be remembered for next
......@@ -637,7 +638,7 @@ FILE *fw;
char errortxt[255], filename[255];
char reflion[255], text[255];
char kyxiin[32], kyxiout[32], kyLxLamp[32], kyRefl[32], kyfwrap[32];
char cparam[12];
char cparam[12], cname[10];
double ide_param[25];
double ear_short[4];
double *time=NULL;
......@@ -697,6 +698,7 @@ int ihorizon, irow, anynul, status = 0;
sprintf(cparam, "%3ld_%02ld_%04ld",
lround(100.*(1.+sqrt(1.-param[0]*param[0]))),
lround(param[1]),lround(10.*param[8]));
sprintf(cname, "kynrefrev");
// am - black hole angular momentum
ide_param[0] = param[0];
am = ide_param[0];
......@@ -1215,7 +1217,7 @@ if(first_h) {
if (status) {
if (status) ffrprt(stderr, status);
ffclos(fptr, &status);
xs_write("\nkynrefionx: set the KYDIR to the directory with the KY tables.", 5);
xs_write("\nkynrefrev: set the KYDIR to the directory with the KY tables.", 5);
for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.;
goto error;
}
......@@ -1561,23 +1563,23 @@ if ((strcmp(kydir, FGMSTR(pname)) != 0) || (rix != rix_old) || first_rix) {
if (status) {
if (status) ffrprt(stderr, status);
ffclos(fptr, &status);
xs_write("\nkynrefionx: set the KYDIR to the directory with the KY tables.", 5);
xs_write("\nkynrefrev: set the KYDIR to the directory with the KY tables.", 5);
for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.;
goto error;
}
if (((strcmp(kydir, FGMSTR(pname)) != 0) || (rix != rix_old)) && !first_rix) {
// Free memory from tmp arrays...
if(abun != NULL); free((void *) abun); abun = NULL;
if(gam != NULL); free((void *) gam); gam = NULL;
if(xi != NULL); free((void *) xi); xi = NULL;
if(logxi != NULL); free((void *) logxi); logxi = NULL;
if(energy0 != NULL); free((void *) energy0); energy0 = NULL;
if(emission != NULL); free((void *) emission); emission = NULL;
if(flux0 != NULL); free((void *) flux0); flux0 = NULL;
if(flux1 != NULL); free((void *) flux1); flux1 = NULL;
if(energy1 != NULL); free((void *) energy1); energy1 = NULL;
if(energy2 != NULL); free((void *) energy2); energy2 = NULL;
if(flx != NULL); free((void *) flx); flx = NULL;
if(abun != NULL){ free((void *) abun); abun = NULL;}
if(gam != NULL){ free((void *) gam); gam = NULL;}
if(xi != NULL){ free((void *) xi); xi = NULL;}
if(logxi != NULL){ free((void *) logxi); logxi = NULL;}
if(energy0 != NULL){ free((void *) energy0); energy0 = NULL;}
if(emission != NULL){ free((void *) emission); emission = NULL;}
if(flux0 != NULL){ free((void *) flux0); flux0 = NULL;}
if(flux1 != NULL){ free((void *) flux1); flux1 = NULL;}
if(energy1 != NULL){ free((void *) energy1); energy1 = NULL;}
if(energy2 != NULL){ free((void *) energy2); energy2 = NULL;}
if(flx != NULL){ free((void *) flx); flx = NULL;}
}
// Let's read tables (binary tables => hdutype=2)
// Move to the first extension ('PARAMETERS')
......@@ -2248,7 +2250,7 @@ if(photar_sw){
// Compute the Fourier transform
fft_reverberation(ear_xspec, ne_xspec, photar, frequency1, frequency2, photar_sw,
time, nn, nt, far, far_prim, flux_bands, flux_bands_prim,
nbands, tot_time_sec, flare_duration_sec, cparam);
nbands, tot_time_sec, flare_duration_sec, cparam, cname);
sprintf(kyfwrap, "%e", fwrap);
FPMSTR(pkyfwrap, kyfwrap);
if(photar_sw == 5 || photar_sw == 7) for(ie=0;ie<ne;ie++)photar[ie] *= ampl_ampl;
......@@ -2328,16 +2330,16 @@ fclose(fw);
error:
//Let's free all allocated arrays
if (time != NULL){free((void *) time); time = NULL;}
if (far != NULL){free((void *) far); far = NULL;}
if (far_final != NULL){free((void *) far_final); far_final = NULL;}
if (flux != NULL){free((void *) flux); flux = NULL;}
if (flux_bands != NULL){free((void *) flux_bands); flux_bands = NULL;}
if (spectrum != NULL){free((void *) spectrum); spectrum = NULL;}
if (far_prim != NULL){free((void *) far_prim); far_prim = NULL;}
if (spectrum_prim != NULL){free((void *) spectrum_prim); spectrum_prim = NULL;}
if (photar1 != NULL){free((void *) photar1); photar1 = NULL;}
if (gf != NULL){free((void *) gf); gf = NULL;}
if (time != NULL){ free((void *) time); time = NULL;}
if (far != NULL){ free((void *) far); far = NULL;}
if (far_final != NULL){ free((void *) far_final); far_final = NULL;}
if (flux != NULL){ free((void *) flux); flux = NULL;}
if (flux_bands != NULL){ free((void *) flux_bands); flux_bands = NULL;}
if (spectrum != NULL){ free((void *) spectrum); spectrum = NULL;}
if (far_prim != NULL){ free((void *) far_prim); far_prim = NULL;}
if (spectrum_prim != NULL){ free((void *) spectrum_prim); spectrum_prim = NULL;}
if (photar1 != NULL){ free((void *) photar1); photar1 = NULL;}
if (gf != NULL){ free((void *) gf); gf = NULL;}
////if (ear1 != NULL){free((void *) ear1); ear1 = NULL;}
return;
......@@ -2392,116 +2394,116 @@ if ((ir0 == 0) || (ir0 >= nradius)) {
// Let's interpolate delay between two radii
delay0 = ttmp * del[ir0] + ttmp1 * del[ir0 - 1];
if((delay0+delay-t0 > nt*dt) || (t0-delay0-delay > flare_duration_rg)){
for(it=0; it < nt; it++) far_loc[ne_loc+(ne_loc+1)*it]=0.;
}else{
// Let's interpolate gfactor between two radii
gfactor = ttmp * gfac[ir0] + ttmp1 * gfac[ir0 - 1];
// Let's interpolate cosmu0 between two radii
// cosmu0 = ttmp * cosin[ir0] + ttmp1 * cosin[ir0 - 1];
// Let's interpolate lensing between two radii
lensing = (ttmp * transf_d[ir0] + ttmp1 * transf_d[ir0 - 1]) *
h * sqrt(1. - 2. * h / (h * h + am2)) / (r * r + h * h) / r;
// Let's compute the emitted flux at the particular radius
if (lensing != 0.) {
if(sw == 1) Finc = Np * LEDD / RG2 / mass * lensing * gfactor;
else Finc = Np * LEDD / RG2 / mass * lensing * pow(gfactor, gamma0 - 1.);
if(nH0 > 0.){
if (qn != 0.) ionisation = log10( Finc / ( nH0 * 1e15 * pow(r, qn) ) );
else ionisation = log10( Finc / ( nH0 * 1e15 ) );
factor1 = 1.;
} else{
ionisation = log10(-nH0) + qn * log10(r);
factor1 = Finc / 1e15 / (-nH0) / pow(r,qn);
}
// given ionisation, find the corresponding index in logXi():
imin = 0;
imax = nxi;
ixi0 = nxi / 2;
while ((imax - imin) > 1) {
if (ionisation >= logxi[ixi0 - 1]) imin = ixi0;
else imax = ixi0;
ixi0 = (imin + imax) / 2;
}
if (ixi0 == 0) ixi0 = 1;
ttmp = (ionisation - logxi[ixi0 - 1]) / (logxi[ixi0] - logxi[ixi0 - 1]);
if (ttmp < 0.) {
ttmp = 0.;
factor1 *= pow(10, ionisation - logxi[0]);
}
if (ttmp > 1.) {
ttmp = 1.;
factor1 *= pow(10, ionisation - logxi[nxi - 1]);
}
ttmp1 = 1. - ttmp;
for (ie = 0; ie < ne_loc; ie++) {
y1 = flux1[ie + ne_loc * (ixi0 - 1)];
y2 = flux1[ie + ne_loc * ixi0];
factor2 = exp( energy1[ie] / EC * (1.-1./gfactor) );
fluxe[ie] = (ttmp1 * y1 + ttmp * y2) * factor * factor1 * factor2;
if (sw == 1) fluxe[ie] *= pow(gfactor, gamma0 - 2.);
if (nH0 > 0. && qn != 0.) fluxe[ie] *= pow(r, qn);
}
// add the thermal component
if(thermalisation != 0.){
// compute thermalised total flux (Ftherm = Finc - Frefl
// or Ftherm = thermalisation * Finc)...
// incident flux in SI units ( i.e. J / s / m^2 )
Finc *= ( 1e-3 / 4. / PI );
if( fabs(thermalisation) <= 1. ) Ftherm = fabs(thermalisation) * Finc;
else{
// reflected flux in SI units
Frefl=0.;
ttmp = fluxe[0] * (*(*ear_loc));
for (ie = 1; ie < ne_loc; ie++){
ttmp1 = fluxe[ie] * (*(*ear_loc+ie));
Frefl += (ttmp + ttmp1) / 2. * ((*(*ear_loc+ie))-(*(*ear_loc+ie-1)));
ttmp = ttmp1;
for(it=0; it < nt; it++) far_loc[ne_loc+(ne_loc+1)*it]=0.;
}else{
// Let's interpolate gfactor between two radii
gfactor = ttmp * gfac[ir0] + ttmp1 * gfac[ir0 - 1];
// Let's interpolate cosmu0 between two radii
// cosmu0 = ttmp * cosin[ir0] + ttmp1 * cosin[ir0 - 1];
// Let's interpolate lensing between two radii
lensing = (ttmp * transf_d[ir0] + ttmp1 * transf_d[ir0 - 1]) *
h * sqrt(1. - 2. * h / (h * h + am2)) / (r * r + h * h) / r;
// Let's compute the emitted flux at the particular radius
if (lensing != 0.) {
if(sw == 1) Finc = Np * LEDD / RG2 / mass * lensing * gfactor;
else Finc = Np * LEDD / RG2 / mass * lensing * pow(gfactor, gamma0 - 1.);
if(nH0 > 0.){
if (qn != 0.) ionisation = log10( Finc / ( nH0 * 1e15 * pow(r, qn) ) );
else ionisation = log10( Finc / ( nH0 * 1e15 ) );
factor1 = 1.;
} else{
ionisation = log10(-nH0) + qn * log10(r);
factor1 = Finc / 1e15 / (-nH0) / pow(r,qn);
}
// given ionisation, find the corresponding index in logXi():
imin = 0;
imax = nxi;
ixi0 = nxi / 2;
while ((imax - imin) > 1) {
if (ionisation >= logxi[ixi0 - 1]) imin = ixi0;
else imax = ixi0;
ixi0 = (imin + imax) / 2;
}
if (ixi0 == 0) ixi0 = 1;
ttmp = (ionisation - logxi[ixi0 - 1]) / (logxi[ixi0] - logxi[ixi0 - 1]);
if (ttmp < 0.) {
ttmp = 0.;
factor1 *= pow(10, ionisation - logxi[0]);
}
if (ttmp > 1.) {
ttmp = 1.;
factor1 *= pow(10, ionisation - logxi[nxi - 1]);
}
ttmp1 = 1. - ttmp;
for (ie = 0; ie < ne_loc; ie++) {
y1 = flux1[ie + ne_loc * (ixi0 - 1)];
y2 = flux1[ie + ne_loc * ixi0];
factor2 = exp( energy1[ie] / EC * (1.-1./gfactor) );
fluxe[ie] = (ttmp1 * y1 + ttmp * y2) * factor * factor1 * factor2;
if (sw == 1) fluxe[ie] *= pow(gfactor, gamma0 - 2.);
if (nH0 > 0. && qn != 0.) fluxe[ie] *= pow(r, qn);
}
// add the thermal component
if(thermalisation != 0.){
// compute thermalised total flux (Ftherm = Finc - Frefl
// or Ftherm = thermalisation * Finc)...
// incident flux in SI units ( i.e. J / s / m^2 )
Finc *= ( 1e-3 / 4. / PI );
if( fabs(thermalisation) <= 1. ) Ftherm = fabs(thermalisation) * Finc;
else{
// reflected flux in SI units
Frefl=0.;
ttmp = fluxe[0] * (*(*ear_loc));
for (ie = 1; ie < ne_loc; ie++){
ttmp1 = fluxe[ie] * (*(*ear_loc+ie));
Frefl += (ttmp + ttmp1) / 2. * ((*(*ear_loc+ie))-(*(*ear_loc+ie-1)));
ttmp = ttmp1;
}
Frefl *= KEV * 1e4;
// fprintf(stdout,"%lg\n", Frefl / Finc);
if(Finc < Frefl){
sprintf(errortext,
"kynrefrev: reflection is larger than illumination at radius %lg GM/c^2!",
r);
xs_write(errortext, 5);
Frefl=Finc;
}
Ftherm = Finc - Frefl;
}
Frefl *= KEV * 1e4;
// fprintf(stdout,"%lg\n", Frefl / Finc);
if(Finc < Frefl){
sprintf(errortext,
"kynrefrev: reflection is larger than illumination at radius %lg GM/c^2!",
r);
xs_write(errortext, 5);
Frefl=Finc;
// compute the temperature of the accretion disc
x = sqrt(r);
Ccal = 1. - 3. / (x*x) + 2. * am / (x * x * x);
if(am < 1.)
Lcal = 1. / x * (x - x0 - 1.5 * am * log(x / x0) -
3. * pow(x1 - am,2.)/x1/(x1 - x2)/(x1 - x3) * log((x - x1)/(x0 - x1))-
3. * pow(x2 - am,2.)/x2/(x2 - x1)/(x2 - x3) * log((x - x2)/(x0 - x2))-
3. * pow(x3 - am,2.)/x3/(x3 - x1)/(x3 - x2) * log((x - x3)/(x0 - x3)));
else Lcal = 1. / x * (x - 1 + 1.5 * log((x + 2.) / 3. / x));
temp = Tnorm * pow(x, -1.5) * pow( arate / mass2 * Lcal / Ccal, 0.25 );
// compute new temperature
Fbb = SIGMA * pow( temp / K_KEVK / f_col, 4. );
temp_new = K_KEVK * pow( ( Fbb + Ftherm ) / SIGMA, 0.25 ) * f_col;
// compute the additional flux due to thermalisation (i.e. "above"
// the stationary thermal flux)
for(ie = 0; ie < ne_loc; ie++){
flux_thermal = flx[ie] * ( 1. / (exp((*(*ear_loc+ie)) / temp_new) - 1.)
- 1. / (exp((*(*ear_loc+ie)) / temp) - 1.) );
if( thermalisation < 0. ) fluxe[ie] = flux_thermal;
else fluxe[ie] += flux_thermal;
}
Ftherm = Finc - Frefl;
}
// compute the temperature of the accretion disc
x = sqrt(r);
Ccal = 1. - 3. / (x*x) + 2. * am / (x * x * x);
if(am < 1.)