/* KYNrefrev - ionised reflection (lamp-post Compton reflection) reverberation * model * * ref. Dovciak M., Karas V., Martocchia A., Matt G., Yaqoob T. (2004) * --> this is very old reference on KY package of models when the reverberation * code has not been developed yet but wasdeveloped from this package * ----------------------------------------------------------------------------- * OTHER REFERENCES: * * Dovciak M., Karas V., Martocchia A., Matt G. & Yaqoob T. (2004). XSPEC model * to explore spectral features from black hole sources. * In Proc. of the workshop on processes in the vicinity of black holes * and neutron stars. S.Hledik & Z.Stuchlik, Opava. In press. [astro-ph/0407330] * * Dovciak M., Karas V. & Yaqoob, T. (2004). An extended scheme * for fitting X-ray data with accretion disk spectra * in the strong gravity regime. ApJS, 153, 205. * * Dovciak M. (2004). Radiation of accretion discs in strong gravity. Faculty of * Mathematics and Physics, Charles University, Prague. PhD thesis. * [astro-ph/0411605] * ----------------------------------------------------------------------------- * * This code computes the emission from an acrretion disc that is * illuminated from the primary power-law source located on the axis above the * central black hole with a flash. All relativistic effects are taken into * account (in all three parts of the light path - from the primary source to * the observer, to disc and from the disc to the observer). The code calls * subroutine ide() for integrating local emission over the disc and uses the * fits file 'KBHtables80.fits' defining the transfer functions needed for the * integration. For details on ide() and the fits file see the subroutine ide(). * This model also uses KBHlamp_qt.fits file with transfer functions for the * light coming from primary source and hitting the accretion disc (see the * description of this file below). The reflection is taken from Ross & Fabian * tables reflionx.mod. * * par1 ... a/M - black hole angular momentum (-1 <= a/M <= 1) * par2 ... theta_o - observer inclination in degrees (0-pole, 90-disc) * par3 ... rin - inner edge of non-zero disc emissivity (in GM/c^2 or in * r_mso) * par4 ... ms - 0 - we integrate from inner edge = par3 * 1 - if the inner edge of the disc is below marginally stable * orbit then we integrate emission above MSO only * 2 - we integrate from inner edge given in units of MSO, i.e. * inner edge = par3 * r_mso (the same applies for outer * edge) * par5 ... rout - outer edge of non-zero disc emissivity (in GM/c^2 or in * r_mso) * par6 ... phi - lower azimuth of non-zero disc emissivity (deg) * par7 ... dphi - (phi + dphi) is upper azimuth of non-zero disc emissivity * 0 <= dphi <= 360 (deg) * par8 ... M/M8 - black hole mass in units of 10^8 solar masses * par9 ... height - height on the axis (measured from the center) at which * the primary source is located (GM/c^2) * par10 ... PhoIndex - power-law energy index of the primary flux * par11 ... Np - dN/dt/dOmega, the intrinsic local (if negative) or the * observed (if positive) primary isotropic flux in the * X-ray energy range 2-10keV in units of Ledd * par12 ... NpNr - ratio of the primary to the reflected normalization * 1 - self-consistent model for isotropic primary source * 0 - only reflection, primary source is hidden * - if negative then Np (par11) means the luminosity towards * the disc * - if positive then Np (par11) means the luminosity towards * the observer * par13 ... nH0 - density profile normalization in 10^15 cm^(-3) * par14 ... q_n - radial power-law density profile * par15 ... abun - Fe abundance (in solar abundance) * par16 ... alpha - position of the cloud centre in GM/c^2 in alpha coordinate * (alpha being the impact parameter in phi direction, * positive for approaching side of the disc) * par17 ... beta - position of the cloud centre in GM/c^2 in beta coordinate * (beta being the impact parameter in theta direction, * positive in up direction, i.e. away from the disc) * par18 ... rcloud - radius of the obscuring cloud * - the meaning of cloud is inverted for negative values of * rcloud, i.e. only the radiation behind the cloud is * computed * par19 ... zshift - overall Doppler shift * par20 ... limb - defines fits file with tables * 0 -> for isotropic emission (flux ~ 1) * 1 -> for Laor's limb darkening (flux ~ (1+2.06*mu)) * 2 -> for Haardt's limb brightening (flux ~ ln (1+1/mu)) * par21 ... tab - which reflion table to use * 1 -> reflion (the old one, lower cut-off energy at 1eV, * not good for PhoIndex > 2) * 2 -> reflionx (the new one, lower cut-off energy at 100eV) * par22 ... sw - swich for the way how to compute the refl. spectra * 1 -> use the computed Xi for the interpolation in reflion, * i.e. use proper total incident intensity * with the shifted cut-offs * 2 -> use the Xi correspondent to the computed normalization * of the incident flux, i.e. do not shift the cut-offs * when computing the total incident intensity * par23 ... ntable - defines fits file with tables (0 <= ntable <= 99) * by default only the tables with ntable=80 are * correct for this model * par24 ... nradius - number of grid points in radius * - if negative than the number of radial grid points is * dependent on height as -nradius / height^0.66 * par25 ... division - type of division in radial integration * 0 -> equidistant radial grid (constant linear step) * 1 -> exponential radial grid (constant logarithmic step) * >1 -> mixed radial grid with a constant logarithmic step * in the inner region and with a constant linear * step in the outer region; the total nradius * (par24) number of points is divided in the 3:2 * ratio in these regions; the value of par25 gives * the transition radius between these regions * (in GM/c^2) * -1 -> mixed radial grid with the transition radius at * 2*height * par26 ... nphi - number of grid points in azimuth * par27 ... deltaT - length of the time bin (GM/c^3) * par28 ... nt - number of time subbins per one time bin * par29 ... t1/f1/E1 - the time to be used in xspec for the spectrum * (0 means average spectrum, i.e. divided by the * flare duration) * - the frequency to be used in XSPEC for the energy * dependent Fourier transform * (0 means average values in the range of 0 to the first * wrapping frequency) * - positive values are in sec or Hz * - negative values are in GM/c^3 or (GM/c^3)^(-1) * - if different than par30, the value gives the lower end * of the time/frequency interval of interest * - if same as par30, then the functions are computed for * this value of the time/frequency of interest * - in case of frequency dependent lags it defines the lower * value of the energy band of interest in keV * par30 ... t2/f2/E2 - used only if different than par29 and if par29 is * nonzero * - its value gives the upper end of the time/frequency * interval of interest * - positive values are in sec or Hz * - negative values are in GM/c^3 or (GM/c^3)^(-1) * - in case of frequency dependent lags it defines the upper * value of the energy band of interest in keV * par31 ... E3 - it defines the lower value of the reference energy band * for lag or amplitude energy dependence as well as in * case of frequency dependent lags and amplitudes * - if zero no reference band is used * - if negative, the whole energy band is used as a * reference band for lag-energy spectra (it must be * non-negative in case of lag-frequency dependence) * always excluding the current energy bin * par32 ... E4 - it defines the upper value of the reference energy band * for lag-energy dependence as well as in case of * frequency dependent lags * par33 ... tshift/Af - lag shift for lag-energy dependence in case of * par35=+6 * - multiplicative factor in case of adding empirical hard * lags Af*f^(qf), used for par35=+16 * par34 ... Amp/qf - multiplicative factor for the amplitude-energy * dependence in case of par35=+5 * - powerlaw index in case of adding empirical hard * lags Af*f^(qf), used for par35=+16 * par35 ... photar_sw - function to be stored in the XSPEC photar array * 0 -> spectrum at time defined by par29 and par30, * the following values correspond to energy * dependent Fourier transform at the frequency band * defined by par29 and par30: * -1 -> real part of FT of the relative reflection * -2 -> imaginary part of FT of the relative reflection * -3 -> amplitude of FT of the relative reflection * -4 -> phase of FT of the relative reflection * -5 -> amplitude for the relative reflection * divided by amplitude in the reference energy band * defined by par31 and par32 * -6 -> lag for the relative reflection with respect * to reference energy band defined by par31 and * par32 * 1 -> real part of FT including primary radiation * 2 -> imaginary part of FT including primary radiation * 3 -> amplitude of FT including primary radiation * 4 -> phase of FT including primary radiation * 5 -> amplitude including the primary radiation * divided by amplitude in the reference energy band * defined by par31 and par32 * 6 -> lag diluted by primary radiation with respect * to reference energy band defined by par31 and * par32 * the following values correspond to frequency dependent * Fourier transform for the energy band of interest * defined by par29 and par30: * -11 -> real part of FT of the relative reflection * -12 -> imaginary part of FT of the relative reflection * -13 -> amplitude of FT of the relative reflection * -14 -> phase of FT of the relative reflection * -15 -> amplitude for the relative reflection * divided by amplitude in the reference energy * band defined by par31 and par32 * -16 -> lag for the relative reflection with respect * to reference energy band defined by par31 and * par32 * 11 -> real part of FT including primary radiation * 12 -> imaginary part of FT including primary radiation * 13 -> amplitude of FT including primary radiation * 14 -> phase of FT including primary radiation * 15 -> amplitude including the primary radiation * divided by amplitude in the reference energy * band defined by par31 and par32 * 16 -> lag diluted by primary radiation with respect * to reference energy band defined by par31 and * par32 * * par36 ... nthreads - how many threads should be used for computations * * NOTES: * -> accuracy vs. speed trade off depends mainly on: nradius, nphi, nt, deltaT * -> the normalization of this model has to be unity! * * ------------------------------------------- * KBHlamp_qt.fits * * This file contains pre-calculated values of the functions needed for the * lamp-post models. It is supposed that a primary source of emission is placed * on the axis at a height h from the centre above the Kerr black hole. * The matter in the disc rotates on stable circular (free) orbits above the * marginally stable orbit and it is freely falling below this orbit * where it has the same energy and angular momentum as the matter * which is on the marginally stable orbit. It is assumed that the corona * between the source and the disc is optically thin, therefore ray-tracing * in the vacuum Kerr space-time could be used for computing the functions. * * There are six functions stored in the KBHlamp_qt.fits file as columns. * They are parametrized by the value of the horizon of the black hole (FITS * table extensions), height of the primary source (rows) and either the * inclination angles of the observer (elements) or the radius in GM/c^2 * at which a photon strikes the disc (elements). All these are defined * as vectors at the beginning of the FITS file. The tables are defined for * r_horizon = 1.00, 1.05, ..., 1.95, 2.00, * h-r_horizon = 0.1 - 100 (100 values with exponentially growing step), * inclination = 0.1, 1, 5, 10, 15, ..., 80, 85, 89 and * r-r_horizon = 0.01 - 1000 (100 values with exponentially growing step). * * The functions included are: * dWadWo - amplification of the primary source flux: * = sin(theta_axis_local) / sin(theta_observer) * * dtheta_axis_local/dtheta_observer * delay_a - delay of photon arrival time from the axis to the observer * q - constant of motion defining the photon angular momentum * pr - the radial component of the photon momentum at the disc * dWadSd - part of the amplification of the incident flux on the disc from the * primary source: * = sin(theta_axis_local) * dtheta_axis_local / dtheta_fake, with the * theta_fake defined as tan(theta_fake) = radius_incident / height * one has to multiply by h/(r^2+h^2) to get the full amplification * delay - delay of photon arrival time from the axis to the disc * * The rest of the functions below are computed from the above ones: * g-factor - the ratio of the energy of a photon hitting the disc to the energy * of the same photon when emitted from a primary source, * cosine of the incident angle - an absolute value of the cosine of the local * incident angle between the incident light ray * and local disc normal * azimuthal incident angle - the angle (in radians) between the projection of * the three-momentum of the incident photon into the * disc (in the local rest frame co-moving with the * disc) and the radial tetrad vector. * * The definition of the file KBHlamp_qt.fits: * 0. All of the extensions defined below are binary. * 1. The first extension contains a vector of the horizon values in GM/c^2 * (1.00 <= r_horizon <= 2.00). * 2. The second extension contains a vector of the values of heights h of * a primary source in GM/c^2 (0.1 <= h-r_horizon <= 100). * 3. The third extension contains a vector of the values of the observer * inclinations (0.1 <= inclination <= 89). * 4. The fourth extension contains a vector of the values of the incident * radius (0.01 <= radius-r_horizon <= 1000). * 5. In the following extensions the functions are defined, each extension is * for a particular value of r_horizon, each row is for a particular value of * height and each element is either for a particular value of inclination * (dWadWo, delay_a) or incident radius (the rest of the functions). * 6. Each of these extensions has six columns. In each column, a particular * function is stored - dWadWo, delay_a, q, pr, dWadSd and delay, * respectively. * *============================================================================== * * non-revreberation code: * 27. 8.2009 changed to work with new tables with azimuthal dependence * 12.11.2009 changed to work with new lamp-post tables, we can fit for height * as well now * * reverberation code: * 21. 7.2012 time variation (box flare) added * 30. 9.2015 code transformed into the C-language * 23. 2.2016 major update so that the code is suitable for fitting lag-energy * as well as lag-frequency dependences inside XSPEC * ******************************************************************************/ /******************************************************************************* *******************************************************************************/ #include #include #include #include #include "fitsio.h" //definition of some parameters #define DURATION -10. // duration of the flare (GM/c^3), // if negative, then duration = -duration * deltaT #define TMAX 256. // maximum time for computation of light curves // and dynamic spectra (for response echo) #define NN 8192 // number of time bins for computing FFT #define SMOOTH 0 /******************************************************************************* *******************************************************************************/ #ifdef OUTSIDE_XSPEC #define IFL 1 #define NPARAM 37 // for the energy dependence #define NE 15 #define E_MIN 0.3 #define E_MAX 10. #define NBANDS 5 /* // for a nice energy dependence #define NE 200 #define E_MIN 0.1 #define E_MAX 80 #define NBANDS 5 */ /* // for the frequency dependence #define NE 200 #define E_MIN 1e-5 #define E_MAX 0.01 #define NBANDS 2 */ /* Let's declare variables that are common for the main and KYNrefrev subroutines */ static double ener_low[NBANDS], ener_high[NBANDS]; int main() { void KYNrefrev(double *ear, int ne, const double *param, int ifl, double *photar, double *photer, const char* init); double ear[NE + 1], param[NPARAM], photar[NE], photer[NE]; char initstr[0] = ""; int ie, ia, iinc, ih; //definition of energy band of interest, reference band is defined as the last //one, usually the whole energy range /* ener_low[0] = 0.3; ener_high[0] = 1.; ener_low[1] = 2.; ener_high[1] = 4.; */ ener_low[0] = 0.3; ener_high[0] = 0.8; ener_low[1] = 1.; ener_high[1] = 3.; ener_low[2] = 3.; ener_high[2] = 9.; ener_low[3] = 12.; ener_high[3] = 40.; ener_low[4] = E_MIN; ener_high[4] = E_MAX; //definition of the KYNrefrev parameters param[ 0] = 1.; // a/M param[ 1] = 30.; // thetaO param[ 2] = 1.; // rin param[ 3] = 1.; // ms param[ 4] = 1000.; // rout param[ 5] = 0.; // phi0 param[ 6] = 360.; // dphi param[ 7] = 0.1; // M/M8 param[ 8] = 3.; // height param[ 9] = 2.; // PhoIndex param[10] = 0.001; // Np param[11] = 1.; // NpNr param[12] = 1.; // nH0 param[13] = 0.; // q_n param[14] = 1.; // abun param[15] = -6.; // alpha param[16] = 0.; // beta param[17] = 0.; // rcloud param[18] = 0.; // zshift param[19] = 0.; // limb param[20] = 2.; // tab param[21] = 2.; // sw param[22] = 80.; // ntable param[23] = -4488.; // nradius param[24] = -1.; // division param[25] = 180.; // nphi param[26] = 1.; // deltaT param[27] = 1.; // nt param[28] = 0.; // time/frequency/energy-lower param[29] = 8.3e-4; // time/frequency/energy-upper param[30] = -1.; // reference energy band-lower param[31] = 3.; // reference energy band-upper // the following should be used only for debugging purposes for the case of // abs(photar_sw) > 10 // the energy bands above should be then re-defined to consist of just 2 bands! // for energy band definitions the param[] values are used, while ener_low[] // and ener_high[] are ignored later on! /* param[28] = ener_low[0]; // time/frequency/energy-lower param[29] = ener_high[0]; // time/frequency/energy-upper param[30] = ener_low[1]; // reference energy band-lower param[31] = ener_high[1]; // reference energy band-upper */ param[32] = 0.; // lag shift or multiplicative factor for hard lags param[33] = 1.; // amplitude multiplicative factor or power-law index for hard lags param[34] = 6.; // photar_sw param[35] = 4.; // nthreads param[36] = 1.; // norm for (ie = 0; ie <= NE; ie++) { // ear[ie] = E_MIN + ie * (E_MAX-E_MIN) / NE; ear[ie] = E_MIN * pow(E_MAX / E_MIN, ((double) ie) / NE); } //for (ia=0;ia<=1;ia++){ for (ia=1;ia<=1;ia++){ // param[0] = (double) ia; // for (iinc=30;iinc<=60;iinc+=30){ for (iinc=30;iinc<=30;iinc+=30){ // param[1] = (double) iinc; // for (ih=1;ih<=5;ih++){ for (ih=2;ih<=2;ih++){ if(ih==1)param[8]=1.5; // else if(ih==2)param[8]=3.; else if(ih==3)param[8]=6.; else if(ih==4)param[8]=15.; else if(ih==5)param[8]=30.; // for (ih=1;ih<=20;ih++){ // param[8] = 1.5 * (100./1.5)**((ih-1.)/19.); if( param[8] >= (1.1+sqrt(1.-param[0]*param[0])) ) KYNrefrev(ear, NE, param, IFL, photar, photer, initstr); } } } return(0); } #endif /******************************************************************************* *******************************************************************************/ #define MPC_2 1.05e-49 #define ERG 6.241509e8 // rg is gravitational unit GM/c^2 for M = 10^8*Msun and rg2=rg^2 #define RG2 2.1819222e26 // if Ledd is changed one need to recalculate also logxi_norm0 below!!! // Ledd is in erg (not W) and multiplied by 10^8 due to (M / (10^8*Msun)) scale #define LEDD 1.26e46 #define HUBBLE 70. #define CLIGHT 299792.458 // sec = G*10^8*Msun/c^3 #define SEC 492.568 // logxi_norm=alog10(Ledd*(M/(10^8*Msun))/(rg*m8)^2/nH)= // alog10(1.26e46/(1.477e5*1e8)^2/1e15)+alog10(mass) // logxi_norm0+alog10(mass) #define LOGXI_NORM0 4.761609554 #define PI 3.14159265359 #define PI2 6.2831853071795865 #define LAMP "KBHlamp_qt.fits" #define REFLION1 "reflion.mod" #define REFLION2 "reflionx.mod" #define REFLIONX_NORM 1.e20 #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 all_energy=0; /* Let's declare variables that are common for the main and emissivity subroutines */ static float *xi=NULL, *logxi=NULL; static double *gfac=NULL, *transf_d=NULL, *energy1=NULL, *flux1=NULL; //static double *cosin=NULL, *phiph=NULL; static double h, gamma0, nH0, qn, mass2, am2, r_plus, Np, logxi_norm, dt, flare_duration_rg, flare_duration_sec; static long nxi; static int sw, limb, nt_ratio; 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 void cutoffpl(double *ear, int ne, double *param, double *photar); void KYNrefrev(double *ear_xspec, int ne_xspec, const double *param, int ifl, double *photar, double *photer, const char* init) { extern int ide(const double *ear, const int ne, const long nt, double *far, double *qar, double *uar, double *var, const double *ide_param, void (*emissivity)(), const int ne_loc); extern void fft_reverberation(const double *ear, int ne, double *photar, double frequency1, double frequency2, int photar_sw, double *time, 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); void emis_KYNrefrev(double** ear_loc, const int ne_loc, const long nt, double *far_loc, double *qar_loc, double *uar_loc, double *var_loc, const double r, const double phi, const double cosmu, const double phiphoton, const double alpha_o, const double beta_o, const double delay, const double g); double incgamma(double a, double x); /* Let's declare static variables whose values should be remembered for next run in XSPEC */ static char kydir[255] = ""; static char pname[128] = "KYDIR", pkyLxLamp[128] = "KYLxLamp"; static char pkyxiin[128] = "KYXIin", pkyxiout[128] = "KYXIout"; static char pkyRefl[128] = "KYRefl", pkyfwrap[128] = "KYfwrap"; static long nrh, nh, nincl, nener; static float *r_horizon, *height, *incl, *dWadWo, *delay_a, *q, *pr, *dWadSd, *delay, *abun, *gam, *emission, *energy0; static long ne_loc, nabun, ngam; static double *energy2, *flux0, transf_o, del_a; static double h_rh_old = -1., gam_old = -1., abun_old = -1., thetaO_old = -1., am_old = -1.; static int rix_old = -1, first_h = 1, first_rix = 1; 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]; double ide_param[25]; double ear_short[4]; double *time=NULL; double *ear=NULL, *far=NULL, *qar=NULL, *uar=NULL, *var=NULL; double *far_final=NULL; //double far_prim[ne_xspec]; //double flux_prim=0.,spectrum_prim[ne_xspec]; double *flux=NULL, *spectrum=NULL, *far_prim=NULL; double flux_prim=0., *spectrum_prim=NULL; #ifndef OUTSIDE_XSPEC #define NBANDS 2 double ener_low[NBANDS], ener_high[NBANDS]; #endif double *flux_bands=NULL, flux_bands_prim[NBANDS]; int ne=0., nbands=0; long nn=NN, ntmax=0; double en1=0., en2=0., en3=0., en4=0.; double lag_shift=0., ampl_ampl=1.; double ttmp, ttmp1, utmp, utmp1, vtmp, vtmp1, y1, y2, y3, y4, y5, y6, y7, y8; double pr_final, pom, pom1, pom2, pom3; double r, r2, delta, ULt, rms, tmp1, Ut, U_r, UrU_r, Lms, Ems; //double q_final, U_phi, Ur; double am, thetaO, rin=0., rout, h_rh, elow, ehigh; double mass, Anorm, Dnorm, g_L, E0, Lx; double flux_prim_tot, flux_refl_tot, refl_ratio, NpNr; double zzshift; double abundance, lensing, gfactor0, ionisation; double time1=0., time1_rg=0., time2=0., time2_rg=0., frequency1=0., frequency2=0., tot_time_sec, deltaT, Af=0., qf=0., fmin, fmax, dt0; double Tmax, Tmax_final, freq_min, freq_max; double sumt, sumf, sumt2, sumtf, qtime, Atime; long nt, iabun, igam, ixi, it; int imin, imax, irh0, ih0, ith0, ir0, iabun0, igam0; int nt0, ms, nspec, polar, stokes, rix, photar_sw, npoints; int i, ie, je, quit, it0, itn, iband; //int itmin; //double photar1_xspec[ne_xspec]; photar1_short[3], *photar1=NULL; //double gf_xspec[ne_xspec+1], gf_short[4], *gf=NULL; double *photar1=NULL, *gf=NULL; // the following are needed for cut-off power-law taken from XSPEC // //double ear1[ne + 1], param1[2]; //double *ear1=NULL, param1[2]; // these are needed to work with a fits file... fitsfile *fptr; char tables_file[255]; int hdutype = 2; int colnum = 1; long frow = 1, felem = 1, nelems, nrow; float float_nulval = 0.; int nelements, nelements1, nelements2; int ihorizon, irow, anynul, status = 0; // Let's initialize parameters for subroutine ide() sprintf(cparam, "%3ld_%02ld_%04ld", lround(100.*(1.+sqrt(1.-param[0]*param[0]))), lround(param[1]),lround(10.*param[8])); // am - black hole angular momentum ide_param[0] = param[0]; am = ide_param[0]; am2 = am * am; pom1 = pow(1. + am, 1. / 3.); pom2 = pow(1. - am, 1. / 3.); pom3 = pow(1. - am2, 1. / 3.); pom = 1. + pom3 * (pom1 + pom2); pom1 = sqrt(3. * am2 + pom * pom); if (am >= 0) rms = 3. + pom1 - sqrt((3. - pom) * (3. + pom + 2. * pom1)); else rms = 3. + pom1 + sqrt((3. - pom) * (3. + pom + 2. * pom1)); r_plus = 1. + sqrt(1. - am2); // thetaO - observer inclination ide_param[1] = param[1]; thetaO = ide_param[1]; // rin - inner edge of non-zero disc emissivity ide_param[2] = param[2]; // ms - whether to integrate from rin or rms ide_param[3] = param[3]; // rout - outer edge of non-zero disc emissivity ide_param[4] = param[4]; rout = param[4]; //set up inner and outer radii: ms = (int) ide_param[3]; if(ms!=0 && ms!=1 && ms!=2){ xs_write("kynrefrev: ms must be 0, 1 or 2",5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } if(ms==1){ if(ide_param[2]= 0.) if (h < r_plus){ h_rh = 0.; h=r_plus;} else h_rh = h - r_plus; else { xs_write("kynrefrev: height has to be positive.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } // PhoIndex - power-law energy index of the lamp emission gamma0 = param[9]; // Np - dN/dt/dOmega primary isotropic flux in Eddington luminosity // as seen by the disc Np = param[10]; // Np:Nr - ratio of the primary to the reflected normalization NpNr = param[11]; if(NpNr > 0.) Np /= NpNr; // nH0 - density profile normalization in 10^15 cm^(-3) nH0 = param[12]; // q_n - radial power-law density profile qn = param[13]; // Fe abundance (in solar abundance) abundance = param[14]; // zshift - overall Doppler shift if (param[18] > 0.) { ide_param[12] = param[18]; Dnorm = pow(HUBBLE / CLIGHT / param[18], 2); }else if (param[18] < 0.) { ide_param[12] = 0.; Dnorm = pow(-HUBBLE / CLIGHT / param[18], 2); }else { ide_param[12] = 0.; Dnorm = 1.; } // zzshift - multiplication factor for gfac from zshift needed for primary zzshift=1.0/(1.0+ide_param[12]); // limb - table model (defines fits file with tables) limb = (int) param[19]; if ((limb < 0) || (limb > 2)) { xs_write("kynrefrev: limb must be >= 0 and <= 2.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } // tab - which reflion table to use rix = (int) param[20]; if (rix == 1) { E0 = 0.001; sprintf(reflion, REFLION1); nener = 500; elow = 0.01; } else{ E0 = 0.1; sprintf(reflion, REFLION2); nener = 375; elow = 0.1; } ehigh = EC; // sw - swich for the way how to compute the refl. spectra sw = (int) param[21]; // edivision - type of division in local energies (0-equidistant, 1-exponential) ide_param[14] = 1.; // variability type ide_param[15]=1.; //photar_sw photar_sw = (int) param[34]; // let's set up the deltaT, ntbin and nn params according to the frequency range // needed dt0 = param[26]; //let's keep the T_tail constant, i.e. nt0 * dt0 = TMAX nt0 = (int) ceil( TMAX / dt0 ); if(param[4] != 1000.){ nt0 = (int) ceil( 1.2 / dt0 * ( rout*(1+sin(thetaO/180.*PI)) + h*cos(thetaO/180.*PI) + h*h/rout/2. ) ); if(nt0 > nn) nn = (int) exp2( ceil( log2( nt0 ) ) ); } Tmax = nt0 * dt0; Tmax_final = nn * dt0; fmin = 1. / (nn * dt0 * SEC * mass); fmax = 1. / (2. * dt0 * SEC * mass); if(photar_sw){ if(abs(photar_sw) <= 10){ freq_min = ( param[28] >=0. ? param[28] : param[28] / ( - SEC * mass)); freq_max = ( param[29] >=0. ? param[29] : param[29] / ( - SEC * mass)); if(freq_max < freq_min)freq_max=freq_min; }else if(abs(photar_sw)>10){ if( ear_xspec[0] <= 0. ){ sprintf(errortxt, "kynrefrev: lower frequency boundary has to be larger than zero."); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } freq_min = ear_xspec[0]; freq_max = ear_xspec[ne_xspec]; } if( freq_min > 0. ){ if( log10( freq_max / freq_min) > 5.){ sprintf(errortxt, "kynrefrev: frequency cannot span more than 5 orders of magnitude!"); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } if( freq_max > fmax ){ dt0 = 1. / freq_max / (2. * SEC * mass); if(dt0 < 0.1){ sprintf(errortxt, "kynrefrev: frequency is too high, much above %lg Hz, please, set it below %lg Hz.", fmax, 1. / (0.2 * SEC * mass)); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } sprintf(errortxt, "kynrefrev: frequency is quite high, above %lg Hz. Do you really want it?", fmax); xs_write(errortxt, 5); sprintf(errortxt, " time step changed to %lg Rg/c.", dt0); xs_write(errortxt, 5); sprintf(errortxt, " computation will be longer!"); xs_write(errortxt, 5); nt0 = (int) ceil(Tmax / dt0); nn = (int) exp2( ceil( log2( Tmax_final / dt0 ) ) ); } if( freq_min < 1. / (nn * dt0 * SEC * mass) ){ if(freq_min < 0.1 * fmin){ sprintf(errortxt, "kynrefrev: frequency is too low, much below %lg Hz, please, set it above %lg Hz.", fmin, 0.1 * fmin); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } nn = (int) exp2( ceil( log2( 1. / freq_min / (dt0 * SEC * mass) ) ) ); sprintf(errortxt, "kynrefrev: frequency is quite low, below %lg Hz. Do you really want it?", fmin); xs_write(errortxt, 5); } } } // duration of the flare (if negative then 10 times the step in time) if( DURATION < 0)flare_duration_rg = - DURATION * dt0; else flare_duration_rg = DURATION; flare_duration_sec = flare_duration_rg * SEC * mass; // nt_ratio (nt) nt_ratio= (int) param[27]; // nt nt = ((long) nt_ratio) * ((long) nt0); if(nt==1){ xs_write("kynrefrev: nt = nt_final * nt_ratio must be larger than 1", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } // dt dt=dt0/nt_ratio; ide_param[16]=dt; // time/frequency if (photar_sw == 0){ time1 = param[28]; if(time1 >= 0.) time1_rg = time1 / SEC / mass; else{ time1_rg = -time1; time1 = -time1 * SEC * mass; } if (time1_rg > (nt-1)*dt) { sprintf(errortxt,"kynrefrev: time has to be smaller than %lg sec or %lg Rg.", (nt-1)*dt*SEC*mass, (nt-1)*dt); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } time2 = param[29]; if(time2 >= 0.) time2_rg = time2 / SEC / mass; else{ time2_rg = -time2; time2 = -time2 * SEC * mass; } if (time2_rg > (nt-1)*dt) { sprintf(errortxt,"kynrefrev: time has to be smaller than %lg sec or %lg Rg.", (nt-1)*dt*SEC*mass, (nt-1)*dt); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } }else if(abs(photar_sw) <= 10){ fmin = 1. / (nn * dt0 * SEC * mass); fmax = 1. / (2. * dt0 * SEC * mass); frequency1 = param[28]; if(frequency1 < 0.) frequency1 /= ( - SEC * mass); if (frequency1 != 0. && frequency1 < fmin ) { sprintf(errortxt,"kynrefrev: frequency has to be larger than %lg Hz or %lg / Rg.", fmin, 1. / (nn * dt0)); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } if (frequency1 != 0. && frequency1 > fmax ) { sprintf(errortxt,"kynrefrev: frequency has to be smaller than %lg Hz or %lg / Rg.", fmax, 1. / (2. * dt0)); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } frequency2 = param[29]; if(frequency2 < 0.) frequency2 /= ( - SEC * mass); if (frequency1 != 0. && frequency2 > fmax ) { sprintf(errortxt,"kynrefrev: frequency has to be smaller than %lg Hz or %lg / Rg.", fmax, 1. / (2. * dt0)); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } }else{ //definition of the energy band of interest for frequency dependent Fourier transform en1 = param[28]; en2 = param[29]; if ( en1 < 0. || en1 >= en2 ) { sprintf(errortxt,"kynrefrev: wrong definition of energy band of interest (0 <= E1 < E2)."); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } //check the frequency range fmin = 1. / (nn * dt0 * SEC * mass); fmax = 1. / (2. * dt0 * SEC * mass); if( ear_xspec[0] < fmin || ear_xspec[0] > fmax || ear_xspec[ne_xspec] < fmin || ear_xspec[ne_xspec] > fmax ){ sprintf(errortxt,"kynrefrev: frequency has to be in the range of %lg to %lg Hz.", fmin, fmax); xs_write(errortxt, 5); } } //definition of the reference energy band for both frequency dependent lags and //amplitudes as well as energy dependent lags and amplitudes en3 = param[30]; en4 = param[31]; if ( ( (abs(photar_sw) == 5 || abs(photar_sw) == 6) && (en3 > 0. && en3 >= en4) ) || ( (abs(photar_sw) == 15 || abs(photar_sw) == 16) && (en3 <= 0. || en3 >= en4) ) ) { sprintf(errortxt,"kynrefrev: wrong definition of the reference energy band."); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } #ifdef OUTSIDE_XSPEC if( abs(photar_sw) > 10 ){ nbands=1; ener_low[0]=en1; ener_high[0]=en2; }else nbands = NBANDS; if(abs(photar_sw) == 15 || abs(photar_sw) == 16){ nbands++; ener_low[nbands-1] = en3; ener_high[nbands-1] = en4; } if(ener_low[nbands-1] == ear_xspec[0] && ener_high[nbands-1] == ear_xspec[ne_xspec])all_energy=1; else all_energy=0; #else if( abs(photar_sw) > 10 ){ nbands=1; ener_low[0]=en1; ener_high[0]=en2; }else nbands = 0; if(abs(photar_sw) == 5 || abs(photar_sw) == 6 || abs(photar_sw) == 15 || abs(photar_sw) == 16){ if(en3 < 0){ nbands++; ener_low[nbands-1] = ear_xspec[0]; ener_high[nbands-1] = ear_xspec[ne_xspec]; all_energy = 1; }else if( en3 > 0 && en3 < en4 ){ nbands++; ener_low[nbands-1] = en3; ener_high[nbands-1] = en4; all_energy = 0; } } #endif //let's define the energy for ide integration for energy dependent lags as well //as for frequency dependent lags of energy bands if(abs(photar_sw) <= 10){ ear=ear_xspec; ne=ne_xspec; }else{ ear=ear_short; //let's define ne and new ear energy array if( (en1==en3 && en2==en4) || abs(photar_sw) < 15){ ne=1; ear[0]=en1; ear[1]=en2; }else if(en1==en3){ ne=2; ear[0]=en1; if(en2 nn ? nt : nn); if ( (time = (double *) malloc( ntmax * sizeof(double)) ) == NULL || (far = (double *) malloc( ne * ntmax * sizeof(double)) ) == NULL || (far_final = (double *) malloc( ne * ntmax * sizeof(double))) == NULL || (flux = (double *) malloc( ntmax * sizeof(double))) == NULL || (flux_bands = (double *) malloc( nbands * nn * sizeof(double))) == NULL || (spectrum = (double *) malloc( ne * sizeof(double))) == NULL || (far_prim = (double *) malloc( ne * sizeof(double))) == NULL || (spectrum_prim = (double *) malloc( ne * sizeof(double))) == NULL || (photar1 = (double *) malloc( ne * sizeof(double))) == NULL || (gf = (double *) malloc( (ne+1) * sizeof(double))) == NULL ) { xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } //if ((ear1 = (double *) malloc( (ne+1) * sizeof(double))) == NULL) { // xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5); // for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; // goto error; //} if(photar_sw == 5 || photar_sw == 6){ lag_shift = param[32]; ampl_ampl = param[33]; } if(photar_sw == 16){ Af = param[32]; qf = param[33]; } // polar - whether we need value of change in polarization angle (0-no,1-yes) //stokes = (int) param[31]; stokes = 0; if ((stokes < 0) || (stokes > 6)) { xs_write("kynrefrev: stokes has to be 0-6.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } polar = 0; if (stokes > 0) polar = 1; ide_param[17] = polar; //delay_r ide_param[18]=-1.; // delay_phi is defined later on //ide_param[19]=del_a // number of threads for multithread computations ide_param[20] = param[35]; // alpha - position of the cloud in alpha impact parameter (in GM/c^2) ide_param[21] = param[15]; // beta - position of the cloud in beta impact parameter (in GM/c^2) ide_param[22] = param[16]; // rcloud - radius of the cloud (in GM/c^2) ide_param[23] = param[17]; //whether the flux defined in emissivity subroutine is local one (0) or the //observed one (1) ide_param[24] = 0.; // check if normalization parameter is equal to 1. if ((param[18] != 0.) && (param[36] != 1.)) { xs_write("kynrefrev: the normalisation parameter par34 should be frozen to unity.", 5); } // Let's clear the flux array for (ie = 0; ie < ne; ie++) for (it = 0; it < nt; it++) far[ie + ne * it] = 0; /******************************************************************************/ // Let's read the lamp post tables if(first_h) { // The status parameter must always be initialized. status = 0; // Open the FITS file for readonly access // - if set try KYDIR directory, otherwise look in the working directory // or in the xspec directory where tables are usually stored... sprintf(kydir, "%s", FGMSTR(pname)); if (strlen(kydir) == 0) sprintf(tables_file, "./%s", LAMP); else if (kydir[strlen(kydir) - 1] == '/') sprintf(tables_file, "%s%s", kydir, LAMP); else sprintf(tables_file, "%s/%s", kydir, LAMP); // Let's read the 'KBHlamp_q' fits file // The status parameter must always be initialized. status = 0; ffopen(&fptr, tables_file, READONLY, &status); if (status) { sprintf(tables_file, "%s%s", FGMODF(), LAMP); status = 0; ffopen(&fptr, tables_file, READONLY, &status); } if (status) { if (status) ffrprt(stderr, status); ffclos(fptr, &status); xs_write("\nkynrefionx: set the KYDIR to the directory with the KY tables.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } // Let's read tables (binary tables => hdutype=2) // Move to the extension 'r_horizon' and read its values ffmrhd(fptr, 1, &hdutype, &status); ffgnrw(fptr, &nrh, &status); //****************************************************************************** // fprintf(stdout,"nrh = %ld\n",nrh); //****************************************************************************** // Allocate memory for r_horizon... if ((r_horizon = (float *) malloc(nrh * sizeof(float))) == NULL) { xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } // Read the data in the 'r_horizon' table nelems = nrh; // FFGCV reads the VALUES from the first column. ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, r_horizon, &anynul, &status); //****************************************************************************** // for ( i=0; i height[nh-1]) { sprintf(errortxt, "kynrefrev: the height must be lower than or equal to %f.", height[nh - 1] + r_plus); xs_write(errortxt, 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } if (h_rh < height[0]) { sprintf(errortxt, "kynrefrev: the height is too low, we set it to %f.", height[0] + r_plus); xs_write(errortxt, 5); h_rh = height[0]; h = h_rh + r_plus; } // nrad - number of grid points in radius if(param[23] > 0) ide_param[7] = param[23]; else ide_param[7] = - param[23] * pow(h,-0.66); // division - type of division in r integration (0-equidistant, 1-exponential) if(param[24]==-1.)ide_param[8]=2.*h; else ide_param[8]=param[24]; // Let's interpolate the tables to desired spin and height if ((am != am_old) || (h_rh != h_rh_old) || (thetaO != thetaO_old)) { // given am->r_plus, find the corresponding index in r_horizon[]: imin = 0; imax = nrh; irh0 = nrh / 2; while ((imax - imin) > 1) { if (r_plus >= r_horizon[irh0 - 1]) imin = irh0; else imax = irh0; irh0 = (imin + imax) / 2; } if (irh0 == 0) irh0 = 1; //if ((imax == nrh) && (r_plus > r_horizon[nrh - 1])) irh0 = nrh; ttmp = (r_plus - r_horizon[irh0 - 1]) / (r_horizon[irh0] - r_horizon[irh0 - 1]); ttmp1 = 1. - ttmp; // given h, find the corresponding index in height[]: imin = 0; imax = nh; ih0 = nh / 2; while ((imax - imin) > 1) { if (h_rh >= height[ih0 - 1]) imin = ih0; else imax = ih0; ih0 = (imin + imax) / 2; } if (ih0 == 0) ih0 = 1; //if ((imax == nh) && (h_rh > height[nh - 1])) ih0 = nh; utmp = (h_rh - height[ih0 - 1]) / (height[ih0] - height[ih0 - 1]); utmp1 = 1. - utmp; // given thetaO, find the corresponding index in incl[]: imin = 0; imax = nincl; ith0 = nincl / 2; while ((imax - imin) > 1) { if (thetaO >= incl[ith0 - 1]) imin = ith0; else imax = ith0; ith0 = (imin + imax) / 2; } if (ith0 == 0) ith0 = 1; //if ((imax == nincl) && (thetaO > incl[nincl - 1])) ith0 = nincl; vtmp = (thetaO - incl[ith0 - 1]) / (incl[ith0] - incl[ith0 - 1]); vtmp1 = 1. - vtmp; // transfer function from the axis to the observer y1 = dWadWo[ith0 - 1 + nincl * (ih0 - 1) + nincl * nh * (irh0 - 1)]; y2 = dWadWo[ith0 - 1 + nincl * (ih0 - 1) + nincl * nh * irh0]; y3 = dWadWo[ith0 - 1 + nincl * ih0 + nincl * nh * irh0]; y4 = dWadWo[ith0 - 1 + nincl * ih0 + nincl * nh * (irh0 - 1)]; y5 = dWadWo[ith0 + nincl * (ih0 - 1) + nincl * nh * (irh0 - 1)]; y6 = dWadWo[ith0 + nincl * (ih0 - 1) + nincl * nh * irh0]; y7 = dWadWo[ith0 + nincl * ih0 + nincl * nh * irh0]; y8 = dWadWo[ith0 + nincl * ih0 + nincl * nh * (irh0 - 1)]; transf_o = (vtmp1 * (utmp1 * (ttmp1 * y1 + ttmp * y2) + utmp * (ttmp * y3 + ttmp1 * y4)) + vtmp * (utmp1 * (ttmp1 * y5 + ttmp * y6) + utmp * (ttmp * y7 + ttmp1 * y8))); // delay from the axis to the observer y1 = delay_a[ith0 - 1 + nincl * (ih0 - 1) + nincl * nh * (irh0 - 1)]; y2 = delay_a[ith0 - 1 + nincl * (ih0 - 1) + nincl * nh * irh0]; y3 = delay_a[ith0 - 1 + nincl * ih0 + nincl * nh * irh0]; y4 = delay_a[ith0 - 1 + nincl * ih0 + nincl * nh * (irh0 - 1)]; y5 = delay_a[ith0 + nincl * (ih0 - 1) + nincl * nh * (irh0 - 1)]; y6 = delay_a[ith0 + nincl * (ih0 - 1) + nincl * nh * irh0]; y7 = delay_a[ith0 + nincl * ih0 + nincl * nh * irh0]; y8 = delay_a[ith0 + nincl * ih0 + nincl * nh * (irh0 - 1)]; del_a = (vtmp1 * (utmp1 * (ttmp1 * y1 + ttmp * y2) + utmp * (ttmp * y3 + ttmp1 * y4)) + vtmp * (utmp1 * (ttmp1 * y5 + ttmp * y6) + utmp * (ttmp * y7 + ttmp1 * y8))); if ((am != am_old) || (h_rh != h_rh_old)) { for (i = 0; i < nradius; i++) { /* // q from the axis to the disc y1 = q[i + nradius * (ih0 - 1) + nradius * nh * (irh0 - 1)]; y2 = q[i + nradius * (ih0 - 1) + nradius * nh * irh0]; y3 = q[i + nradius * ih0 + nradius * nh * irh0]; y4 = q[i + nradius * ih0 + nradius * nh * (irh0 - 1)]; q_final = utmp1 * (ttmp1 * y1 + ttmp * y2) + utmp * (ttmp * y3 + ttmp1 * y4); */ // pr at the disc y1 = pr[i + nradius * (ih0 - 1) + nradius * nh * (irh0 - 1)]; y2 = pr[i + nradius * (ih0 - 1) + nradius * nh * irh0]; y3 = pr[i + nradius * ih0 + nradius * nh * irh0]; y4 = pr[i + nradius * ih0 + nradius * nh * (irh0 - 1)]; pr_final = utmp1 * (ttmp1 * y1 + ttmp * y2) + utmp * (ttmp * y3 + ttmp1 * y4); // temporary variables r = r_plus + radius[i]; r2 = r * r; delta = r2 - 2. * r + am2; ULt = sqrt((h * h + am2) / (h * h - 2. * h + am2)); if (r >= rms) { tmp1 = sqrt(r2 - 3. * r + 2. * am * sqrt(r)); Ut = (r2 + am * sqrt(r)) / r / tmp1; // U_phi = (r2 + am2 - 2. * am * sqrt(r)) / sqrt(r) / tmp1; U_r = 0.; // Ur = 0.; UrU_r = 0.; } else { tmp1 = sqrt(rms * (rms - 3.) + 2. * am * sqrt(rms)); Lms = (rms * rms + am2 - 2. * am * sqrt(rms)) / sqrt(rms) / tmp1; Ems = (rms * (rms - 2.) + am * sqrt(rms)) / rms / tmp1; Ut = (Ems * (r * (r2 + am2) + 2. * am2) - 2. * am * Lms) / r / delta; // U_phi = Lms; UrU_r = -1. + ((r2 + am2 + 2. * am2 / r) * Ems * Ems - 4. * am / r * Ems * Lms - (1. - 2. / r) * Lms * Lms) / delta; if (UrU_r < 0.) UrU_r = 0.; U_r = -sqrt(UrU_r / delta) * r; // Ur = -sqrt(delta * UrU_r) / r; } tmp1 = Ut - pr_final * U_r; // gfactor from the axis to the disc gfac[i] = tmp1 / ULt; // cosin at the disc // cosin[i] = q_final / r / tmp1; // phip_i at the disc // phiph[i] = atan2(-U_phi, r * (pr_final - Ur *tmp1)); // dWadSd from the axis to the disc y1 = dWadSd[i + nradius * (ih0 - 1) + nradius * nh * (irh0 - 1)]; y2 = dWadSd[i + nradius * (ih0 - 1) + nradius * nh * irh0]; y3 = dWadSd[i + nradius * ih0 + nradius * nh * irh0]; y4 = dWadSd[i + nradius * ih0 + nradius * nh * (irh0 - 1)]; transf_d[i] = utmp1 * (ttmp1 * y1 + ttmp * y2) + utmp * (ttmp * y3 + ttmp1 * y4); // delay from the axis to the disc y1 = delay[i + nradius * (ih0 - 1) + nradius * nh * (irh0 - 1)]; y2 = delay[i + nradius * (ih0 - 1) + nradius * nh * irh0]; y3 = delay[i + nradius * ih0 + nradius * nh * irh0]; y4 = delay[i + nradius * ih0 + nradius * nh * (irh0 - 1)]; del[i] = utmp1 * (ttmp1 * y1 + ttmp * y2) + utmp * (ttmp * y3 + ttmp1 * y4); } } //****************************************************************************** // fprintf(stdout,"%f %f\n", thetaO, transf_o); // for(i = 0; i < nradius; i++) // fprintf(stdout,"%d %f %f %f %f %f\n", i, radius[i], gfac[i], cosin[i], // phiph[i], transf_d[i]); //****************************************************************************** } /******************************************************************************/ // Let's read the reflionx tables if ((strcmp(kydir, FGMSTR(pname)) != 0) || (rix != rix_old) || first_rix) { sprintf(text, "kynrefrev: initializing %s tables...", reflion); xs_write(text, 5); xs_write("Ross & Fabian (2005), MNRAS, 358, 211",5); // The status parameter must always be initialized. status = 0; // Open the FITS file for readonly access // - if set try KYDIR directory, otherwise look in the working directory // or in the xspec directory where tables are usually stored... sprintf(kydir, "%s", FGMSTR(pname)); if (strlen(kydir) == 0) sprintf(tables_file, "./%s", reflion); else if (kydir[strlen(kydir) - 1] == '/') sprintf(tables_file, "%s%s", kydir, reflion); else sprintf(tables_file, "%s/%s", kydir, reflion); // Let's read the reflion(x) tables // The status parameter must always be initialized. status = 0; ffopen(&fptr, tables_file, READONLY, &status); if (status) { sprintf(tables_file, "%s%s", FGMODF(), reflion); status = 0; ffopen(&fptr, tables_file, READONLY, &status); } if (status) { if (status) ffrprt(stderr, status); ffclos(fptr, &status); xs_write("\nkynrefionx: 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; } // Let's read tables (binary tables => hdutype=2) // Move to the first extension ('PARAMETERS') ffmrhd(fptr, 1, &hdutype, &status); // Read the values of parameters --> abundance, photon index and ionisation nelems = 1; colnum = 9; frow = 1; ffgcv(fptr, TLONG, colnum, frow, felem, nelems, &float_nulval, &nabun, &anynul, &status); frow = 2; ffgcv(fptr, TLONG, colnum, frow, felem, nelems, &float_nulval, &ngam, &anynul, &status); frow = 3; ffgcv(fptr, TLONG, colnum, frow, felem, nelems, &float_nulval, &nxi, &anynul, &status); nspec = nabun * ngam * nxi; //****************************************************************************** // fprintf(stdout,"%d %d %d\n", nabun, ngam, nxi); //****************************************************************************** // Allocate memory for abun, gam and xi... if( (abun = (float *) malloc(nabun * sizeof(float))) == NULL || (gam = (float *) malloc( ngam * sizeof(float))) == NULL || (xi = (float *) malloc( nxi * sizeof(float))) == NULL || (logxi = (float *) malloc( nxi * sizeof(float))) == NULL ) { xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } colnum = 10; frow = 1; nelems = nabun; ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, abun, &anynul, &status); frow = 2; nelems = ngam; ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, gam, &anynul, &status); frow = 3; nelems = nxi; ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, xi, &anynul, &status); for (i = 0; i < nxi; i++) logxi[i] = log10(xi[i]); //****************************************************************************** // for(i = 0; i < 20; i++) // fprintf(stdout,"%4d %12.1f %12.1f %12.1f %12.1f\n", i, abun[i], gam[i], // xi[i], logxi[i]); //****************************************************************************** // Move to the second extension 'ENERGIES' and read energy values ffmrhd(fptr, 1, &hdutype, &status); ffgnrw(fptr, &ne_loc, &status); ne_loc++; //****************************************************************************** // fprintf(stdout,"%d\n", ne_loc); //****************************************************************************** // Allocate memory for energy... if ((energy0 = (float *) malloc(ne_loc * sizeof(float))) == NULL) { xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } nelems = ne_loc - 1; colnum = 1; frow = 1; ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, energy0, &anynul, &status); nelems = 1; colnum = 2; frow = ne_loc - 1; ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, &energy0[ne_loc - 1], &anynul, &status); //****************************************************************************** // for(i = 0; i < ne_loc; i++) fprintf(stdout,"%d %f\n", i, energy0[i]); //****************************************************************************** // Allocate memory for emission... if ((emission = (float *) malloc(nspec * (ne_loc - 1) * sizeof(float))) == NULL) { xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } // Let's read the tables for emission ffmrhd(fptr, 1, &hdutype, &status); // to read the file only once we have to read in blocks (all columns // from the extension are put to buffer together) // let's find out how many rows are going to be read into the buffer ffgrsz(fptr, &nrow, &status); nelements = nrow * (ne_loc - 1); for (irow = 0; irow < nspec; irow += nrow) { // the last block to read may be smaller: if ((nspec - irow) < nrow) nelements = (nspec - irow) * (ne_loc - 1); iabun = irow / (ngam * nxi) + 1; igam = (irow - (iabun - 1) * ngam * nxi) / nxi + 1; ixi = irow + 1 - (iabun - 1) * ngam * nxi - (igam - 1) * nxi; ffgcv(fptr, TFLOAT, 2, irow + 1, 1, nelements, &float_nulval, &emission[(ne_loc - 1) * (ixi - 1) + (ne_loc - 1) * nxi * (igam - 1) + (ne_loc - 1) * nxi * ngam * (iabun - 1)], &anynul, &status); } // The FITS file must always be closed before exiting the program. ffclos(fptr, &status); /******************************************************************************* iabun = 7; igam = 5; ixi = 3; fprintf(stdout, "%f %f %f %f\n", abun[iabun - 1], gam[igam - 1], pow(10, xi[ixi - 1]), xi[ixi - 1]); for(ie = 0; ie < ne_loc - 1; ie++) fprintf(stdout, "%d %f %f\n", ie, energy0[ie], emission[ie + (ne_loc - 1) * (ixi - 1) + (ne_loc - 1) * nxi * (igam - 1) + (ne_loc - 1 ) * nxi * ngam * (iabun - 1)]); *******************************************************************************/ // Allocate memory for local flux... if( (flux0 = (double *) malloc(ne_loc * nxi * sizeof(double))) == NULL || (flux1 = (double *) malloc( nener * nxi * sizeof(double))) == NULL || (energy1 = (double *) malloc( (nener + 1) * sizeof(double))) == NULL || (energy2 = (double *) malloc( (nener + 1) * sizeof(double))) == NULL ) { xs_write("kynrefrev: Failed to allocate memory for tmp arrays.", 5); for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } first_rix = 0; xs_write("------------------------------------------------", 5); } // end of reading reflionx fits file............................................ /******************************************************************************/ // Let's interpolate the emission to get rid of the non-radial dependences, i.e. // interpolate in abundance and gamma0 if ((rix != rix_old) || (strcmp(kydir, FGMSTR(pname)) != 0) || ((abun_old == -1.) || (abundance != abun_old)) || ((gam_old == -1) || (gamma0 != gam_old))) { // given gamma0, find the corresponding index in gam(): imin = 0; imax = ngam; igam0 = ngam / 2; while ((imax - imin) > 1) { if (gamma0 >= gam[igam0 - 1]) imin = igam0; else imax = igam0; igam0 = (imin + imax) / 2; } if (igam0 == 0) igam0 = 1; ttmp = (gamma0 - gam[igam0 - 1]) / (gam[igam0] - gam[igam0 - 1]); if (ttmp < 0.) ttmp = 0.; if (ttmp > 1.) ttmp = 1.; ttmp1 = 1. - ttmp; // given abundance, find the corresponding index in abun[]: imin = 0; imax = nabun; iabun0 = nabun / 2; while ((imax - imin) > 1) { if (abundance >= abun[iabun0 - 1]) imin = iabun0; else imax = iabun0; iabun0 = (imin + imax) / 2; } if (iabun0 == 0) iabun0 = 1; utmp = (abundance - abun[iabun0 - 1]) / (abun[iabun0] - abun[iabun0 - 1]); if (utmp < 0.) utmp = 0.; if (utmp > 1.) utmp = 1.; utmp1 = 1. - utmp; // In the following the last energy index is not used // (we have one less flux indices then in energy) for (ixi = 0; ixi < nxi; ixi++) for (ie = 0; ie < ne_loc - 1; ie++) { y1 = emission[ie + (ne_loc - 1) * ixi + (ne_loc - 1) * nxi * (igam0 - 1) + (ne_loc - 1) * nxi * ngam * (iabun0 - 1)]; y2 = emission[ie + (ne_loc - 1) * ixi + (ne_loc - 1) * nxi * igam0 + (ne_loc - 1) * nxi * ngam * (iabun0 - 1)]; y3 = emission[ie + (ne_loc - 1) * ixi + (ne_loc - 1) * nxi * igam0 + (ne_loc - 1) * nxi * ngam * iabun0]; y4 = emission[ie + (ne_loc - 1) * ixi + (ne_loc - 1) * nxi * (igam0 - 1) + (ne_loc - 1) * nxi * ngam * iabun0]; flux0[ie + ne_loc * ixi] = (utmp1 * (ttmp1 * y1 + ttmp * y2) + utmp * (ttmp * y3 + ttmp1 * y4)); } /******************************************************************************* ixi = 8; fprintf(stdout, "%d %d %d\n", iabun0, igam0, ixi); fprintf(stdout, "%f %f %f %f\n", abun[iabun0 - 1], gam[igam0 - 1], pow(10, logxi[ixi - 1]), xi[ixi - 1]); for(ie = 0; ie < ne_loc - 1; ie++) fprintf(stdout, "%d %f %f\n", ie, energy0[ie], flux0[ie + ne_loc * (ixi - 1)] / (energy0[i + 1] - energy0[i])); *******************************************************************************/ // Let's rebin the spectra to the evenly log spaced energies if ((rix != rix_old) || (strcmp(kydir, FGMSTR(pname)) != 0)) { energy2[0] = elow / (1. + pow(ehigh / elow, 1. / (nener - 1.))) * 2.; for (ie = 1; ie <= nener; ie++) { energy2[ie] = energy2[0] * pow(ehigh / elow, ie / (nener - 1.)); energy1[ie - 1] = elow * pow(ehigh / elow, (ie - 1.) / (nener - 1.)); } } for (ixi = 0; ixi < nxi; ixi++) for (ie = 0; ie < nener; ie++) flux1[ie + nener * ixi] = 0.; ie = 1; while (energy0[ie] <= energy2[0]) ie++; je = 1; quit = 0; while ((ie <= (ne_loc - 1)) && (energy0[ie - 1] < energy2[nener])) { ttmp = energy0[ie] - energy0[ie - 1]; while (energy2[je - 1] < energy0[ie]) { if (energy2[je - 1] < energy0[ie - 1]) utmp = energy0[ie - 1]; else utmp = energy2[je - 1]; if (energy2[je] < energy0[ie]) utmp1 = energy2[je]; else utmp1 = energy0[ie]; if (utmp1 > utmp) for (ixi = 0; ixi < nxi; ixi++) flux1[je - 1 + nener * ixi] += flux0[ie - 1 + ne_loc * ixi] * (utmp1 - utmp) / ttmp; if (je < nener) je++; else { quit = 1; break; } } if (quit) break; if (energy2[je - 1] > energy0[ie]) je--; ie++; } for (ixi = 0; ixi < nxi; ixi++) for (ie = 0; ie < nener; ie++) flux1[ie + nener * ixi] /= (energy2[ie + 1] - energy2[ie]); /******************************************************************************* ixi = 8; fprintf(stdout, "%d %d %d\n", iabun0, igam0, ixi); fprintf(stdout, "%f %f %f %f\n", abun[iabun0 - 1], gam[igam0 - 1], pow(10, logxi[ixi - 1]), xi[ixi - 1]); for(ie = 0; ie < nener; ie++) fprintf(stdout, "%d %f %f\n", ie, energy1[ie], flux1[ie + nener * (ixi - 1)]); *******************************************************************************/ abun_old = abundance; gam_old = gamma0; } sprintf(kydir, "%s", FGMSTR(pname)); rix_old = rix; am_old = am; thetaO_old = thetaO; h_rh_old = h_rh; /******************************************************************************/ //check if energy ranges defined are reasonable for(iband=0;iband energy1[nener-1]){ sprintf(errortxt,"kynrefrev: defined energy band reaches beyond the energy of local flux definition!"); xs_write(errortxt, 5); } } // delay_phi ide_param[19]=del_a; /******************************************************************************/ //Let's compute the Np according to its definition, i.e. transform from //intrinsic or observed in 2-10keV to total intrinsic photon flux //Lx is intrinsic photon flux in 2-10keV g_L = sqrt(1. - 2. * h / (am2 + h * h)); if( Np < 0. ){ Lx = -Np; Np *= - incgamma(2. - gamma0, E0 / EC) / ( incgamma(2. - gamma0, 2. / EC) - incgamma(2. - gamma0, 10. / EC)); }else{ Lx = Np / g_L / g_L / transf_o / ( incgamma(2. - gamma0, 2. / g_L / EC) - incgamma(2. - gamma0, 10. / g_L / EC)); Np = Lx * incgamma(2. - gamma0, E0 / EC); Lx *= ( incgamma(2. - gamma0, 2. / EC) - incgamma(2. - gamma0, 10. / EC)); } //Let's write the intrinsic photon flux in 2-10keV into the xset XSPEC variable //KYLxLamp sprintf(kyLxLamp, "%e", Lx); FPMSTR(pkyLxLamp, kyLxLamp); // Let's compute the ionisation at rin and rout for (i = 0; i < 2; i++) { if (i == 0) if (rin <= (radius[0] + r_plus)) r = radius[1] + r_plus; else r = rin; else r = rout; imin = 0; imax = nradius; ir0 = nradius / 2; while ((imax - imin) > 1) { if (r >= (radius[ir0 - 1] + r_plus)) imin = ir0; else imax = ir0; ir0 = (imin + imax) / 2; } if ((ir0 == nradius) && (r == (radius[nradius - 1] + r_plus))) ir0 = ir0 - 1; ttmp = (r - radius[ir0 - 1] - r_plus) / (radius[ir0] - radius[ir0 - 1]); ttmp1 = 1. - ttmp; // Let's interpolate gfactor between two radii gfactor0 = ttmp * gfac[ir0] + ttmp1 * gfac[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; if (lensing != 0.) { if (qn != 0.) { if (sw == 1) ionisation = logxi_norm + log10(pow(r, -qn) * lensing * gfactor0 * Np / mass2 / nH0); //sw==2 else ionisation = logxi_norm + log10(pow(r, -qn) * lensing * pow(gfactor0, gamma0 - 1.) * Np / mass2 / nH0); } else { if (sw == 1) ionisation = logxi_norm + log10(lensing * gfactor0 * Np / mass2 / nH0); //sw==2 else ionisation = logxi_norm + log10(lensing * pow(gfactor0, gamma0 - 1.) * Np / mass2 / nH0); } if (i == 0) { sprintf(kyxiin, "%e", pow(10, ionisation)); FPMSTR(pkyxiin, kyxiin); } else { sprintf(kyxiout, "%e", pow(10, ionisation)); FPMSTR(pkyxiout, kyxiout); } } } //****************************************************************************** //****************************************************************************** //****************************************************************************** // Let's integrate local emission over the accretion disc if (ide(ear, ne, nt, far, qar, uar, var, ide_param, emis_KYNrefrev, nener)) { for (ie = 0; ie < ne_xspec; ie++) photar[ie] = 0.; goto error; } //****************************************************************************** //****************************************************************************** // Let's normalize the reflected flux properly for (ie = 0; ie < ne; ie++) for (it = 0; it g_L*zzshift*E0){ far_prim[ie] = Anorm * photar1[ie]; flux_prim_tot += far_prim[ie]; } } flux_prim_tot *= (flare_duration_rg / dt); refl_ratio = flux_refl_tot / flux_prim_tot; } sprintf(kyRefl, "%e", refl_ratio); FPMSTR(pkyRefl, kyRefl); /******************************************************************************/ #ifdef OUTSIDE_XSPEC // let's write the input parameters to a file sprintf(text,"kynrefrev_%s.txt",cparam); fw = fopen(text, "w"); fprintf(fw, "a/m %12.6f\n", param[0]); fprintf(fw, "thetaO %12.6f\n", param[1]); fprintf(fw, "rin %12.6f\n", param[2]); fprintf(fw, "ms %12d\n", (int) param[3]); fprintf(fw, "rout %12.6f\n", param[4]); fprintf(fw, "phi %12.6f\n", param[5]); fprintf(fw, "dphi %12.6f\n", param[6]); fprintf(fw, "M/M8 %12.6f\n", param[7]); fprintf(fw, "height %12.6f\n", param[8]); fprintf(fw, "PhoIndex %12.6f\n", param[9]); fprintf(fw, "L/Ledd %12.6f\n", param[10]); fprintf(fw, "Np:Nr %12.6f\n", param[11]); fprintf(fw, "nH0 %12.6f\n", param[12]); fprintf(fw, "qn %12.6f\n", param[13]); fprintf(fw, "abun %12.6f\n", param[14]); fprintf(fw, "alpha %12.6f\n", ide_param[21]); fprintf(fw, "beta %12.6f\n", ide_param[22]); fprintf(fw, "rcloud %12.6f\n", ide_param[23]); fprintf(fw, "zshift %12.6f\n", param[18]); fprintf(fw, "limb %12d\n", (int) param[19]); fprintf(fw, "tab %12d\n", (int) param[20]); fprintf(fw, "sw %12d\n", (int) param[21]); fprintf(fw, "duration_rg %12.6f\n", DURATION); fprintf(fw, "deltaT %12.6f\n", dt0); fprintf(fw, "ntbin %12d\n", nt0); fprintf(fw, "ntable %12d\n", (int) param[22]); fprintf(fw, "nradius %12d\n", (int) param[23]); fprintf(fw, "division %12.6f\n", ide_param[8]); fprintf(fw, "nphi %12d\n", (int) param[25]); fprintf(fw, "nt %12d\n", nt_ratio); fprintf(fw, "t1/f1/E1 %12.4g\n", param[28]); fprintf(fw, "t2/f2/E2 %12.4g\n", param[29]); fprintf(fw, "E3 %12.4g\n", param[30]); fprintf(fw, "E4 %12.4g\n", param[31]); fprintf(fw, "Af %12.4g\n", param[32]); fprintf(fw, "qf %12.4g\n", param[33]); fprintf(fw, "photar_sw %12d\n", (int) param[34]); fprintf(fw, "smooth %12d\n", SMOOTH); fprintf(fw, "Stokes %12d\n", stokes); fprintf(fw, "polar %12d\n", polar); fprintf(fw, "r_horizon %12.6f\n", r_plus); fprintf(fw, "r_ms %12.6f\n", rms); fprintf(fw, "edivision %12d\n", (int) ide_param[14]); fprintf(fw, "nener %12ld\n", nener); fprintf(fw, "norm %12.6f\n", param[36]); fprintf(fw, "nthreads %12d\n", (int) param[35]); fprintf(fw, "Xi_in %s\n", kyxiin); fprintf(fw, "Xi_out %s\n", kyxiout); fprintf(fw, "refl_ratio %12.4g\n", refl_ratio); fclose(fw); #endif /******************************************************************************/ // let's increase the time bin due to finate time resolution of the detector for(it=0; it= ener_low[iband] && ear[ie+1] <= ener_high[iband])ttmp=1.; // first bin else if(ear[ie] < ener_low[iband] && ear[ie+1] > ener_low[iband]){ if(ear[ie+1] > ener_high[iband])ttmp=ener_high[iband]-ener_low[iband]; else ttmp=ear[ie+1]-ener_low[iband]; ttmp /= (ear[ie+1]-ear[ie]); // last bin }else if(ear[ie+1] > ener_high[iband] && ear[ie] < ener_high[iband] && ear[ie] >= ener_low[iband]) ttmp = (ener_high[iband]-ear[ie]) / (ear[ie+1]-ear[ie]); else ttmp=0; for(it=0;it= time2){ // given time1, find the corresponding index in time[]: it0 = (int) ceil( (time1 - time[0]) / deltaT + 1 ); if(it0 < 1) it0 = 1; else if(it0 > nt-1) it0 = nt-1; ttmp = (time1 - time[it0 - 1]) / (time[it0] - time[it0 - 1]); ttmp1 = 1. - ttmp; for(ie=0;ie nt-1) it0 = nt; itn = (int) ceil( (time2 - time[0]) / deltaT ); if(itn < 1) itn = 0; else if(itn > nt-1) itn = nt-1; // the first and last partial bin ttmp = (time1 - time[it0 - 1]) / (time[it0] - time[it0 - 1]); ttmp1 = 1. - ttmp; utmp = (time2 - time[itn - 1]) / (time[itn] - time[itn - 1]); utmp1 = 1. - utmp; for(ie=0;ie far_loc(:) array // local energy array --> ear_loc() // ne_loc --> number of points in local energy where the local photon flux // density in keV is defined; // disc surface in polar coords r, phi; // cosine of local emission angle --> cosmu double factor, factor1, gfactor, lensing, ionisation, fluxe[ne_loc]; double time, delay0; double ttmp, ttmp1, y1, y2; long ixi0, it; int ie, imin, imax, ir0; *ear_loc = energy1; // Normalization due to an imposed emissivity law if (limb == 0) factor = 1. / PI; if (limb == 1) factor = 1. / PI / (1. + 2.06 * 2. / 3.) * (1. + 2.06 * cosmu); if (limb == 2) factor = 1. / PI * log(1. + 1. / cosmu); factor *= nH0; // given r, find corresponding indices in radius: imin = 0; imax = nradius; ir0 = nradius / 2; while ((imax - imin) > 1) { if (r >= (radius[ir0 - 1] + r_plus)) imin = ir0; else imax = ir0; ir0 = (imin + imax) / 2; } //if (ir0 == 0) ir0 = 1; //if ((imax == nradius) && (r > (radius[nradius - 1] + r_plus)) ir0 = nradius; if ((ir0 == nradius) && (r == (radius[nradius - 1] + r_plus))) ir0 -= 1; if ((ir0 == 0) || (ir0 >= nradius)) { for(it=0;it 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 (qn != 0.) { if (sw == 1) ionisation = logxi_norm + log10(pow(r, -qn) * lensing * gfactor * Np / mass2 / nH0); // sw==2 else ionisation = logxi_norm + log10(pow(r, -qn) * lensing * pow(gfactor, gamma0 - 1.) * Np / mass2 / nH0); } else { if (sw == 1) ionisation = logxi_norm + log10(lensing * gfactor * Np / mass2 / nH0); // sw==2 else ionisation = logxi_norm + log10(lensing * pow(gfactor, gamma0 - 1.) * Np / mass2 / nH0); } // 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]); factor1 = 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]; fluxe[ie] = (ttmp1 * y1 + ttmp * y2) * factor * factor1; if (sw == 1) fluxe[ie] *= pow(gfactor, gamma0 - 2.); if (qn != 0.) fluxe[ie] *= pow(r, qn); } }else{ for(it=0;it=0. && time<=flare_duration_rg){ far_loc[ne_loc+(ne_loc+1)*it]=1.; for(ie=0;ie=1;i--) gfunc=(a-i)*i/(2.*i+1.+x-a+gfunc); gfunc = pow(x,a) * exp(-x) / (1.+x-a+gfunc); return(gfunc); } /******************************************************************************* *******************************************************************************/