Planet simulation real time and size

    ##################################################################################
    #
    # Kinematiikan harjoituksia 2
    # Ohjelmoi kuu-maa-animaatio tasaisella ympyräliikkeellä xy-tasossa
    # Kohdat 1., 2. ja 3.
    # Alkuarvot http://fi.wikipedia.org/wiki/Kuu)
    # Määritä alkuarvot siten, että
    # - Maa sijaitsee origossa paikoillaan
    # - Kuu kiertää maata ympyräradalla
    # Etsi sopiva arvot ympyräliikkeen liikeyhtälöön (t:n kertoimeksi) niin,
    # että kuun kiertoajaksi tulee näytöllä n. 28 vrk (tark Wikipediassa)
    # Simulaation kannattaa tietenkin esittää nopeutettuna muuttamalla arvoa rate()
    # Koordinaatiston ajan kuluminen ilmoitetaan tekstikehyksellä grafiikkaikkunassa
    #
    # Jori Niskanen
    #
    # 19.9.2013
    #
    ##################################################################################

from visual import *
from math import *

#Ohjelman alustus tehtävän mukainen/laaja
Mode = 2
while Mode != 1 and Mode != 2:
    Mode = input("Missä tilassa sovellus käynnistetään?\n1) Tehtävä     2) Laaja\nValinta: ")

#Infonäyttö
sceneinfo = display(title="Kinematiikka 2. Aurinko, Maa ja Kuu: Info",x=0, y=0, width=1000, height=220, userzoom = False, userspin = False)
sceneinfo.range = (10,10,10)
timeLabel1 = label(pos=(6.1,-0.5,0))
timeLabel2 = label(pos=(6.1,0.5,0))
timeLabel3 = label(pos=(6.1,1.5,0))
timeLabel4 = label(pos=(-7.5,0.5,0))
timeLabel5 = label(pos=(-7.5,1.3,0))
timeLabel7 = label(pos=(-1.4,-1.5,0))
timeLabel6 = label(pos=(-1.4,1.3,0))
timeLabel8 = label(pos=(-1.4,0,0))
timeLabel5.text = "Ohjaa simulaation nopeutta nuolilla\nNuolet 0.1  ctrl+nuolet 1"
timeLabel6.text = "Nuolet sivuille vaihtavat seurattavaa kappaletta\nEnd pysayttaa seurattavan kappaleen liikkeen"
timeLabel8.text = "Suurilla dt-arvoilla kiertoajat saattavat olla viallisia\nKaytetaan dt-arvoa laskentamarginaalina"
timeLabel7.text = "Maa liike= True"
if Mode == 1: timeLabel7.text = "Maa liike= False"

#Simulaationäyttö
scene = display(title="Kinematiikka 2. Aurinko, Maa ja Kuu: Simulaatio",x=0, y=220, width=1000, height=800,forward=(-1,-0.08,-0.1))



    # Alla netistä otettu hyödyllinen funktio make_time_string(t), joka ilmoittaa parametreina annetun
    # ajan t sopivina aikayksikköinä
def make_time_string(t):
    "Accept a number of seconds, return a relative string."
    if t < 60: return "%02i seconds"%t
    mins,secs = divmod(t, 60)
    if mins < 60: return "%02i minutes %02i seconds"%(mins,secs)
    hours, mins = divmod(mins,60)
    if hours < 24: return "%02i hours %02i minutes %02i seconds"%(hours,mins, secs)
    days,hours = divmod(hours,24)
    return "%02i days %02i hours %02i minutes %02i seconds"%(days,hours,mins,secs)



'Simulate' # Simuloidaanko kappaleen liike (True/False)
'Orbiting' # Mitä kappaletta kyseinen kappale kiertää
'Semi-major_axis' # Kiertoradan keskietäisyys kierrettävästä esineestä. http://fi.wikipedia.org/wiki/Radan_isoakselin_puolikas
        #Isoakselin puolikas on ellipsin pisimmän halkaisijan puolikas. Tätä rataelementtiä käytetään tyypillisesti kuvaamaan ellipsiradalla kiertävän kappaleen etäisyyttä keskuskappaleesta, joka sijaitsee ellipsin yhdessä polttopisteessä.
        #Isoakselin puolikasta nimitetään usein "keskietäisyydeksi". Nimitys on siinä mielessä perusteltu, että isoakselin puolikas on kappaleen suurimman ja pienimmän etäisyyden keskiarvo
'Eccentricity' # Eksentrisyys. Lasketaan SMA:n avulla kiertoradan elliptisyys
'Orbital_period' #Kiertoaika kierrettävän kappaleen ympäri sekunneissa
'Inclination' # Kiertoradan inklinaatio muodossa radians(inc) tai ilman inklinaatiota 0. Ecliptic inclination
'Radius' # Objektin säde
'Sidereal_rotation' #0=ei pyörimistä. Muuten pyörähdysaika x sekuntia
'Axial_tilt' # akselin kaltevuus muodossa radians(tilt°)
'Material' # Ei toimi vielä! Määritettävä "käsin" objektin alustuksessa
'Ascend_Longitude' #Longitude of the Ascending Node
'Ascend_Longitude_Progress_W' # 2*pi/s Ascend Longituden Progressio s. 2*pi/-s Regressiossa -s
'Argument_of_periapsis'
'Argument_of_periapsis_Progress_W'
#lightyear = 9.4607e15 #valovuosi
#sirius_et = 8.60*lightyear
#solarRadius = 695500000 #      1.711
#rSirius = 1.711*solarRadius*10000




Objects = ['Aurinko', 'Merkurius', 'Venus', 'Maa', 'Kuu', 'Mars', 'Phobos', 'Deimos'] #Simuloitavat kappaleet. Kierrettävän kappaleen oltava ennen kiertävää kappaletta.
Item = {}
Data = {}
Data['Aurinko'] = {'Simulate': False,'Orbiting':vector(0,0,0),'Semi-major_axis':0,'Eccentricity':0,'Orbital_period':0,'Inclination':0,'Radius':696342000,'Sidereal_rotation':0,'Axial_tilt':0,'Material':'emissive'}
Data['Maa'] = {'Simulate':True,'Orbiting':'Aurinko','Semi-major_axis':1496.0e8,'Eccentricity':0.01671123,'Orbital_period':31558150,'Inclination':0,'Radius':6378100,'Sidereal_rotation':86164.1,'Axial_tilt':radians(23.4392810833),'Material':'BlueMarble','Ascend_Longitude':radians(348.73936)}
Data['Kuu'] = {'Simulate':True,'Orbiting':'Maa','Semi-major_axis':384400000,'Eccentricity':0.0549,'Orbital_period':2360622,'Inclination':radians(5.1454),'Radius':1737100,'Sidereal_rotation':2360584.68,'Axial_tilt':radians(1.5424),'Material':'wood','Ascend_Longitude':radians(0),'Ascend_Longitude_Progress_W':2*pi/-586958823,'Argument_of_periapsis_Progress_W':2*pi/279278795} #Ascend Longitude regress

Data['Merkurius'] = {'Simulate':True,'Orbiting':'Aurinko','Semi-major_axis':57909100000,'Eccentricity':0.205630,'Orbital_period':7.60053e6,'Inclination':radians(7.005),'Radius':2439700,'Sidereal_rotation':5.06701e6,'Axial_tilt':radians(2.11),'Material':'BlueMarble','Ascend_Longitude':radians(48.331)}
Data['Venus'] = {'Simulate':True,'Orbiting':'Aurinko','Semi-major_axis':108208000000,'Eccentricity':0.0067,'Orbital_period':1.94142e7,'Inclination':radians(3.39458),'Radius':6051800,'Sidereal_rotation':-2.09968e7,'Axial_tilt':radians(177.36),'Material':'BlueMarble','Ascend_Longitude':radians(76.678)}
Data['Mars'] = {'Simulate':True,'Orbiting':'Aurinko','Semi-major_axis':227939100000,'Eccentricity':0.093315,'Orbital_period':5.93543e7,'Inclination':radians(1.850),'Radius':3396200,'Sidereal_rotation':88642.7,'Axial_tilt':radians(25.19),'Material':'BlueMarble','Ascend_Longitude':radians(49.562)}

Data['Phobos'] = {'Simulate':True,'Orbiting':'Mars','Semi-major_axis':9377200,'Eccentricity':0.0151,'Orbital_period':27553.8,'Inclination':radians(26.04),'Radius':11100,'Sidereal_rotation':27553.8,'Axial_tilt':radians(0),'Material':'BlueMarble','Ascend_Longitude':radians(0)}
Data['Deimos'] = {'Simulate':True,'Orbiting':'Mars','Semi-major_axis':23460000,'Eccentricity':0.0002,'Orbital_period':109074.816,'Inclination':radians(27.58),'Radius':6200,'Sidereal_rotation':109074.816,'Axial_tilt':radians(0),'Material':'BlueMarble','Ascend_Longitude':radians(0)}

Data['Kuu']['Eccentricity'] = 0.5
#Data['Kuu']['Ascend_Longitude_Progress_W'] = Data['Kuu']['Ascend_Longitude_Progress_W']*100
#Data['Kuu']['Argument_of_periapsis_Progress_W'] = 0

Asteroid = {}
Population = {}
i = 0
Asteroid_tot = 0;
while i < Asteroid_tot:
    #Longitude of the ascending node = Ascend_Longitude
    Asteroid[i] = {'Semi-major_axis':random.randint(0, 200000)*1000000+239000000000,'Eccentricity':random.randint(0,1000)/10000,'Orbital_period':0,'Radius':random.randint(1000,339620),'Inclination':radians(random.randint(0,10000)/1000),'Rand':random.randint(0,100000)*10e6,'Ascend_Longitude':radians(random.randint(0,360))}
    Orb = Asteroid[i]['Semi-major_axis'] /227939100000
    Asteroid[i]['Orbital_period'] = 5.93543e7*Asteroid[i]['Semi-major_axis']/227939100000
    spos = rotate(vector(0,0,0) + Asteroid[i]['Semi-major_axis']*vector((1-Asteroid[i]['Eccentricity'])*cos(Asteroid[i]['Rand'])*cos(Asteroid[i]['Inclination']),-(1+Asteroid[i]['Eccentricity'])*cos(Asteroid[i]['Rand'])*sin(Asteroid[i]['Inclination']),sin(Asteroid[i]['Rand'])),Asteroid[i]['Ascend_Longitude'],vector(0,1,0))
    Population[i] = sphere( pos= spos, radius= Asteroid[i]['Radius'], make_trail=True)
    Asteroid[i]['W_orbit'] = -2*pi/Asteroid[i]['Orbital_period']
    #Population[i].trail_object.color=Color("#999999")
    i += 1
    # 1. Ohjelmoi

if Mode == 1:
    Objects = ['Maa', 'Kuu']
    Data['Maa']['Simulate'] = False

    # 2. Ohjelmoi :
    # Kuvaa maa ja kuu palloina (todelliset säteet voi suurentaa visualisoinnin vuoksi
    # aseta kuulle attribuutiksi make_trail= True, jolloin kiertorata piirtyy ikkunaan
else: #Ei tehtävämoodissa tehdään aurinko
    spos = Data['Aurinko']['Orbiting']
    sradius = Data['Aurinko']['Radius']
    smaterial = Data['Aurinko']['Material']
    Item['Aurinko'] = sphere (pos=spos, radius=sradius, material=materials.emissive) #Aurinko origossa. Kokona auringon koko. Materiaalina hohtava valoa varten
    local_light(pos=spos, color=color.white) #Tehdään auringosta valonlähde
    scene.lights = [] #Valo tulee vain auringosta


s = 'Maa'
if (Mode == 2):spos = Item[Data[s]['Orbiting']].pos + Data[s]['Semi-major_axis']*vector((1-Data[s]['Eccentricity'])*cos(0)*cos(Data[s]['Inclination']),-(1+Data[s]['Eccentricity'])*cos(0)*sin(Data[s]['Inclination']),sin(0))
else: spos = vector(0,0,0)
sradius = Data[s]['Radius']
#Item[s] = sphere( pos= rotate(spos,Data[s]['Ascend_Longitude'],vector(0,1,0)), radius= sradius,material=materials.BlueMarble,make_trail=True)

loop_len = len(Objects) #Tarkistetaan kappaleiden lukumäärä muuttujaan käyttöä varten
if Mode == 1: i = 1
else: i = 1
while i < loop_len:
    s = Objects[i]
    spos = Item[Data[s]['Orbiting']].pos + Data[s]['Semi-major_axis']*vector((1-Data[s]['Eccentricity'])*cos(0)*cos(Data[s]['Inclination']),-(1+Data[s]['Eccentricity'])*cos(0)*sin(Data[s]['Inclination']),sin(0))
    sradius = Data[s]['Radius']
    Item[s] = sphere( pos= rotate(spos,Data[s]['Ascend_Longitude'],vector(0,1,0)), radius= sradius,material=materials.earth ,make_trail=True)
    i += 1



Item['Maa'].trail_object.color=color.green  ################
Item['Kuu'].trail_object.color=color.orange #Seurannan väri#

if Mode == 2:
    Item['Mars'].trail_object.color=color.red #Seurannan väri#
    Item['Phobos'].trail_object.color=color.green #Seurannan väri#
    Item['Deimos'].trail_object.color=color.orange #Seurannan väri#





    ##################################################################################
    # simulaatioluuppi
    ##################################################################################

scene.forward = -norm(Item['Maa'].pos)+vector(0,-0.2,0) #Käännetään kuva oikeaan suuntaan
Size_multiplier = 1 #Kappaleiden halkaisijoiden kerroin
t=0L
dt=0.001 # simulaation nopeutta voi muuttaa ohjelmassa
timeLabel4.text = "dt= "+str(dt)

finished=True

 # Oletuksena seurataan maata
stepper = 0 #Asteroid simulations
i = 0
while i < loop_len: #Viimeistellään kappaleet ja data ennen simulointia
    s = Objects[i]
    if ( s == 'Maa' ): cam_follow = i #Seurataan alussa maata
    if (Data[s]['Sidereal_rotation']): #Jos kappaleelle on määritetty pyöriminen, pyöräytetään z-akselin suhteen akselin kallistumisen verran
        Item[s].rotate(angle=Data[s]['Axial_tilt'], axis=vector (0,0,1) , origin=Item[s].pos)
    if (Data[s]['Orbital_period']):
        Data[s]['W_orbit'] = -2*pi/Data[s]['Orbital_period'] #Maan kiertokulman laskeminen
        Data[s]['W_axis'] = 2*pi/Data[s]['Sidereal_rotation']    #Maan pyörimisen kiertokulman laskeminen
        Data[s]['W_axis_vect'] = vector(0,1,0).rotate(angle=Data[s]['Axial_tilt'], axis=vector (0,0,1)) #Lasketaan pyöräytys
    if (Size_multiplier != 1):
        Item[s].radius *= Size_multiplier #Muutetaan kappaleiden kokoa kertoimen mukaan


    i += 1
viiva = curve(pos=[(300000000,0,0),(0,0,0),(0,300000000,0),(0,0,0),(0,0,300000000),(0,0,0)])

#debugball = sphere(pos = (0,0,0), radius= 5000000, color=color.white, material=materials.wood)


while finished :
    rate(1000)
    scene.center = Item[Objects[cam_follow]].pos #Pidetään katse tiukasti planeetassamme

    if dt <> 0:
        # 3. Ohjelmoi kuun kiertoliike tähän
        i = 0
        while i < loop_len: #Käydään kaikki kappaleet läpi
            s = Objects[i]
            if (Data[s]['Simulate']):#Kappaleiden liike.
                if s == 'Kuu':



                    #pos1 = Item[Data[s]['Orbiting']].pos + rotate(vector(ecc1.x*cos(Data[s]['W_orbit']*t)*cos(Data[s]['Inclination']),ecc1.y*-cos(Data[s]['W_orbit']*t)*sin(Data[s]['Inclination']),Data[s]['Semi-major_axis']*sin(Data[s]['W_orbit']*t)),t*Data[s]['Ascend_Longitude_Progress_W'],vector(0,1,0))
                    pos2 =rotate( Item[Data[s]['Orbiting']].pos + Data[s]['Semi-major_axis']*vector((1-Data[s]['Eccentricity'])*cos(radians(90))*cos(Data[s]['Inclination']),(1+Data[s]['Eccentricity'])*-cos(radians(90))*sin(Data[s]['Inclination']),sin(radians(90))),t*Data[s]['Ascend_Longitude_Progress_W'],vector(0,1,0))
                    pos3 =rotate( Item[Data[s]['Orbiting']].pos + Data[s]['Semi-major_axis']*vector((1-Data[s]['Eccentricity'])*cos(radians(0))*cos(Data[s]['Inclination']),(1+Data[s]['Eccentricity'])*-cos(radians(0))*sin(Data[s]['Inclination']),sin(radians(0))),t*Data[s]['Ascend_Longitude_Progress_W'],vector(0,1,0))
                    pos1 =rotate( Item[Data[s]['Orbiting']].pos + Data[s]['Semi-major_axis']*vector((1-Data[s]['Eccentricity'])*cos(Data[s]['W_orbit']*t)*cos(Data[s]['Inclination']),(1+Data[s]['Eccentricity'])*-cos(Data[s]['W_orbit']*t)*sin(Data[s]['Inclination']),sin(Data[s]['W_orbit']*t)),t*Data[s]['Ascend_Longitude_Progress_W'],vector(0,1,0))
                    #debugball.pos = rotate( pos2,t*Data[s]['Ascend_Longitude_Progress_W']-t*Data[s]['Argument_of_periapsis_Progress_W'],norm(cross(pos3, pos2)))#-t*Data[s]['Ascend_Longitude_Progress_W']

                    #if pos1.x < 0: pos3 = -pos1
                    #elif pos1.x >= 0: pos3 = pos1
                    #pos2 = norm(pos2)
                    #viiva.pos = [(0,0,0),-norm(cross(pos3, pos2))*300000000,(0,0,0),rotate(-norm(cross(pos3, pos2))*500000000,radians(90),vector(0,0,1)),(0,0,0),rotate(-norm(cross(pos3, pos2))*500000000,radians(90),vector(1,0,0))]
                    Item[s].pos = rotate( pos1,t*Data[s]['Ascend_Longitude_Progress_W']-t*Data[s]['Argument_of_periapsis_Progress_W'],norm(cross(pos3, pos2)))#-t*Data[s]['Ascend_Longitude_Progress_W']

                else:
                    Item[s].pos =rotate( Item[Data[s]['Orbiting']].pos + Data[s]['Semi-major_axis']*vector((1-Data[s]['Eccentricity'])*cos(Data[s]['W_orbit']*t)*cos(Data[s]['Inclination']),(1+Data[s]['Eccentricity'])*-cos(Data[s]['W_orbit']*t)*sin(Data[s]['Inclination']),sin(Data[s]['W_orbit']*t)),Data[s]['Ascend_Longitude'],vector(0,1,0))

                #Item[s].pos =                  #Muutetaan kappaleen sijainti
                #Item[Data[s]['Orbiting']].pos+ #Kierrettävän kappaleen sijainti +
                #Data[s]['Semi-major_axis']*    #SMA * vektori
                #x=(1-Data[s]['Eccentricity'])*cos(Data[s]['W_orbit']*t)*cos(Data[s]['Inclination'])
                #y=(1+Data[s]['Eccentricity'])*-cos(Data[s]['W_orbit']*t)*sin(Data[s]['Inclination'])
                #z=sin(Data[s]['W_orbit']*t)



            if (Data[s]['Sidereal_rotation']): #Pyöritetään kappaleita lasketun vektorin ympäri
                Item[s].rotate(angle=Data[s]['W_axis']*dt, axis=Data[s]['W_axis_vect'] , origin=Item[s].pos)

            i += 1 #
        stepper += 1
        if stepper == 10 :
            i = 0
            while i < Asteroid_tot:
                Population[i].pos = rotate(vector(0,0,0) + Asteroid[i]['Semi-major_axis']*vector((1-Asteroid[i]['Eccentricity'])*cos(Asteroid[i]['W_orbit']*t+Asteroid[i]['Rand'])*cos(Asteroid[i]['Inclination']),(1+Asteroid[i]['Eccentricity'])*-cos(Asteroid[i]['W_orbit']*t+Asteroid[i]['Rand'])*sin(Asteroid[i]['Inclination']),sin(Asteroid[i]['W_orbit']*t+Asteroid[i]['Rand'])),Asteroid[i]['Ascend_Longitude'],vector(0,1,0))
                i += 1
            stepper = 0



    # tekstikehyksien päivitys
    text = "Kulunut aika: " + make_time_string(t)
    timeLabel3.text=text
    if ((t) % 2360622 < dt): #Tulostetaan kiertoaika sen täyttyessä. Suurella dt-arvolla voi poiketa oikeasta ajasta. Otetaan jakojäännös kuluneesta ajasta kiertokulmalla.
        text="Kuunkierto: " + make_time_string(t)
        print text
        timeLabel2.text=text
    if ((t) % 31558150 < dt):
        text="Maankierto: " + make_time_string(t)
        print text
        timeLabel1.text=text

    t+=dt

    if scene.kb.keys or sceneinfo.kb.keys: # onko painettu?
        if scene.kb.keys: key = scene.kb.getkey() #Kummassa ikkunassa painettu
        if sceneinfo.kb.keys: key = sceneinfo.kb.getkey()
        if key == "down":
            dt -= 0.1
        elif key == "ctrl+down":
            dt -= 1
        elif key == "shift+down":
            dt = 0
        elif key == "up":
            dt += 0.1
        elif key == "ctrl+up":
            dt += 1
        elif key == "shift+up":
            dt += 10
        elif key == "left":
            cam_follow -= 1
            print degrees(diff_angle(Item['Kuu'].pos, vector(0,1,0)))%90
        elif key == "right":
            cam_follow += 1
            print cross(pos1,Item['Maa'].pos)
            print Item['Maa'].pos
            print norm(pos1)
            print norm(cross(pos3, pos2))
            Item['Kuu'].pos += norm(cross(pos3, pos2))*500000
            Item['Kuu'].trail_object.pos = ()
        if key == "end": #Hallitaan kappaleen liikkumista end-painikkeella. Oletuksena simuloimattomien kappaleiden simulointi kaataa ohjelman
            if Data[Objects[cam_follow]]['Simulate']: Data[Objects[cam_follow]]['Simulate'] = False
            else: Data[Objects[cam_follow]]['Simulate'] = True

        if (cam_follow < 0): cam_follow = 0 #Pidetään seurattavan kappaleen raja-arvot hallinnassa
        elif (cam_follow >= loop_len): cam_follow -= 1
        timeLabel7.text = Objects[cam_follow]+" liike= " + str(Data[Objects[cam_follow]]['Simulate'])


        timeLabel4.text = "dt= "+str(dt)


            #if Item['Maa'].material == materials.earth: Item['Maa'].material = materials.BlueMarble
            #else: Item['Maa'].material = materials.earth