--+--------------------------------++--------------------------------+-- ==+================================++================================+== | PROGRAM COMPTON | | version 2 | | external specification | | procedure specifications | | | | written by Marek Gierlinski | | Marek.Gierlinski@gmail.com | | | | Version history: | | - started: 2 February 1996 | | - last update: 15 September 2008 | ==+================================++================================+== --+--------------------------------++--------------------------------+-- (not finished) ======================================================================== FORMAT OF THE INPUT FILE ======================================================================== NumberTrajectories = N WeightLimit = w Absorption = YES | NO AngularBins = cos1 cos2 ... cosn InputEnergyRange = Min Max ElectronEnergyRange = Min Max SpectrumEnergyRange = Min Max N logflag PathRange = Min Max N logflag PhotonTemperature = theta ElectronTemperature = T PhotonEnergyIndex = alpha ElectronEnergyIndex = s ThomsonThickness = tau Geometry = SLAB | SPHERE OutputPhotons = ALL | SCATTERED_ONLY Distribution = HOMOGENEOUS | CENTRAL | S_T | CONST_IRRADIATION | THIN_SOURCE | THICK_SOURCE InitialSpectrum = BLACK_BODY | POWER_LAW | FILE ElectronSpectrum = MAXWELL | POWER_LAW IrradiationAngle = eta0 OutputUnits = mc2 | keV InitialFile = filename TransmittedFile = filename ReflectedFile = filename PathFile = filename LogFile = filename Notes: - The order of lines is not important, names of the left side of '=' sign are important. - All text strings are case sensitive. - 'a | b | c' means: a or b or c (only one of them) - 'N' in the energy and path range denotes number of discrete points: the first point is 'Min', the last one is 'Max'. The number of subintervals is equal N - 1. - logflag might be LIN or LOG. - Photon energies are given in keV, both in linear and logarithmic scale. - Electron energies are given in kT units: E = (gamma - 1) m c^2 / (k T), T - electron temperature - Temperatures are given in keV. - Irradiation angle is cosine angle from 'z' axis. - Angular bins are defined by cos angle from the z axis (slab geometry only). For spherical geometry this is always assumed to be 0 1 (i.e. one bin) - The maximum number of trajectories is 2^64 - 1 = 1,844,674,407,370,955,161 - The input seed photon file (case of InitialSpectrum = FILE) is a photon spectrum output from xspec, i.e. a QDP file (first three lines are ignored) containing three columns: energy, delta_energy and photon_flux. Example: NumberTrajectories = 1000000 WeightLimit = 1e-7 InputEnergyRange = 0.001 2 ElectronEnergyRange = 0.00316 40 SpectrumEnergyRange = 0.001 2000 100 LOG ElectronSpectrum = MAXWELL ElectronTemperature = 170 ElectronEnergyIndex = 2 ThomsonThickness = 0.80 InitialSpectrum = BLACK_BODY PhotonTemperature = 0.1 PhotonEnergyIndex = 3.0 Geometry = SLAB OutputPhotons = SCATTERED_ONLY Distribution = THICK_SOURCE IrradiationAngle = 0.45 OutputUnits = keV InitialFile = bb.dat TransmittedFile = trans.dat ReflectedFile = refl.dat LogFile = compton.log ======================================================================== GLOBAL DATA STRUCTURES defined in 'defs.h' file ======================================================================== ------------------------------------------------------------------------ Distribution ------------------------------------------------------------------------ struct Distribution { double Tab[N_DISTR]; int N; double x; double logx; double L0; double x1L0; }; This structure contains tabulated INVERSE cumulative probability distribution. The distribution is tabulated in the array Tab, and indices run from 0 to 2 * N. Let's think about this distribution as of the function F(x). The range for x is of course [0, 1], because it is inversion of the cumulative distribution. But nodes on the x-axis (array elements) are NOT distributed homogeneously over the range [0,1]. The arrangement of the nodes is geometrical, as follows: 0.0 0.5 1.0 |--|---|----|-----|--...--|----------|----------|--...--|----|---|--| 0 1 2 3 4 ... N-1 N N+1 ... 2N-1 2N The length of the first interval, between indices 0 and 1, is L0 (this is also the field of the structure Distribution). The next one is of the length L0 * x, where x is some constant value a little greater than 1. The next is L0 * x * x, and so on: Li = L{i-1} * x. The length of the interval between indices i and i+1 is L0 * x^i. The position of the node 'i' on the [0, 1] range is then: r(i) = L0 + L1 + ... + L{i-1} = L0 * (1 + x + x^2 + ... + x^{i-1}) = = L0 * S{i-1}, where Si = (x^(i+1) - 1) / (x - 1) so: r(i) = L0 * (x^i - 1) / (x - 1) The constant value (x - 1) / L0 is stored as the field 'x1L0' of this structure. In order to find an index 'i' for a given randomly drawn value 'R' one has to invert the above formula: i = log(R * (x - 1) / L0 + 1) / log x The constant value log x is stored in the field 'logx' of this structure. The proper interpolation (e.g. quadratic) is required to obtain a precise value of F(R). ------------------------------------------------------------------------ SourceDescription ------------------------------------------------------------------------ struct SourceDescription { int Geometry; int Distribution; int InitialSpectrum; int OutputPhotons; double tau; double eta0; }; where: Geometry = SLAB | SPHERE Distribution = HOMOGENEOUS | CENTRAL | S_T | CONST_IRRADIATION | THIN_SOURCE | THICK_SOURCE InitialSpectrum = BLACK_BODY | POWER_LAW OutputPhotons = ALL | SCATTERED_ONLY tau - Thomson thermal thickness of the cloud eta0 - angle at which slab is irradiated (CONST_IRRADIATION case) ------------------------------------------------------------------------ Photon ------------------------------------------------------------------------ struct Photon { double x; double eta; double tau; double w; double x_init; double w_init; double Path; int Scatter; }; where: x = (h nu) / (m c^2) - dimensionless energy (m - electron mass) eta = cos(theta), theta - angle between 'z' axis and photon direction tau - optical depth in the cloud (one-dimensional position) w - statistical weight x_init - initial energy w_init - initial weight Path - total optical path travelled by the photon Scatter - scatter counter ------------------------------------------------------------------------ Electron ------------------------------------------------------------------------ struct Electron { double beta, mu, phi; }; where: beta = v / c - electron velocity mu - cosine theta, measured from the 'z' axis phi - other angle ------------------------------------------------------------------------ Range ------------------------------------------------------------------------ struct Range { double Min, Max; }; where: Min, Max - lower and higher limits This structure stores a range of a varying quantity. ------------------------------------------------------------------------ Scale ------------------------------------------------------------------------ struct Scale { double Min, Max; int N; double Delta; int LogFlag; }; where: Min, Max - lower and higher limits N - number of discrete steps Delta - the step, Delta = (Max - Min) / N LogFlag - =0 for linear scale, =1 for logarithmic scale This structure is used to convert between continuous physical quantities and discrete, tabulated values used in program. When LogFlag = LOG, the limits and 'Delta' are stored as logarithms. ------------------------------------------------------------------------ InputData ------------------------------------------------------------------------ struct InputData { unsigned long long int NumberTrajectories; double WeightLimit; int WritePhotons; int Absorption; struct Intervals Angles; struct Range RInput; struct Range RElectron; struct Scale SSpectrum; struct Scale SPathDistr; double PhotonTemp; double ElectronTemp; double PhotonEnergyIndex; double ElectronEnergyIndex; int ElectronSpectrum; struct SourceDescription Source; int OutputUnits; char SeedFileName[MAX_FILE_NAME]; char InitialFile[MAX_FILE_NAME]; char TransmittedFile[MAX_FILE_NAME]; char ReflectedFile[MAX_FILE_NAME]; char PathFile[MAX_FILE_NAME]; char LogFile[MAX_FILE_NAME]; }; ======================================================================== GLOBAL VARIABLES defined in 'compton.c' file ======================================================================== struct Distribution BlackBodyDistr; Contains tabulated black body inverse cumulative distribution, prepared by the function TabulateBBDistr. Tabulated values are photon energies, h nu / (m c^2), in the range from RInput.Min to RInput.Max. struct Distribution ElVelocityDistr[MAX_PHOTON_ENERGY]; Tabulated inverse cumulative distributions of electron velocities for a range of photon energies. Indices are from 0 to S2Photon.N - 1. Photon energies, h nu / (m c^2), are from S2Photon.Min to S2Photon.Max. Tabulated values are electron Lorentz gamma factors, in the range from RElectron.Min to RElectron.Max. This array is initialized by the function 'TabulateEVDistr', and used to generate an electron velocity. struct Distribution SunyaevTitarchukDistr; Tabulated inverse cumulative Sunyaev-Titarchuk distribution. Tabulated values are optical depths, in the range from 0 to the Thomson thickness of the source. struct PhotonSpectrum SeedFile; Energy spectrum of the seed photons in case of InitialSpectrum = FILE. The spectrum (photon spectrum) is normalized to its integral. double CrossSectionTab[MAX_PHOTON_ENERGY]; Contains tabulated mean cross section; prepared by the function TabulateCrossSection. The photon energy is h nu / (m c^2). double TransmittedSpectrum[MAX_PHOTON_ENERGY]; Counts the transmitted photons during the simulation. This is the case of the spherical geometry, or the side of the slab opposite to the irradiating source. Each time a photon escapes the cloud, the element of this array corresponding to a photon energy is increased by the photon statistical weight. The photon energy is h nu / (m c^2). double ReflectedSpectrum[MAX_PHOTON_ENERGY]; Counts the reflected photons during the simulation. In this case we count photons on the same side of the slab, where the irradiating source is located. Each time a photon escapes the cloud, the element of this array corresponding to a photon energy is increased by the photon statistical weight. The photon energy is h nu / (m c^2). double TransmittedError[MAX_PHOTON_ENERGY][MAX_ANGLE]; double ReflectedError[MAX_PHOTON_ENERGY][MAX_ANGLE]; Errors for output spectra. They are calculated as a running standard deviation. double PathDistr[MAX_PATH_DISTR][MAX_ANGLE]; Distribution of total optical paths (distance) travelled by the emerging photons. struct Intervals Angles; Angular bins for output spectra, defined by cos angle from the z-axis. it is [0, 1] by default for sphere. struct Range RInput, RSpectrum, RElectron; Ranges of quantities: input photon energy, admissible photon energy, electron gamma Lorentz factor. struct Scale SCross, S2Photon, SSpectrum, SPathDistr; Scales for tabulated values: cross section, photon energies for electron velocities, emerging spectra and total optical path distribution. double ElectronTemp; Dimensionless electron temperature = (k T) / (m c^2), where T - thermal temperature, m - electron mass. double PhotonTemp; Dimensionless photon temperature = (k T) / (m c^2), where T - radiation temperature, m - electron mass. double PhotonEnergyIndex; Energy index 'alpha' of the initial power-law spectrum. double ElectronEnergyIndex; Energy index 's' of the electron power-law distribution. int ElectronSpectrum; Electron spectrum: MAXWELL or POWER_LAW. int OutputUnits; Units of the output spectrum: mc2 or keV double LogMin, LogMax, PowerLawNorm; Auxiliary variables for power-law generation. ======================================================================== LIST OF THE FUNCTIONS ======================================================================== compton.c: CrossSection GenerateInitialPhoton NextPosition ElectronVelocity Scattering init.c Usage ConvertScale InitData TabulateBlackBodyDistr TabulateSunyaevTitarchukDistr TabulateElVelocityDistr TabulateCrossSection io.c: Error Warning Message OpenLogFile CloseLogFile ReadScale ReadRange RecognizeDescriptor Descriptor ReadData CopyData WriteInitialSpectrum WriteData auxil.c: Integral Idx Val TabulateInverseDistribution Draw DrawRnd NormalizeSpectra func.c: F InvF MaxwellDistr PowerLawDistr BlackBodySpectrum ElectronDistr PhotonDistr Ro1Prime SunyaevTitarchuk ======================================================================== FUNCTIONS defined in 'init.c' file ======================================================================== ------------------------------------------------------------------------ Usage ------------------------------------------------------------------------ PURPOSE: Writes short information about how to use program. DECLARATION: void Usage() DESCRIPTION: This function does what it should do. ------------------------------------------------------------------------ ConvertScale ------------------------------------------------------------------------ PURPOSE: Converts a scale to mc^2 units and do logarithms if necessary. DECLARATION: void ConvertScale(S); struct Scale *S; INPUT: S - a scale OUTPUT: S - converted scale DESCRIPTION: Convert 'Min' and 'Max' to mc^2 units, performs 10-logarithms when 'LogFlag' is set to LOG, and calculates 'Delta'. ------------------------------------------------------------------------ InitData ------------------------------------------------------------------------ PURPOSE: Initializes the converting scales, and some global variables. DECLARATION: void InitData(Data); struct InputData Data; INPUT: Data - structure containing input data OUTPUT: SInput, SS_T, S2Photon, S2Electron, SCross, SSpectrum, RSpectrum, PhotonTemp, ElectronTemp, ElectronEnergyIndex, OutputUnits, LogMin, LogMax, PowerLawNorm - global variables. DESCRIPTION: This function initializes the values in all converting scales, plus photon and electron temperatures, spectral index and some auxiliary variables. COMMENT: 'Data' structure should be initialized with 'ReadData' function. ------------------------------------------------------------------------ TabulateBlackBodyDistr ------------------------------------------------------------------------ PURPOSE: Tabulates the black body spectrum cumulative distribution. DECLARATION: void TabulateBlackBodyDistr() INPUT: PhotonTemp = (k T) / (m c^2) - dimensionless photon temperature (T - radiation temperature). This is a global variable. OUTPUT: BlackBodyDistr - tabulated distribution (global structure). DESCRIPTION: Tabulates the invention of cumulative distribution of the black body spectrum from Min to Max (as given by the global 'RInput' structure). COMMENTS: Global structure 'BlackBodyDistr' must exist. Global variable 'RInput' must be initialized. ------------------------------------------------------------------------ TabulateElVelocityDistr ------------------------------------------------------------------------ PURPOSE: Tabulates electron velocity distributions. DECLARATION: void TabulateElVelocityDistr() INPUT: OUTPUT: ElVelocityDistr - tabulated distribution (global array). DESCRIPTION: Tabulates a set of electron velocity inverse cumulative distributions for a range of photon energies. The distribution is denoted as \ro_{1}' and described in a separate documentation. When those distributions are tabulated, structure 'ElVelocityDistr' is used to draw electron velocity for a given photon energy. COMMENTS: Global structure 'ElVelocityDistr' must exist. Global variables 'S2Photon' and 'RElectron' must be initialized. ------------------------------------------------------------------------ TabulateCrossSection ------------------------------------------------------------------------ PURPOSE: Tabulates the cross section. DECLARATION: void TabulateCrossSection() INPUT: OUTPUT: CrossSectionTab - tabulated cross section (global array). DESCRIPTION: Tabulates the mean cross section for Compton scattering in the photon energy range [SCross.Min, SCross.Max] (SCross is a global structure). The cross section is approximated, and denoted as \Sigma_i in the separate documentation. COMMENTS: Global array 'CrossSeciotnTab' must exist. Global variable 'SCross' must be initialized. ------------------------------------------------------------------------ TabulateSunyaevTitarchukDistr ------------------------------------------------------------------------ PURPOSE: Tabulates Sunyaev-Titarchuk cumulative distribution. DECLARATION: void TabulateSunyaevTitarchukDistr() INPUT: OUTPUT: SunyaevTitarchukDistr - tabulated distribution (global array). DESCRIPTION: Tabulates Sunyaev-Titarchuk cumulative distribution in the Thomson depth range [SS_T.Min, SS_T.Max] (i.e.: from 0 to Thomson thickness of the source). COMMENTS: Global array 'SunyaevTitarchukDistr' must exist. Global variable 'SS_T' must be initialized. ======================================================================== FUNCTIONS defined in 'io.c' file ======================================================================== ------------------------------------------------------------------------ Error ------------------------------------------------------------------------ PURPOSE: The error handling function. DECLARATION: void Error(Msg) char *Msg; INPUT: Msg - a string containing an error message. OUTPUT: Nothing. DESCRIPTION: This instance of the error handling function writes the message to the log file. Program is terminated. ------------------------------------------------------------------------ Warning ------------------------------------------------------------------------ PURPOSE: Writes the warning message to the log file. DECLARATION: void Warning(Msg) char *Msg; INPUT: Msg - a string containing a warning message. OUTPUT: Nothing. DESCRIPTION: Writes the message to the log file. Program is not terminated. ------------------------------------------------------------------------ Message ------------------------------------------------------------------------ PURPOSE: Writes the message to the log file. DECLARATION: void Message(Msg) char *Msg; INPUT: Msg - a string containing a message. OUTPUT: Nothing. DESCRIPTION: Writes the message to the log file. Program is not terminated. ------------------------------------------------------------------------ OpenLogFile ------------------------------------------------------------------------ PURPOSE: Opens the log file. DECLARATION: int OpenLogFile(FileName) char *FileName; INPUT: FileName - name of the file. DESCRIPTION: Opens the log file (global variable 'LogFile'). Returns NULL if error encountered. ------------------------------------------------------------------------ CloseLogFile ------------------------------------------------------------------------ PURPOSE: Closes the log file. DECLARATION: void CloseLogFile(); DESCRIPTION: Well, it's really simple. ------------------------------------------------------------------------ ReadScale ------------------------------------------------------------------------ PURPOSE: Reads the scale from the input stream. DECLARATION: int ReadScale(F, R) FILE *F; struct Scale *R; INPUT: F - input stream, must be open OUTPUT: R - structure containing read scale ReadScale - value passed from the fscanf function. DESCRIPTION: Reads three numbers and one string from the input stream and stores it into 'R' structure. The numbers are: lower and higher limits and number of steps. The text is "LIN" or "LOG". When other text is read the message is written to the console and program is terminated. ------------------------------------------------------------------------ RecognizeDescriptor ------------------------------------------------------------------------ PURPOSE: Finds the text descriptor in a converting table and returns its ID. DECLARATION: int RecognizeDescriptor(s) char *s; INPUT: s - text descriptor OUTPUT: RecognizeDescriptor - descriptor ID DESCRIPTION: Searches in the converting table for the text 's'. When found, returns the ID of the descriptor (from the table). When the text is not found, error message is written to the console and program is terminated. COMMENTS: Function 'Descriptor' performs action opposite to this one. ------------------------------------------------------------------------ Descriptor ------------------------------------------------------------------------ PURPOSE: Finds the descriptor identifier in a converting table and returns the text of the descriptor. DECLARATION: char *Descriptor(ID) int ID; INPUT: ID - descriptor identifier OUTPUT: Descriptor - text descriptor DESCRIPTION: Searches for the identifier 'ID' in the converting table. When found, returns the text of the descriptor. When the identifier is not found, returns "UNKNOWN" string. ------------------------------------------------------------------------ ReadData ------------------------------------------------------------------------ PURPOSE: Reads the input data. DECLARATION: int ReadData(FileName, Data) char *FileName; struct InputData *Data; INPUT: FileName - name of the input file OUTPUT: Data - structure containing all the input data DESCRIPTION: Format of the input file is given above. 'ReadData' reads all the data from the input file and returns it in the structure 'Data'. When error encountered while reading, program is terminated and the message written to the console. 'Error' function cannot be called, since the log file is not known. The input file is divided into lines. Each line contains the text descriptor, '=' sign and the value(s) assigned to the descriptor, like in the example: SpectrumEnergyRange = 0.001 25000 300 LIN The function 'ReadData' reads the descriptor, parses it with the function "RecognizeDescriptor', and reads the data on the right-hand side of the '=' sign, storing it in the appropriate place in 'Data' structure. ------------------------------------------------------------------------ CopyData ------------------------------------------------------------------------ PURPOSE: Copies the input data to the log file. DECLARATION: void CopyData(Data) struct InputData Data; INPUT: Data - structure containing the input data OUTPUT: Output is written to the log file; some global variables are initialized. DESCRIPTION: The contents of 'Data' structure is decoded and written to the log file (this is good for verification). ------------------------------------------------------------------------ WriteInitialSPectrum ------------------------------------------------------------------------ PURPOSE: Writes the initial spectrum to the file. DECLARATION: void WriteInitialSpectrum(Distr, FileName) int Distr; char *FileName; INPUT: Distr - distribution of the spectrum (BLACK_BODY or POWER_LAW) FileName - name of the output file DESCRIPTION: Writes the initial spectrum to the file given by the 'FileName'. ------------------------------------------------------------------------ WriteData ------------------------------------------------------------------------ PURPOSE: Writes the output data. DECLARATION: void WriteData(TransFileName, ReflFileName) char *TransFileName, *ReflFileName; INPUT: TransFileName - name of the output file with transmitted spectrum ReflFileName - name of the output file with reflected spectrum DESCRIPTION: Writes the resulting spectra to the files given by names. ======================================================================== FUNCTIONS defined in 'auxil.c' file ======================================================================== ------------------------------------------------------------------------ Integral ------------------------------------------------------------------------ PURPOSE: Calculates the integral of the given function. DECLARATION: double Integral(f, a, b, eps) double (*f)(); double a, b; double eps; INPUT: f - a function to be integrated over a - left limit of integration b - right limit of integration eps - fractional accuracy of calculation OUTPUT: Integral - integral of 'f' from 'a' to 'b'. DESCRIPTION: Integrates the function 'f' using Simpson's rule. The additional define: #define MAX_ITER 20 before function declaration is required. MAX_ITER is the maximum number of iterations (2^(MAX_ITER-1) is the maximum number of steps). ------------------------------------------------------------------------ Idx ------------------------------------------------------------------------ PURPOSE: Index and remainder for a given continuous value and a scale. DECLARATION: int Idx(Val, S, R) double Val; struct Scale S; double *R; INPUT: Val - continuous value S - a converting scale OUTPUT: Idx - discrete index for a value 'Val' and scale 'S' R - the remainder DESCRIPTION: The scale converts between continuous values and discrete indices of an array. This function returns an index for a value 'Val' and a given scale 'S'. The returned index plus the remainder targets exactly the value of 'Val'. This means that: Val(Idx) + R * Delta = Val COMMENTS: When the value 'Val' is outside range [S.Min, S.Max], 'Error' function is called, and program is terminated. ------------------------------------------------------------------------ Val ------------------------------------------------------------------------ PURPOSE: Continuous value for a given index and a scale. DECLARATION: double Val(Idx, S) int Idx; struct Scale S; INPUT: Idx - discrete index S - a converting scale OUTPUT: Val - continuous value for an index 'Idx' and scale 'S' DESCRIPTION: The scale converts between continuous values and discrete indices of an array. This function returns a value for an index 'Idx' and a given scale 'S'. COMMENTS: When the index 'Idx' is outside range [0, S.N-1], 'Error' function is called, and program is terminated. ------------------------------------------------------------------------ Draw ------------------------------------------------------------------------ PURPOSE: Draw a random number according to a given distribution DECLARATION: double Draw(D) struct Distribution *D; INPUT: D - tabulated inverse cumulative distribution OUTPUT: Draw - generated random value DESCRIPTION: Structure 'D' contains tabulated inverse distribution. The drawn value is generated according to it. COMMENTS: This function generates a uniformly distributed random number, and calls the 'DrawRnd' function. ------------------------------------------------------------------------ DrawRnd ------------------------------------------------------------------------ PURPOSE: Produces a random number according to a given distribution from a given random number with the homogeneous distribution over [0, 1] DECLARATION: double DrawRnd(D, Rand) struct Distribution *D; double Rand; INPUT: D - tabulated inverse cumulative distribution Rand - a random number, uniformly distributed over [0, 1] range OUTPUT: DrawRnd - generated random value DESCRIPTION: Structure 'D' contains tabulated inverse distribution. The drawn value is generated according to it. COMMENTS: Apart from some special cases one should not call this function, but rather the 'Draw' function. ------------------------------------------------------------------------ TabulateInverseDistribution ------------------------------------------------------------------------ PURPOSE: Tabulates inversion of the cumulative probability distribution DECLARATION: void TabulateInverseDistribution(f, Rng, D) double (*f)(); struct Range Rng; struct Distribution *D; INPUT: f - probability density distribution function Rng - range over which distribution is calculated OUTPUT: D - tabulated distribution DESCRIPTION: The inversion of the cumulative probability distribution is tabulated in the structure 'D'. COMMENTS: ------------------------------------------------------------------------ NormalizeSpectra ------------------------------------------------------------------------ PURPOSE: Normalizes the output spectra before printing out. DECLARATION: void NormalizeSpectra(NumPhoton) int NumPhoton; INPUT: NumPhoton - number of photons generated during simulation OUTPUT: TransmittedSpectrum, ReflectedSpectrum - global arrays DESCRIPTION: Spectra in global arrays 'TransmittedSpectrum' and 'ReflectedSpectrum' are normalized. COMMENTS: ======================================================================== FUNCTIONS defined in 'func.c' file ======================================================================== ------------------------------------------------------------------------ F ------------------------------------------------------------------------ PURPOSE: Returns the function F. DECLARATION: double F(x) double x; INPUT: x = h nu / (m c^2) - dimensionless photon energy OUTPUT: F - function estimated. DESCRIPTION: Function F estimates the integral over cross section: 1 / (2 Pi r0^2) Integral[y sigma[y], {y, 0, x}] where sigma - Compton cross section. More detailed description is given in the separate documentation. ------------------------------------------------------------------------ InvF ------------------------------------------------------------------------ PURPOSE: Returns the inverse function F. DECLARATION: double InvF(F) double F; INPUT: F - value of the function to be inverted OUTPUT: InvF - estimated F^{-1}. DESCRIPTION: Function F estimates the integral over cross section: 1 / (2 Pi r0^2) Integral[y sigma[y], {y, 0, x}] where sigma - Compton cross section. This function estimated the inverse of the function F. More detailed description is given in the separate documentation. ------------------------------------------------------------------------ ElectronDistr ------------------------------------------------------------------------ PURPOSE: Returns the energetic distribution of electrons. DECLARATION: double ElectronDistr(gamma) double gamma; INPUT: gamma - electron Lorentz factor (=1 / Sqrt(1 - (v/c)^2) OUTPUT: ElectronDistr - distribution of electrons. DESCRIPTION: Returns the velocity distribution of electrons (in gamma Lorentz factors). This function is actually used by program. COMMENTS: This function just makes a call to the appropriate distribution function (e.g. Maxwellian or power-law). It also passes additional arguments to the distribution, like electron temperature or electron energy index. This function must have only one argument, since it is integrated using function 'Integrate' ------------------------------------------------------------------------ PhotonDistr ------------------------------------------------------------------------ PURPOSE: Returns the number distribution of initial photons. DECLARATION: double PhotonDistr(x) double x; INPUT: x = (h nu) / (m c^2) - photon energy in m c^2 units, OUTPUT: PhotonDistr - number distribution of photons. DESCRIPTION: Returns the number distribution of photons. This function is actually used by program. COMMENTS: This function returns the number distribution of initial black body spectrum. This means that the spectrum described by this function is a photon spectrum, i.e. number of photons per keV. The original black body spectrum, as defined by 'BlackBodySpectrum' function is an energetic distribution, and therefore it should be divided by energy in order to get the required photon spectrum. So, the function 'PhotonDistr' makes a call of 'BlackBodySpectrum' and divides the result by energy. Another reason of making this function is the fact that it must have only one parameter, because it is to be integrated by the function 'Integrate'. Original 'BlackBodySpectrum' function has two parameters and cannot be directly integrated. ------------------------------------------------------------------------ Ro1Prime ------------------------------------------------------------------------ PURPOSE: Probability distribution for the velocity of electron colliding with the photon. DECLARATION: double Ro1Prime(gamma) double gamma; INPUT: gamma - electron Lorentz factor (=1 / Sqrt(1 - (v/c)^2) [PhotonX - photon energy (= (h nu) / (m c^2), global variable] OUTPUT: Ro1Prime - velocity probability distribution DESCRIPTION: This is the distribution \ro_1', as denoted in the separate documentation. This distribution is NOT normalized. COMMENTS: In fact, this is two-dimensional distribution, and this function has TWO input parameters. However, for technical reasons, the second parameter is located as a global variable in this file. This is because the function is to be integrated by "Integrate" function. That's nasty, but this is not an object oriented code... ------------------------------------------------------------------------ SunyaevTitarchuk ------------------------------------------------------------------------ PURPOSE: Returns sin-like Sunyaev-Titarchuk probability density distribution. DECLARATION: double SunyaevTitarchuk(x) double x; INPUT: x - variable [0..Pi] OUTPUT: SunyaevTitarchuk - S-T distribution, normalized to [0, 1] integral DESCRIPTION: Function returns Sunyaev-Titarchuk distribution f(x) = sin x. COMMENT: This function is used by 'TabulateSunyaevTitarchukDist' function in order to tabulate this distribution. ------------------------------------------------------------------------ MaxwellDistr ------------------------------------------------------------------------ PURPOSE: Gives the electron Maxwellian distribution. DECLARATION: double MaxwellDistr(gamma, Te) double gamma; double Te; INPUT: gamma - electron Lorentz factor (=1 / Sqrt(1 - (v/c)^2) Te = (k T) / (m c^2) - dimensionless electron temperature. OUTPUT: MaxwellDistr - Maxwellian distribution for a given 'gamma' and 'Te'. DESCRIPTION: Calculates the relativistic Maxwellian energy distribution. Due to changing variables to gamma Lorentz factor, the factor gamma * sqrt(gamma^2 - 1) appeared: p^2 dp = gamma sqrt(gamma^2 - 1) dgamma ------------------------------------------------------------------------ PowerLawDistr ------------------------------------------------------------------------ PURPOSE: Gives the electron power-law distribution. DECLARATION: double PowerLawDistr(gamma, s) double gamma; double s; INPUT: gamma - electron Lorentz factor (=1 / Sqrt(1 - (v/c)^2) s = index of the power-law distribution OUTPUT: PowerLawDistr - power-law distribution of electrons for the given 'gamma' and 's' DESCRIPTION: Calculates the number density of electrons in the range (gamma, gamma + dgamma) in the power-law distribution: N(gamma) = const gamma^(-s) ------------------------------------------------------------------------ BlackBodySpectrum ------------------------------------------------------------------------ PURPOSE: Gives the black body spectrum. DECLARATION: double BlackBodySpectrum(x, Theta) double x, Theta; INPUT: x = (h nu) / (m c^2) - photon energy in m c^2 units, theta = (k T) / (m c^2) - dimensionless photon temperature (T - radiation temperature). OUTPUT: BlackBodySpectrum - the spectrum for a given energy and temperature. DESCRIPTION: Function returns the black body spectrum density distribution. The spectrum is normalized in the range [0, Infinity). COMMENTS: Used by 'PhotonDistr' function. ======================================================================== FUNCTIONS defined in 'compton.c' file ======================================================================== ------------------------------------------------------------------------ CrossSection ------------------------------------------------------------------------ PURPOSE: Returns the approximated mean cross section. DECLARTION: double CrossSection(x) double x; INPUT: x = h nu / (m c^2) - dimensionless photon energy. OUTPUT: CrossSection - mean cross section. DESCRIPTION: The mean cross section for different photon energies is tabulated in the global array 'CrossSectionTab'. This function interpolates linearly in the array. COMMENTS: The global array 'CrossSectionTab' must be prepared by then by the function 'TabulateCrossSection'. Global variable 'SCross' must be initialized. ------------------------------------------------------------------------ GenerateInitialPhoton ------------------------------------------------------------------------ PURPOSE: Generates the initial photon according to appropriate energy and spatial distributions. DECLARATION: void GenerateInitialPhoton(S, Ph) struct SourceDescription S; struct Photon *Ph; INPUT: S - description of the source geometry and initial distribution. OUTPUT: Ph - generated photon. DESCRIPTION: This function generates initial photon energy 'x', direction 'eta', location 'tau' and statistical weight 'w'. According to the source description given by 'S', the photon parameters are generated for black body or power law distributions, for slab or sphere geometry, and for different initial space and directional distributions. Energy for black body distribution is generated by inverting the cumulative distribution tabulated in 'BlackBodyDistr' array, statistical weight is set to 1. Power law is a special case (...to be explained later...). For the slab geometry we have four initial configurations: homogeneous distribution over the slab; slab irradiated at a given angle 'eta0;; slab irradiated by the optically thin source; slab irradiated by the optically thick source. Details are discussed in a separate document. COMMENTS: ------------------------------------------------------------------------ NextPosition ------------------------------------------------------------------------ PURPOSE: Calculates the position of the next scattering and probability of escape. DECLARATION: void NextPosition(S, Ph) struct SourceDescription S; struct Photon *Ph; INPUT: S - description of the source geometry and initial distribution. Ph - input photon. OUTPUT: Ph - photon with the new position. DESCRIPTION: For a given photon this function determines the distance to a cloud edge, according to geometry given by 'S'. The probability of scattering (and the probability of escape) is calculated along the path of the photon. The statistical weight and the spectrum of emerging photons (global arrays 'TransmittedSpectrum' and 'ReflectedSpectrum') is updated. Then the photon free path is generated according to the appropriate probability distribution. Hence the position of the next scattering can be determined. The transmitted spectrum is updated for a spherical geometry. When the slab is considered, the transmitted spectrum is updated for photons emerging from the side of the slab opposite to the initial source, and the reflected spectrum gathers those photons that leave the slab on the side of the source. The details are described in the separate documentation. ------------------------------------------------------------------------ ElectronVelocity ------------------------------------------------------------------------ PURPOSE: Generates the electron velocity for a given photon energy. DECLARATION: void ElectronVelocity(x, El) double x; struct Electron *El; INPUT: x = h nu / (m c^2) - dimensionless photon energy OUTPUT: El - generated electron velocity DESCRIPTION: The electron velocity is determined with respect to the probability distribution \ro_1(v) (as denoted in the separate documentation). The velocity is given in the frame in which photon goes along the 'z' axis. First the angle 'phi' is generated with the uniform distribution. Then, velocity 'beta' is found by inverting the two- dimensional distribution \ro_1'(v), tabulated in ElVelocityDistr. Finally, the photon energy 'x' (in the electron rest frame) is generated with the distribution D_1"(x|v), and the electron cosine theta 'mu' is found, by transforming to the plasma frame. More details can be found in the separate documentation. ------------------------------------------------------------------------ Scattering ------------------------------------------------------------------------ PURPOSE: Determine the direction and the energy of the photon after scattering. DECLARATION: void Scattering(El, Ph); struct Electron El; struct Photon *Ph; INPUT: El - electron velocity Ph - photon before scattering OUTPUT: Ph - photon after scattering DESCRIPTION: The direction of the photon after scattering is generated with the probability distribution \ro_2(mu, phi) (as denoted in the separate documentation). During this procedure we change the frame several times. In the original frame (plasma frame) the incident photon moves in the x-z plane (we can rotate it around z axis freely thanks to symmetry). Than we rotate the original frame around 'y' axis by the angle Theta, hence the photon is moving along 'z' axis. In this frame the electron direction (theta, phi) is generated by procedure 'ElectronVelocity'. After this we rotate the frame around 'z' axis by angle phi and around 'y' axis by angle theta, so the electron moves along 'z' axis. This is our frame for scattering (in fact, the scattering is calculated in the electron rest frame and than moved to this frame) in which we find the direction (theta', phi') and the energy x' of the scattered photon. Than we move back to the original frame in order to find the angle Theta' (cos(Theta') = eta') of the scattered photon and its energy. We determine theta' and phi' by the rejection method. More detailed description along with pictures is given in the separate documentation.