174 lines
3.8 KiB
Python
174 lines
3.8 KiB
Python
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)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|