xskynrefrev.c

/* 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 was developed 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 ... L/Ledd   - 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 Ledd
 * 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 positive then L/Ledd (par11) means the luminosity 
 *                    towards the observer
 *                  - 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
 *                     = 0 - only the reverberation of reflected radiation is 
 *                           computed
 *                     < 0 - only the reverberation of thermal radiation is 
 *                           computed
 *                     > 0 - both the thermal and reflection reverberation is
 *                           included
 *                     abs(par16) > 1 - the fraction of thermalisation is 
 *                                      computed from difference between the 
 *                                      incident and reflected fluxes
 * par17 ... arate  - accretion rate in units of Ledd if positive or in 
 *                    Solar mass per Julian year (365.25days) if negative
 * par18 ... f_col  - spectral hardening factor
 * par19 ... 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)
 * par20 ... 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)
 * par21 ... 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
 * par22 ... zshift - overall Doppler shift
 * par23 ... limb   - limb darkening/brightening law (emission directionality)
 *                  - if = 0 the local emisivity is not multiplied by anything
 *                  - if = 1 the local emisivity is multiplied by 1+2.06*mu
 *                    (limb darkening)
 *                  - if = 2 the local emisivity is multiplied by ln(1+1/mu)
 *                    (limb brightening)
 * par24 ... 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)
 * par25 ... 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
 * par26 ... ntable - table of relativistic transfer functions used in the model
 *                    (defines fits file with tables), 0<= ntable <= 99
 * par27 ... nrad   - number of grid points in radius
 *                  - if negative than the number of radial grid points is 
 *                    dependent on height as -nradius / height^0.66
 * par28 ... 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 
 *                            (par27) number of points is divided in the 3:2 
 *                            ratio in these regions; the value of par28 gives 
 *                            the transition radius between these regions 
 *                            (in GM/c^2)
 *                      -1 -> mixed radial grid with the transition radius at 
 *                            2*height
 * par29 ... nphi     - number of grid points in azimuth
 * par30 ... deltaT   - length of the time bin (GM/c^3)
 * par31 ... nt       - number of time subbins per one time bin
 * par32 ... 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 par33, the value gives the lower end
 *                      of the time/frequency interval of interest
 *                    - if same as par33, 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
 * par33 ... t2/f2/E2 - used only if different than par32 and if par32 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
 * par34 ... Eref1   - 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:
 *                     * for lag-energy spectra, the whole energy band is used 
 *                       as a reference band, always excluding the current 
 *                       energy bin
 *                     * for lag-frequency dependence, the energy reference band 
 *                       is abs(par34) to abs(par35) excluding overlaping part 
 *                       with energy band of interest abs(par32) to abs(par33)
 * par35 ... Eref2   - it defines the upper value of the reference energy band
 *                     for lag-energy dependence as well as in case of 
 *                     frequency dependent lags
 * par36 ... dt/Af   - lag shift for lag-energy dependence in case of 
 *                     par38=+6
 *                   - multiplicative factor in case of adding empirical hard
 *                     lags Af*f^(qf), used for par38=+16 and par38=+18; 
 *                     if par36=-1 then the following hard lags prescription
 *                     is used (see Epitropakis & Papadakis, 2017):
 *                     100 * log10(Eref/E) * (f/1e-4)^(-1) s
 *                     with Eref being middle of the reference energy band
 *                     and E middle of the energy band of interest
 * par37 ... Amp/qf  - multiplicative factor for the amplitude-energy 
 *                     dependence in case of par38=+5
 *                   - powerlaw index in case of adding empirical hard 
 *                     lags Af*f^(qf), used for par38=+16 and par38=+18
 * par38 ... xsw - function to be stored in the XSPEC photar array
 *                  0 -> spectrum at time defined by par32 and par33,
 *                 the following values correspond to energy
 *                 dependent Fourier transform at the frequency band 
 *                 defined by par32 and par33:
 *                 -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 par34 and par35 (integration in frequencies
 *                       is done in real and imaginary parts first and then 
 *                       the amplitudes are computed)
 *                 -6 -> lag for the relative reflection with respect
 *                       to reference energy band defined by par34 and 
 *                       par35 (integration in frequencies is done in real and
 *                       imaginary parts first and then the lags are computed
 *                       with frequency at half of the wrapping frequency or 
 *                       middle of the frequency band)
 *                 -7 -> amplitude  for the relative reflection
 *                       divided by amplitude in the reference energy band
 *                       defined by par34 and par35 (integration in frequencies
 *                       here is done in amplitudes directly)
 *                 -8 -> lag for the relative reflection with respect
 *                       to reference energy band defined by par34 and 
 *                       par35 (integration in frequencies here is done in 
 *                       lags directly)
 *                  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 par34 and par35 (integration in frequencies
 *                       is done in real and imaginary parts first and then 
 *                       the amplitudes are computed)
 *                  6 -> lag diluted by primary radiation with respect
 *                       to reference energy band defined by par34 and 
 *                       par35 (integration in frequencies is done in real and
 *                       imaginary parts first and then the lags are computed  
 *                       with frequency at half of the wrapping frequency or 
 *                       middle of the frequency band)
 *                  7 -> amplitude including the primary radiation
 *                       divided by amplitude in the reference energy band
 *                       defined by par34 and par35 (integration in frequencies
 *                       here is done in amplitudes directly)
 *                  8 -> lag diluted by primary radiation with respect
 *                       to reference energy band defined by par34 and 
 *                       par35 (integration in frequencies here is done in 
 *                       lags directly)
 *                 the following values correspond to frequency dependent
 *                 Fourier transform for the energy band of interest
 *                 defined by par32 and par33:
 *                 -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 par34 and par35
 *                        (rebinning here is done in real and imaginary parts 
 *                         first and then the amplitudes are computed)
 *                 -16 -> lag for the relative reflection with respect
 *                        to reference energy band defined by par34 and 
 *                        par35 (rebinning here is done in real and imaginary
 *                        parts first and then the lags are computed)
 *                 -17 -> amplitude  for the relative reflection
 *                        divided by amplitude in the reference energy
 *                        band defined by par34 and par35
 *                        (rebinning here is done in amplitudes directly)
 *                 -18 -> lag for the relative reflection with respect
 *                        to reference energy band defined by par34 and 
 *                        par35 (rebinning here is done in lags directly)
 *                  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 par34 and par35
 *                        (rebinning here is done in real and imaginary parts 
 *                         first and then the amplitudes are computed)
 *                  16 -> lag diluted by primary radiation with respect
 *                        to reference energy band defined by par34 and 
 *                        par35 (rebinning here is done in real and imaginary
 *                        parts first and then the lags are computed)
 *                  17 -> amplitude including the primary radiation
 *                        divided by amplitude in the reference energy
 *                        band defined by par34 and par35
 *                        (rebinning here is done in amplitudes directly)
 *                  18 -> lag diluted by primary radiation with respect
 *                        to reference energy band defined by par34 and 
 *                        par35 (rebinning here is done in lags directly)
 * 
 * par39 ... 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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 40

/*
// 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     100
#define E_MIN  1e-4
#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] = 0.8;
ener_low[1] = 1.;
ener_high[1] = 3.;
/*
ener_low[0] = 0.3;
ener_high[0] = 0.9;
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;    // L/Ledd
param[11] = 1.;       // Np:Nr
param[12] = 1.;       // density
param[13] = 0.;       // den_prof
param[14] = 1.;       // abun
param[15] = 0.;       // thermalisation
param[16] = 0.1;      // arate
param[17] = 2.4;      // f_col
param[18] = -6.;      // alpha
param[19] = 0.;       // beta
param[20] = 0.;       // rcloud
param[21] = 0.;       // zshift
param[22] = 0.;       // limb
param[23] = 2.;       // tab
param[24] = 1.;       // sw
param[25] = 80.;      // ntable
param[26] = -4488.;   // nrad
param[27] = -1.;      // division
param[28] = 180.;     // nphi
param[29] = 1.;       // deltaT
param[30] = 1.;       // nt
/*
param[31] = 0.;       // t1/f1/E1
param[32] = 8.3e-4;   // t2/f2/E2
param[33] = -1.;      // Eref1
param[34] = 3.;       // Eref2
*/
// 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[31] = ener_low[0];       // t1/f1/E1
param[32] = ener_high[0];      // t2/f2/E2
param[33] = ener_low[1];       // Eref1
param[34] = ener_high[1];      // Eref2

param[35] = 0.;       // dt/Af
param[36] = 1.;       // Amp/qf
param[37] = 16.;      // xsw
param[38] = 4.;       // nthreads
param[39] = 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
/*******************************************************************************
*******************************************************************************/
//Mpc in cm
#define MPC      3.0856776e+24
#define ERG      6.241509e8
// rg is gravitational unit GM/c^2 in cm^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.
//speed of light in km/s
#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(Ledd/RG2/1e15)+alog10(mass)=
//            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
//kev in SI units (not in ergs!)
#define KEV      1.6022e-16
#define H_KEVS   4.13566743e-18
#define K_KEVK   8.6174e-8
//speed of light in cm/s
#define C_CMS    2.99792458e10
#define MSOLAR   1.989e+30
#define SIGMA    5.6704e-8
#define YEAR     31557600.0
#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       exclude_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 *flx=NULL;
static double h, gamma0, nH0, qn, mass, mass2, am2, r_plus, Np, dt,
              flare_duration_rg, flare_duration_sec;
static double thermalisation, arate, x0, x1, x2, x3, Tnorm, am, f_col;
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, Ut_rms;
//double q_final, U_phi, Ur;
double thetaO, rin=0., rout, h_rh, elow, ehigh;
double Anorm, Dnorm, g_L, E0, Lx;
double flux_prim_tot, flux_refl_tot, refl_ratio, NpNr;
double zzshift;
double abundance, lensing, gfactor0, ionisation;
double arcosa3;
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);
//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];
// 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]<rms)rin=rms-r_plus;
 else rin=ide_param[2]-r_plus;
 rout=ide_param[4]-r_plus;
}else if(ms==2){
 rin=ide_param[2]*rms-r_plus;
 rout=ide_param[4]*rms-r_plus;
}else if(ms==0){
 rin=ide_param[2]-r_plus;
 rout=ide_param[4]-r_plus;
}
if(rin<0.)rin=0.;
if(rout<0.)rout=0.;
// phi  - lower azimuth of non-zero disc emissivity (deg)
ide_param[5] = param[5];
// dphi - (phi+dphi) is upper azimuth of non-zero disc emissivity (deg)
ide_param[6] = param[6];
// nphi - number of grid points in azimuth
ide_param[9] = param[28];
// smooth - whether to smooth the resulting spectrum (0-no, 1-yes)
ide_param[10] = SMOOTH;
// normal - how to normalize the final spectrum
ide_param[11] = -1.;
// ntable - table model (defines fits file with tables)
ide_param[13] = param[25];
// M/M8 - black hole mass in units of 10^8 solar masses
mass = param[7];
mass2 = mass * mass;
// height - height on the axis (measured from the center) at which the primary
//          source is located (GM/c^2)
h = param[8];
if (h >= 0.)
  if (h < r_plus + 0.1){
    h_rh = 0.1; h=r_plus+0.1;
      xs_write("kynrefrev: too low height, we set it to 0.1 above horizon", 5);
  }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];
// L/Ledd - dE/dt 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/ionisation profile normalization in 10^15 cm^(-3)
nH0 = param[12];
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];
// zshift - overall Doppler shift
if (param[21] > 0.) {
  ide_param[12] = param[21];
  Dnorm = pow( HUBBLE / MPC / CLIGHT / param[21], 2 );
}else if (param[21] < 0.) {
  ide_param[12] = 0.;
  Dnorm = pow( - HUBBLE / MPC / CLIGHT / param[21], 2 );
}else {
  ide_param[12] = 0.;
  Dnorm = 1. / ( MPC * MPC );
}
// 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[22];
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[23];
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[24];
// 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[37];
// let's set up the deltaT, ntbin and nn params according to the frequency range 
// needed
dt0 = param[29];
//let's keep the T_tail constant, i.e. nt0 * dt0 = TMAX
nt0 = (int) ceil( TMAX / dt0 );
if(param[4] != 1000.){
//here the last term is for inner echo, specially needed for low rout and low h
  nt0 = (int) ceil( 1.2 / dt0 * 
        ( rout*(1+sin(thetaO/180.*PI)) + h*cos(thetaO/180.*PI) + h*h/rout/2. +
          321./(h+1.) ) );
  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[31] >=0. ? param[31] : param[31] / ( - SEC * mass));
    freq_max = ( param[32] >=0. ? param[32] : param[32] / ( - 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[30];
// 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[31];
  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[32];
  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[31];
  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[32];
  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[31];
  en2 = param[32];
  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[33];
en4 = param[34];
if ( ( ( abs(photar_sw) == 5 || abs(photar_sw) == 6 ||
         abs(photar_sw) == 7 || abs(photar_sw) == 8 ) && 
       ( en3 > 0. && en3 >= en4 ) ) ||
     ( ( abs(photar_sw) == 15 || abs(photar_sw) == 16 || 
         abs(photar_sw) == 17 || abs(photar_sw) == 18 ) && ( fabs(en3) >= fabs(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;
}
exclude_energy=0;
#ifdef OUTSIDE_XSPEC
if( abs(photar_sw) > 10 ){
  nbands=1;
  ener_low[0]=en1;
  ener_high[0]=en2;
}else{
  nbands = NBANDS;
  if(ener_low[nbands-1] == ear_xspec[0] &&
     ener_high[nbands-1] == ear_xspec[ne_xspec])exclude_energy=1;
}
if( abs(photar_sw) == 15 || abs(photar_sw) == 16 || 
    abs(photar_sw) == 17 || abs(photar_sw) == 18 ){
  nbands++;
  ener_low[nbands-1] = fabs(en3);
  ener_high[nbands-1] = fabs(en4);
  if(en3<0.)exclude_energy=1.;
}
#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) == 7 || abs(photar_sw) == 8 ){
  if(en3 < 0){
    nbands++;
    ener_low[nbands-1] = ear_xspec[0];
    ener_high[nbands-1] = ear_xspec[ne_xspec];
    exclude_energy = 1;
  }else if( en3 > 0 && en3 < en4 ){
    nbands++;
    ener_low[nbands-1] = en3;
    ener_high[nbands-1] = en4;
    exclude_energy = 0;
  }
}
if( abs(photar_sw) == 15 || abs(photar_sw) == 16 ||
    abs(photar_sw) == 17 || abs(photar_sw) == 18 ){
  if( fabs(en3) < fabs(en4) ){
    nbands++;
    ener_low[nbands-1] = fabs(en3);
    ener_high[nbands-1] = fabs(en4);
    if( en3 < 0 )exclude_energy = 1;
    else exclude_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==fabs(en3) && en2==fabs(en4)) || abs(photar_sw) < 15){
    ne=1;
    ear[0]=en1;
    ear[1]=en2;
  }else if(en1==fabs(en3)){
    ne=2;
    ear[0]=en1;
    if(en2<fabs(en4)) {ear[1]=en2;ear[2]=fabs(en4);}
    else {ear[1]=fabs(en4);ear[2]=en2;}
  }else if(en2==fabs(en4)){
    ne=2;
    ear[2]=en2;
    if(en1<fabs(en3)) {ear[0]=en1;ear[1]=fabs(en3);}
    else {ear[0]=fabs(en3);ear[1]=en1;}
  }else if(en2==fabs(en3)){
    ne=2;
    ear[0]=en1;
    ear[1]=en2;
    ear[2]=fabs(en4);
  }else if(en1==fabs(en4)){
    ne=2;
    ear[0]=fabs(en3);
    ear[1]=en1;
    ear[2]=en2;
  }else if(en1<fabs(en3)){
    ne=3;
    ear[0]=en1;
    if(en2<fabs(en3)){ear[1]=en2;ear[2]=fabs(en3);ear[3]=fabs(en4);}
    else if(en2<fabs(en4)){ear[1]=fabs(en3);ear[2]=en2;ear[3]=fabs(en4);}
    else {ear[1]=fabs(en3);ear[2]=fabs(en4);ear[3]=en2;}
  }else if(fabs(en3) < en1){
    ne=3;
    ear[0]=fabs(en3);
    if(fabs(en4)<en1){ear[1]=fabs(en4);ear[2]=en1;ear[3]=en2;}
    else if(fabs(en4)<en2){ear[1]=en1;ear[2]=fabs(en4);ear[3]=en2;}
    else {ear[1]=en1;ear[2]=en2;ear[3]=fabs(en4);}
  }
}
//Create arrays of different dimensions according to photar_sw, i.e. different
//dimensions if ear_xspec contains energy bins or frequency bins!
ntmax=( nt > 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 || photar_sw == 7 || photar_sw == 8 ){
  lag_shift = param[35];
  ampl_ampl = param[36];
}
if(photar_sw == 16 || photar_sw == 18){
  if(param[35]==-1.){
    Af = 350.e-4*log10((fabs(en3)+fabs(en4))/(en1+en2));
    qf = 1.;
  }else{
    Af = param[35];
    qf = param[36];
  }
}
// polar - whether we need value of change in polarization angle (0-no,1-yes)
//stokes = (int) param[39];
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[38];
// thermalisation - thermalisation fraction
thermalisation = param[15];
// arate - accretion rate (note that LEDD is in erg/s not W and CLIGHT in km/s
//          not m/s and 1e-13 = erg/s/W / (km/m)^2 = 1e-7 / (1e3)^2)
if(param[16] < 0.) arate = fabs(param[16]);
else arate = param[16]*LEDD*1e-13*mass/MSOLAR*YEAR/CLIGHT/CLIGHT/(1-Ut_rms);
// f_col - colour correction factor
f_col = param[17];
if(thermalisation != 0.){
  Tnorm = f_col * K_KEVK * sqrt( C_CMS ) *
          pow( 3. * MSOLAR / (8. * PI * RG2 * YEAR * SIGMA), 0.25);
  x0 = sqrt(rms);
  arcosa3 = acos(am) / 3.;
  x1 = 2 * cos(arcosa3 - PI / 3.);
  x2 = 2 * cos(arcosa3 + PI / 3.);
  x3 = -2 * cos(arcosa3);
}
// alpha - position of the cloud in alpha impact parameter (in GM/c^2)
ide_param[21] = param[18];
// beta - position of the cloud in beta impact parameter (in GM/c^2)
ide_param[22] = param[19];
// rcloud - radius of the cloud (in GM/c^2)
ide_param[23] = param[20];
//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[21] != 0.) && (param[39] != 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<nrh; i++)fprintf(stdout,"%f\n",r_horizon[i]);
//******************************************************************************   
// Move to the extension 'height' and read its values
  ffmrhd(fptr, 1, &hdutype, &status);
  ffgnrw(fptr, &nh, &status);
//******************************************************************************
//  fprintf(stdout,"nh = %ld\n",nh);
//******************************************************************************   
// Allocate memory for height...
  if ((height = (float *) malloc(nh * 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 'height' table
  nelems = nh;
// FFGCV reads the VALUES from the first column.
  ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, height,
        &anynul, &status);
//******************************************************************************
//  for ( i=0; i<nh; i++)fprintf(stdout,"%f\n",height[i]);
//******************************************************************************   
// Move to the extension 'inclination' and read its values
  ffmrhd(fptr, 1, &hdutype, &status);
  ffgnrw(fptr, &nincl, &status);
//******************************************************************************
//  fprintf(stdout,"nincl = %ld\n",nincl);
//******************************************************************************   
// Allocate memory for inclination...
  if ((incl = (float *) malloc(nincl * 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 'inclination' table
  nelems = nincl;
// FFGCV reads the VALUES from the first column.
  ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, incl,
        &anynul, &status);
//******************************************************************************
//  for ( i=0; i<nincl; i++)fprintf(stdout,"%f\n",incl[i]);
//******************************************************************************   
// Move to the extension 'r_rh' and read its values
  ffmrhd(fptr, 1, &hdutype, &status);
  ffgnrw(fptr, &nradius, &status);
//******************************************************************************
//  fprintf(stdout,"nradius = %ld\n",nradius);
//******************************************************************************   
// Allocate memory for r_rh...
  if ((radius = (float *) malloc(nradius * 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_rh' table
  nelems = nradius;
// FFGCV reads the VALUES from the first column.
  ffgcv(fptr, TFLOAT, colnum, frow, felem, nelems, &float_nulval, radius,
        &anynul, &status);
//******************************************************************************
//  for ( i=0; i<nradius; i++)fprintf(stdout,"%f\n",radius[i]);
//******************************************************************************   
// Let's read the tables for dWadWo, q, p^r and dWadSd
// allocate memory for the arrays
  if ((dWadWo  = (float *) malloc(nincl * nh * nrh * sizeof(float))) == NULL ||
      (delay_a = (float *) malloc(nincl * nh * nrh * sizeof(float))) == NULL ||
      (q     = (float *) malloc(nradius * nh * nrh * sizeof(float))) == NULL ||
      (pr    = (float *) malloc(nradius * nh * nrh * sizeof(float))) == NULL ||
      (dWadSd= (float *) malloc(nradius * nh * nrh * sizeof(float))) == NULL ||
      (delay = (float *) malloc(nradius * nh * 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 tables
  for (ihorizon = 0; ihorizon < nrh; ihorizon++) {
    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);
    nelements1 = nrow * nincl;
    nelements2 = nrow * nradius;
    for (irow = 0; irow < nh; irow += nrow) {
//    the last block to read may be smaller:
      if ((nradius - irow) < nrow) {
        nelements1 = (nradius - irow) * nincl;
        nelements2 = (nradius - irow) * nradius;
      }
      ffgcv(fptr, TFLOAT, 1, irow + 1, 1, nelements1, &float_nulval, 
            &dWadWo[irow * nincl + nh * nincl * ihorizon],
            &anynul, &status);
      ffgcv(fptr, TFLOAT, 2, irow + 1, 1, nelements1, &float_nulval, 
            &delay_a[irow * nincl + nh * nincl * ihorizon],
            &anynul, &status);
      ffgcv(fptr, TFLOAT, 3, irow + 1, 1, nelements2, &float_nulval, 
            &q[irow * nradius + nh * nradius * ihorizon],
            &anynul, &status);
      ffgcv(fptr, TFLOAT, 4, irow + 1, 1, nelements2, &float_nulval, 
            &pr[irow * nradius + nh * nradius * ihorizon],
            &anynul, &status);
      ffgcv(fptr, TFLOAT, 5, irow + 1, 1, nelements2, &float_nulval, 
            &dWadSd[irow * nradius + nh * nradius * ihorizon],
            &anynul, &status);
      ffgcv(fptr, TFLOAT, 6, irow + 1, 1, nelements2, &float_nulval, 
            &delay[irow * nradius + nh * nradius * ihorizon],
            &anynul, &status);
    }
  }
// The FITS file must always be closed before exiting the program.
  ffclos(fptr, &status);
/*******************************************************************************
  irh=20;
  ih=99;
  for ( i=0; i<nincl; i++ ) fprintf(stdout,"%d\t%f\t%f\n",i,incl[i],
    dWadWo[i+nincl*ih+nincl*nh*irh]);
  for ( i=0; i<nradius; i++) fprintf(stdout,"%d\t%f\t%f\t%f\t%f\n",i,radius[i],
    q[i+nradius*ih+nradius*nh*irh],pr[i+nradius*ih+nradius*nh*irh],
    dWadSd[i+nradius*ih+nradius*nh*irh]);
*******************************************************************************/
// Firstly we have to free allocated memory for the arrays gfac,
// cosin, phiph, transf_d
  if((gfac    = (double *) malloc(nradius * sizeof(double))) == NULL ||
//     (cosin   = (double *) malloc(nradius * sizeof(double))) == NULL ||
//     (phiph   = (double *) malloc(nradius * sizeof(double))) == NULL ||
     (transf_d= (double *) malloc(nradius * sizeof(double))) == NULL ||
     (del     = (double *) malloc(nradius * 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_h = 0;
}
/******************************************************************************/
if (h_rh > 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[26] > 0) ide_param[7] = param[26];
else ide_param[7] = - param[26] * pow(h,-0.66);
// division - type of division in r integration (0-equidistant, 1-exponential)
if(param[27]==-1.)ide_param[8]=2.*h;
else ide_param[8]=param[27];
// 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(