A major part of this work was the computer implementation of the simulation. The resulting program units shall now be described here. The program was completely written in IDL (Interactive Data Language), version 5.4. This language provides easy methods and functions to manipulate, query and evaluate large data arrays, which avoid the necessity to program loops. Several useful mathematic functions are implemented, including numerical integration and random number generation. The algorithms used for those features are mostly taken from the famous ”Numerical Recipes“ (Press et al. (1988)). A large number of useful libraries are available on the Internet. NASA provides the astrolib (NASA (2001)), containing lots of procedures covering astrophysical problems. This library also includes routines to read and write the fits-format, being an interface to a lot of data reduction programs. Finally IDL is able to write plots directly to Postscript format. The major flaw of IDL is, that it is a commercial program and thus not easily available to anyone. It is sold by Research Systems.
For programming convenience the program was divided into three parts. cosmology.pro provides for functions related to cosmology, mainly distance calculations. It uses the equations derived in appendix A. The involved integrations are done numerically, using the Romberg method from Press et al. (1988). utility.pro includes the SED model described in section 2.1 and some modified routines for reading table files and plotting data with asymmetrical errors. Finally, sim.pro contains the Monte Carlo simulation itself, which is described in section B.2. It must be said that this part is affected most by the development during the work, undergoing several major changes. So, the procedures found there are laden with historical ballast and would benefit from a general cleanup and rearrangement. This may already have happened while you read this lines.
In the following sections the routines of the program units are described. The syntax and usage is listed, and also known issues are discussed. The routines are available on request for interested users. Routines written as Name (Argument) are functions, those written as Name, Parameter are procedures. Not all functions and procedures contained in the units are described, as some routines are only for the internal computation or may be unfinished or have unfixed errors putting them beyond usefulness.
The final modeling is done through Monte Carlo simulations. This method builds a realization of a distribution using random numbers. The distribution is described by a function p(x). The function is then normalized to be 1 at its maximum:
|
| (B.1) |
To x a random value chosen from all allowed values for x is assigned. A second random value is generated from an uniform distribution of all numbers between 0 and 1. If this value is less or equal than P(x), the distribution is regarded as ’hit’ and the value of x is made part of the realization of p(x). This procedure is repeated until the desired number of values is generated. One possible realization of p(x) is found.
The random numbers are generated by the IDL built-in random number generator. It is initialized with the system time of the computer, measured in seconds. To generate uniformly distributed random numbers, IDL uses a method from ”Random Number Generators: Good Ones are Hard to Find“, Park and Miller, Communications of the ACM, Oct 1988, Vol 31, No. 10, p. 1192, which is modified by a Bays-Durham shuffle, resulting in a random generator similar to ran1() in Press et al. (1988), section 7.1. The formula used is:
 | (B.2) |
where i1 and i2 are integer random deviates in the range [1...imax] and imax = 231 - 2 is the largest possible integer random deviate. yu is in the range 0 < yu < 1. For normally distributed (Gaussian) random numbers the Box-Muller method is used. For details see Press et al. (1988), section 7.2. The generated numbers have a mean of zero and a standard deviation of one. Such a random number in the following is designated by yn.
For a given solid angle of sky a range in redshift to be populated with starburst
galaxies is defined, since the simulation out to z = 
(or out to z 
1000, where the
CMB radiation indicates a hot dense state of the universe without galaxies) would take to
much computing time. The assigned range is divided into equal sized redshift intervals
(redshift bins). For each interval the number of galaxies in it is calculated. This is done by
integrating the luminosity function to a luminosity cutoff Lmin and multiplying
the result by the comoving volume corresponding to the redshift bin and solid
angle.
To assign luminosities to the galaxies of each redshift bin, a Monte Carlo realization of the luminosity distribution is conducted until the required number of galaxies is found. To minimize the number of drawn blanks, the drawing is done in logarithmic space. Further, an upper luminosity limit Lmax is set to reduce the number of drawn blanks. To each galaxy a SED is assigned. This is done by Monte Carlo realizations of the parameter distributions given in sections 2.1.3 and 2.2.2.
The position assignment for the generation of artificial maps is done by a Monte Carlo realization of the spatial galaxy distribution, derived from a biased dark matter distribution. The latter is computed from the numerical simulations of the Hubble Volume project.
This procedure initializes and provides the common block Cosmology:
COMMON
Cosmology,
OmegaM,
OmegaK,
OmegaL,
Omega0,
H0,
h,
c,
k,
Dh,
tH
where OmegaM, OmegaL, OmegaK, H0 and k are 
m, 
, k,H0 and k as defined in
appendix A. c is the speed of light in km s-1. h is h
100, Dh and tH are the Hubble
distance DH and the Hubble time tH, also defined in appendix A. The latter given
in years:
 | (B.3) |
The values of the cosmological parameters are stored in the procedure. To change them, one has to edit this procedure. Later versions may make use of IDL keywords to override the settings in this procedure.
This four functions do what they are named for, returning the various distance measures. The argument of the functions is redshift z. See appendix A for details.
This function returns the age of universe at a given redshift z. As the IDL
integration with open integral from z to does not work well with the needed
integrand, the integration was changed from z to upper limit 10000. This inflicts a
neglectable error. However, this function should only be used out to a redshift of
about a thousand.
This function returns the lookback time according to redshift z.
This function returns the comoving volume, all sky out to redshift z, using eq. (A.33).
This function hands over the comoving volume corresponding to an area of sky (given in steradian), enclosed between redshift z and z + dz. It uses AllSkyComovingVolume (z) and the method described in section A.2.5.
This function contains the comoving volume element defined by eq. (A.31).
The comoving volume is calculated here between solid angle (given in radians)
![[[h, h + dh],[f,f + df]]](DiplWeb138x.png)
and redshift interval [z,z + dz]. The angles are defined as
usual in polar coordinates. The final integration then is:
 | (B.4) |
which is evaluated by numerical integration using IDL’s INT3D function.
This function returns the absorption coefficient of dust at (a) given
wavelength(s) as described in section 2.1.2. Beta is 
as defined there.
This is the greybody making use of Kappa(). Distance and z seem redundant. They are both to be given to speed up computing. In the simulation of this work the luminosity distances are evaluated at start and so the integration has not to be rerun here. The function returns the flux density in Jy. If Distance and z are set to zero, the function returns a differential luminosity in erg s-1 Hz-1.
This functions interpolate the SED given in Wavelength and Luminosity to return a differential luminosity at wavelength AtWavelength. The unit of the returned value is the same as of Luminosity. Spectrum interpolates linearly, SpectrumLog logarithmically. Again, ValueCount may be obsolete in later versions.
This procedure reads a file into the arrays Wavelength, FluxDensity and
Errors. The first column in the file is assumed to be wavelength in 
m, the
second flux density in Jy and the third is the error of flux density in Jy. If the
file does not provide a third column, the errors are set to zero. ValueCount
contains the number of data points and may be obsolete in later versions.
This function evaluates the Lorentz function as defined by eq. (2.36).
This function returns the PAH line features described by Lorentz functions.
This function returns the flux density generated from the SED model as described in section 2.1. For a given set of parameters the differential luminosity in erg s-1 Hz-1 is returned for a set of (restframe) wavelengths.
This function returns the bolometric luminosity of the model SED, described
by a given parameter set, in erg s-1. The bolometric luminosity is calculated
between 10 and 3.105 m. Attention: This function at the moment only works
properly for = 2!
This function returns a set of parameters describing a model SED. This set of parameters results in a SED having L60 = L60. Seed is used by the random number generator. If subsequent calls shall be made or BuildSED is part of a routine using random numbers, this variable should always be passed between all involved routines in order to get non repeating random number series.
This procedure should be used to start any procedure in sim.pro. It loads the SED of an object at redshift z , stored in SpectrumFile into memory, converts it into differential luminosities and stores the data in the common block SepctralData. Another common block initialized in this procedure is PlotData, containing the variables PlotColor and SavePS. The latter is a flag. If it is set to one, all plots are saved to Postscript files. Zero indicates output on screen. Some procedures in this unit overrides this flag on request. PlotColor is the color used for b/w-plots. It is white on screen, and black on Postscript. The keyword AUTO is used to run the simulation without interactive parts. If the keyword is set, SavePS is forced to be one and a flag is handed to the simulation in order to suppress input requests. This mode was implemented to run the simulation on the linux cluster. At the end of Main the procedure(s) to be executed must be given, in order to branch to them. This routine is one of those that badly needs to be streamlined.
This function returns the luminosity function used in the simulation. The z-dependence is achieved through a common block, as it is necessary to integrate the luminosity function numerically. The routines doing this in IDL only accept functions depending on one parameter for a single dimension integration.
This is the Monte Carlo simulation. AutoMode is the flag handed over by Main, suppressing interactive mode. As the simulation depends on a lot of parameters, they are not handed over to this procedure but written directly into the procedure. That part of the procedure is emphasized as parameter block. It contains the following parameters and flags:
This is the startup flag. Is it set to one, the parameters are loaded from Parameters.dat and instead of building the universe by Monte Carlo, it is loaded from Universe.fits. be sure those files belong to each other! In the future, the Monte Carlo simulation and the following plots and consistency checks will be divided into two procedures, making this flag obsolete. If this flag is set to one, the following parameters are ignored.
Those variables define at which redshift the simulation starts, where it ends and how broad the bins are.
ViewArea is the solid angle the observation looks into, given in sr. AreaText is a string that may contain a comprehensive version of the solid angle, e.g. ”1 deg2“ or ”100 arcmin“. This string is used in plots. It may contain IDL font and format codes as e.g. ”1 deg!U2!N“.
These two arrays must be of the same size. ObservingFrequency contains the frequencies the simulation shall take a close look to, e.g. by evaluating the number counts, the total flux or by generating a map, To be given in GHz. FluxCutoff defines the sensitivity (in Jy) of the imaginary detector at each frequency.
This variable contains the number of flux bins the number counts shall be calculated in.
These variables define the faintest and the brightest objects that shall be placed in the simulated part of the universe. They have to be defined in erg s-1.
This flag indicates whether the lower luminosity cutoff may be subjected to the k correction and thus vary with redshift. This method is hardly tested, handle with care!
This flag tells if to use GHz in plot texts (set to one), or to use m (set
to two).
This parameter tells how to convert the 60 m luminosity function into
number counts at other wavelengths. The following values and options may be
used:
Three methods to derive the object positions in the volume are possible, set by this variable:
If this flag is set to one, a map for each observing frequency is generated. Zero suppresses this, which may save a vast amount of memory.
This flag is only used in auto-mode. If set to one, the universe is saved to Parameters.dat and Universe.fits . If it is set to zero, the simulation is discarded after doing some summary output. This may save lots of disk space. In interactive mode this setting is ignored and the user is asked if the universe shall be saved (your chance to be a hero!).
It is only necessary to define this array if DoMap is set to one. This array must be of the same size as ObservingFrequency and gives the beam size in arc seconds for each frequency.
This parameter states how many pixels per dimension shall be used for one beam size. The number of pixels per map NMap is then:
After the parameter block the simulation follows, and after that the results may be saved, some consistency checks are done and various plots can be viewed. The latter part will eventually migrated into a single procedure.
This function returns the extragalactic far infrared background as described by Fixsen et al. (1998). The function returns the intensity in erg s-1 m-2 sr-1. The function is only valid between ~ 100 and 3000 GHz! A warning message is displayed if this range is exceeded.