StarPopulations

Stellar populations

Above, a view from the middle of 31 000 stars generated at random (in Pov-Ray) and confined to a narrow disk. Note
the similarity to the Milky Way as seen in the sky from Earth - the stars aren't simply confined to a band far away, as it
seems, in fact on average the stars are of the same density near by as they are far away, but as you are looking from
within the disk the stars appear to get denser further away as you look through the disk. This article looks at how we
can model stellar populations using physical principles. This is a useful exercise to facilitate one's understanding of
astrophysics, by drawing together many different results and ideas, and also useful for computer games in which we
might want a realistic region of space to explore!

Stars can usually be grouped into star clusters and also into populations of stars of a similar type. The spectral
classification of stars is discussed under the section on Main Sequence stars. Basically, more massive stars are hotter
and hotter stars are bluer than cooler stars which are redder. The colour sequence in order of increasing temperature
is: red, orange, yellow, white and blue. The very hottest stars emit most of their light in the ultraviolet or at even
shorter wavelengths. (Recall that for a particle of light, or photon, those with a shorter wavelength or higher frequency
have a higher energy and are perceived as bluer - hence a hotter object emits more blue photons). As more massive
Main Sequence (MS) stars (stars that are mature but not old) are hotter and thus bluer, we speak of white and blue
giants, but red, orange and yellow dwarfs. The Sun, Sol, is a yellow dwarf. (Red giants are cooler stars that have left
the Main Sequence and swollen in size and are discussed in the section on supergiants also white dwarfs are old
stellar remnants and not MS stars).

Now, giant stars may have more fuel, but they burn it so much faster and are consequently relatively short-lived. A
large blue giant has a MS lifespan of only about 2 million years, whereas a small red dwarf may live on the MS for
some 6 000 000 000 000 (6 million million years) which greatly exceeds the current age of the Universe (about 14 000
000 000 or 14 thousand million years, let us say 14 billion by taking one billion to be one thousand million). Thus,
what we find is that blue giants near are relatively young stars whereas a red dwarf could be new or it could be one of
the first stars that formed in the Universe!).

Two more views from within our model disk of stars - occasionally a star appears
quite far above or below the plane of the disk, but most are close to it. This is a
model of a galactic disk, though our model is not to scale!

#include "rand.inc"

camera { location <0,0,-60-10> look_at <0,0,0> }
background {color White*0.0}

//Red dwarfs
#declare stars = 10000;
#declare star = 1;
#while(star <= stars)
sphere
{ 0, 3 hollow
pigment { rgbt 1 }
interior
{ media { emission 0.5*RRand(0.1,2,1)
density
{ spherical density_map
   {
     [ 0.0 rgbt <0.0, 0, 0, 1> ]
     [ 0.25 rgb <.5, 0, 0> ]
     [ 0.5 rgb <1, 0, 0> ]
     [ 0.75 rgb <.5, 0, 0> ]
     [ 1.0 rgb <0, 0, 0> ]
   }
   turbulence 0.001
   frequency 0.5
}
}
}
scale 1*5
translate <1,1,1>*VRand_In_Box( <-10000,-100,-10000>, <10000,100,10000>, 1)
}
#declare star = star + 1;
#end

//Blue giants
#declare stars = 10000;
#declare star = 1;
#while(star <= stars)
sphere
{ 0, 3 hollow
pigment { rgbt 1 }
interior
{ media { emission 0.5*RRand(0.5,1.5,1)
density
{ spherical density_map
   {
     [ 0 rgbt <0, 0, 0, 1> ]
     [ 0.25 rgb <0, 0.4, 0.8> ]
     [ 0.5 rgb <0.2, 0.6, 1> ]
     [ 0.75 rgb <0.6, 0.8, 1> ]
     [ 1.0 rgb <0.6, 1, 1> ]
   }
   turbulence 0.001
   frequency 0.5
}
}
}
scale 1*10
translate <1,1,1>*VRand_In_Box( <-10000,-100,-20000>, <10000,100,20000>, 1)
}
#declare star = star + 1;
#end

//White giants
#declare stars = 1000;
#declare star = 1;
#while(star <= stars)
sphere
{ 0, 3 hollow
pigment { rgbt 1 }
interior
{ media { emission 0.5*RRand(0.1,2,1)
density
{ spherical density_map
   {
     [ 0 rgbt <0, 0, 0, 1> ]
     [ 0.25 rgb <0.8, 0.8, 0.8> ]
     [ 0.5 rgb <1, 1, 1> ]
     [ 0.75 rgb <0.6, 0.8, 1> ]
     [ 1.0 rgb <0.6, 1, 1> ]
   }
   turbulence 0.00
   frequency 0.5
}
}
}
scale 1*10
translate <1,1,1>*VRand_In_Box( <-10000,-100,-10000>, <10000,100,20000>, 1)
}
#declare star = star + 1;
#end

//Yellow dwarfs
#declare stars = 5000;
#declare star = 1;
#while(star <= stars)
sphere
{ 0, 3 hollow
pigment { rgbt 1 }
interior
{ media { emission 0.5*RRand(0.1,1,2)
density
{ spherical density_map
   {
       [ 0.0 rgbt <0, 0, 0, 1> ]
       [ 0.3 rgb <1, 1, 0.4> ]
       [ 0.5 rgb <1, 0.5, 0> ] //yellow
       [ 0.75 rgb <.5, 0, 0> ]
       [ 1.0 rgb <0, 0, 0> ]
   }
   turbulence 0.001
   frequency 0.5
}
}
}
scale 1*5
translate <1,1,1>*VRand_In_Box( <-10000,-100,-10000>, <10000,100,20000>, 1)
}
#declare star = star + 1;
#end

//Orange dwarfs
#declare stars = 5000;
#declare star = 1;
#while(star <= stars)
sphere
{ 0, 3 hollow
pigment { rgbt 1 }
interior
{ media { emission 0.5*RRand(0.1,2,1)
density
{ spherical density_map
   {
       [ 0.0 rgbt <0, 0.0, 0, 1> ]
       [ 0.3 rgb <1, 0.2, 0> ]
       [ 0.5 rgb <1, 0.5, 0> ] //yellow
       [ 0.75 rgb <.5, 0, 0> ]
       [ 1.0 rgb <0, 0, 0> ]
   }
   turbulence 0.001
   frequency 0.5
}
}
}
scale 1*5
translate <1,1,1>*VRand_In_Box( <-10000,-100,-10000>, <10000,100,20000>, 1)
}
#declare star = star + 1;
#end

Left: the Pov-Ray code for our population of stars. The numbers of
each star type can be set as desired and the computer then
generates random coordinates for each one; generating a total of
31 000 stars. Each star is also given a randomly varying luminosity
(emission). This is all very good and useful for artistic effect, but
what if we want a statistically realistic star population? How many of
each spectral class should there be? How large and how luminous
should each star be? How many should be single or binary or
multiple stars?

What I have done is use Java to generate a more detailed and more
realistic population of about 30 000 stars (many more can be
generated, but for the purposes of using Excel to plot graphs, Excel
cannot handle more than some 32000 data points, though in house
software could be developed to plot more points). Details of the
program code will be provided soon.

First we decide on the mean stellar density and the volume of space
we would like to populate. Typical values for a galaxy range from
0.003 to 1.5 stars per cubic light-year, with a value of 0.006 for the
region of space near to the Sun. I chose a value of 0.0081 star
systems per cubic light-year and a cubic volume of length 136.92
light-years, simply because this generates the right number of stars
at the given volume (the number of stars being given as the nearest
whole number value of density x volume).

Multiplicity

Star systems can be single, double or multiple. Based on actual
statistics for the Milky Way Galaxy the following Normal statistical
distribution was used:

   multiplicity = (int)(Math.abs(generator.nextGaussian()*2));

In a typical run this gave 31,258 stars grouped into 20,791 separate
star systems with 14,292 of these systems single (45.7% of stars or
68.7% of star systems), 3,778 double (24% of stars or 18.2% of star
systems), 1,787 triple (17.2% of stars or 8.6% of star systems) and
934 (4.5% of star systems) with more than three stars. This is about
right, since about half of observed stars belong to systems with more
than one star.

Mass and Spectral Class

The summary statistics for our stellar population are given in the
extract from the output below:

Total stars = 31258
Number of star systems = 20791
Mean Multiplicity = 1.503438988023664
Single star systems = 14292
Double star systems = 3778
Triple star systems = 1787
Star systems with > 3 stars = 934
Percentage of stars in single systems = 45.72269499008254
Mean stellar mass = 0.9733209783803813
Maximum mass = 19.92247551646834
Minimum mass = 0.010010464049533319

//Spectral class of all stars in their main sequence
Number BD stars = 2192
Number M stars = 11034
Number K stars = 7402
Number G stars = 4231
Number F stars = 3888
Number A stars = 2495
Number B stars = 16
Number O stars = 0

//Stellar remnants
Number RG stars = 4
Number ASG stars = 4
Number WD stars = 4010
Number NS stars = 36
Number BH stars = 0

Where BD are brown dwarf stars. It is difficult to be certain how to
distribute stars among the spectral classes. Brown dwarfs are dim
and so hard to detect and some theorists predict that the majority of
stars may be brown dwarfs, however, the proportions generated
here correspond roughly to know proportions in the Sun's
neighbourhood. M-class red dwarfs are the most numerous and the
giant hot blue O-stars are so rare that only occasionally does one
appear in a sample of 30 000 stars! The relevant code for
generating the star types, based on mass (the initial mass
function, IMF) is:

double mean = -1;
double stdDev = 1.05;
while (mass[i] < 0.01)
{
mass[i] = generator.nextGaussian()*stdDev + mean;
mass[i] = mass[i]*Math.exp(mass[i]*0.35) - Math.exp(-mass[i]*0.1);
}

Initially a gaussian (Normal) distribution is used to generate a spread
of stars that gives us the desired proportions of dwarf stars, whilst
an exponential distribution superimposed upon this gives the
proportions of the smallest and largest stars. The results
approximately agree with the empirical result for stars (other than
brown dwarfs) given by the Salpeter equation:

   IMF ~ M^-2.35

   (^ means 'raised to the power of' or superscript)

This means that for every 5000 G-class stars (like the Sun) there is,
on average, 1 O-star and 25 B-stars, though many observations
give lower counts than this. Our choice adopts more conservative
frequencies of UV-stars (O and B-stars are collectively called
UV-stars since most of their 'light' is emitted as ultraviolet because
they are so hot!). The function could be tweaked if more UV-stars
were desired and, in any case, the number generated is quite
variable with each run. The spectral class can be determined by
setting a minimum mass (taken from empirical data) for each
spectral class.

Having generated the mass and, from that, the spectral class of the
star when on the main sequence, we need to generate the
luminosity, radius and effective surface temperature of each star.

Radius

The stellar radius is given by a theoretical approximation like:

   R is proportional to M^index

where the radius is proportional to the mass raised by some index,
0.4 for smaller stars, 0.8 for larger stars. The value of this index
depends upon the main cycles of thermonuclear reactions occurring
within the star's core. The approximation I adopted was:

       //Calculate radius
       double index = 0.4 + (Math.log10(mass[i])/3);
       if(index < 0.4)
       { index = 0.4; }
       radius[i] = Math.pow(mass[i], index);

Using a minimum index of 0.4 which increases gradually to about 0.8
for very massive stars.

Luminosity and Metallicity

A star's luminosity depends mainly on its mass and its metallicity
(or mean 'molecular mass'). These terms are slightly misleading
since in astrophysics any element heavier than helium is called a
'metal' for convenience and stars do not consist of molecules
(though molecules may form in the upper atmospheres of the cooler
stars) but of plasma, so we really mean 'mean particle mass' with the
particles being electrons, ions or neutral atoms. In other words,
metallicity gives us a measure of the chemical composition of the
star. First generation stars, formed first after the Big Bang, contain
only light elements, mainly hydrogen and helium, but some lithium
and beryllium. However, during their life stars synthesise heavier
elements, which are ultimately spread into space when the star dies.
stars which are formed from several generations of star dust, like the
Sun, have higher metallicities, and as the Universe ages, average
metallicity increases.

Population I stars are young metal-rich and highly luminous stars
found in the spiral arms of galaxies, and so concentrated in the
galactic plane, and associated with gas and dust. These are later
generation stars as the spiral arms of galaxies are regions of active
and periodic or continuous star formation.

Population II stars are old, red and metal-poor dwarfs,
concentrated in elliptical and lenticular galaxies and in the galactic
centres of spiral galaxies and the galactic halo, and also in a thicker
disc through the galactic plane. Smaller, cooler stars have longer
life-spans, and the smaller red dwarfs have life-spans in excess of
the age of the Universe and so these first-born stars still burn bright
to this day!

This scheme is a simplification, and intermediate types between
populations one and two exist.

With metallicity, m, set equal to 0.60 (a typical value for population 1
stars) the luminosity, L, can be estimated from the theoretical
approximation:

   L = 4.6 x M^3 x m^4

where the factor of 4.6 is arbitrary (as this is an approximation) but
chosen to agree with the observed luminosity for the Sun for M = 1
(in units of solar mass).

Effective Surface Temperature

The surface temperature of our model stars can be calculated once
we know the radius and the luminosity, using teh standard stellar
equation:

Using the Java code:

temperature[i] = Math.pow( (luminosity[i]*solarLum) / (4*Math.PI *
Math.pow(radius[i]*solarRad, 2) * sigma), 0.25);

Now we can calculate the life-span of each star on the main
sequence:

//Calculate MS life-span in billions of years
MSLifespan[i] = 10*mass[i]*(1/luminosity[i]);

We can also calculate the time spent by a forming proto-star before
it enters the main sequence (MS) but this is negligible compared to
the length of the MS life-span for all but the largest stars:

//Calculate protostellar life span in billions of years = thermal
timescale
PSLifespan[i] =
(1E15/60/60/24/365/1E9)*Math.pow(mass[i],2)*
(1/radius[i])*(1/luminosity[i]);

Now we can determine an arbitrary age for each star, according to
whatever model we chose, and then compare this to the set time
elapsed since the star was born and so determine what stage of its
life it is in - main sequence, giant branch or a post-stellar remnant
(such as a white dwarf or neutron star). If a star has passed the
end of the MS and become a red giant, or a later red supergiant,
white dwarf, neutron star, pulsar or black hole, then we can
recalculate its properties according to other models.

Red Giants

For red giants:

if (LS == Star.lifeStage.RG)
{
remnantType = "Red Giant";
remnantSurfaceT = 5000 - (mass*100 * generator.nextDouble());
remnantR = 10+(mass/10 * (generator.nextDouble()*2)) *
MSRadius;
remnantL = MSLuminosity * (1 +
Math.abs(generator.nextGaussian()));
  remnantM = mass;
}

We have simply used empirical ballpark figures with some
randomisation to generate the temperature, radius and luminosity
and we have assumed no mass loss, so the mass of the giant is
the same as in its MS youth. In reality, some mass will be lost as
stellar wind, especially in the giant stage and this could be factored
in. Red giants typically have a surface temperature of 5000K or
below and a radius of ten times that of the Sun, but the stars
luminosity remains approxiametely constant as it expands and
cools (though cooler, it is larger and so just as luminous).

Red Supergiants

A similar approach has been used for supergiants:

if (LS == Star.lifeStage.SG)
{
  remnantType = "Supergiant";
  remnantSurfaceT = 4500 - (mass*100 * generator.nextDouble());
  remnantR = 100+(mass/10 * (generator.nextDouble()*2)) *
MSRadius;
  remnantL = MSLuminosity * (1 +
Math.abs(generator.nextGaussian()));
  remnantM = mass;
}

With a slightly cooler temperature than a red giant and a radius of
about 100 times the Sun's radius! red supergiants are often
variable stars, undergoing oscillations in luminosity, which our
model does not incorporate.

White Dwarfs

For white dwarfs, the process is a bit more involved:

if(LS == Star.lifeStage.WD)
{
remnantType = "White Dwarf";
//Generate WD mass
if(mass <= 1)
{
      remnantM = (0.2+(0.2*mass)) * (generator.nextDouble()+0.5);
}
else
{
      //upper limit = (1.44 - mass/330)
      remnantM = Math.abs(generator.nextGaussian()) + 0.6;
      if(remnantM >= 1.11)
      {
           remnantM = (1.44 - mass/330);
       }
   }
   //Generate WD radius
   //Nauenberg empirical formula
   //remnantR = 7.83E06 * Math.pow( Math.pow(1.44/mass, 2/3) -
Math.pow(mass/1.44, 2/3), 0.5 );
   remnantR = (solarRad/74) * Math.pow(solarMass/mass, 1/3) *
0.001;
   lengthUnits = "km";
   //Generate WD surface temperature
   double cooling = age/1000 * (1/3.4);
   //Gaussian cooling
   remnantSurfaceT = ( 10000 * (generator.nextDouble() + 0.8 ) ) -
cooling;
   //Generate WD luminosity
   remnantL = 4*Math.PI*Math.pow((remnantR*1000),
2)*sigma*Math.pow(remnantSurfaceT, 4);
   remnantL = remnantL/solarLum;
   remnantL = remnantL + ( (generator.nextDouble()+0.01) *0.02);
}

The upper mass limit for white dwarfs varies between about 1.11 to
1.44 solar masses, depending upon the chemical composition
(which in turn depends upon the mass of the parent main
sequence star). This limit is called the Chandrasekhar mass. If
the remnant has more mass than this, then it will become a neutron
star / pulsar. Two different formula for calculating the radius of a
white dwarf, based upon its mass, are given, the second one being
used in this instance. The surface temperature of a white dwarf is
around 10 000 K when the remnant is newly formed. To this is
added some statistical variation and cooling - white dwarfs undergo
slow radiative cooling, so that older white dwarfs tend to be slightly
cooler. (Eventually, if enough time elapsed, they would cool into
black dwarfs). Using the temperature and radius, we can calculate
luminosity in the usual way by using the standard formula given in
the box above.

Neutron Stars and Pulsars

The Java code for generating a neutron star is:

if(LS == Star.lifeStage.NS)
{
   remnantType = "Neutron Star";
   remnantR = (10.0 + mass/10) *
(Math.abs(generator.nextGaussian())+1);
   lengthUnits = "km";
   //remnant mass: 1.11 to 3 solar masses
   remnantM = generator.nextDouble()*2 + (1.44 - mass/330);
   remnantSurfaceT = 10000 * (generator.nextDouble()*9 + 1);
   remnantL =
4*Math.PI*Math.pow(remnantR*1000,2)*sigma*Math.pow(remnantS
urfaceT,4);
   remnantL = remnantL/solarLum;
}

the radius is calculated using the main sequence star mass and a
minimum radius of 10 km and a certain randomness is introduced
(to account for the uncertain amount of mass the star may have
lost during its life and when undergoing a supernova explosion).
The mass is calculated using a minimum of 1.11 to 1.44 (the
Chandrasekhar limit) and an upper limit of about 3 solar masses
(above which a black hole is expected to form instead of a neutron
star). In generating the surface temperature we have assumed a
minimum of 10 000 K (really just a guess!) with a maximum of
about 100 000 K having been observed. Luminosity is then
calculated in the usual way.

It is assumed that 90% of neutron stars will be pulsars (rotating
neutron stars that beam radio pulses as it spins). These are given
one extra parameter - a pulse frequency or spin rate 9rotation
rate):

double spinRate = generator.nextDouble()*9999 + 1;

This generates a pulse frequency between 1 and 10 000 ms, in
accordance with observation.

Above: our Java program results for luminosity (on a log scale) and temperature for
each star plotted in Excel. The log (base 10) of luminosity in terms of solar units
(which set the luminosity of the Sun = 1) is plotted on the vertical axis and effective
surface temperature, in Kelvins, on the horizontal axis, starting with the highest
value on the left to follow convention. A diagram of this type is called a
Hertzsprung-Russell diagram (in this case for a main sequence population).

Java Code

The Java computer code for our star population generator is given below.

Design issues (for computer programmers): The main class drives the program and
prints out the results and save s them to a number of text files (MSHRData.txt for the
main sequence H-R data, HRData.txt for the evolved H-R diagram data and Stars.txt
for all the details about every star, main sequence and evolved). The Main class
generates an ArrayList of star objects and each star object contains all the stars
within its star system, since each Star object can be single, binary or multiple. If one
or more of these stars have left the main sequence, then the characteristics of their
stellar remnants are generated and stored in an ArrayList of Remnant class objects,
contained within the Star class. Originally I thought that the Remnant class might be
a subclass of Star, then i abandoned this idea and made it an external class.
However, I got tired writing accessor (getter and setter) methods when a Remnant
object should have access to the data of its parent Star. Thus I left the fields of
Remnant with default access (which in Java means that any class in the same
package can access them). This is generally considered bad practise. It would have
been better to make the Remnant class an inner class of the Star class, since every
star has a remnant, even at some time in the future, and this would have made
access to the required data secure.

Notice that each Star also contains an ArrayList of Planet objects. The Planet class
is still being developed and will be described later. This allows us to generate a
range of star systems with detailed descriptions of their planets. This is useful for
gaining insight into the possible variety of planetary systems, and also useful for
computer games, such as an adventure game (which we have planned for the near
future).

Above: the generated MS with approximate spectral classes based on surface
temperature. True stars enter the main sequence several hundreds of thousands to
tens of millions of years after their initial formation as protostars. The lowest
temperature that they can enter at, and be stable, is around 3000K 9the starting
point of the main sequence on a standard H-R diagram). Cooler stars are brown
dwarfs. Brown dwarfs only undergo one brief episode of nuclear fusion, as they fuse
tritium heavy hydrogen) an possibly one or two other uncommon isotopes.
Therefore, the brown dwarf branch of the 'main sequence' is not really the main
sequence proper. Our model does not make this distinction. According to theory,
brown dwarfs with masses down to about 0.01 solar masses form by stellar
processes, that is by the gravitational collapse of fragments of dark, cold clouds of
gas/dust and so we have modeled brown dwarfs as the tail-end of the main
sequence continuum (however, the final stages of their formation do differ, with
brown dwarfs failing to make it onto the main sequence proper, so our model is a
simplification).

The HR diagram generated for a later evolutionary time-point, above, and below -
part of the data seen closer up.

Below: the life-stages indicated: ASG asymptotic giant branch (supergiants); G:
giant branch (red giants); MS: stars still on the main sequence; NS: neutron stars,
most of which are pulsars; and WD: white dwarfs.

Neutron stars are not usually shown on H-R diagrams, because they are so faint,
due to their small size, that they disappear off the bottom! One feature of our model
that is lacking in realism is the lack of statistical variation on the main sequence - all
the stars pretty much fit a neat curve. In reality, there will be a spread of stars at
each temperature band either side of this curve.

Star Clusters and Stellar Populations

Stars generally form in groups from collapsing clouds of gas, which shrink under
their own gravitational force (they are said to be self-gravitating). The gravitational
force for a spherical cloud or solid object acts toward the centre of the mass, with all
matter closer to the centre contributing to the pull on matter outside, so given
enough mass the cloud will shrink. Heat energy opposes this collapse, however,
since hot molecules move around faster and this motion tends to disperse the cloud.
For a given mass there is a minimum temperature, above which the cloud will not
collapse and will instead spread outwards. Therefore, we say that stars form from
the gravitational collapse of cold clouds of matter. As the cloud collapses, it
fragments, as the matter is not evenly spread and local centres of high density form
and collapse independently. Each of these fragments can become a star system
containing one or more stars orbiting their common centre of gravity - it may
fragment again with each fragment forming an individual star. Thus, each
star-forming nebula will form a cluster of stars.

Open Clusters

These contain loose aggregates of less than 20 to over a thousand stars. These are
Population I stars and these clusters or situated in or close to the galactic plane of
spiral galaxies. Theory predicts that these clusters only survive intact for one or two
galactic revolutions, and so all are young. As these clusters rotate around the
galactic centre, they spread out and intermingle, so that a population may result that
contains a stars of mixed ages. (Our model assumes an even mixture of stars aged
between

Within Population I stars, there are two population sub-types, depending upon the
age of the open cluster. Extreme Population I stars belong to young open
clusters. These are HII regions (see energetic processes) as they still contain
much of the gas and dust of the original star-forming nebula. These clusters follow
circular orbits and occur in the spiral arms of spiral galaxies. They contain many T
Tauri stars, OB associations (see main sequence stars), supergiants and classical
Cepheids.

Older open clusters contain older Population I stars with strong spectral lines
(since the stars have had time to synthesise more heavy elements). They contain
more Me dwarfs (red dwarfs of spectral class M with strong emission lines due to
dusty shells around them), A stars and giants. They tend to be more concentrated
toward the galactic centre and have approximately circular orbits around the galactic
centre. The Sun belongs to such an old open cluster.

Globular Clusters

These are compact clusters of between about ten thousand and one million stars.
The stars are more densely packed towards the centre of the cluster, reaching
densities of 1000 stars per cubic parsec (about 29 stars per cubic light-year!). Most
occur in giant and elliptical galaxies, but about 150 occur in the Milky-Way Galaxy
(MWG). They have very eccentric orbits and are distributed in a rough sphere
around the galactic centre - the Galactic Halo (and are not concentrated toward the
galactic plane as are open clusters). The youngest occur in the galactic disc (about
20%) where the orbits are more circular. The stars in globular clusters are
Population II stars and all are older than the Sun. Due to their great age, being
born nearer the beginning of the universe, they have very low metallicities. Many are
red giants and bluer horizontal branch giants, X-ray binaries, Cataclysmic Variables
and RR Lyrae stars.

Hertzsprung-Russel (H-R) Diagrams

package randomspace;

import java.util.*;

public class Star
{
//Constants
public static final double sigma = 5.671E-08; //Jm^-2K^-4s^-1
public static final double solarMass = 1.99E+30; // kg
public static final double solarLum = 3.83E+26; // J/s
public static final double solarRad = 6.96E+08; // m
public static final double metal = 0.60;

int multiplicity;
String[] designation;
enum spectralClass
{ BD, M, K, G, F, A, B, O };
spectralClass SC[];
enum lifeStage
{ PS, MS, RG, SG, PN, WD, NS, P, BH };
lifeStage LS[];
double mass[];
double maxAge = 10; //billions of years
double minAge = 0.01; //billions of years
double age;
double MSLifespan[];
double PSLifespan[];
int spaceVolume = 0;
double XCoord;
double YCoord;
double ZCoord;
double luminosity[];
double temperature[];
double radius[];
int numberPlanets;
static int starSystemCount;
Planet planets[];
Random generator = new Random();
Remnant remnants[];

/** Creates a new instance of Star */
public Star()
{
    starSystemCount++;
    //Generate number of stars in system
    multiplicity = (int)(Math.abs(generator.nextGaussian()*2));
    if(multiplicity <= 0)
    { multiplicity = 1; }

    designation = new String[multiplicity+1];
    mass = new double[multiplicity];
    SC = new spectralClass[multiplicity];
    LS = new lifeStage[multiplicity];
    planets = new Planet[multiplicity];
    radius = new double[multiplicity];
    luminosity = new double[multiplicity];
    temperature = new double[multiplicity];
    MSLifespan = new double[multiplicity];
    PSLifespan = new double[multiplicity];
    remnants = new Remnant[multiplicity];
    double mean = -1;
    double stdDev = 1.05;

    //Generate age of system in billions of years
    age = generator.nextDouble() * maxAge + minAge;

    //Star system designation
    //designation[0] = "(" + Double.toString(XCoord) + "," + Double.toString(YCoord) + "," + Double.toString(ZCoord) + ")";
    designation[0] = Integer.toString(starSystemCount);

    //Generate attributes of each star in system
    for(int i = 0; i < multiplicity; i++)
    {
        designation[i+1] = "" + Integer.toString(i+1);
        //Generate mass
        mass[i] = 0;
        while (mass[i] < 0.01)
        {
            mass[i] = generator.nextGaussian()*stdDev + mean;
            mass[i] = mass[i]*Math.exp(mass[i]*0.35)-Math.exp(-mass[i]*0.1);
        }
        //Generate spectral class from mass
        if(mass[i] >= 0.01 && mass[i] < 0.075)
        {
            //Brown dwarf
            SC[i] = spectralClass.BD;
        }
        if(mass[i] >= 0.075 && mass[i] < 0.5)
        {
            //M
            SC[i] = spectralClass.M;
        }
        if(mass[i] >= 0.5 && mass[i] < 1.0)
        {
            //K
            SC[i] = spectralClass.K;
        }
        if(mass[i] >= 1.0 && mass[i] < 1.5)
        {
            //G
            SC[i] = spectralClass.G;
        }
        if(mass[i] >= 1.5 && mass[i] < 2.5)
        {
            //F
            SC[i] = spectralClass.F;
        }
        if(mass[i] >= 2.5 && mass[i] < 10)
        {
            //A
            SC[i] = spectralClass.A;
        }
        if(mass[i] >= 10 && mass[i] < 40)
        {
            //B
            SC[i] = spectralClass.B;
        }
        if(mass[i] >= 40)
        {
            //O
            SC[i] = spectralClass.O;
        }

        //Generate number of planets
        numberPlanets = (int)(generator.nextGaussian()*10);
        if (numberPlanets < 0)
        { numberPlanets = 0; }

        //Calculate radius
        double index = 0.4 + (Math.log10(mass[i])/3);
        if(index < 0.4)
        { index = 0.4; }
        //if(index > 0.8)
        //{ index = 0.8; }
        radius[i] = Math.pow(mass[i], index);

        //Calculate luminosities
        luminosity[i] = 4.6*Math.pow(mass[i], 3)*Math.pow(metal, 4);
        //Add luminosity variation
        //luminosity[i] = luminosity[i] + ( 0 + ( (generator.nextDouble() - 0.5 ) * 0.01) );

        //Calculate surface temperatures
        temperature[i] = Math.pow( (luminosity[i]*solarLum) / (4*Math.PI * Math.pow(radius[i]*solarRad, 2) * sigma), 0.25);

        //Calculate MS life-span in billions of years
        MSLifespan[i] = 10*mass[i]*(1/luminosity[i]);

        //Calculate protostellar life span in billions of years = thermal timescale
        PSLifespan[i] = (1E15/60/60/24/365/1E9)*Math.pow(mass[i],2)*(1/radius[i])*(1/luminosity[i]);

         //Generate life stage from age and mass
        if(MSLifespan[i] <= PSLifespan[i])
        {
            LS[i] = lifeStage.PS;
        }
        else
        {
            LS[i] = lifeStage.MS;
        }
        double remnantAge = age % MSLifespan[i];
        if(age > MSLifespan[i])
        {
            //double remnantAge = age % MSLifespan[i];
            if(remnantAge <= 0.001) { LS[i] = lifeStage.RG; }
            if(remnantAge >= 0.001 && remnantAge <= 0.0015 ) { LS[i] = lifeStage.SG; }
            if(remnantAge >= 0.0015 && mass[i] >= 0.7 && mass[i] < 5 ) { LS[i] = lifeStage.WD; }
            double Rand = generator.nextDouble();
            if(remnantAge >= 0.0015 && mass[i] >= 5 && mass[i] < 60 && Rand <= 0.1) { LS[i] = lifeStage.NS; }
            if(remnantAge >= 0.0015 && mass[i] >= 5 && mass[i] < 60 && Rand > 0.1) { LS[i] = lifeStage.P; }
            if(remnantAge >= 0.0015 && mass[i] >= 60) { LS[i] = lifeStage.BH; }
        }
        else
        {
            LS[i] = lifeStage.MS;
        }

        remnants[i] = new Remnant(LS[i], mass[i], radius[i], luminosity[i], temperature[i], remnantAge);
    }

}

public Star(int spaceSize)
{
    this();
    spaceVolume = spaceSize;
    this.setXCoord();
    this.setYCoord();
    this.setZCoord();
}

public String toString()
{
    StringBuffer masses = new StringBuffer();
    for(int i = 0; i < multiplicity; i++)
    {
        masses.append(mass[i]);
        masses.append(", ");
    }
    StringBuffer radii = new StringBuffer();
    for(int i = 0; i < multiplicity; i++)
    {
        radii.append(radius[i]);
        radii.append(", ");
    }
    StringBuffer lumin = new StringBuffer();
    for(int i = 0; i < multiplicity; i++)
    {
        lumin.append(luminosity[i]);
        lumin.append(", ");
    }
    StringBuffer temp = new StringBuffer();
    for(int i = 0; i < multiplicity; i++)
    {
        temp.append(temperature[i]);
        temp.append(", ");
    }
    StringBuffer spectral = new StringBuffer();
    for(int i = 0; i < multiplicity; i++)
    {
        spectral.append(SC[i]);
        spectral.append(", ");
    }
    StringBuffer lifespans = new StringBuffer();
    for(int i = 0; i < multiplicity; i++)
    {
        lifespans.append(MSLifespan[i]);
        lifespans.append(", ");
    }
    StringBuffer lifestages = new StringBuffer();
    for(int i = 0; i < multiplicity; i++)
    {
        lifestages.append(LS[i]);
        lifestages.append(", ");
    }
    StringBuffer catDes = new StringBuffer();
    for(int i = 0; i < multiplicity+1; i++)
    {
        catDes.append(designation[i]);
        catDes.append(", ");
    }
    StringBuffer stellarRemnants = new StringBuffer();
    for(int i = 0; i < multiplicity; i++)
    {
        stellarRemnants.append(remnants[i]);
        stellarRemnants.append(", ");
    }
    String thisStar = ("Designation = " + catDes + "Coords = (" + (int)XCoord + "," + (int)YCoord + "," + (int)ZCoord + "); " + "Multiplicity = " + multiplicity
            + "; Mass(es) = " + masses.toString() + "; Planets = " + numberPlanets + "; R = " + radii + "; L = " + lumin + "; T = " + temp + "; age = " + age
            + "; SC = " + spectral + "; Life-span(s) = " + lifespans + "; Life stage(s) = " + lifestages + "\nRemnants: " + stellarRemnants);
    return thisStar;
}

//Getters and Setters
//Multiplicity
public void setMultiplicity(int multiplicity)
{
    this.multiplicity = multiplicity;
}
public int getMultiplicity()
{
    return multiplicity;
}
//SC
public void setSpectralClass(spectralClass[] SC)
{
    this.SC = SC;
}
public spectralClass[] getSpectralClass()
{
    return SC;
}
//Masses
public double[] getMasses()
{
    return mass;
}
//Life stages
public lifeStage[] getLifeStages()
{
    return LS;
}
//Space Volume for coordinate generation
public void setSpaceSize(int spaceSize)
{
    //sets the length of a cuboidal region of space
    spaceVolume = spaceSize;
}
public int getSpaceSize()
{
    //Returns the length of the local region of space
    return spaceVolume;
}
//Get surface temperatures
public double[] getTemperatures()
{
    return temperature;
}
//Get luminosities
public double[] getLuminosities()
{
    return luminosity;
}
//Get coords
public double getXCoord()
{ return XCoord; }
public double getYCoord()
{ return YCoord; }
public double getZCoord()
{ return ZCoord; }
//Set coords
public void setXCoord()
{ XCoord = generator.nextDouble() * spaceVolume; }
public void setYCoord()
{ YCoord = generator.nextDouble() * spaceVolume; }
public void setZCoord()
{ ZCoord = generator.nextDouble() * spaceVolume; }
//Get remnants
public Remnant[] getRemnants()
{
    return remnants;
}
}

package randomspace;

import java.util.*;

public class Remnant
{
public static final double sigma = 5.671E-08; //Jm^-2K^-4s^-1
public static final double solarMass = 1.99E+30; // kg
public static final double solarLum = 3.83E+26; // J/s
public static final double solarRad = 6.96E+08; // m

double remnantM;
double remnantR;
double remnantL;
double remnantSurfaceT;
double remnantDensity;
String remnantType = "Not remnant";
String lengthUnits = "";

private Random generator = new Random();

/** Creates a new instance of Remnant */
public Remnant()
{
    super();
}

public Remnant(Star.lifeStage LS, double mass, double MSRadius, double MSLuminosity, double MSTemp, double age)
{
    if(LS == Star.lifeStage.MS)
    {
        remnantType = "Main Sequence";
        remnantSurfaceT = MSTemp;
        remnantR = MSRadius;
        remnantM = mass;
        remnantL = MSLuminosity;
    }
    if (LS == Star.lifeStage.RG)
    {
        remnantType = "Red Giant";
        remnantSurfaceT = 5000 - (mass*100 * generator.nextDouble());
        remnantR = 10+(mass/10 * (generator.nextDouble()*2)) * MSRadius;
        remnantL = MSLuminosity * (1 + Math.abs(generator.nextGaussian()));
        remnantM = mass;
    }
    if (LS == Star.lifeStage.SG)
    {
        remnantType = "Supergiant";
        remnantSurfaceT = 4500 - (mass*100 * generator.nextDouble());
        remnantR = 100+(mass/10 * (generator.nextDouble()*2)) * MSRadius;
        remnantL = MSLuminosity * (1 + Math.abs(generator.nextGaussian()));
        remnantM = mass;
    }
    if(LS == Star.lifeStage.WD)
    {
        remnantType = "White Dwarf";
        //Generate WD mass
        if(mass <= 1)
        {
            remnantM = (0.2+(0.2*mass)) * (generator.nextDouble()+0.5);
        }
        else
        {
            //upper limit = (1.44 - mass/330)
            remnantM = Math.abs(generator.nextGaussian()) + 0.6;
            if(remnantM >= 1.11)
            {
                remnantM = (1.44 - mass/330);
            }
        }
        //Generate WD radius
        //Nauenberg empirical formula
        //remnantR = 7.83E06 * Math.pow( Math.pow(1.44/mass, 2/3) - Math.pow(mass/1.44, 2/3), 0.5 );
        remnantR = (solarRad/74) * Math.pow(solarMass/mass, 1/3) * 0.001;
        lengthUnits = "km";
        //Generate WD surface temperature
        double cooling = age/1000 * (1/3.4);
        //Gaussian cooling
        remnantSurfaceT = ( 10000 * (generator.nextDouble() + 0.8 ) ) - cooling;
        //Generate WD luminosity
        remnantL = 4*Math.PI*Math.pow((remnantR*1000), 2)*sigma*Math.pow(remnantSurfaceT, 4);
        remnantL = remnantL/solarLum;
        remnantL = remnantL + ( (generator.nextDouble()+0.01) *0.02);
    }
    if(LS == Star.lifeStage.NS)
    {
        remnantType = "Neutron Star";
        remnantR = (10.0 + mass/10) * (Math.abs(generator.nextGaussian())+1);
        lengthUnits = "km";
        //remnant mass: 1.11 to 3 solar masses
        remnantM = generator.nextDouble()*2 + (1.44 - mass/330);
        remnantSurfaceT = 10000 * (generator.nextDouble()*9 + 1);
        remnantL = 4*Math.PI*Math.pow(remnantR*1000,2)*sigma*Math.pow(remnantSurfaceT,4);
        remnantL = remnantL/solarLum;
    }
    if(LS == Star.lifeStage.P)
    {
        remnantType = "Pulsar";
        double spinRate = generator.nextDouble()*9999 + 1; // 1 to 10 000 ms
        remnantType = "Pulsar " + spinRate + "ms";
        remnantR = (10.0 + mass/10) * (generator.nextGaussian()+1);
        remnantM = generator.nextDouble()*2 + (1.44 - mass/330);
        remnantSurfaceT = 10000 * (generator.nextDouble()*9 + 1);
        remnantL = 4*Math.PI*Math.pow(remnantR*1000,2)*sigma*Math.pow(remnantSurfaceT,4);
        remnantL = remnantL/solarLum;
    }

}

public String toString()
{
    String myDescription = "";
    myDescription = "\nRemnant = " + remnantType + "; T = " + Double.toString(remnantSurfaceT)
    + "; R = " + Double.toString(remnantR) + lengthUnits + "; L = " + Double.toString(remnantL) + "; M = " + Double.toString(remnantM);
    return myDescription;
}

}

package randomspace;

import java.util.*;
import java.io.*;

public class Main
{

private Random generator = new Random(1);
private static int number = 0;
private static double totalMass = 0;
private static double maxMass = 0;
private static double minMass = 1;
private static int BDCount = 0;
private static int MCount = 0;
private static int KCount = 0;
private static int GCount = 0;
private static int FCount = 0;
private static int ACount = 0;
private static int BCount = 0;
private static int OCount = 0;

private static int PSCount = 0;
private static int RGCount = 0;
private static int ASGCount = 0;
private static int WDCount = 0;
private static int NSCount = 0;
private static int PCount = 0;
private static int BHCount = 0;

private static final double length = 32.6*3*1.4; // light-years
private static final double density = 0.0081; //0.003 to 1.5 per cubic light-year
private static double volume = Math.pow(length, 3);

private static ArrayList<Star> starSystems = new ArrayList<Star>();
private static int numStars = (int)(density*volume);
private static int singles = 0;
private static int binaries = 0;
private static int triples = 0;
private static int others = 0;
private static int starCount = 0;

/** Creates a new instance of Main */
public Main()
{
}

public static void main(String[] args)
{
    generateStars();
    displayResults();
    try
    {
        saveResults();
    }
    catch ( IOException io )
    {
        System.out.println("Error accessing Stars file " + io.getMessage());
    }
    //testPlanets();
    HRData();
    MSHRData();
}

public static void generateStars()
{
    //Galaxy
    for(int i = 0; i < numStars; i++)
    {
        Star aStar = new Star((int)volume);
        starSystems.add(aStar);
    }

    //Generate star stats
    for(int i = 0; i < numStars; i++)
    {
        //Multiplicity stats
        number = starSystems.get(i).getMultiplicity();
        starCount = starCount + number;
        String myStar = (starSystems.get(i).toString());
        System.out.println(myStar);
        if (number == 1)
        { singles = singles + 1; }
        else if (number == 2)
        { binaries++; }
        else if (number == 3)
        { triples++; }
        else
        { others++; }

        //Mass stats
        double[] myMasses = starSystems.get(i).getMasses();
        for(int j = 0; j < number; j++)
        {
            totalMass = totalMass + myMasses[j];
            //Max mass
            if (myMasses[j] > maxMass)
            { maxMass = myMasses[j]; }
            //Min mass
            if (myMasses[j] < minMass)
            { minMass = myMasses[j]; }
            //Spectral distributions
            if(myMasses[j] >= 0.01 && myMasses[j] < 0.075)
            { BDCount++; }
            if(myMasses[j] >= 0.075 && myMasses[j] < 0.5)
            { MCount++; }
            if(myMasses[j] >= 0.5 && myMasses[j] < 1.0)
            { KCount++; }
            if(myMasses[j] >= 1.0 && myMasses[j] < 1.5)
            { GCount++; }
            if(myMasses[j] >= 1.5 && myMasses[j] < 2.5)
            { FCount++; }
            if(myMasses[j] >= 2.5 && myMasses[j] < 10)
            { ACount++; }
            if(myMasses[j] >= 10 && myMasses[j] < 40)
            { BCount++; }
            if(myMasses[j] >= 40)
            { OCount++; }
        }

        //Life stage stats
        Star.lifeStage[] myLifeStages;
        myLifeStages = starSystems.get(i).getLifeStages();
        for(int j = 0; j < number; j++)
        {
            if(myLifeStages[j] == Star.lifeStage.RG)
            { RGCount++; }
            if(myLifeStages[j] == Star.lifeStage.SG)
            { ASGCount++; }
            if(myLifeStages[j] == Star.lifeStage.WD)
            { WDCount++; }
            if(myLifeStages[j] == Star.lifeStage.NS)
            { NSCount++; }
            if(myLifeStages[j] == Star.lifeStage.P)
            { PCount++; }
            if(myLifeStages[j] == Star.lifeStage.BH)
            { BHCount++; }
        }

    }
}

public static void displayResults()
{
    //display multiplicity stats
    System.out.println("\n" + "Total stars = " + starCount);
    System.out.println("Number of star systems = " + numStars);
    double averageMultiplicity = (double)(starCount)/numStars;
    System.out.println("Multiplicity = " + Double.toString(averageMultiplicity));
    System.out.println("Single star systems = " + singles);
    System.out.println("Double star systems = " + binaries);
    System.out.println("Triple star systems = " + triples);
    System.out.println("Star systems with > 3 stars = " + others);
    System.out.println("Percentage of stars in single systems = " + ((double)singles/starCount)*100);

    //display mass stats
    System.out.println("Mean stellar mass = " + totalMass/starCount);
    System.out.println("Maximum mass = " + maxMass);
    System.out.println("Minimum mass = " + minMass);

    //display spectral type stats
    System.out.println("Number BD stars = " + BDCount);
    System.out.println("Number M stars = " + MCount);
    System.out.println("Number K stars = " + KCount);
    System.out.println("Number G stars = " + GCount);
    System.out.println("Number F stars = " + FCount);
    System.out.println("Number A stars = " + ACount);
    System.out.println("Number B stars = " + BCount);
    System.out.println("Number O stars = " + OCount);

    //display life stage stats
    System.out.println("Number RG stars = " + RGCount);
    System.out.println("Number ASG stars = " + ASGCount);
    System.out.println("Number WD stars = " + WDCount);
    System.out.println("Number NS stars = " + NSCount);
    System.out.println("Number Pulsars = " + PCount);
    System.out.println("Number BH stars = " + BHCount);
}

public static void saveResults() throws IOException
{
    PrintWriter pw = new PrintWriter("Stars.txt");

    for(int i = 0; i < numStars; i++)
    {
        //individual stars
        number = starSystems.get(i).getMultiplicity();
        String myStar = (starSystems.get(i).toString());
        pw.println(myStar);
    }

    //save multiplicity stats
    pw.println("\n" + "Total stars = " + starCount);
    pw.println("Number of star systems = " + numStars);
    double averageMultiplicity = (double)(starCount)/numStars;
    pw.println("Multiplicity = " + Double.toString(averageMultiplicity));
    pw.println("Single star systems = " + singles);
    pw.println("Double star systems = " + binaries);
    pw.println("Triple star systems = " + triples);
    pw.println("Star systems with > 3 stars = " + others);
    pw.println("Percentage of stars in single systems = " + ((double)singles/starCount)*100);

    //save mass stats
    pw.println("Mean stellar mass = " + totalMass/starCount);
    pw.println("Maximum mass = " + maxMass);
    pw.println("Minimum mass = " + minMass);

    //save spectral type stats
    pw.println("Number BD stars = " + BDCount);
    pw.println("Number M stars = " + MCount);
    pw.println("Number K stars = " + KCount);
    pw.println("Number G stars = " + GCount);
    pw.println("Number F stars = " + FCount);
    pw.println("Number A stars = " + ACount);
    pw.println("Number B stars = " + BCount);
    pw.println("Number O stars = " + OCount);

    //save life stage stats
    pw.println("Number RG stars = " + RGCount);
    pw.println("Number ASG stars = " + ASGCount);
    pw.println("Number WD stars = " + WDCount);
    pw.println("Number NS stars = " + NSCount);
    pw.println("Number BH stars = " + BHCount);

    pw.close();
}

public static void HRData() //throws IOException
{
    try
    {
        PrintWriter pw = new PrintWriter("HRData.txt");
        pw.println("T, L");
        for(int i = 0; i < numStars; i++)
        {
            //individual stars
            number = starSystems.get(i).getMultiplicity();
            double[] temps = starSystems.get(i).getTemperatures();
            double[] lumins = starSystems.get(i).getLuminosities();
            Remnant[] myRemnants = starSystems.get(i).getRemnants();
            for(int j = 0; j < number; j++)
            {
                pw.print(myRemnants[j].remnantSurfaceT + ", " + myRemnants[j].remnantL + '\n');
            }
            pw.println();
        }
        pw.close();
    }
    catch (IOException io)
    {
        System.out.println("Error accessing HRData file! " + io.getMessage());
    }
}

public static void MSHRData()
{
    try
    {
        PrintWriter pw = new PrintWriter("MSHRData.txt");
        pw.println("T, L");
        for(int i = 0; i < numStars; i++)
        {
            //individual stars
            number = starSystems.get(i).getMultiplicity();
            double[] temps = starSystems.get(i).getTemperatures();
            double[] lumins = starSystems.get(i).getLuminosities();
            for(int j = 0; j < number; j++)
            {
                    pw.print(temps[j] + "," + lumins[j] + '\n');
            }
        }
        pw.close();
    }
    catch (IOException io)
    {
            System.out.println("Error accessing MSHRData file " + io.getMessage());
    }
}

}

Above: part of a main sequence group of stars generated by our Java program (I added in a function to
generate a Pov-Ray file and rendered it) shown in orthographic camera view and below in perspective view.
Click the images to enlarge them. See if you can spot a brown dwarf (they are small and faint!). Somewhere in
the group there are 16 B giants, but none are visible in this view!

Analyse a sample output of some 200 star systems with the Astro Navigator.