Small corrections in documentation.

README.md

@@ -177,16 +177,16 @@ source code. Summary of the parameters:
   * **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~)
+    - inner edge of non-zero disc emissivity (in GM/c^2 or in r<sub>mso</sub>)
   * **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 &times; r~mso~ (the same applies for outer edge)
+         edge = par3 &times; r<sub>mso</sub> (the same applies for outer edge)
   * **par5  ... rout**
-    - outer edge of non-zero disc emissivity (in GM/c^2 or in r~mso~)
+    - outer edge of non-zero disc emissivity (in GM/c^2 or in r<sub>mso</sub>)
   * **par6  ... phi**
     - lower azimuth of non-zero disc emissivity (degrees)
   * **par7  ... dphi**
@@ -199,17 +199,17 @@ source code. Summary of the parameters:
       source is located (GM/c^(2))
   * **par10 ... PhoIndex**
     - power-law energy index of the primary flux
-  * **par11 ... L/L~Edd~**
+  * **par11 ... L/L<sub>Edd</sub>**
     - dE/dt, the intrinsic local (if negative) or the observed 
       (if positive) primary isotropic flux in the X-ray energy range 2-10keV 
-      in units of L~Edd~
+      in units of L<sub>Edd</sub>
   * **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/L~Edd~ (par11) means the luminosity towards the 
+    - if positive then L/L<sub>Edd</sub> (par11) means the luminosity towards the 
       observer
-    - if negative then L/L~Edd~ (par11) means the luminosity towards the disc
+    - if negative then L/L<sub>Edd</sub> (par11) means the luminosity towards the disc
   * **par13 ... density/ionisation**
     - density profile normalization in 10^15 cm^(-3) if positive
     - ionisation profile normalisation if it is negative
@@ -230,7 +230,7 @@ source code. Summary of the parameters:
     - abs(par16) > 1: the fraction of thermalisation is computed from 
                       difference between the incident and reflected fluxes
   * **par17 ... arate**
-    - accretion rate in units of L~Edd~ if positive or in Solar mass per  
+    - accretion rate in units of L<sub>Edd</sub> if positive or in Solar mass per  
       Julian year (365.25 days) if negative
   * **par18 ... f_col**  
     - spectral hardening factor
@@ -330,8 +330,8 @@ source code. Summary of the parameters:
       Af&times;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(E~ref~/E) * (f/1e-4)^(-1) s
-      with E~ref~ being middle of the reference energy band and E middle of 
+      100 * log10(E<sub>ref</sub>/E) * (f/1e-4)^(-1) s
+      with E<sub>ref</sub> 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 
@@ -527,16 +527,16 @@ source code. Summary of the parameters:
   * **par2  ... theta_o**
     - observer inclination in degrees (0&deg;-pole, 90&deg;-disc)
   * **par3  ... rin**
-    - inner edge of non-zero disc emissivity (in GM/c^2 or in r~mso~)
+    - inner edge of non-zero disc emissivity (in GM/c^2 or in r<sub>mso</sub>)
   * **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 &times; r~mso~ (the same applies for outer edge)
+         edge = par3 &times; r<sub>mso</sub> (the same applies for outer edge)
   * **par5  ... rout**
-    - outer edge of non-zero disc emissivity (in GM/c^2 or in r~mso~)
+    - outer edge of non-zero disc emissivity (in GM/c^2 or in r<sub>mso</sub>)
   * **par6  ... phi**
     - lower azimuth of non-zero disc emissivity (degrees)
   * **par7  ... dphi**
@@ -549,17 +549,17 @@ source code. Summary of the parameters:
       source is located (GM/c^(2))
   * **par10 ... PhoIndex**
     - power-law energy index of the primary flux
-  * **par11 ... L/L~Edd~**
+  * **par11 ... L/L<sub>Edd</sub>**
     - dE/dt, the intrinsic local (if negative) or the observed 
       (if positive) primary isotropic flux in the X-ray energy range 2-10keV 
-      in units of L~Edd~
+      in units of L<sub>Edd</sub>
   * **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/L~Edd~ (par11) means the luminosity towards the 
+    - if positive then L/L<sub>Edd</sub> (par11) means the luminosity towards the 
       observer
-    - if negative then L/L~Edd~ (par11) means the luminosity towards the disc
+    - if negative then L/L<sub>Edd</sub> (par11) means the luminosity towards the disc
   * **par13 ... density/ionisation**
     - density profile normalization in 10^15 cm^(-3) if positive, 
       i.e. n = par13 &times; r^(par14)
@@ -588,7 +588,7 @@ source code. Summary of the parameters:
     - abs(par17) > 1: the fraction of thermalisation is computed from 
                       difference between the incident and reflected fluxes
   * **par18 ... arate**
-    - accretion rate in units of L~Edd~ if positive or in Solar mass per  
+    - accretion rate in units of L<sub>Edd</sub> if positive or in Solar mass per  
       Julian year (365.25 days) if negative
   * **par19 ... f_col**
     - spectral hardening factor
@@ -695,8 +695,8 @@ source code. Summary of the parameters:
       Af&times;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(E~ref~/E) * (f/1e-4)^(-1) s
-      with E~ref~ being middle of the reference energy band and E middle of 
+      100 * log10(E<sub>ref</sub>/E) * (f/1e-4)^(-1) s
+      with E<sub>ref</sub> 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 

lmodel-kynreverb.dat

@@ -26,7 +26,7 @@ $tab     2.
 $sw      1.
 ntable   " "      80.     0.      0.      99.     99.     -1.
 nrad     " "      -4488.  -10000. -10000. 10000.  10000.  -100.
-division " "      -1.     0.      0.      1.      1.      -1.
+division " "      -1.     -1.     -1.      1.      1.      -1.
 nphi     " "      180.    1.      1.      20000.  20000.  -100.
 deltaT   GM/c^3   1.      1e-3    1e-3    10.     10.     -0.05
 nt_ratio " "      1.      1.      1.      10.     10.     -1.
@@ -67,7 +67,7 @@ limb     " "      0.      0.      0.      2.      2.      -1.
 $tab     11.
 ntable   " "      80.     0.      0.      99.     99.     -1.
 nrad     " "      -4488.  -10000. -10000. 10000.  10000.  -100.
-division " "      -1.     0.      0.      1.      1.      -1.
+division " "      -1.     -1.     -1.      1.      1.      -1.
 nphi     " "      180.    1.      1.      20000.  20000.  -100.
 deltaT   GM/c^3   1.      1e-3    1e-3    10.     10.     -0.05
 nt_ratio " "      1.      1.      1.      10.     10.     -1.

xskynrefrev.c

@@ -457,8 +457,8 @@ 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[12] = 1.;       // density/ionisation
+param[13] = 0.;       // den_prof/ion_prof
 param[14] = 1.;       // abun
 param[15] = 0.;       // thermalisation
 param[16] = 0.1;      // arate
@@ -2060,8 +2060,8 @@ for(ie=0;ie<ne;ie++){
   spectrum[ie]=0.;
   for(it=0;it<nt;it++) spectrum[ie] += far[ie+ne*it];
 // we have to divide by duration of the flare, flare_duration_rg, however,
-// we did not multiply by deltaT, thus we have to divide by
-// flare_duration_rg / (deltaT)
+// we did not multiply by deltaT (i.e. far is per second), thus we have to 
+// divide by flare_duration_rg / (deltaT)
   spectrum[ie] *= dt / flare_duration_rg;
   if(NpNr != 0.)spectrum_prim[ie] = far_prim[ie];
 }
@@ -2268,6 +2268,7 @@ if(photar_sw){
 //      photar[ie] /= flare_duration_rg;  <-- we have already divided by flare duration!
     }
   }else if(time1 >= time2){
+//  we compute flux per second at time time1
 //  given time1, find the corresponding index in time[]:
     it0 = (int) ceil( (time1 - time[0]) / deltaT + 1 );
     if(it0 < 1) it0 = 1;
@@ -2310,7 +2311,10 @@ if(photar_sw){
 //    all the whole bins
       for(it=it0+1;it<itn;it++) photar[ie] += (far[ie+ne*it]+far[ie+ne*(it-1)])/2.;
       photar[ie] *= dt / flare_duration_rg;
-      if(NpNr != 0.)photar[ie] += far_prim[ie];
+      if(NpNr != 0. && time1_rg <= flare_duration_rg ){
+        if(time2_rg < flare_duration_rg)photar[ie] += far_prim[ie] * (time2-time1);
+        else photar[ie] += far_prim[ie] * (flare_duration_sec-time1);
+      }
     }
   }
 }

xskynxilrev.c

@@ -475,8 +475,8 @@ 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[12] = 1.;       // density/ionisation
+param[13] = 0.;       // den_prof/ion_prof/density
 param[14] = 1.;       // abun
 param[15] = 300.;     // E_cut
 param[16] = 0.;       // thermalisation
@@ -3064,8 +3064,8 @@ for(ie=0;ie<ne;ie++){
   spectrum[ie]=0.;
   for(it=0;it<nt;it++) spectrum[ie] += far[ie+ne*it];
 // we have to divide by duration of the flare, flare_duration_rg, however,
-// we did not multiply by deltaT, thus we have to divide by
-// flare_duration_rg / (deltaT)
+// we did not multiply by deltaT (i.e. far is per second), thus we have to 
+// divide by flare_duration_rg / (deltaT)
   spectrum[ie] *= dt / flare_duration_rg;
   if(NpNr != 0.)spectrum_prim[ie] = far_prim[ie];
 }
@@ -3272,6 +3272,7 @@ if(photar_sw){
 //      photar[ie] /= flare_duration_rg;  <-- we have already divided by flare duration!
     }
   }else if(time1 >= time2){
+//  we compute flux per second at time time1
 //  given time1, find the corresponding index in time[]:
     it0 = (int) ceil( (time1 - time[0]) / deltaT + 1 );
     if(it0 < 1) it0 = 1;
@@ -3314,7 +3315,10 @@ if(photar_sw){
 //    all the whole bins
       for(it=it0+1;it<itn;it++) photar[ie] += (far[ie+ne*it]+far[ie+ne*(it-1)])/2.;
       photar[ie] *= dt / flare_duration_rg;
-      if(NpNr != 0.)photar[ie] += far_prim[ie];
+      if(NpNr != 0. && time1_rg <= flare_duration_rg ){
+        if(time2_rg < flare_duration_rg)photar[ie] += far_prim[ie] * (time2-time1);
+        else photar[ie] += far_prim[ie] * (flare_duration_sec-time1);
+      }
     }
   }
 }