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)