aboutsummaryrefslogtreecommitdiff
path: root/noao/artdata/lists/stmix.x
blob: 4254e8eb6312738315d02355ef9eef85db168ba0 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
include	<math.h>
include "starlist.h"

# ST_ESMIX -- Compute the percentage of elliptical galaxies.

procedure st_esmix (egal, nstar, esmix, seed)

int	egal[ARB]		# array of types
int	nstar			# number of objects
real	esmix			# fraction of elliptical galaxies
long	seed			# seed for random number generator

int	i
real	fraction
real	urand()

begin
	do i = 1, nstar {
	    fraction = urand (seed)
	    if (fraction <= esmix)
	        egal[i] = ST_DEVAUC
	    else
	        egal[i] = ST_EXP
	}
end


# ST_ROUND -- Compute an array of roundness parameters.
# For ellipical models use a uniform distribution of axial ratios.
# For spiral models use a uniform inclination distribution.  Compute
# the axial ratio with a .1 flatness parameter from Holmberg (Stars
# and Stellar Systems).  Then add cosecant absorption factor to the
# the magnitudes with a limit of 10 in the cosecant.

procedure st_round (egal, mag, round, nstars, ar, alpha, seed)

int	egal[ARB]		# array of galaxy types
real	mag[ARB]		# magnitudes
real	round[ARB]		# array of roundness values
int	nstars			# number of stars
real	ar			# minimum roundness value
real	alpha			# absorption coefficent
long	seed			# seed for the random number generator

int	i
real	dr, s, urand()

begin
	dr = (1. - ar)
	do i = 1, nstars {
	    if (egal[i] == ST_DEVAUC)
	        round[i] = ar + dr * urand (seed)
	    else {
		s = sin (HALFPI * urand (seed))
		round[i] = sqrt (s**2 * .99 + .01)
		mag[i] = mag[i] + alpha * (1 / max (0.1, s) - 1) / 9.
	    }
	}

end


# ST_PHI -- Compute an array of position angles.

procedure st_phi (phi, nstars, seed)

real	phi[ARB]		# array of position angles
int	nstars			# number of stars
long	seed			# seed for random number generator

int	i
real	urand()

begin
	do i = 1, nstars
	    phi[i] = urand (seed)
	call amapr (phi, phi, nstars, 0.0, 1.0, 0.0, 360.0)
end


# ST_DIAMETERS -- Compute the effective diameters of the galaxies based
# on their magnitudes.  The relation used is from Holmberg (Stars and
# Stellar Systems).  A uniform dispersion of 20% is added.

procedure st_diameters (mag, egal, axis, nstars, minmag, maxmag, eradius,
	sradius, seed)

real	mag[ARB]		# input array of magnitudes
int	egal[ARB]		# array of galaxy types
real	axis[ARB]		# output array of diameters
int	nstars			# number of stars
real	minmag			# minimum magnitude
real	maxmag			# maximum magnitude
real	eradius			# maximum elliptical radius
real	sradius			# maximum spiral radius
long	seed			# seed for random number generator

int	i
real	urand()

begin
	do i = 1, nstars {
	    if (egal[i] == ST_DEVAUC)
	        axis[i] = eradius * 10.0 ** ((minmag - mag[i]) / 6.)
	    else
	        axis[i] = sradius * eradius * 10.0 ** ((minmag - mag[i]) / 6.)
	    axis[i] = axis[i] * (0.8 + 0.4 * urand (seed))
	}
end


# ST_ZDIAMETERS -- Compute the effective diameters of the galaxies based
# on their magnitudes and redshift. The relation used is that the redshift
# is proportional to the luminousity and the diameters include the
# factor of (1+z)**2.  A uniform dispersion of 50% is added.

procedure st_zdiameters (mag, egal, axis, nstars, minmag, maxmag,
	z, eradius, sradius, seed)

real	mag[ARB]		# input array of magnitudes
int	egal[ARB]		# array of galaxy types
real	axis[ARB]		# output array of diameters
int	nstars			# number of stars
real	minmag			# minimum magnitude
real	maxmag			# maximum magnitude
real	z			# minumum redshift
real	eradius			# maximum elliptical radius
real	sradius			# maximum spiral radius
long	seed			# seed for random number generator

int	i
real	z0, z1, urand()

begin
	z0 = z / (1 + z) ** 2
	do i = 1, nstars {
	    z1 = z * 10.0 ** (-0.2 * (minmag - mag[i]))
	    if (egal[i] == ST_DEVAUC)
	        axis[i] = eradius * z0 * (1 + z1) ** 2 / z1
	    else
	        axis[i] = sradius * eradius * z0 * (1 + z1) ** 2 / z1
	    axis[i] = axis[i] * (0.5 + urand (seed))
	}
end