Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Models orbits using keplers equations, then tries to makes sounds from them.
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()
"""
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment