initial import

--HG--
branch : stargen
This commit is contained in:
cecilkorik 2010-04-27 03:43:05 +00:00
commit 3066899c07
3 changed files with 372 additions and 0 deletions

22
rnd.py Normal file
View file

@ -0,0 +1,22 @@
import random as random_lib
import math
def random():
return random_lib.random()
def gaussian(mean=0.0, stddev=1.0):
while True:
x1 = (2.0 * random()) - 1.0
x2 = (2.0 * random()) - 1.0
w = (x1 * x1) + (x2 * x2)
if w < 1.0:
break
w = math.sqrt((-2.0 * math.log(w)) / w)
return mean + (x1 * w * stddev)
def randrange(min, max):
return min + (random() * (max - min))

176
stargen.py Normal file
View file

@ -0,0 +1,176 @@
from rnd import *
import math
systemsize = 6e9
stars = (1.0, 0.5)
planets = (8.0, 4.0)
star_spectra = ['O', 'B', 'A', 'F', 'G', 'K', 'M']
star_lumin = [(500000.0, 30000.0), (30000.0, 25.0), (25.0, 5.0), (5.0, 1.5), (1.5, 0.6), (0.6, 0.008), (0.08, 0.0008)]
star_masses = [(40.0, 16.0), (16.0, 2.1), (2.1, 1.4), (1.4, 1.04), (1.04, 0.8), (0.8, 0.45), (0.45, 0.1)]
star_radii = [(22.0, 8.0), (7.0, 1.8), (1.8, 1.4), (1.4, 1.15), (1.15, 0.96), (0.96, 0.7), (0.7, 0.25)]
star_prob = [0.0, 0.05, 0.35, 0.6, 0.8, 0.97, 1.0]
num_stars = round(gaussian(*stars))
if num_stars < 1.0:
num_stars = 1.0
num_stars = int(num_stars)
num_planets = round(gaussian(*planets))
if num_planets < 0.0:
num_planets = 0.0
num_planets = int(num_planets)
stars = []
x = y = z = 0.0
while len(stars) < num_stars:
r = random()
for i, e in enumerate(star_prob):
if r <= e:
si = i
break
r = random()
slumin = (star_lumin[si][1] * 0.8) + (star_masses[si][0] * 1.2 * r)
r2 = r + (random() * 0.3) - 0.15
smass = (star_masses[si][1] * 0.8) + (star_masses[si][0] * 1.2 * r)
r2 = r + (random() * 0.3) - 0.15
sradius = (star_radii[si][1] * 0.8) + (star_radii[si][0] * 1.2 * r)
print "Star lumin %s mass %s class %s" % (slumin, smass, star_spectra[si])
stars += [(slumin, smass, star_spectra[si], sradius, x, y, z)]
# first star is always at 0.0, 0.0, 0.0, next star gets random coords
x = gaussian(0.0, systemsize)
y = gaussian(0.0, systemsize)
if len(stars) > 1:
z = gaussian(0.0, systemsize / 10e3)
4
planets = []
sys_size = math.sqrt(sys_lumin)
while len(planets) < num_planets:
def dist(x, y, z):
return math.sqrt(x**2 + y**2 + z**2)
x = gaussian(0.0, systemsize * sys_size)
y = gaussian(0.0, systemsize * sys_size)
z = gaussian(0.0, systemsize * sys_size / 10e3)
reject = False
for star in stars:
slumin, smass, sclass, sradius, sx, sy, sz
d = dist(x-sx,y-sy,z-sz)
deathzone = sradius + (sradius * 0.25 * math.sqrt(slumin))
if d < deathzone:
reject = True
if reject:
continue
r = random()
rocky = False
gas_giant = False
if r < 0.6:
# rock planet
rocky = True
massmin = 1.0e19
massmax = 8.8e25
mass = 0.0
while mass < massmin or mass > massmax:
massexp = gaussian(24.5, 1.25)
mass = 10.0 ** massexp
density = 0.0
while density < 0.75:
density = gaussian(4.5, 1.0)
mag = -1.0
while mag < 0.0:
mag = gaussian(density, 3.0)
else:
# gas giant
gas_giant = True
massmin = 1e25
massmax = 4e28
mass = 0.0
while mass < massmin or mass > massmax:
massexp = gaussian(26.5, 1.0)
mass = 10.0 ** massexp
density = 0.0
while density < 0.3:
density = gaussian(1.5, 0.75)
mag = -1.0
while mag < 2.0:
mag = gaussian(7.0, 3.0)
volume = mass / density / 1e12
radius = pow(volume * 0.75 / math.pi, 1.0/3.0)
circsurfarea = math.pi * radius**2
surfarea = circsurfarea * 4.0
insol = 0.0
for star in stars:
slumin, smass, sclass, sradius, sx, sy, sz
d = dist(x-sx,y-sy,z-sz)
spower = 3.8e26 * slumin
distsphere = 4.0 * math.pi * d**2
insol += (circsurfarea / distsphere) * spower
areainsol = insol / surfarea / 1e6
if areainsol > 75000:
continue
atmo = None
while atmo == None or atmo > density:
if gas_giant:
atmo_exp = randrange(-3, 0.5)
else:
atmo_exp = randrange(-9.2, -0.9)
atmo = 10 ** atmo_exp
albedo = -1.0
while albedo < 0.0 or albedo > 1.0:
albedo = gaussian(0.35, 0.05)
emissivity = 1.0 - albedo
emissivity = randrange(0.00001, 1.0)
emiss_diff = (((10.0 - pow(atmo, 1.0/4.0)) / 10.0) ** 100.0) / 2.0
emissivity_min = emissivity - (emissivity * emiss_diff)
emissivity_max = emissivity + (emissivity * emiss_diff)
avg_temp = pow(((1.0 - albedo) * areainsol) / (4.0 * emissivity * 5.67e-8), 1.0/4.0)
min_temp = pow(((1.0 - albedo) * areainsol) / (4.0 * emissivity_min * 5.67e-8), 1.0/4.0)
max_temp = pow(((1.0 - albedo) * areainsol) / (4.0 * emissivity_max * 5.67e-8), 1.0/4.0)

174
stargen.py.bak Normal file
View file

@ -0,0 +1,174 @@
from rnd import *
import math
systemsize = 6e9
stars = (1.0, 0.5)
planets = (8.0, 4.0)
star_spectra = ['O', 'B', 'A', 'F', 'G', 'K', 'M']
star_lumin = [(500000.0, 30000.0), (30000.0, 25.0), (25.0, 5.0), (5.0, 1.5), (1.5, 0.6), (0.6, 0.008), (0.08, 0.0008)]
star_masses = [(40.0, 16.0), (16.0, 2.1), (2.1, 1.4), (1.4, 1.04), (1.04, 0.8), (0.8, 0.45), (0.45, 0.1)]
star_radii = [(22.0, 8.0), (7.0, 1.8), (1.8, 1.4), (1.4, 1.15), (1.15, 0.96), (0.96, 0.7), (0.7, 0.25)]
star_prob = [0.0, 0.05, 0.35, 0.6, 0.8, 0.97, 1.0]
num_stars = round(gaussian(*stars))
if num_stars < 1.0:
num_stars = 1.0
num_stars = int(num_stars)
num_planets = round(gaussian(*planets))
if num_planets < 0.0:
num_planets = 0.0
num_planets = int(num_planets)
stars = []
x = y = z = 0.0
while len(stars) < num_stars:
r = random()
for i, e in enumerate(star_prob):
if r <= e:
si = i
break
r = random()
slumin = (star_lumin[si][1] * 0.8) + (star_masses[si][0] * 1.2 * r)
r2 = r + (random() * 0.3) - 0.15
smass = (star_masses[si][1] * 0.8) + (star_masses[si][0] * 1.2 * r)
r2 = r + (random() * 0.3) - 0.15
sradius = (star_radii[si][1] * 0.8) + (star_radii[si][0] * 1.2 * r)
print "Star lumin %s mass %s class %s" % (slumin, smass, star_spectra[si])
stars += [(slumin, smass, star_spectra[si], sradius, x, y, z)]
# first star is always at 0.0, 0.0, 0.0, next star gets random coords
x = gaussian(0.0, systemsize)
y = gaussian(0.0, systemsize)
if len(stars) > 1:
z = gaussian(0.0, systemsize / 10e3)
planets = []
sys_size = math.sqrt(sys_lumin)
while len(planets) < num_planets:
def dist(x, y, z):
return math.sqrt(x**2 + y**2 + z**2)
x = gaussian(0.0, systemsize * sys_size)
y = gaussian(0.0, systemsize * sys_size)
z = gaussian(0.0, systemsize * sys_size / 10e3)
reject = False
for star in stars:
slumin, smass, sclass, sradius, sx, sy, sz
d = dist(x-sx,y-sy,z-sz)
deathzone = sradius + (sradius * 0.25 * math.sqrt(slumin))
if d < deathzone:
reject = True
if reject:
continue
r = random()
rocky = False
gas_giant = False
if r < 0.6:
# rock planet
rocky = True
massmin = 1.0e19
massmax = 8.8e25
mass = 0.0
while mass < massmin or mass > massmax:
massexp = gaussian(24.5, 1.25)
mass = 10.0 ** massexp
density = 0.0
while density < 0.75:
density = gaussian(4.5, 1.0)
mag = -1.0
while mag < 0.0:
mag = gaussian(density, 3.0)
else:
# gas giant
gas_giant = True
massmin = 1e25
massmax = 4e28
mass = 0.0
while mass < massmin or mass > massmax:
massexp = gaussian(26.5, 1.0)
mass = 10.0 ** massexp
density = 0.0
while density < 0.3:
density = gaussian(1.5, 0.75)
mag = -1.0
while mag < 2.0:
mag = gaussian(7.0, 3.0)
volume = mass / density / 1e12
radius = pow(volume * 0.75 / math.pi, 1.0/3.0)
circsurfarea = math.pi * radius**2
surfarea = circsurfarea * 4.0
insol = 0.0
for star in stars:
slumin, smass, sclass, sradius, sx, sy, sz
d = dist(x-sx,y-sy,z-sz)
spower = 3.8e26 * slumin
distsphere = 4.0 * math.pi * d**2
insol += (circsurfarea / distsphere) * spower
areainsol = insol / surfarea / 1e6
if areainsol > 75000:
continue
atmo = None
while atmo == None or atmo > density:
if gas_giant:
atmo_exp = randrange(-3, 0.5)
else:
atmo_exp = randrange(-9.2, -0.9)
atmo = 10 ** atmo_exp
albedo = -1.0
while albedo < 0.0 or albedo > 1.0:
albedo = gaussian(0.35, 0.05)
emissivity = 1.0 - albedo
emissivity = randrange(0.00001, 1.0)
emiss_diff = (((10.0 - pow(atmo, 1.0/4.0)) / 10.0) ** 100.0) / 2.0
emissivity_min = emissivity - (emissivity * emiss_diff)
emissivity_max = emissivity + (emissivity * emiss_diff)
avg_temp = pow(((1.0 - albedo) * areainsol) / (4.0 * emissivity * 5.67e-8), 1.0/4.0)
min_temp = pow(((1.0 - albedo) * areainsol) / (4.0 * emissivity_min * 5.67e-8), 1.0/4.0)
max_temp = pow(((1.0 - albedo) * areainsol) / (4.0 * emissivity_max * 5.67e-8), 1.0/4.0)