bjet_core
bj_core02 namespace
-
namespace bj_core02
Functions
-
int ifexist(char *fname)
Checks if a file exists.
This function checks if the specified file exists on the system.
- Parameters:
fname – A pointer to the name of the file to check.
- Returns:
int Returns 1 if the file exists, 0 otherwise.
-
double N_e_BknPowLaw(double gamma)
Calculates the electron spectrum using a broken power-law model.
- Parameters:
gamma – The gamma value for which the electron spectrum is calculated
- Returns:
The electron spectrum value for the given gamma
-
double N_e_SmoothBknPowLaw(double gamma)
Calculates the electron energy distribution function using a smooth broken power law with an exponential cutoff.
The function is based on the equation described in the paper by Tavecchio et al. (2001), with an additional exponential cutoff at high energies.
- Parameters:
gamma – The electron energy
- Returns:
The value of the electron energy distribution function at the given energy
-
double N_e_PileUp(double gamma)
Calculates the Pile Up factor (N_e) based on the given gamma value.
The Pile Up factor is computed using the following formula: N_e = K1 * gamma^2 * exp(-2.0 * gamma / GAMMA_BRK)
- Parameters:
gamma – The value of gamma.
- Returns:
The calculated Pile Up factor.
-
double N_e(double gamma)
Calculates the number of electrons in the given electron spectrum shape.
This function calculates the number of electrons in a given electron spectrum shape. It supports three different shapes: PileUp, SmoothBknPowLaw, and BknPowLaw.
See also
N_e_PileUp
See also
N_e_SmoothBknPowLaw
See also
N_e_BknPowLaw
- Parameters:
gamma – The gamma value representing the energy of the electron.
- Returns:
The number of electrons in the given electron spectrum shape.
-
double N_e_PowLaw(double gamma)
Calculate the number density of electrons using a power law function.
This function calculates the number density of electrons based on the provided gamma value using a power law function. The power law function is defined as N(e) = N_VAL[SL_CUR] * gamma^(-n_n).
- Parameters:
gamma – The gamma value for which to calculate the number density of electrons.
- Returns:
The number density of electrons calculated using the power law function.
-
double N_e_PowLawExpcut(double gamma)
Calculates the value of a power law distribution with an exponential cutoff.
- Parameters:
gamma – The input parameter for the power law distribution.
- Returns:
The calculated value of the power law distribution with exponential cutoff.
-
double N_e_Jet(double gamma)
Switches between different electron spectrum shapes for the jet.
- Parameters:
gamma – The value of gamma.
- Returns:
The calculated value of foo.
-
double ftr(double F)
Calculates the value of F multiplied by the power of DOP_JET raised to the third power, and then multiplied by (1.0 + z).
Note
It is important to ensure that DOP_JET and z are properly defined before calling this function.
Note
The result of the calculation may vary depending on the values of DOP_JET and z.
- Parameters:
F – The input value to be multiplied.
- Returns:
The result of the calculations described above.
-
double int_spec_j(double N_0, int a)
Calculate the integral spectrum using a power law (jet).
- Parameters:
N_0 – The normalization factor
a – The power law parameter (use 1 for particle density, 2 for energy density)
- Returns:
The calculated integral spectrum
-
double int_spec_b(int a)
Calculates the integral of a power-law function for a given input parameter for the blob.
- Parameters:
a – The input parameter for the power-law function. Use 1 for particle density, 2 for energy density.
- Returns:
The calculated integral of the power-law function.
-
void description()
This method provides a description of the code and the processes it simulates.
This code computes the radiative output from a homogeneous blob and a stratified jet given their particle energy distribution. The processes accounted for are:
synchrotron radiation (blob and jet)
SSC radiation (blob and jet)
2nd order SSC radiation (blob)
external inverse Compton & absoption on a disk, hot corona & BLR radiation field
external inverse Compton on the jet synchrotoron radiation field
See also
-
int load_params(char *name)
Load parameters from a parameter file.
This function loads parameters from a parameter file specified by the
nameparameter. The parameter file should be formatted in a specific way, with each parameter followed by its value. If any error occurs while reading the parameter file or if any of the parameter values are out of range, the function will return 0, indicating failure.- Parameters:
name – The name of the parameter file.
- Returns:
Returns 1 if the parameters were successfully loaded, 0 otherwise.
-
int load_params_from_list(char **list, int model_type, int starting_index, int list_length)
Loads parameters from a string list.
The function loads various parameters from a string list and assigns them to corresponding variables. It first reads transformation parameters such as redshift, Hubble constant, and angle. Then it reads blob parameters such as Doppler factor, particle density, slopes, gamma min and max values, gamma break, magnetic field strength, radius, emitting length, absorption level, and blob distance. Finally, if the model type is 1, it reads nucleus parameters such as disc black body temperature, tore black body temperature, nucleus luminosity, fraction of luminosity scattered/reprocessed isotropically, and tore luminosity.
The function validates the loaded parameters and returns 0 if there is an error in any of the parameters.
- Parameters:
list – The list of parameters as strings.
model_type – The type of the model.
starting_index – The starting index in the list.
list_length – The length of the list.
- Returns:
0 if there is an error in the parameters, otherwise 1.
-
int run_models()
Executes the run_models function.
The run_models function is responsible for running the models in the software. This function should be called to execute the models and obtain the desired output.
- Returns:
void
Variables
-
int INPUT_MODE = 0
Sarah Youngquist added; 0 = normal 1 = no prev files
-
char DATA_DIRECTORY[254] = "output"
-
int CASE_JET
-
int CASE_X
-
int CASE_EIC
-
const int EBLFLAG = 3
0: Kneiske et al. (2002,2004), 1: Franceschini (2008), 2: Finke (2010), 3: Franceschini (2017)
-
time_t T_START
-
time_t T_END
-
char PARA_FN[256]
-
int IIR_level = 0
-
int NU_DIM = 50
CURRENT NUMBER OF SPECTRAL POINTS.
-
double NU_STR = 0.0
START FREQUENCY Obs frame.
-
double NU_END = 0.0
END FREQUENCY Obs frame.
-
double NU_STR_B = 0.0
START FREQUENCY Blob frame.
-
double NU_END_B = 0.0
END FREQUENCY Blob frame.
-
double Utot_e = 0.0
-
double Utot_B = 0.0
-
double Utot_p = 0.0
-
double I_BASE = 0.0
-
double I_JET = 0.0
-
double z
-
double H_0
-
double THETA
-
double D_L
-
double DOP_B
-
double LOR_B
-
double V_B
-
double V_B_APP
DOP_BB,.
-
double R_src
-
double R_blr
-
double L_src
-
double L_nuc
-
double B
-
double K1
-
double K2
-
double N1
-
double N2
-
double GAMMA_MIN
-
double GAMMA_BRK
-
double GAMMA_MAX
-
double T_BB
-
double tau
fraction of L_nuc scattered/rerocessed isotropically (EIC)
-
double tau_tor
fraction of L_tor scattered/rerocessed isotropically (EIC)
-
double B_0
-
double N_0
-
double n_n
-
double n_N
-
double n_G
-
double D_b
-
double D_b_src_j
-
double D_BJ
-
double U_B
-
double U_e
-
double jj1
-
double kk1
-
double L_jet_eic
-
double L_tor
-
double T_BB_tor
-
double Tcool_synch
-
double Tcool_vhe
-
double Tcool_radio
-
double Tcool_synch_gbreak
-
double Tcool_SSC_gbreak
-
double Tcool_EIC_gbreak
-
double DOP_JET
-
double LOR_JET
-
double V_JET
-
double V_JET_APP
-
double DOP_B_J
-
double LOR_B_J
-
double V_B_J
-
int SL_DIM
-
int jet_slice_blob
-
double PHI
-
double PHI_src
-
double THETA_src_j
-
double J_LEN
-
double J_LEN_src
-
double A
-
double X_MIN
-
double X_MAX
-
double X_OUT
-
double Y_MIN
-
double Y_MAX
-
double R_MOY
-
double Xcut
-
double antisym
-
double Sprev
-
double Sum_S
-
int SL_CUR
-
double GAMMA_MIN1
-
double GAMMA_MAX_0
-
int null0
-
double V_exp
-
double DEL_Tph
-
double DEL_Tj
-
char prefix[256]
-
double NU[NU_DIM_MAX + 1]
-
double I_rad[NU_DIM_MAX + 1]
-
double I_rad1st[NU_DIM_MAX + 1]
for 2nd order SSC
-
double I_CMB[NU_DIM_MAX + 1]
-
double X_VAL[SL_DIM_MAX + 2]
-
double Y_VAL[SL_DIM_MAX + 2]
-
double DEL_X[SL_DIM_MAX + 2]
-
double Sum_DEL[SL_DIM_MAX + 2]
-
double Stot[SL_DIM_MAX]
-
double Sr_TOT[SL_DIM_MAX]
-
double B_VAL[SL_DIM_MAX + 1]
-
double N_VAL[SL_DIM_MAX + 1]
-
double G_VAL[SL_DIM_MAX + 1]
-
double UJ_SLICE[SL_DIM_MAX + 1]
-
double PJ_SLICE[SL_DIM_MAX + 1]
-
double I_SYN_EXT[NU_DIM_MAX + 1]
-
double L_BB_nuc[NU_DIM_MAX + 1]
-
double I_eic_jet[NU_DIM_MAX + 1]
-
double I_rad_syn[NU_DIM_MAX + 1]
-
double I_rad_syn2[NU_DIM_MAX + 1]
-
double I_rad_com[NU_DIM_MAX + 1]
-
double I_rad2nd[NU_DIM_MAX + 1]
-
double I_rad_ext[NU_DIM_MAX + 1]
-
double I_rad_ext_Int[NU_DIM_MAX + 1]
-
double I_rad_ext1[NU_DIM_MAX + 1]
-
double I_com_ext[NU_DIM_MAX + 1]
-
double I_com_ext1[NU_DIM_MAX + 1]
-
double I_com_disc[NU_DIM_MAX + 1]
-
double I_rad_ext_s[NU_DIM_MAX + 1]
-
double I_rad_ext_D[200 + 1]
-
double D[200 + 1]
-
double C_e[NU_DIM_MAX + 1]
-
double C_a[NU_DIM_MAX + 1]
-
double C_e1[NU_DIM_MAX + 1]
-
double C_a1[NU_DIM_MAX + 1]
-
double F_IC_disk[NU_DIM_MAX + 1]
-
double F_IC_tot[NU_DIM_MAX + 1]
-
double NU_IC_disk[NU_DIM_MAX + 1]
-
double I_SYN_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double I_SYN_JET_EDGE[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double I_SYN_JET_BASE[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double J_SYN_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double K_ESA_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double F_SYN_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double I_COM_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double I_COM_JET_EDGE[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double I_COM_JET_BASE[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double J_COM_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double J_COM_EXT_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double K_ABS_SSC_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double K_ABS_COM_EXT_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double F_COM_JET[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double I_SYN_TOT[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double I_SYN_BLOB_TOT[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double I_COM_TOT[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double F_SYN_BLOB[SL_DIM_MAX + 1][NU_DIM_MAX + 1]
-
double Sl[SL_DIM_MAX + 1][SL_DIM_MAX + 1]
-
double Sr[SL_DIM_MAX + 1][SL_DIM_MAX + 1]
-
double S[SL_DIM_MAX + 1][SL_DIM_MAX + 1]
-
double Dcut[SL_DIM_MAX + 1][SL_DIM_MAX + 1]
-
int ifexist(char *fname)
Functions and Variables in processes_supp_core
Slightly modified version of processes_supp_core.cpp for use in MCMC code; modified by Sarah Youngquist, Feb 2022 (see below) Only change is switching to using namespace bj_core02 and including bj_core.h instead of bj02 and bj02.h
Functions
-
double N_p_PowLawExpCutoff(double energy)
Calculates the value of N_p_PowLawExpCutoff function for a given energy.
This function takes the energy in electron volts (eV) as input and calculates the value of N_p_PowLawExpCutoff function. The N_p_PowLawExpCutoff function is defined as fNorm_prot multiplied by energy raised to the power of -NP and multiplied by the exponential of -energy divided by EPCUTOFF.
- Parameters:
energy – The energy in eV for which to calculate the N_p_PowLawExpCutoff.
- Returns:
The calculated value of N_p_PowLawExpCutoff for the given energy.
-
double N_p(double energy)
Calculates the value of N_p function for a given energy.
This function takes the energy in electron volts (eV) as input and calculates the value of the N_p function. The N_p function is defined as the result of the N_p_PowLawExpCutoff function.
- Parameters:
energy – The energy in eV for which to calculate the N_p value.
- Returns:
The calculated value of N_p for the given energy.
-
double J_p_PowLawExpCutoff(double energy)
Calculates the value of J_p_PowLawExpCutoff for a given energy.
This function calculates the value of J_p_PowLawExpCutoff using the provided energy value. The energy is expected to be in electron volts (eV).
- Parameters:
energy – The energy value for which to calculate the J_p_PowLawExpCutoff.
- Returns:
The calculated value of J_p_PowLawExpCutoff.
- Pre:
The variables Kp, NP, and EPCUTOFF must be defined and have valid values.
-
double J_p(double energy)
Calculates the value of J_p for a given energy.
This function calculates the value of J_p using the provided energy value. The energy is expected to be in electron volts (eV). The function internally calls J_p_PowLawExpCutoff to perform the calculation.
See also
- Parameters:
energy – The energy value for which to calculate the J_p.
- Returns:
The calculated value of J_p.
- Pre:
The variables Kp, NP, and EPCUTOFF must be defined and have valid values.
-
double Simpson(double func[], double ics[], int res, int start, int end)
Modified Simpson Integration Routine.
This function calculates the integral of a given function using the Modified Simpson integration method.
- Parameters:
func – Pointer to the array storing the function values
ics – Pointer to the array storing the x-values of the function
res – Length of the arrays func[] and ics[]
start – Index of the starting point in the arrays func[] and ics[]
end – Index of the ending point in the arrays func[] and ics[]
- Returns:
The integral value calculated by the Modified Simpson integration method
-
void gauleg(double x1, double x2, double x[], double w[], int n)
This function calculates the weights for a Gauss-Legendre integrating routine.
The gauleg function calculates the weights for a Gauss-Legendre integrating routine using the specified lower and upper bounds, number of points, and arrays to store the calculated values. The weights are calculated using the Gauss-Legendre method.
Remark
This function assumes that the arrays x and w have already been allocated with enough memory to store the calculated values.
Remark
The calculated points and weights are stored in the arrays x and w respectively, in the range from index 1 to n. The values at index 0 and index n+1 are not used.
Remark
The tolerance for convergence is set to EPS = 3.0e-11. If the maximum absolute difference between z and z1 is smaller than EPS, the iteration is considered to have converged and the loop is exited.
- Parameters:
x1 – The lower bound of the interval.
x2 – The upper bound of the interval.
x – Pointer to the array to store the calculated points.
w – Pointer to the array to store the calculated weights.
n – The number of points to calculate.
- Returns:
void
-
double intgl(double (*func)(double), double a, double b, int n, int m)
This function calculates the integral of a given function using a combined integrating routine. The routine uses both the trapezoid rule and Gauss-Legendre integration.
The function first checks if the lower limit ‘a’ is greater than or equal to the upper limit ‘b’. If so, it returns 0 and prints an error message. Otherwise, it calculates the step size based on the number of steps ‘n’ and the logarithms of ‘a’ and ‘b’.
Then, it iterates from 1 to ‘n-1’ and calculates the lower and upper limits for each subinterval. The gauleg function is called to calculate the Gauss-Legendre points and weights for each subinterval. The function ‘func’ is evaluated at each point and multiplied by the corresponding weight. The evaluations are summed up and multiplied by the step size to calculate the integral for each subinterval. Finally, all the subinterval integrals are summed up to obtain the total integral, which is returned as the result.
See also
Note
The function assumes that the provided function ‘func’ is continuous within the given limits ‘a’ and ‘b’. The function ‘gauleg’ is used internally to calculate the Gauss-Legendre points and weights for each subinterval.
- Parameters:
func – A pointer to the function to be integrated.
a – The lower limit of integration.
b – The upper limit of integration.
n – The number of steps in the trapezoid rule.
m – The number of points in the Gauss-Legendre quadrature.
- Returns:
The value of the integral.
-
double linint(double x, double xvec[], int xdim, double x_min, double x_max)
Perform linear interpolation on a log-log array.
This function performs linear interpolation on a log-log array. The x value provided is used to determine the interpolated value between two adjacent values in the xvec array. The x values in the xvec array must be sorted in ascending order. If the given x value is less than or equal to zero, it is replaced with the x_min value. If any input parameters are negative or if x_min is equal to x_max, an error message is displayed and the program exits. The interpolated value is calculated based on the provided x and xvec values by determining the position of x in the xvec array and using the surrounding values for interpolation. The interpolation is performed in logarithmic space if the loglog flag is set, otherwise it is performed in linear space. The resulting interpolated value is returned.
- Parameters:
x – The value to interpolate at.
xvec – The array of x values to interpolate between.
xdim – The size of the xvec array.
x_min – The minimum x value in the xvec array.
x_max – The maximum x value in the xvec array.
- Returns:
The interpolated value at x.
-
double j_syn(double (*elec_spec)(double), double gamma_min, double gamma_max, double nu, double B, int prec1, int prec2)
Calculates the synchrotron emission coefficient.
This function calculates the synchrotron emission coefficient based on the given parameters. It integrates over the electron energy spectrum defined by elec_spec and returns the result.
See also
elec_spec
Note
The function prints error messages and returns 0 if any of the input values are out of range.
- Parameters:
elec_spec – Function pointer to the electron energy spectrum function.
gamma_min – The minimal energy of electrons (Lorentz factor).
gamma_max – The maximal energy of electrons (Lorentz factor).
nu – The frequency for which the spectrum will be calculated.
B – The value of the magnetic field.
prec1 – The integration precision for trapezoid integration.
prec2 – The integration precision for Gauss-Legendre.
- Returns:
The synchrotron emission coefficient.
-
double k_esa(double (*elec_spec)(double), double gamma_min, double gamma_max, double nu, double B, int prec1, int prec2)
Calculate the electrons self-absorption coefficient.
This function calculates the self-absorption coefficient of electrons for a given energy spectrum, frequency, and magnetic field. It performs numerical integration to compute the coefficient.
Note
Make sure to provide accurate and valid values for the parameters as specified.
Note
The integration precisions (prec1, prec2) should be between 3 and 1.0e+4.
Note
The gamma_min parameter should be between 1.0 and 1.0e+10.
Note
The gamma_max parameter should be between 1.0 and 1.0e+10 and should be greater than gamma_min.
Note
The nu parameter should be between MIN_FREQ and MAX_FREQ.
Note
The B parameter should be between 1.0e-10 and 1.0e+5.
Note
The elec_spec function should be defined to provide the electrons energy spectrum.
Note
The function uses several constants: MIN_FREQ, MAX_FREQ, e, m_e, c, C2, sig_T, C1, C3. Make sure these constants are defined in the code.
- Parameters:
elec_spec – Pointer to the function that defines the electrons energy spectrum
gamma_min – Minimal energy of electrons (Lorentz factor)
gamma_max – Maximal energy of electrons (Lorentz factor)
nu – Frequency for which the spectrum will be calculated
B – Value of the magnetic field
prec1 – Integration precision for the trapezoid integration
prec2 – Integration precision for the Gauss-Legendre method
- Returns:
The self-absorption coefficient
-
double Compton_kernel(double ec, double gg, double es)
Calculates the inverse-Compton emissivity.
This function calculates the inverse-Compton emissivity based on the given parameters.
- Parameters:
ec – Energy of the electron
gg – Lorentz factor of the electron
es – Energy of the soft photon
- Returns:
Computed inverse-Compton emissivity
-
double j_com(double (*elec_spec)(double), double gamma_min, double gamma_max, double I_rad[], double nu_rad_min, double nu_rad_max, int nu_rad_dim, double nu, int prec1, int prec2)
Calculates the inverse Compton emission coefficient.
This function calculates the inverse Compton emission coefficient using the Inoue & Takahara 1996 formula. Given the input parameters and the electron energy spectrum, it calculates the emission coefficient for a given frequency.
See also
The Inoue & Takahara 1996 paper for detailed information on the formula.
Note
Before calling this function, make sure to define the following constants:
MIN_FREQ: Minimum frequency value (double)
MAX_FREQ: Maximum frequency value (double)
h: Planck’s constant (double)
m_e: Electron mass (double)
c: Speed of light (double)
Note
The function uses Gauss-Legendre quadrature for numerical integration.
- Parameters:
elec_spec – Function pointer to the electron energy spectrum function. The function should take a double as an argument and return a double.
gamma_min – Minimal energy of electrons (Lorentz factor).
gamma_max – Maximal energy of electrons (Lorentz factor).
I_rad – Array containing the intensity of the radiation field.
nu_rad_min – Minimum frequency of the radiation field photons.
nu_rad_max – Maximum frequency of the radiation field photons.
nu_rad_dim – Dimension of the matrix I_rad.
nu – Frequency for which the emission spectrum will be calculated.
prec1 – Integration precision for the trapezoid integration.
prec2 – Integration precision for the Gauss-Legendre.
- Returns:
The inverse Compton emission coefficient as a double value.
-
double Interpolate1(double ee, double zz, int print_flag)
Interpolates a value using bilinear interpolation.
This function computes the interpolated value of a given point (ee, zz) using bilinear interpolation. It checks if the given point falls within the range of the provided data arrays E1 and Z1. If the point is within the range, it performs bilinear interpolation using the formula: t = a*ee*zz + b*zz + c*ee + d
- Parameters:
ee – The x-coordinate of the point to be interpolated.
zz – The y-coordinate of the point to be interpolated.
print_flag – A flag indicating whether to print debug information. If set to 1, debug information will be printed to stderr. If set to 0, no debug information will be printed.
- Returns:
The interpolated value of the given point (ee, zz). If the point is outside the range of the data arrays, it returns 0.
-
double Interpolate2(double ee, double zz, int print_flag)
Interpolates a value using a bilinear interpolation method.
This function takes an coordinate (ee, zz) and returns the interpolated value based on a bilinear interpolation method. The function checks if the coordinate lies within the valid range of the given datasets. If the coordinate is within the range, the interpolation is performed. Otherwise, the function returns 0.
- Parameters:
ee – The E coordinate.
zz – The Z coordinate.
print_flag – Flag to indicate whether to print debug information.
- Returns:
The interpolated value.
-
double tau_IRA_Kneiske(double nu, double zz, int range)
Calculates the optical depth tau based on the frequency nu, redshift zz, and range parameter.
- Parameters:
nu – - the frequency in nu
zz – - the redshift in zz
range – - the range parameter
- Returns:
the calculated optical depth tau
-
double Interpolate3(double ee, double zz, int print_flag)
Interpolate3 function calculates the interpolated value for a given input E and z.
The Interpolate3 function takes in two parameters ee (the input value for E) and zz (the input value for z), and calculates the interpolated value using the function: t = a * ee * zz + b * zz + c * ee + d The coefficients a, b, c, and d are calculated based on the surrounding data points and used to determine the interpolated value. If the input values ee and zz fall within the defined range of E and z values, the function will return the interpolated value. Otherwise, it returns 0.
The function uses the following data arrays:
E3 - An array of size EDIM_Franceschini containing the E values. EDIM_Franceschini - The size of the E3 array.
Z3 - An array of size ZDIM_Franceschini containing the z values. ZDIM_Franceschini - The size of the Z3 array.
Note: The values of E3, EDIM_Franceschini, Z3, and ZDIM_Franceschini are not shown in the generated documentation.
- Parameters:
ee – - The input value for E.
zz – - The input value for z.
print_flag – - A flag indicating whether to print debug information.
- Returns:
The interpolated value.
-
double tau_IRA_Franceschini(double nu, double zz)
Calculates the optical depth using the IRA Franseschini method.
This function calculates the optical depth using the IRA Franseschini method. It takes in two parameters:
nu: The frequency in hertz.
zz: The redshift.
It performs the following steps:
Convert the frequency from hertz to keV.
Convert the frequency from keV to GeV and TeV.
Interpolate the value of the optical depth using the Interpolate3 function.
The function returns the calculated optical depth.
Note: The function relies on the Interpolate3 function to perform the interpolation. The Interpolate3 function takes in two parameters ee (the input value for E) and zz (the input value for z), and calculates the interpolated value using the function: t = a * ee * zz + b * zz + c * ee + d The coefficients a, b, c, and d are calculated based on the surrounding data points and used to determine the interpolated value. If the input values ee and zz fall within the defined range of E and z values, the function will return the interpolated value. Otherwise, it returns 0.
The Interpolate3 function uses the following data arrays:
E3: An array of size EDIM_Franceschini containing the E values.
EDIM_Franceschini: The size of the E3 array.
Z3: An array of size ZDIM_Franceschini containing the z values.
ZDIM_Franceschini: The size of the Z3 array.
Note: The values of E3, EDIM_Franceschini, Z3, and ZDIM_Franceschini are not shown in the generated documentation.
- Parameters:
nu – - The frequency in hertz.
zz – - The redshift.
- Returns:
The optical depth.
-
double Interpolate4(double ee, double zz, int print_flag)
Interpolates a value using a 4-point interpolation algorithm.
- Parameters:
ee – The input value for E.
zz – The input value for z.
print_flag – A flag indicating whether to print debug information.
- Returns:
The interpolated value.
-
double tau_IRA_Finke(double nu, double zz)
Calculates the optical depth, tau, using the IRA_Finke interpolation method.
- Parameters:
nu – The frequency in Hz
zz – The redshift value
- Returns:
The calculated optical depth
-
double Interpolate5(double ee, double zz, int print_flag)
This function performs interpolation to estimate the value of a function at a given point (ee, zz). The function uses bilinear interpolation to estimate the value based on a 2D table of values. The interpolation is performed within a specific range defined by the input parameters.
Example usage:
double result = Interpolate5(0.01, 0.5, 1);
This function requires the following external variables to be defined:
E5 : An array representing the range of possible values for the first input parameter. Defined as double E5[EDIM_Franceschini17].
EDIM_Franceschini17 : The number of elements in the array E5. Defined as an integer constant.
Z5 : An array representing the range of possible values for the second input parameter. Defined as double Z5[ZDIM_Franceschini17].
ZDIM_Franceschini17 : The number of elements in the array Z5. Defined as an integer constant.
T5 : A 2D array representing the values of the function to be interpolated. Defined as double T5[EDIM_Franceschini17 * ZDIM_Franceschini17]. The values in the array represent the function values at specific (ee, zz) coordinate pairs. The indexing of the array is as follows:
T5[ei][zi] represents the value at E5[ei] and Z5[zi].
For correct usage, make sure to define the external variables before calling this function.
- Parameters:
ee – The value of the first input parameter.
zz – The value of the second input parameter.
print_flag – Determines if additional debug information should be printed. Set to non-zero to enable printing, zero to disable printing.
- Returns:
The estimated value of the function at the given point (ee, zz). If the inputs are outside the defined range, the function returns 0.
-
double tau_IRA_Franceschini17(double nu, double zz)
Calculates the absorption coefficient (tau) for a given frequency (nu) and redshift (zz) using the IRA_Franceschini17 model. The absorption coefficient represents the amount of radiation that is absorbed as it travels through the universe. The function performs interpolation to estimate the value of tau based on a 2D table of function values.
Example usage:
double tau = tau_IRA_Franceschini17(4.8e+14, 0.5);
- Parameters:
nu – The frequency in Hz.
zz – The redshift.
- Returns:
The absorption coefficient tau.
-
double a_i(int i_i, double z_i, int level)
Calculates the absorption coefficient using the IIR model (Stecker & Jager 1998).
- Parameters:
i_i – The index value to determine the row in the absorption coefficient matrix.
z_i – The value for which the absorption coefficient is calculated.
level – The radiation field level. Set to non-zero for a higher IR radiation field, and zero for a lower IR radiation field.
- Returns:
The absorption coefficient for the given parameters.
-
double tau_IIR(double nu, double z, int level)
Calculates the absorption coefficient that describes the absorption of gamma-rays by the intergalactic infrared medium.
- Parameters:
nu – The observed frequency of the gamma-ray.
z – The redshift value.
level – The absorption level. Use 1 for high absorption level and 0 for low absorption level.
- Returns:
The absorption coefficient for the given parameters.
-
double sigma_gg(double nu_s, double nu_c)
Calculates the complete photon-photon cross-section according to the formula given by Aharonian et al. (2008)
- Parameters:
nu_s – The target photon frequency
nu_c – The High-Energy Photon Frequency
- Returns:
The cross-section value
-
double gg_abs(double nu_c, double I_c[], int nu_dim, double nu_min, double nu_max, int prec1, int prec2)
Calculates the absorption coefficient for pair production between VHE gamma rays and synchrotron radiation.
- Parameters:
nu_c – Frequency for which the coefficient will be calculated.
I_c – Array which contains intensity of radiation field.
nu_dim – Dimension of matrix I_rad.
nu_min – Minimum frequency of radiation field photons.
nu_max – Maximum frequency of radiation field photons.
prec1 – Integration precision for trapezoid integration.
prec2 – Integration precision for Gauss-Legendre.
- Returns:
The calculated absorption coefficient.
-
double Planck(double nu_BB, double T_BB)
Calculate the Planck function at a given frequency and temperature. Planck law, to take into account the external inverse compton. Radiation of the infrared photons reprocessed by the broad line regions
- Parameters:
nu_BB – Frequency for which the Planck function will be calculated.
T_BB – Temperature of the black body.
- Returns:
The intensity of the black body radiation in erg/s/sr/m²/Hz. Returns 0.0 if the input values are outside the valid range.
-
double CylTransfEquat(double I_inp, double jj, double kk, double ll)
Analytical solution of the transfer equation for cylindrical geometry with constant emission and absorption coefficients along the emitting region.
This function calculates the combined intensity after absorption and radiation for a given input intensity, emission coefficient, absorption coefficient, and length of the emitting region. It first calculates the optical depth using the length of the emitting region and the absorption coefficient. Based on the value of the optical depth, it calculates the absorbed intensity using the input intensity and the exponential decay factor. If the emission coefficient is greater than a very small threshold, it also calculates the radiated intensity using the emission and absorption coefficients.
The implementation of this function accounts for three different cases depending on the value of the optical depth:
If the optical depth is below a very small threshold, the absorbed intensity is equal to the input intensity. In this case, the radiated intensity is not considered.
If the optical depth is above a large threshold, both the absorbed and radiated intensities are set to zero.
If the optical depth is between the two thresholds, the absorbed intensity is calculated by multiplying the input intensity with the exponential decay factor, and the radiated intensity is calculated using the emission and absorption coefficients.
- Parameters:
I_inp – The input intensity.
jj – The emission coefficient.
kk – The absorption coefficient.
ll – The length of the emitting region.
- Returns:
The combined intensity after absorption and radiation.
-
double SphTransfEquat(double jj, double kk, double ll)
Analytical solution of the transfer equation for spherical geometry.
This function calculates the intensity of radiation given the emission and absorption coefficients along with the length of the emitting region. The solution assumes constant emission and absorption along the emitting region.
Remark
The intensity of radiation is computed using the following conditions:
If the emission coefficient is less than 1.0e-300 or has a value of “-inf”, “inf”, or “nan”, the intensity is set to 0.0.
Otherwise, the intensity is calculated based on the length of the emitting region and the absorption coefficient. If the length is very large (tau > 7.00e+2), the intensity is set to jj divided by kk. If the length is very small (tau < 1.0e-4), the intensity is set to (4/3) times the emission coefficient times the length. Otherwise, the intensity is calculated using a formula that considers the relation between emission and absorption.
See also
SphTransfEquat2()
Note
The emission coefficient and absorption coefficient must be positive numbers.
Warning
The function does not handle cases where the emission coefficient is negative or zero, or where the absorption coefficient is zero.
- Parameters:
jj – The emission coefficient.
kk – The absorption coefficient.
ll – The length of the emitting region.
- Returns:
The intensity of radiation.
-
double Intens2Flux(double Intens, double Radius, double Doppler, double z, double Hubble)
Calculates the flux from the intensity based on given parameters.
The following constants are referenced in the calculation:
c: The speed of light in centimeters per second. Defined as: const double c = 2.997924 * 1.0e+10;
The function calculates the Hubble constant H_0 by dividing ‘Hubble’ by 3.086 and multiplying by 1.0e-19. It then calculates the luminosity distance D_L using the cosmological formula and the Hubble constant. Finally, it calculates the flux using the given formula and returns the result.
- Parameters:
Intens – The intensity.
Radius – The radius.
Doppler – The Doppler shift.
z – The redshift.
Hubble – The Hubble constant.
- Returns:
The calculated flux.
- Pre:
The value of ‘Radius’ must be between 1.0 and 1.0e+100 (inclusive).
- Pre:
The value of ‘Doppler’ must be between 1.0 and 1.0e+2 (inclusive).
- Pre:
The value of ‘z’ must be greater than 0.0 and less than or equal to 5.0.
- Pre:
The value of ‘Hubble’ must be between 50.0 and 100.0 (inclusive).
- Post:
If the input parameters are outside the valid ranges, the function prints an error message and returns 0.0.
-
double CylIntens2Flux(double Intens1, double Intens2, double Radius, double Length, double Doppler, double z, double Hubble, double Theta)
Converts light intensity values to flux.
This function takes the following parameters:
Intens1: Light intensity value 1.
Intens2: Light intensity value 2.
Radius: Radius of the cylinder.
Length: Length of the cylinder.
Doppler: Doppler value.
z: Redshift value.
Hubble: Hubble constant value.
Theta: Angle value in degrees.
The function first validates the input parameters to ensure they fall within the expected ranges. It validates the Radius should be greater than or equal to 1.0 and less than or equal to 1.0e+100. It validates the Doppler should be greater than or equal to 1.0 and less than or equal to 1.0e+2. It validates the z should be greater than 0.0 and less than or equal to 5.0. It validates the Hubble should be greater than or equal to 50.0 and less than or equal to 100.0. It validates the Theta should be greater than or equal to 0.0 and less than or equal to 90.0.
If any of the parameters fail the validation, an error message will be printed to the standard error output, and the function will return 0.0.
If the parameters pass the validation, the function calculates the value of H_0 based on the Hubble constant, and D_L based on the redshift value. Then it calculates and returns the flux value using the input parameters and intermediate calculations.
Note
The value of c is a constant equal to 2.997924 * 1.0e+10.
Note
The function assumes that the M_PI constant is defined and represents the value of pi.
- Parameters:
Intens1 – Light intensity value 1.
Intens2 – Light intensity value 2.
Radius – Radius of the cylinder.
Length – Length of the cylinder.
Doppler – Doppler value.
z – Redshift value.
Hubble – Hubble constant value.
Theta – Angle value in degrees.
- Returns:
The calculated flux based on the input parameters.
-
double RingIntens2Flux(double Intens, double InnRadius, double OutRadius, double Doppler, double z, double Hubble, int check)
Convert ring intensity to flux.
This function calculates the flux from a ring intensity based on various parameters.
See also
c
Note
The function checks if the given values for InnRadius, OutRadius, Doppler, z, and Hubble are within valid ranges and outputs an error message if any of them are invalid.
- Parameters:
Intens – The intensity of the ring.
InnRadius – The inner radius of the ring.
OutRadius – The outer radius of the ring.
Doppler – The Doppler value.
z – The redshift value.
Hubble – The Hubble value.
check – A flag indicating if error checking should be performed.
- Returns:
The flux calculated from the given ring intensity.
-
double FreqTransS2O(double nu, double Doppler, double z)
Perform frequency transformation from the source frame to the observer frame.
This function takes the source frequency, Doppler factor, and redshift as input, and calculates the transformed frequency. The transformed frequency is obtained by multiplying the Doppler factor with the source frequency and dividing it by the sum of 1 and the redshift. If the input values do not fall within the specified range, an error message is printed and 0.0 is returned.
Note
The function assumes that the source frequency is in the range of 1.0e+5 to 1.0e+35.
Note
The function assumes that the Doppler factor is in the range of 0.0 to 1.0e+2.
Note
The function assumes that the redshift is in the range of 0.0 to 5.0.
- Parameters:
nu – The source frequency.
Doppler – The Doppler factor.
z – The redshift.
- Returns:
The transformed frequency.
-
double FreqTransO2S(double nu, double Doppler, double z)
This function performs a frequency transformation from the observer frame to the source frame.
The function takes in three parameters:
nurepresents the source frequency.Doppleris the Doppler factor.zis the redshift.
Note
The function expects
nuto be in the range between 1.0e+5 and 1.0e+35,Dopplerto be in the range between 1.0 and 1.0e+2, andzto be in the range between 0.0 and 5.0. If any of the parameters fall outside these ranges, an error message is printed to stderr and 0.0 is returned.- Parameters:
nu – The source frequency.
Doppler – The Doppler factor.
z – The redshift.
- Returns:
The transformed frequency in the source frame.
-
double LuminDist(const double z)
Calculates the luminosity distance based on the redshift value.
The calculation is based on the following formula:
dl = c/H0 * ( z + (z*z*(1.-q0))/(1.+q0*z+sqrt(1.+2.*q0*z)) )
where:
dl is the luminosity distance
c is the speed of light constant
H0 is the Hubble constant multiplied by the relevant conversion factors for distance units
z is the redshift value
q0 is the deceleration parameter
- Parameters:
z – The redshift value
- Returns:
The calculated luminosity distance
-
double Distance_Luminosity(double z, double H0, double WM)
Calculates the distance luminosity in cm.
This function calculates the distance luminosity in centimeters based on the redshift (z), Hubble constant (H0), and Omega matter (WM). The function is adapted from a script by Ned Wright. https://www.astro.ucla.edu/~wright/CosmoCalc.html
- Parameters:
z – The redshift.
H0 – The Hubble constant.
WM – The Omega matter.
- Returns:
The distance luminosity in cm.
-
void setHadronicParameters()
This function sets the parameters for hadronic interactions, specifically for pion decay.
The function calculates and sets various parameters needed for the hadronic interactions’ calculations. It computes the energetics of electrons to normalize the proton spectrum, calculates the normalization of the proton spectrum from ENERGY50, and sets other related parameters.
The function does not take any input parameters and does not return any values. It modifies the global variables ENERGY50, fEnergy_prot, fNorm_prot, and fNb_prot.
-
double piondecay(double energy)
Calculates the flux of gamma ray spectra from the decay of pions. The calculations are based on equations from Kelner, Aharonian & Bugayov, PhRvD, 74, 034018, 2006. The function takes the energy of the gamma ray as input and returns the calculated flux.
- Parameters:
energy – The energy of the gamma ray in eV.
- Returns:
The flux of gamma ray spectra in erg cm-2 s-1.
-
void setHadronicParametersTest()
Calculates the hadronic parameters for testing purposes.
This function calculates the normalization of the proton spectrum and sets the values of the fNorm_prot and Kp variables. It performs an integration from 1 TeV to infinity to calculate the normalization. The integration is done using a step size of 1/(per_decade) and the result is stored in the fNorm_prot variable. The Kp variable is also set to the same value. The integration is done using the N_p function, which calculates the value of the N_p function for a given energy.
- Returns:
None.
-
double piondecaytest(double energy)
Calculates the decay test for a given energy.
This function calculates the decay test for a given energy using equations from the paper “Kelner, Aharonian & Bugayov, PhRvD, 74, 034018, 2006: 2006PhRvD..74c4018K”. The function takes the energy in electron volts (eV) as input and returns the result of the decay test.
- Parameters:
energy – The energy in eV for which to calculate the decay test.
- Returns:
The result of the decay test for the given energy.
-
double piondecay_ahaprecise_electron(double energy)
Calculates the energy spectrum for electron emission from pion decay.
This function calculates the energy spectrum for electron emission resulting from pion decay. The input energy parameter represents the electron energy in electron volts (eV). The output is the calculated energy spectrum for the given electron energy.
- Parameters:
energy – The energy of the electrons in eV for which to calculate the energy spectrum.
- Returns:
The calculated energy spectrum for the given electron energy.
-
void energetics()
Calculate the energetics of the system.
This function calculates the energetics of the electrons and protons in the system. It computes various parameters such as the number of particles, energy density, equipartition parameter, mean energy, and total energy. It also prints out the derived parameters if the PRINT flag is set to 1.
- Returns:
void
Variables
-
const double C1 = 0.78
-
const double C2 = 0.25
-
const double C3 = 2.175
-
const double MIN_FREQ = 1.0e+4
-
const double MAX_FREQ = 1.0e+40
-
double fDistance
in pc
-
double fBtesla
in Tesla
-
double Ksync
Facteur Ksync*E*E pour le calcul des pertes d’energie synchrotron.
-
double Kcompt
Facteur Kcompt*E*E pour le calcul des pertes d’energie ic.
-
double Kloss1
Facteur Kloss1*E*E pour le calcul des pertes totales en E*E.
-
double Kbremss
Facteur Kbremss*E pour le calcul des pertes d’energie bremsstrahlung.
-
double Rcar
-
double Tcar
-
double Vcar
-
double fVelocity
-
double fEmin
in eV
-
double fEmax
in eV
-
double fEcut
exp. cutoff in eV
-
double fSlope
slope of the electronic distribution
-
double fEnergy50
in 10^50 ergs
-
double fNorm
spectrum normalization
-
double fEnergy
average electron energy
-
double fNbelec
no. of electrons
-
double fEnergy_prot
av. proton energy
-
double fNb_prot
no. of protons
-
double fFluxfact
per m^2 1/(4*PI*fDistance^2)
-
int NBSECONDARYELECTRONS
fElectrons; nb of secondary electrons for the synchrotron process
-
double fEnergy_electrons[10000]
energy of the secondary electrons produced by pp interaction
-
double fFlux_electrons[10000]
Flux of the secondary electrons produced by pp interaction.
-
double fPerDecade
Numbers of steps per decade used for the integration in energy for radiative process.
-
double Kp
Normalization of the proton spectrum, in terms of Lorentz factor [cm^-3].
-
double ENERGY50
total energy in protons [10^50 erg]
-
double EPMIN
minimum proton energy [eV]
-
double EPMAX
maximum proton energy [eV]
-
double EPCUTOFF
minimum proton energy [eV]
-
double GAMMA_P_MAX
maximum proton energy in terms of Lorentz factor [E/mp c^2]
-
double SLOPEP
slope of proton energy spectrum
-
double NP
slope of proton energy spectrum
-
double ATOMICDENSITY
atomic density of target protons [cm^-3]
-
double SIZE
radius of emitting zone [pc]
-
double DISTANCE
distance of the source [pc]
-
double AGE
Age of the source [s].
-
double fNorm_prot
-
double RATIO
ratio proton over electron number (ratio of total number of particles)
-
int CASE_PION
-
const int LDIM = 51
-
const int MDIM = 2601
51x51;
-
double Z1[LDIM] = {0.000, 0.002, 0.004, 0.006, 0.008, 0.010, 0.012, 0.014, 0.016, 0.018, 0.020, 0.022, 0.024, 0.026, 0.028, 0.030, 0.032, 0.034, 0.036, 0.038, 0.040, 0.042, 0.044, 0.046, 0.048, 0.050, 0.052, 0.054, 0.056, 0.058, 0.060, 0.062, 0.064, 0.066, 0.068, 0.070, 0.072, 0.074, 0.076, 0.078, 0.080, 0.082, 0.084, 0.086, 0.088, 0.090, 0.092, 0.094, 0.096, 0.098, 0.100}
-
double E1[LDIM] = {-1.000, -0.820, -0.640, -0.460, -0.280, -0.100, 0.080, 0.260, 0.440, 0.620, 0.800, 0.980, 1.160, 1.340, 1.520, 1.700, 1.880, 2.060, 2.240, 2.420, 2.600, 2.780, 2.960, 3.140, 3.320, 3.500, 3.680, 3.860, 4.040, 4.220, 4.400, 4.580, 4.760, 4.940, 5.120, 5.300, 5.480, 5.660, 5.840, 6.020, 6.200, 6.380, 6.560, 6.740, 6.920, 7.100, 7.280, 7.460, 7.640, 7.820, 8.000}
-
double Z2[LDIM] = {0.000, 0.100, 0.200, 0.300, 0.400, 0.500, 0.600, 0.700, 0.800, 0.900, 1.000, 1.100, 1.200, 1.300, 1.400, 1.500, 1.600, 1.700, 1.800, 1.900, 2.000, 2.100, 2.200, 2.300, 2.400, 2.500, 2.600, 2.700, 2.800, 2.900, 3.000, 3.100, 3.200, 3.300, 3.400, 3.500, 3.600, 3.700, 3.800, 3.900, 4.000, 4.100, 4.200, 4.300, 4.400, 4.500, 4.600, 4.700, 4.800, 4.900, 5.000}
-
double E2[LDIM] = {-1.000, -0.820, -0.640, -0.460, -0.280, -0.100, 0.080, 0.260, 0.440, 0.620, 0.800, 0.980, 1.160, 1.340, 1.520, 1.700, 1.880, 2.060, 2.240, 2.420, 2.600, 2.780, 2.960, 3.140, 3.320, 3.500, 3.680, 3.860, 4.040, 4.220, 4.400, 4.580, 4.760, 4.940, 5.120, 5.300, 5.480, 5.660, 5.840, 6.020, 6.200, 6.380, 6.560, 6.740, 6.920, 7.100, 7.280, 7.460, 7.640, 7.820, 8.000}
-
double Z_MIN = 0.0001
-
double Z_MAX = 0.0999
-
double E_MIN = log10(500)
500 GeV
-
double E_MAX = log10(5.0e+4)
50 TeV
-
const int IDIM = 100
-
const int ZDIM_Franceschini = 9
-
const int EDIM_Franceschini = 50
-
const int MDIM_Franceschini = 450
-
double Z3[ZDIM_Franceschini] = {0.01, 0.03, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0}
-
double E3[EDIM_Franceschini] = {0.0200, 0.0240, 0.0289, 0.0347, 0.0417, 0.0502, 0.0603, 0.0726, 0.0873, 0.104, 0.126, 0.151, 0.182, 0.219, 0.263, 0.316, 0.381, 0.458, 0.550, 0.662, 0.796, 0.957, 1.15, 1.38, 1.66, 2.000, 2.40, 2.89, 3.47, 4.17, 5.02, 6.03, 7.26, 8.73, 10.5, 12.6, 15.2, 18.2, 21.9, 26.4, 31.7, 38.1, 45.8, 55.1, 66.2, 79.6, 95.7, 115., 138., 166.}
-
double T3[MDIM_Franceschini]
-
const int ZDIM_Finke = 500
-
const int EDIM_Finke = 99
-
const int MDIM_Finke = 49500
-
double Z4[ZDIM_Finke] = {0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2, 1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3, 1.31, 1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4, 1.41, 1.42, 1.43, 1.44, 1.45, 1.46, 1.47, 1.48, 1.49, 1.5, 1.51, 1.52, 1.53, 1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.6, 1.61, 1.62, 1.63, 1.64, 1.65, 1.66, 1.67, 1.68, 1.69, 1.7, 1.71, 1.72, 1.73, 1.74, 1.75, 1.76, 1.77, 1.78, 1.79, 1.8, 1.81, 1.82, 1.83, 1.84, 1.85, 1.86, 1.87, 1.88, 1.89, 1.9, 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97, 1.98, 1.99, 2.0, 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08, 2.09, 2.1, 2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.18, 2.19, 2.2, 2.21, 2.22, 2.23, 2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.3, 2.31, 2.32, 2.33, 2.34, 2.35, 2.36, 2.37, 2.38, 2.39, 2.4, 2.41, 2.42, 2.43, 2.44, 2.45, 2.46, 2.47, 2.48, 2.49, 2.5, 2.51, 2.52, 2.53, 2.54, 2.55, 2.56, 2.57, 2.58, 2.59, 2.6, 2.61, 2.62, 2.63, 2.64, 2.65, 2.66, 2.67, 2.68, 2.69, 2.7, 2.71, 2.72, 2.73, 2.74, 2.75, 2.76, 2.77, 2.78, 2.79, 2.8, 2.81, 2.82, 2.83, 2.84, 2.85, 2.86, 2.87, 2.88, 2.89, 2.9, 2.91, 2.92, 2.93, 2.94, 2.95, 2.96, 2.97, 2.98, 2.99, 3.0, 3.01, 3.02, 3.03, 3.04, 3.05, 3.06, 3.07, 3.08, 3.09, 3.1, 3.11, 3.12, 3.13, 3.14, 3.15, 3.16, 3.17, 3.18, 3.19, 3.2, 3.21, 3.22, 3.23, 3.24, 3.25, 3.26, 3.27, 3.28, 3.29, 3.3, 3.31, 3.32, 3.33, 3.34, 3.35, 3.36, 3.37, 3.38, 3.39, 3.4, 3.41, 3.42, 3.43, 3.44, 3.45, 3.46, 3.47, 3.48, 3.49, 3.5, 3.51, 3.52, 3.53, 3.54, 3.55, 3.56, 3.57, 3.58, 3.59, 3.6, 3.61, 3.62, 3.63, 3.64, 3.65, 3.66, 3.67, 3.68, 3.69, 3.7, 3.71, 3.72, 3.73, 3.74, 3.75, 3.76, 3.77, 3.78, 3.79, 3.8, 3.81, 3.82, 3.83, 3.84, 3.85, 3.86, 3.87, 3.88, 3.89, 3.9, 3.91, 3.92, 3.93, 3.94, 3.95, 3.96, 3.97, 3.98, 3.99, 4.0, 4.01, 4.02, 4.03, 4.04, 4.05, 4.06, 4.07, 4.08, 4.09, 4.1, 4.11, 4.12, 4.13, 4.14, 4.15, 4.16, 4.17, 4.18, 4.19, 4.2, 4.21, 4.22, 4.23, 4.24, 4.25, 4.26, 4.27, 4.28, 4.29, 4.3, 4.31, 4.32, 4.33, 4.34, 4.35, 4.36, 4.37, 4.38, 4.39, 4.4, 4.41, 4.42, 4.43, 4.44, 4.45, 4.46, 4.47, 4.48, 4.49, 4.5, 4.51, 4.52, 4.53, 4.54, 4.55, 4.56, 4.57, 4.58, 4.59, 4.6, 4.61, 4.62, 4.63, 4.64, 4.65, 4.66, 4.67, 4.68, 4.69, 4.7, 4.71, 4.72, 4.73, 4.74, 4.75, 4.76, 4.77, 4.78, 4.79, 4.8, 4.81, 4.82, 4.83, 4.84, 4.85, 4.86, 4.87, 4.88, 4.89, 4.9, 4.91, 4.92, 4.93, 4.94, 4.95, 4.96, 4.97, 4.98, 4.99}
-
double E4[EDIM_Finke] = {0.001, 0.001122018, 0.001258925, 0.001412538, 0.001584893, 0.001778279, 0.001995262, 0.002238721, 0.002511886, 0.002818383, 0.003162278, 0.003548134, 0.003981072, 0.004466836, 0.005011872, 0.005623413, 0.006309573, 0.007079458, 0.007943282, 0.008912509, 0.01, 0.01122018, 0.01258925, 0.01412538, 0.01584893, 0.01778279, 0.01995262, 0.02238721, 0.02511886, 0.02818383, 0.03162278, 0.03548134, 0.03981072, 0.04466836, 0.05011872, 0.05623413, 0.06309573, 0.07079458, 0.07943282, 0.08912509, 0.1, 0.1122018, 0.1258925, 0.1412538, 0.1584893, 0.1778279, 0.1995262, 0.2238721, 0.2511886, 0.2818383, 0.3162278, 0.3548134, 0.3981072, 0.4466836, 0.5011872, 0.5623413, 0.6309573, 0.7079458, 0.7943282, 0.8912509, 1.0, 1.122018, 1.258925, 1.412538, 1.584893, 1.778279, 1.995262, 2.238721, 2.511886, 2.818383, 3.162278, 3.548134, 3.981072, 4.466836, 5.011872, 5.623413, 6.309573, 7.079458, 7.943282, 8.912509, 10.0, 11.22018, 12.58925, 14.12538, 15.84893, 17.78279, 19.95262, 22.38721, 25.11886, 28.18383, 31.62278, 35.48134, 39.81072, 44.66836, 50.11872, 56.23413, 63.09573, 70.79458, 79.43282}
-
double T4[MDIM_Finke]
-
const int ZDIM_Franceschini17 = 10
-
const int EDIM_Franceschini17 = 56
-
const int MDIM_Franceschini17 = 560
-
double Z5[ZDIM_Franceschini17] = {0.01, 0.03, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0}
-
double E5[EDIM_Franceschini17] = {0.0052, 0.00631, 0.00767, 0.00932, 0.01132, 0.01375, 0.01671, 0.0203, 0.02466, 0.02997, 0.03641, 0.04423, 0.05374, 0.06529, 0.07932, 0.09636, 0.11708, 0.14224, 0.17281, 0.20995, 0.25507, 0.30989, 0.3765, 0.45742, 0.55573, 0.67516, 0.82027, 0.99657, 1.21076, 1.47098, 1.78712, 2.17122, 2.63787, 3.20481, 3.8936, 4.73042, 5.7471, 6.9823, 8.48296, 10.3061, 12.5211, 15.2122, 18.4817, 22.4539, 27.2798, 33.1429, 40.2661, 48.9203, 59.4344, 72.2083, 87.7276, 106.582, 129.489, 157.319, 191.131, 232.21}
-
double T5[MDIM_Franceschini17]