Galaxies

from sage.all import Tachyon
from math import sin,cos,tan,log,sqrt
from random import random

pi = 3.14159

# How long the arms should be
th_max = 3.75*pi

A = 20
N = 8
B = 0.04

star_size_max = 0.02

if(N < th_max/2):
    print("N has to be greater than th_max/2")

# Defines how the star spread changes along the arms
def star_spread(dist):
    return (cos(pi/2*dist)+0.8)**1*2
# Defines how the vertical star spread changes along the radius
def star_vert_spread(rad):
    return  (cos(pi*rad)+1)**2/4
# Defines how the number of stars changes along the arms
def star_count(dist):
    return 200*(1-th/th_max)**2

# Galactic arm equation
def r(theta):
    return -A/log(B*tan(theta/(2*N)))

def X(theta):
    return cos(theta)*r(theta)
def Y(theta):
    return sin(theta)*r(theta)

t = Tachyon(xres=1000,yres=1000, camera_center=(17,0,0))
t.light((4,3,2), 100.2, (1,1,1))
t.texture('star', color=(0,0,0), ambient = .4)

# Sample the path every 0.5 degree
th_delta = pi/360
old_px = 0
old_py = 0
old_sx = 0
old_sy = 0
path_delta = 0.1

# Move along the path until we've moved a distance of path_delta
th = th_delta
while th < int(th_max):
    # How far have we moved in units of path_delta?
    # I.e. how many steps do we need to make this loop
    path_x = X(th)
    path_y = Y(th)
    dx = path_x-old_px
    dy = path_y-old_py
    delta = sqrt(dx**2 + dy**2)/path_delta

    # Normalize to a distance of path_delta
    dx = dx/delta
    dy = dy/delta
    print(str(th/th_max*100)+"% done")

    # Always make at least one star cluster
    delta = max(int(delta), 1)

    for step in range(0, delta):
        step_x = old_px + dx*step
        step_y = old_py + dy*step

        # Ensure that the distance between drawing clusters is always >= path_delta
        if (step_x - old_sx)**2 + (step_y - old_sy)**2 < path_delta**2:
            continue
        old_sx = step_x
        old_sy = step_y

        # How many stars, and how spread out should they be?
        count   = int(star_count(th/th_max))
        max_rad = star_spread(th/th_max)

        for n in range(0, count):

            # Move the star a random distance in a random direction
            # Square the distance to make smaller radii more likely
            alpha = random()*2*pi
            beta  = (random()-0.5)*pi
            displacement = random()**2*max_rad

            # Move away from the path
            x = step_x + cos(alpha)*cos(beta)*displacement
            y = step_y + sin(alpha)*cos(beta)*displacement
            z = 0      + sin(beta)*star_vert_spread(r(th)/r(th_max))*displacement

            # Render two spheres, one on each arm
            t.sphere((x,y,z),random()*star_size_max, 'star')
            t.sphere((-x,-y,-z),random()*star_size_max, 'star')
    old_px = path_x
    old_py = path_y

    th += th_delta

t.save()