Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
StrongGravity
KYNreverb
Commits
e6da4aa0
Commit
e6da4aa0
authored
Aug 04, 2016
by
Michal Dovčiak
Browse files
Update introductory comments in the source files and add draft documentation.
parent
1a680f9d
Changes
2
Hide whitespace changes
Inline
Sidebyside
README.md
View file @
e6da4aa0
**To use the KYNrefrev model in XSPEC you need:**
*
[
source files
](
https://projects.asu.cas.cz/files/note/345/KYNrefrevv1.3.4.tar.gz
)
,
*
KY tables:
[
KBHlamp_qt.fits
](
http://www.astro.cas.cz/dovciak/pub/KY/KBHlamp_qt.fits
)
and [KBHtables80.fits]
(http://www.astro.cas.cz/dovciak/pub/KY/KBHtables80.fits),
*
[
REFLION(X)
](
https://heasarc.gsfc.nasa.gov/xanadu/xspec/models/reflion.html
)
tables:

[
reflion.mod
](
https://heasarc.gsfc.nasa.gov/xanadu/xspec/models/reflion.mod.gz
)
(
old
)
,

[
reflionx.mod
](
https://heasarc.gsfc.nasa.gov/xanadu/xspec/models/reflionx.mod.gz
)
,
or in case the links are not available

[
reflion.mod
](
http://www.astro.cas.cz/dovciak/pub/KYexternal/reflion.mod
)
(
old
)
,

[
reflionx.mod
](
http://www.astro.cas.cz/dovciak/pub/KYexternal/reflionx.mod
)
.

The code is compiled inside XSPEC with the following command (assuming all the source files and FITS tables are in the directory /path/to/KYNrefrev):
*
`initpackage kynrefrev lmodel.dat /path/to/KYNrefrev`
.
To use the KYNrefrev model inside XSPEC, first the package needs to be loaded and directory with KYNrefrev set:
*
`lmod kynrefrev /path/to/KYNrefrev`
,
*
`xset KYDIR /path/to/KYNrefrev`
.
Then the model may be used:
*
`mo kynrefrev`
.
_Note_
:
In case of segmentation fault, one may need to increase the stack size, e.g. with the command
`ulimit s unlimited`
or
`ulimit s 65532`
.

**For use outside of XSPEC:**
*
One also needs the Makefile and libxspec library available
[
here
](
https://projects.asu.cas.cz/files/note/344/outside.tar.gz
)
.
*
The library to work with FITS files (libcfitsio.so) is needed, thus one needs to
define the name of the library and path to it in the provided Makefile.
*
The model parameters have to be changed inside the source file.
*
Compile with the make command:
*
`make kynrefrev`
for everything, i.e. light curves, spectra and Fourier transforms (real,
imaginary parts, amplitudes and phases).
*
Run the code:
*
`./kynrefrev`
.
*
The model creates various files described below.
_Note_
:
In case of segmentation fault, one may need to increase the stack size, e.g. with the command
`ulimit s unlimited`
or
`ulimit s 65532`
.

#### **1. Parameters of the model**
*
The meaning of the input parameters are explained at the beginning of the kynrefrev.c file.
The parameters when the code runs under XSPEC are defined in the usual way as for other XSPEC
models. The parameter definitions when run outside of XSPEC must be changed directly inside the source
code. Summary of the parameters:
* **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 nonzero disc emissivity (in GM/c^2 or in r~mso~)
* **par4 ... ms**
 switch for inner edge
 0: we integrate from inner edge = par3
 1: if the inner edge of the disc is below marginally stable orbit (MSO) 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 nonzero disc emissivity (in GM/c^2 or in r~mso~)
* **par6 ... phi**
 lower azimuth of nonzero disc emissivity (degrees)
* **par7 ... dphi**
 (phi + dphi) is upper azimuth of nonzero disc emissivity 0° ≤ dphi ≤ 360°
* **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**
 powerlaw energy index of the primary flux
* **par11 ... Np**
 dN/dt/dΩ, the intrinsic local (if negative) or the observed (if positive)
primary isotropic flux in the Xray energy range 210keV in units of L~Edd~
* **par12 ... NpNr**
 ratio of the primary to the reflected normalization
 1: selfconsistent model for isotropic primary source
 0: only reflection, primary source is hidden
 if positive then Np (par11) means the luminosity towards the observer
 if negative then Np (par11) means the luminosity towards the disc
* **par13 ... nH0**
 density profile normalization in 10^15 cm^(3)
* **par14 ... q_n**
 radial powerlaw 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 φ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 θ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**
 0: for isotropic emission (flux ~ 1)
 1: for Laor's limb darkening (flux ~ 1+2.06μ)
 2: for Haardt's limb brightening (flux ~ ln (1+1/μ))
* **par21 ... tab**
 which reflion table to use
 1: reflion (the old one, lower cutoff energy at 1eV, not good for PhoIndex > 2)
 2: reflionx (the newer one, lower cutoff energy at 100eV)
* **par22 ... sw**
 switch for the way how to compute the refl. spectra
 1: use the computed ionisation parameter, ξ, for the interpolation in reflion, i.e. use proper
total incident intensity with the shifted cutoffs
 2: use the ionisation parameter, ξ, correspondent to the computed normalization of the incident flux, i.e. do not shift
the cutoffs when computing the total incident intensity
* **par23 ... ntable**
 defines fits file with tables (0 ≤ ntable ≤ 99), currently 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 lagenergy spectra, always excluding
the current energy bin; it must be nonnegative in case of lagfrequency dependence
* **par32 ... E4**
 it defines the upper value of the reference energy band for lagenergy dependence as well as in case
of frequency dependent lags
* **par33 ... tshift/Af**
 lag shift for lagenergy 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 amplitudeenergy 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**
 defines output in the XSPEC (photar array)
 0: spectrum for time interval 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
* **par37 ... norm**
 **has to be set to unity!**
*
Parameters that need to be defined inside the
**kynrefrev.c**
code when run outside of XSPEC:
 _energy_ in the following lines:
#define NE 30
#define E_MIN 0.3
#define E_MAX 80.
 choose the _energy bands of interest_ in the following lines:
#define NBANDS 5
.
.
.
//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] = 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;
 _all basic parameters_ of the model (physical ones as well as those defining resolution grid for
computations) are defined in the following lines:
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] = 2.e4; // time/frequency/energylower
param[29] = 8.e4; // time/frequency/energyupper
param[30] = 1.; // reference energy bandlower
param[31] = 3.; // reference energy bandupper
param[32] = 0.; // lag shift or multiplicative factor for hard lags
param[33] = 1.; // amplitude multiplicative factor or powerlaw index for hard lags
param[34] = 6.; // photar_sw
param[35] = 4.; // nthreads
param[36] = 1.; // norm
 some parameters are later changed in the loops for convenience (to create files for grid of
parameters), see lines as:
for (ia=0;ia<=1;ia++){
param[0] = (double) ia;
for (iinc=20;iinc<=80;iinc+=20){
param[1] = (double) iinc;
for (ih=1;ih<=20;ih++){
param[8] = 1.5 * (100./1.5)**((ih1.)/19.);
#### **2. Output files created**
*
The output files are created only when the code is run outside XSPEC.
*
The following naming scheme is used for the output files:

**AAA**
is 100
×
the horizon value (thus 100 means a=1 and 200 means a=0),

**BB**
is the inclination in degrees,

**CCCC**
is 10
×
the height (e.g. 0030 means h=3),

**u1**
is used for phase unwrapped in frequency dependence,

**u2**
is used for phase unwrapped in energy dependence.
*
All spectra and light curves are always computed in photon numbers and per keV and per second, time is in
seconds and frequency in Hz.
*
The output files created by
**kynrefrev.c**
code:

below by relative reflection it is meant the disc response divided by the total primary flux in the
flash.

**kynrefrev_photar.dat**
→
the values as would be given inside XSPEC,

**kynrefrev_AAA_BB_CCCC.txt**
→
the summary of parameter values,

**kynrefrev_AAA_BB_CCCC_far.dat**
→
the time evolving observed reflection spectrum where the
energy changes with rows and the time changes with columns,

**kynrefrev_AAA_BB_CCCC_flux_prim.dat**
→
the total observed primary flux (i.e. integrated in
energy) per second, it is constant during the duration of the flare,

**kynrefrev_AAA_BB_CCCC_lc.dat**
→
the light curve of the observed reflection (2^nd column)
where we integrated over the whole energy range, the 1^st column contains the time,

**kynrefrev_AAA_BB_CCCC_spectrum.dat**
→
the time integrated spectrum of the observed
reflection (2^nd column) and the observed primary (3^rd column), both are divided by the flare
duration, the 1^st column contains the central value of the energy bins in keV.

**kynrefrev_AAA_BB_CCCC_bands_lc.dat**
→
the light curves for the observed reflection (2^nd
and higher columns) where in each column the light curve is integrated over different energy band
(defined in the code), the 1^st column contains the time,

**kynrefrev_AAA_BB_CCCC_bands_prim.dat**
→
the observed primary flux per second, it is
constant during the duration of the flare where in each column the flux is
integrated over different energy band (defined in the code).

**kynrefrev_AAA_BB_CCCC_real.dat**
→
the real part of the FFT of the relative reflection with frequency
changing with rows and energy with columns,

**kynrefrev_AAA_BB_CCCC_imag.dat**
→
the imaginary part of the FFT of the relative reflection with frequency
changing with rows and energy with columns,

**kynrefrev_AAA_BB_CCCC_ampl.dat**
→
the amplitude of the FFT of the relative reflection with frequency
changing with rows and energy with columns,

**kynrefrev_AAA_BB_CCCC_phase.dat**
,
**kynrefrev_AAA_BB_CCCC_phase_u1.dat**
,
**kynrefrev_AAA_BB_CCCC_phase_u2.dat**
→
the phase of the FFT of the relative reflection with
frequency changing with rows and energy with columns,

**kynrefrev_AAA_BB_CCCC_real_tot.dat**
→
the real part of the FFT of the total signal
(reflection response plus primary flash) with frequency changing with rows and energy with columns,

**kynrefrev_AAA_BB_CCCC_imag_tot.dat**
→
the imaginary part of the FFT of the total signal
(reflection response plus primary flash) with frequency changing with rows and energy with columns,

**kynrefrev_AAA_BB_CCCC_ampl_tot.dat**
→
the amplitude of the FFT of the total signal
(reflection response plus primary flash) with frequency changing with rows and energy with columns,

**kynrefrev_AAA_BB_CCCC_phase_tot.dat**
,
**kynrefrev_AAA_BB_CCCC_phase_tot_u1.dat**
,
**kynrefrev_AAA_BB_CCCC_phase_tot_u2.dat**
→
the phase of the FFT of the total signal (reflection
response plus primary flash) with frequency changing with rows and energy with columns,

**kynrefrev_AAA_BB_CCCC_bands_real.dat**
→
the real part, as a function of frequency, of the
FFT of the relative reflection integrated in different energy bands as defined in the code (2^nd and
higher columns), the 1^st column contains the frequency,

**kynrefrev_AAA_BB_CCCC_bands_imag.dat**
→
the imaginary part, as a function of frequency, of
the FFT of the relative reflection integrated in different energy bands as defined in the code (2^nd
and higher columns), the 1^st column contains the frequency,

**kynrefrev_AAA_BB_CCCC_bands_ampl.dat**
→
the amplitude, as a function of frequency, of the
FFT of the relative reflection integrated in different energy bands as defined in the code (2^nd and
higher columns), the 1^st column contains the frequency,

**kynrefrev_AAA_BB_CCCC_bands_phase.dat**
,
**kynrefrev_AAA_BB_CCCC_bands_phase_u1.dat**
→
the phase,
as a function of frequency, of the FFT
of the relative reflection integrated in different energy bands as defined in the code (2^nd and higher
columns), the 1^st column contains the frequency,

**kynrefrev_AAA_BB_CCCC_bands_real_tot.dat**
→
the real part, as a function of frequency,
of the FFT of the total signal (reflection response plus primary flash) integrated in different energy
bands as defined in the code (2^nd and higher columns), the 1^st column contains the frequency,

**kynrefrev_AAA_BB_CCCC_bands_imag_tot.dat**
→
the imaginary part, as a function of frequency,
of the FFT of the total signal (reflection response plus primary flash) integrated in different energy
bands as defined in the code (2^nd and higher columns), the 1^st column contains the frequency,

**kynrefrev_AAA_BB_CCCC_bands_ampl_tot.dat**
→
the amplitude, as a function of frequency,
of the FFT of the total signal (reflection response plus primary flash) integrated in different energy
bands as defined in the code (2^nd and higher columns), the 1^st column contains the frequency,

**kynrefrev_AAA_BB_CCCC_bands_phase_tot.dat**
,
**kynrefrev_AAA_BB_CCCC_bands_phase_tot_u1.dat**
→
the phase, as a function of frequency, of the
FFT of the total signal (reflection response plus primary flash) integrated in different energy bands as
defined in the code (2^nd and higher columns), the 1^st column contains the frequency,

**kynrefrev_AAA_BB_CCCC_freq_wrap.dat**
→
the first wrapping frequency for the phase computed
for the relative reflection (the lowest one from all energy bins)

**kynrefrev_AAA_BB_CCCC_fft_tot_int.dat**
→
energy dependent average values of the FFT of the total
signal (real part, imaginary part, amplitude, phase, unwrapped phase, delay, ratio of the amplitudes for
the energy band of interest and reference energy band, delay between the energy band of interest and
reference energy band computed from wrapped and unwrapped phases and directly averaged delay between the
two energy bands as well as the ratio of the amplitudes and delay difference between the energy band of
interest and reference energy band where reference energy band always excludes the current energy bin),
the average is computed in the range of 0 to the first wrapping frequency, the 1^st column contains
the central value of energy bins in keV; the FFT here is averaged over frequencies for real and imaginary
parts first and then, from the result, all the rest quantities are computed except for the
directly averaged delay, where the delay is computed first from real and imaginary parts of the FFT for
each frequency and only then it is averaged (just for comparison).

**kynrefrev_AAA_BB_CCCC_fft_tot_fband.dat**
→
energy dependent average values of the FFT of the total
signal (real part, imaginary part, amplitude, phase, unwrapped phase, delay, ratio of the amplitudes for
the energy band of interest and reference energy band, delay between the energy band of interest and
reference energy band computed from wrapped and unwrapped phases and directly averaged delay between the
two energy bands as well as the ratio of the amplitudes and delay difference between the energy band of
interest and reference energy band where reference energy band always excludes the current energy bin),
the average is computed in the frequency range given by param[28] and param[29], the 1^st column
contains the central value of energy bins in keV; the FFT here is averaged over frequencies for real and
imaginary parts first and then, from the result, all the rest quantities are computed except for the
directly averaged delay, where the delay is computed first from real and imaginary parts of the FFT for
each frequency and only then it is averaged (just for comparison).
xside.c
View file @
e6da4aa0
...
...
@@ 35,9 +35,10 @@
* f) azimuthal emission angle (only needed when emission depends on this angle)
*
* These functions differ for different spin, a/M, and inclination, theta_o.
* They are read and interpolated for particular a/M and theta_o from the fits
* file called 'KBHtables80.fits'. The format of this fits file is described in
* detail below.
* Some are read and interpolated for particular a/M and theta_o from the fits
* file called 'KBHtables80.fits', others are computed from other functions
* stored in this FITS file. The format of this fits file is described in detail
* below.
*
* By 'local' it is meant 'with respect to the local inertial frame
* connected with the fluid in the accretion disc' everywhere in this code
...
...
@@ 49,14 +50,14 @@
* ear  array of energy bins (same as 'ear' for local XSPEC models)
* ne  number of energy bins (same as 'ne' for local XSPEC models)
* nt  number of grid points in time (nt=1 means stationary model)
* far
(
ne,nt
)
 array of photon number density flux per bin
* far
[
ne,nt
]
 array of photon number density flux per bin
* (same as 'photar' for local XSPEC models but with time
* resolution)
* qar
(
ne,nt
)
 array of Stokes parameter Q divided by energy
* qar
[
ne,nt
]
 array of Stokes parameter Q divided by energy
* ("photon number density flux per bin for Stokes parameter Q")
*  it characterizes a linear polarization perpendicular or
* parallel to the disc
* uar
(
ne,nt
)
 array of Stokes parameter U divided by energy
* uar
[
ne,nt
]
 array of Stokes parameter U divided by energy
* ("photon number density flux per bin for Stokes parameter U")
*  it characterizes a linear polarization in the direction +45
* (if positive) or 45 degrees (if negative) with the direction
...
...
@@ 64,7 +65,7 @@
* looking towards the coming emitted light beam, the direction of
* polarization is 45 degrees counterclockwise from the "up"
* direction...)
* var
(
ne,nt
)
 array of Stokes parameter V divided by energy
* var
[
ne,nt
]
 array of Stokes parameter V divided by energy
* ("photon number density flux per bin for Stokes parameter V")
*  it characterizes circular polarization  counterclockwise
* (righthanded) if positive and clockwise (lefthanded) if
...
...
@@ 114,250 +115,248 @@
* (int) ide_param[9]: nphi  number of grid points in azimuth
* (int) ide_param[10]: smooth  whether to smooth the resulting spectrum
* (0no, 1yes)
ide_param[11]: normal  how to normalize the final spectrum
= 0.  normalization to unity (total flux=1.)
(e.g. for line)
> 0.  normalization to the flux at 'normal' keV
(usually used for continuum)
= 1.  final spectrum is not normalized in ide
= 2.  final spectrum is normalized to have maximum
photon number density flux equal unity
ide_param[12]: zshift  overall Doppler shift
(int) ide_param[13]: ntable  table of transfer functions used in the model
(defines fits file with tables), 0<= ntable <= 99
(int) ide_param[14]: edivision  type of division in local energies
(0equidistant, 1exponential = more points
with lower energy)
(int) ide_param[15]: variability_type  variable emission is for
reverberation (1)
 need not to be set if nt=1
ide_param[16]: dt  time step (in GM/c^2)
 need not to be set if nt=1
(int) ide_param[17]: polar  whether polarization from fits tables is needed
1  yes, 0  no
(we usually do not need polarization in XSPEC and
so we need not to read this information from FITS
file  the tables will occupy less space in
memory)
ide_param[18]: delay_r  (rr_plus)  delay at this coord. will be
subtracted in nonstationary calculations
ide_param[19]: delay_phi  phi  delay in degrees at this BL. coord. will be
subtracted in nonstationary calculations
> the last two parameters are used if we want to set the beginning of
time axis, e.g. for the moment when photons emitted from the spot
that was at the closest approach to the observer come to the observer
ide_param[20]: nthreads  how many threads should be used for computations
ide_param[21]: alpha  position of the cloud in alpha impact parameter (in GM/c^2)
ide_param[22]: beta  position of the cloud in beta impact parameter (in GM/c^2)
ide_param[23]: rcloud  radius of the cloud (in GM/c^2)
(int) ide_param[24]: observed_flux  whether the flux defined in emissivity
subroutine is local one (0) or the observed
one (1)
NOTE: accuracy vs. speed trade off depends mainly on: nrad, nphi

emissivity subroutine:
emissivity(ear_loc,ne_loc,nt,far_loc,qar_loc,uar_loc,var_loc,
$ rnow,pnow,cosine,phiphoton,alpha_o,beta_o,delay,g)
External emissivity subroutine has 8 parameters:
ear_loc(0:ne_loc)  array of local energies where local photon number density
flux far_loc is defined
 if ear_loc(0)>0 then the local emissivity consists of two
energy regions where flux is nonzero, local flux between
these regions is zero, this applies only for stationary
models, i.e. nt=1
 ear_loc(0) defines number of points in local energy where
the local flux is zero (in the middle region), this
applies only for stationary models, i.e. nt=1
ne_loc  number of points where local photon number density flux is defined
(in energies)
nt  number of grid points in time (nt=1 means stationary model)
far_loc(0:ne_loc,nt)  array of local photon number density flux (per keV)
for 'each time'
 if the local emissivity consists of two separate
nonzero regions (i.e. ear_loc > 0.) then far_loc(0,1)
is the index of the last point of the first nonzero
local energy region (used only in stationary case with
nt=1)
 if the local emissivity is zero for a certain time
(jt) for all energies (je) we can put far_loc(ne_loc,jt)=0
and then we exclude adding zeros to the
total flux (speeding up the computations),
if local flux for jt is NOT zero for all energies
je=1..ne_loc then far_loc(ne_loc,jt) MUST BE NONZERO!!!

NOTE: in the following dener means:
dener  width of the equidistant interval between
neighbouring points in local energy (for edivision=0) or in
logarithm of local energy (for edivision=1)
Example: ear_loc(0)=5
far_loc(0,1)=12 and nt=1
=> this means that:
 ear_loc(1)..ear_loc(12) is the first local energy region with nonzero
local flux far_loc(1,1)..far_loc(12,1)
 then there should be 5 points in local energy where local flux is zero
 ear_loc(13)..ear_loc(ne_loc) is the second local energy region with
nonzero local flux far_loc(13,1)..far_loc(ne_loc,1)
 while ear_loc(i)ear_loc(i1) = dener if edivision = 0 or
log10(ear_loc(i))log10(ear_loc(i1) = dener if edivision = 1
for i!=13
for i=13 it is:
ear_loc(13)ear_loc(12) = (ear_loc(0)+1)*dener if edivision = 0 or
log10(ear_loc(13))log10(ear_loc(12)) = (ear_loc(0)+1)*dener
(i.e. we are skipPIng the zerolocalemissivity region)

qar_loc(ne_loc,nt)  array of local Stokes parameter Q divided by local
energy for 'each time' ("local photon number density
flux per keV for Stokes parameter Q")
 it characterizes a linear polarization perpendicular
or parallel to the disc in local frame
uar_loc(ne_loc,nt)  array of local Stokes parameter U divided by local
energy for 'each time' ("local photon number density
flux per keV for Stokes parameter U")
 it characterizes a linear polarization in the direction
+45 (if positive) or 45 (if negative) degrees with the
direction perpendicular to the disc in local frame (the
angle +45 means that when we are looking toward coming
emitted light beam, the direction of polarization is 45
degrees counterclockwise from the "up" direction...)
var_loc(ne_loc,nt)  array of local Stokes parameter V divided by local
energy for 'each time' ("local photon number density
flux per keV for Stokes parameter V")
 it characterizes circular polarization 
counterclockwise (righthanded) if positive and
clockwise (lefthanded) if negative in local frame
rnow  radius where the local photon flux far_loc (or Stokes parameters
qar_loc, uar_loc, var_loc) at local energy ear_loc is wanted
pnow  azimuth where the local photon flux far_loc (or Stokes parameters
qar_loc, uar_loc, var_loc)at local energy ear_loc is wanted
cosine  cosine of the local angle between emitted ray and local disc normal
phiphoton  angle between emitted ray projected onto the plane of the disc
(in the local frame of the moving disc) and radial component of
the local tetrade (in rad)
IMPORTANT NOTE: while integrated arrays far, qar and var are evaluated per
energy bin (i.e. these quantities are integrated over energy
bin), local arrays far_loc, qar_loc and var_loc are evaluated
per keV (i.e. they are NOT integrated over energy bin)!!!

Transfer functions in the fits file KBHtablesNN.fits
Transfer functions are defined here as tables for different values of horizon
(r_horizon, not a/M!) and observer inclination angle theta_o. Each table
consists of values of particular transfer function for different r (rows) and
phi (columns) coordinates. Particular r_horizon, theta_o, r and phi where the
functions are given are defined at the beginning of the fits file as vectors.
Definition of KBHtablesNN.fits:
0. All of the extensions defined below are binary.
1. The first extension contains 6 integers defining which of the transfer
functions are present in the tables. The integers correspond to delay,
gfactor, cosine, lensing, polarization angle and azimuthal emission angle,
respectively. Value 0 means that the function is not present in the tables,
value 1 means it is.
2. The second extension contains vector of horizon values.
(1.00 <= r_horizon <= 2.00)
3. The third extension contains vector of values of the observer's
inclination angle in degrees. (0 <= theta_o <= 90, 0pole, 90disc)
4. The fourth extension contains vector of (rr_horizon) values.
Note it is not rcoordinate itself but 'coordinate distance' from the
horizon!!! (0 <= r < infinity, 0  horizon)
5. The fifth extension contains vector of phi values. Subroutine ide expects
phi to be Kerr ingoing coordinate, not BoyerLindquist one!
(0 <= phi <= 2*PI)
6. All previous vectors have to have values sorted in increasing order!
7. In following extensions the transfer functions are defined, each extension
is for particular value of r_horizon and theta_o. The values of r_horizon
and theta_o are changing with each extension in the following order:
r_horizon[1] x theta_o[1],
r_horizon[1] x theta_o[2],
r_horizon[1] x theta_o[3],
...
...
r_horizon[2] x theta_o[1],
r_horizon[2] x theta_o[2],
r_horizon[2] x theta_o[3],
...
...
8. Each of these extensions has the same number of columns (up to 6). In each
column a particular transfer function is defined  delay, gfactor, cosine,
lensing, polarization angle and azimuthal emission angle, repsectively.
The order of the functions is important but some of the functions may be
missing as defined in the first extension (see 1. above).
delay  BoyerLindquist time that elapses between emission of a photon from
the disc and absorption of the photon by the observers eye at
infinity minus a particular constant (so that delay = zero close to
the xaxis in the equatorial plane, observer is at infinity in the
direction of yaxis with x=0, the black hole is at x=y=z=0)
gfactor  ratio of the energy of a photon received by the observer at
infinity to the local energy of the same photon when emitted
from an accretion disc
cosine  cosine of the local angle between the emitted light ray and local
disc normal
lensing  ratio of the area at infinity perpendicular to the light rays
through which photons come to the area on the disc from which
these photons are emitted
polarization angle
 if the light emitted from the disc is linearly polarized then the
direction of polarization will be changed by this angle in infinity 
counterclockwise if positive, clockwise if negative (we are looking
toward the coming emitted beam), on the disc we measure the angle of
polarization with respect to "up" direction perpendicular to the disc
with respect to the local rest frame, in infinity we also measure the
angle of polarization with respect to "up" direction perpendicular to
the disc  "polarization angle" is the difference between these two
angles
phi_photon  angle between emitted ray projected onto the plane of the disc
(in the local frame of the moving disc) and radial component of
the local tetrade (in rad)
9. Each row corresponds to a particular value of r (see 4. above).