Models orbits using keplers equations, then tries to makes sounds from them.

kepler.py

import math
#import sys
#import sunau
#import array
"""
ecc=0.5
P=1
p=1
if len(sys.argv)>1:
    ecc=float(sys.argv[1])
if len(sys.argv)>2:
    P=float(sys.argv[2])
if len(sys.argv)>3:
    p=float(sys.argv[3])

orbit = []

if ecc>0.999:
   ecc=0.999
"""


def m_anomaly(t, P):
    return 6.28319 * t / P


def e_anomaly(M, ecc):
    """Calculate the eccentric anomaly

    This uses newtons method to solve the Kepler equation"""
    if ecc > 0.8:
        E = 3.14
    else:
        E = M

    new_E = E - (E - ecc * math.sin(E) - M) / (1 - ecc * math.cos(E))
    while math.fabs(E - new_E) > 0.01:
        E = new_E
        new_E = E - (E - ecc * math.sin(E) - M) / (1 - ecc * math.cos(E))
    return new_E


def true_anomaly(E, ecc):
    return 2 * math.atan(math.sqrt((1 + ecc) / (1 - ecc)) * math.tan(E / 2))


def distance(th, ecc, p):
    return p / (1 + ecc * math.cos(th))

"""
def amplify(value):
    maxvalue=max(value)
    minvalue=min(value)
    return int(60000*(value-min)/(max-min)-30000)

step=0.01
t=0
max=0
min=p

while t<=P:
    samp=distance(true_anomaly(e_anomaly(m_anomaly(t,P),ecc),ecc),p)
    if samp<min:
        min=samp
    elif samp>max:
        max=samp
    orbit.append(samp)
    t+=step


orbit = map(amplify,orbit)
print orbit
sampdata=array.array('h', orbit).tostring()


sound=sunau.open("out.au","w")

sound.setnchannels(1)
sound.setframerate(44100)
sound.setsampwidth(2)

for i in range(5000):
    sound.writeframes(sampdata)

sound.close()
"""