Chapter 6
Random Number Generation

Monte Carlo method: A method of using random sampling of numbers in order to estimate the solution to a numerical problem (D. Thompson, The Concise Oxford English Dictionary ninth edition, 1995, Clarendon Press)

Random numbers are used within the program as a source of decision making. The generation of these numbers is done using a computer code written by Dr. Jos Thijssen. A good random number generator must reasonably represent a known probability distribution function (usually uniform over some finite domain). There should be no built in trends or periodicites, and there should be no correlation in anyway between previously generated values. Originally, a “NAG library” function was used to generate the random numbers. However, it was found that this function was responsible for making the program operate incorrectly. It is believed that the function incorporated periodicity into its random number generation. This Monte-Carlo method of radiative transfer follows individual photons as they scatter through the dust until they escape or are absorbed.

The following derivation will aid in the understanding of the equations derived later in this section:

If is the probability of being in the range and is continuous in the interval . When we pick a random number the probability of its being in is just , e.g. .

We want a function such that putting gives us values distributed according to in .

If then:

probability of getting between

The first random number chosen is to help determine a starting direction, in the co-ordinate system of the emerging photon, where is the component of the photon. Random number is used to select the initial azimuthal direction of travel , the angle at which a photon departs the star in the plane. & are calculated thus:

The component of the emerging photon must vary between the positive and negative hemispheres, the North and South Poles of the star, to ensure a truly random distribution of photons. The simple equation used to convert the random number produced by the computer, which will be in the range , is:

Numbers between will produce a negative , while numbers between will produce a positive .

We can now say that, , and calculate it thus:

The second random number is used to calculate the mean free path of the photon, , how far it will travel before it interacts with a dust grain. is calculated by means of the equation:

since the random numbers are uniformly distributed between 0 and 1 we can change this, without loss of generality, to:

The third random number needed is to decide whether a photon is scattered or absorbed by an interaction with a dust particle. The decision is a simple one, if the albedo of the dust grain is greater than then the photon is deemed to have been scattered, if not the photon is said to have been absorbed.

Assuming the photon is scattered, then it requires a new direction to travel in. This new direction is chosen by two random numbers, & . is the azimuthal angle through which the photon is scattered. This can now be calculated in exactly the same way as that of, . Therefore , becomes:

The fifth random number, , is the angle declination through which the photon is scattered, and is calculated using the scattering phase function , we can say that .

using the subsitution

Rearranging we get:

6.1 The lifetime of a photon: birth to circumstellar disc escape

Once a photon has been released from the central star it travels out through the void between the star and the circumstellar disc, on entry to the disc its passage is impeded by the dust comprising the disc. The interactions that take place attenuate the photon. The large majority of light which finally emerges from the disc is, therefore not in the condition it was, when discharged from the star.

---

On release from the star, the photon travels in a direction arbitrarily chosen by random numbers & as explained above:

From the diagram above and trigonometry we can see that the initial & directions are

From these calculations we get the direction that the photon will initially travel. The photon travels in this direction out through the void between the star and the dust disc, un-attenuated by physical interaction. The void exists due to the extreme temperatures within this region, which vaporises all matter.

The photon is advanced, from the inner edge of the dust cloud, through the dust disc by means of an adaptive integration step. The distance that the photon will move by is determined by means of the random number . This distance is measured in mean-free-paths , however accurate integration requires that is much less than a mean free path and therefore the photon only actually travels a fraction of a full mean free path. This is because the dust distribution is not uniform throughout one mean free path. Visualising a straight line, from infinity for analytic modelling or , the outer radius of the dust shell, for numerical computations, the closer to the star along this line one is, the denser the dust will be. The distance travelled is therefore calculated thus:

but

Where n is the number-density of particles and is the average cross section presented by a single grain. is the monochromatic opacity coefficient, which is the total absorption cross-section presented by all matter in a unit volume to photons of frequency , and is equivalent to the inverse of a mean-free-path. The larger the opacity of a system the smaller the distance a photon can travel without interaction with a grain. We can therefore say:

Where is calculated in accordance with the step-function that determines the degree of flatness the disc has. The derivation of can be seen in appendix I

The position of the photon is calculated at each step of its travel through the disc. In this way, the photon is ray-traced throughout its travel. The photon’s position and distance moved is, at this first stage calculated by:

Once the position of the photon has been calculated, the next stage of computation is to check if it has travelled a straight-line distance greater than that of the depth of the disc. If this condition is satisfied, the photon has escaped from the circumstellar disc and our interest with it is concluded.

Should it not yet have finished it journey, which is most likely at this stage of events, the remaining segment of the optical depth left to be travel by the photon is calculated thus

This loop of calculating position, travelling a distance, and checking for circumstellar escape, in accordance with the dust distribution, continues until the photon travels the full length of the computed mean free path.

The position of the photon being re-calculated at each stage along the mean-free-path using the modified equations;

The dust particle that the photon must encounter, at the end of its calculated mean-free-path, has its interaction properties chosen by random number . If the random number is greater than the albedo, which is input by the user, then the photon is absorbed and not scattered. Should this be the case then the simulation returns to the star to inject another photon into the circumstellar dust.

---

However should the photon scatter, a new direction through which the photon is to be scattered must be computed, this is done using & . So that the program can be written in a loop structure, the calculation of this new direction must occur in the photon’s co-ordinate system. The working frame is always the frame in which the photon is situated, and the origin of this working frame is the grain the photon interacts with. This means transforming back to the initial co-ordinate system once the calculation has been done in the photon’s co-ordinate system; so that the path the photon takes can be kept track of. The transformation from the photon’s co-ordinate system, the double-prime system, so called because the transformation is done in two stages, can be seen in appendix I

As explained above (EQUATION 49):

but from FIGURE 13:

and therefore in the photons co-ordinate system:

Therefore, the direction of the photon is calculated using the equation:

and as before,

Where it is & that are chosen by the random numbers. The transformation back to the initial frame results in the equations taking the following form:

At this point, the entire loop is reiterated until an escape condition is satisfied; this could be either the photon being absorbed or the photon escaping the boundaries of the circumstellar dust.

Should the photon escape, and it is on a line of sight that will influence the observer’s field of view, it will do so along . The plane of the sky, for this observer is orthogonal to . The observer’s co-ordinate system, on the plane of the sky (his sky) consists of a horizontal co-ordinate which is perpendicular to & plus a vertical co-ordinate which is perpendicular to &

The circumstellar disc is viewed through this new co-ordinate system and the escape points of the photons are translated on to this planer view. The photons & co-ordinates are translated to the plane by use of the equations:

The derivation of which can be seen in appendix I

These values, along with , are then stored in an array for future use by the program.

To generate an isophotal map, a grid of and values must be defined; this is done using the computer code:

where:

• = the outer edge of the dust disc,
• = the side dimension of each grid square,
• = the number of grid squares in the isophotal map,

With the isophotal is square; this helps in visual verification of the map’s symmetry. is multiplied by 1.1 to allow for photons which appear to come from the extreme edge of the circumstellar disc. These photons will have an influence beyond the disc edge, setting the co-ordinate system slightly larger than the disc allows one to see more fully the extent of the photons influence.

The stored values of & must now be examined to see which of them falls within a pre-designated limit, of the observing angle which is entered by the user. If the escaping photon can be seen from the observing angle, the position is saved to an array, if not it is discarded.

This new array of values is then subjected to examination by the & values, in the same way as they were by the search. The designated values of and defining a sub-grid within which each photon contributes to the over all brightness of the grid squares. A weighting function is used to ensure that the photon contributes less light to a point further away from its escaping point, than to an element closer to its escaping point. A grid square can have several photons adding to its overall brightness level, so that the value of one grid square is the culmination of many photons.

where,

• = The co-ordinate of the emerging photon.
• = The co-ordinate of the grid.
• = The co-ordinate of the emerging photon.
• = The angle of the emerging photon.
• = The angle of the observer.
• = The radius of the circle of influence.

The outputted values from these computations are then stored in a file in matrix format:

In this way an isophotal map is merely a colour indication, or schematic representation, of a matrix, a colour intensity graph in a simple 2D co-ordinate system.

6.2 Diagnostics

The computer program contains some very simple diagnostics to ensure that the random number generation is indeed random. This coding takes the following form; the number of times an individual random number is generated is recorded . The sum of each individual set of random numbers is also recorded , as is the sum of the squares . The diagnostic routine then checks the mean of each set of random numbers:

which should be equal to, or very close to . The second diagnostic performed is to check the standard deviation of the random numbers:

This should result in an answer equal to or very close to 2.94. The results of the operations performed with the random numbers were also tested. If the answers to any of the tests were too far away from the standard answer then there was a fault with the program, most likely the random number generator. The program was then scrutinised until this problem was found and solved.

The values currently produced from these diagnostics are of the order:

---

As can be seen from the results of the program (see TABLE 2), the numbers that are produced are random.

There was statistical noise produced within the program which was outputted on the isophotal maps, this noise was calculated and deemed using the following conjecture:

• = Number of pixels.
• = Number of photons.

Therefore, on average the number of photons per pixel and the noise . Therefore increasing the number of pixels decreases noise. The program used a pixel grid and photons, therefore:

The average number of photons per pixel and the noise . As a percentage the noise is:

However, this is only an average result for the noise as, by the nature of program, there will be more photons at he centre of the isophotal map, as this is where the star is.

6.3 Virgin Radiation

The will be a certain number of photons that are able to escape the confines of the circumstellar disc without interacting with a dust grain, this radiation is deemed virgin radiation. Since the photons will be coming directly from the star, they will contribute to the central area of the isophotal map. The brighter the central area of the map, the smaller the dynamic range available for the outer regions of the map will be, thus hiding some of the detail of the modelled circumstellar disc. For this reason a number of runs were made where any of the virgin radiation was removed. The results can be seen in appendix II.