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

Add direct definition of ionisation of the disc

* Add direct definition of ionisation of the disc, i.e. instead of indirect
  definition by specifying the density
* Add better definition of Ut(r_ms) in thermal reverberation that is well
  defined for a=1 as well
* Update definition of shifted Ecut
parent 8fb964bd
......@@ -175,7 +175,7 @@ Definition in XSPEC
* **par10 ... PhoIndex**
- power-law energy index of the primary flux
* **par11 ... L/L~Edd~**
- dE/dt;, the intrinsic local (if negative) or the observed
- dE/dt, the intrinsic local (if negative) or the observed
(if positive) primary isotropic flux in the X-ray energy range 2-10keV
in units of L~Edd~
* **par12 ... Np:Nr**
......@@ -185,10 +185,14 @@ Definition in XSPEC
- if positive then L/L~Edd~ (par11) means the luminosity towards the
observer
- if negative then L/L~Edd~ (par11) means the luminosity towards the disc
* **par13 ... density**
- density profile normalization in 10^15 cm^(-3)
* **par14 ... den_prof**
- radial power-law density profile
* **par13 ... density/ionisation**
- density profile normalization in 10^15 cm^(-3) if positive
- ionisation profile normalisation if it is negative
- this parameter cannot be zero
* **par14 ... den_prof/ion_prof**
- radial power-law density profile if par13 is positive
- radial ionisation profile if par13 is negative
- the radial profiles in both cases are given by abs(par13)*r^(par14)
* **par15 ... abun**
- Fe abundance (in solar abundance)
* **par16 ... therm**
......
......@@ -23,10 +23,10 @@ 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 2.
$sw 2.
$sw 1.
ntable " " 80. 0. 0. 99. 99. -1.
nrad " " -4488. -10000. -10000. 10000. 10000. -100.
division " " -1. -1. -1. 1. 1. -1.
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.
......@@ -35,6 +35,6 @@ 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.
Amp/qf " " 1. -1e8 -1e8 1e8 1e8 -1.
k/qf " " 1. -1e8 -1e8 1e8 1e8 -1.
$xsw 16.
nthreads " " 4. 1. 1. 100. 100. -1.
......@@ -59,12 +59,15 @@
* par12 ... Np:Nr - 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 L/Ledd (par11) means the luminosity
* towards the disc
* - if positive then L/Ledd (par11) means the luminosity
* towards the observer
* par13 ... density - density profile normalization in 10^15 cm^(-3)
* par14 ... den_prof - radial power-law density profile
* - if negative then L/Ledd (par11) means the luminosity
* towards the disc
* par13 ... density - density profile normalization in 10^15 cm^(-3)
* if positive
* - ionisation profile normalisation if it is negative
* par14 ... den_prof - radial power-law density profile if par13 is positive
* - radial ionisation profile if par13 is negative
* par15 ... abun - Fe abundance (in solar abundance)
* par16 ... therm - fraction of thermalised flux from the overal incident
* flux illuminating the disc
......@@ -385,6 +388,7 @@
#define IFL 1
#define NPARAM 40
/*
// for the energy dependence
#define NE 15
......@@ -401,9 +405,9 @@
*/
// for the frequency dependence
#define NE 80
#define E_MIN 3e-5
#define E_MAX 0.004
#define NE 100
#define E_MIN 1e-4
#define E_MAX 0.01
#define NBANDS 2
......@@ -423,9 +427,9 @@ 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] = 1.2;
ener_high[1] = 4.;
ener_high[0] = 0.8;
ener_low[1] = 1.;
ener_high[1] = 3.;
/*
ener_low[0] = 0.3;
ener_high[0] = 0.9;
......@@ -463,7 +467,7 @@ param[20] = 0.; // rcloud
param[21] = 0.; // zshift
param[22] = 0.; // limb
param[23] = 2.; // tab
param[24] = 2.; // sw
param[24] = 1.; // sw
param[25] = 80.; // ntable
param[26] = -4488.; // nrad
param[27] = -1.; // division
......@@ -703,7 +707,8 @@ 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);
Ut_rms=(rms*rms-2.*rms+am*sqrt(rms))/rms/sqrt(rms*rms-3.*rms+2.*am*sqrt(rms));
//Ut_rms=(rms*rms-2.*rms+am*sqrt(rms))/rms/sqrt(rms*rms-3.*rms+2.*am*sqrt(rms));
Ut_rms = ( 4 * ( sqrt(rms) - am ) + am ) / sqrt(3.) / rms;
// thetaO - observer inclination
ide_param[1] = param[1];
thetaO = ide_param[1];
......@@ -770,9 +775,14 @@ 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 - density/ionisation profile normalization in 10^15 cm^(-3)
nH0 = param[12];
// q_n - radial power-law density profile
if (nH0 == 0.) {
xs_write("kynrefrev: density/ionisation must be non-zero!", 5);
for (ie = 0; ie < ne; ie++) photar[ie] = 0.;
return;
}
// q_n - radial power-law density/ionisation profile
qn = param[13];
// Fe abundance (in solar abundance)
abundance = param[14];
......@@ -1866,20 +1876,21 @@ for (i = 0; i < 2; i++) {
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_NORM0 + log10(pow(r, -qn) * lensing *
gfactor0 * Np / mass / nH0);
//sw==2
else ionisation = LOGXI_NORM0 + log10(pow(r, -qn) * lensing *
pow(gfactor0, gamma0 - 1.) * Np / mass / nH0);
}
else {
if (sw == 1) ionisation = LOGXI_NORM0 +
log10(lensing * gfactor0 * Np / mass / nH0);
//sw==2
else ionisation = LOGXI_NORM0 +
log10(lensing * pow(gfactor0, gamma0 - 1.) * Np / mass / nH0);
}
if(nH0 > 0.){
if (qn != 0.) {
if (sw == 1) ionisation = LOGXI_NORM0 + log10(pow(r, -qn) * lensing *
gfactor0 * Np / mass / nH0);
//sw==2
else ionisation = LOGXI_NORM0 + log10(pow(r, -qn) * lensing *
pow(gfactor0, gamma0 - 1.) * Np / mass / nH0);
} else {
if (sw == 1) ionisation = LOGXI_NORM0 +
log10(lensing * gfactor0 * Np / mass / nH0);
//sw==2
else ionisation = LOGXI_NORM0 +
log10(lensing * pow(gfactor0, gamma0 - 1.) * Np / mass / nH0);
}
} else ionisation = log10(-nH0) + qn * log10(r);
if (i == 0) {
sprintf(kyxiin, "%e", pow(10, ionisation));
FPMSTR(pkyxiin, kyxiin);
......@@ -2350,7 +2361,7 @@ void emis_KYNrefrev(double** ear_loc, const int ne_loc, const long nt,
// disc surface in polar coords r, phi;
// cosine of local emission angle --> cosmu
double factor, factor1, gfactor, lensing, ionisation, fluxe[ne_loc];
double factor, factor1, factor2, gfactor, lensing, ionisation, fluxe[ne_loc];
double temp, x, Ccal, Lcal, Fbb, Finc, Frefl, Ftherm, temp_new, flux_thermal;
double time, delay0;
double ttmp, ttmp1, y1, y2;
......@@ -2363,7 +2374,8 @@ char errortext[255];
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 * REFLIONX_NORM;
factor *= REFLIONX_NORM;
if(nH0 > 0.) factor *= nH0;
// given r, find corresponding indices in radius:
imin = 0;
imax = nradius;
......@@ -2397,8 +2409,14 @@ if ((ir0 == 0) || (ir0 >= nradius)) {
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 (qn != 0.) ionisation = log10( Finc / ( nH0 * 1e15 * pow(r, qn) ) );
else ionisation = log10( Finc / ( nH0 * 1e15 ) );
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;
......@@ -2410,22 +2428,22 @@ if ((ir0 == 0) || (ir0 >= nradius)) {
}
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]);
factor1 *= pow(10, ionisation - logxi[0]);
}
if (ttmp > 1.) {
ttmp = 1.;
factor1 = pow(10, ionisation - logxi[nxi - 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;
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 (qn != 0.) fluxe[ie] *= pow(r, qn);
if (nH0 > 0. && qn != 0.) fluxe[ie] *= pow(r, qn);
}
// add the thermal component
if(thermalisation != 0.){
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment