habe hier einen code geschrieben zur Planetensimulation.habe mich dabei nach einem matlabcode gehalten, der auf folgender seite zu finden ist:
http://www.am.uni-duesseldorf.de/de/Leh ... eminar.pdf
Code: Alles auswählen
########################################
from numpy import *
from visual import *
n=3 # anzahl der planeten
Q=array([(0,0,0),(2,0.2,0),(0,5,0)]) # ortsMatrix[sonne, erde, mond]
V=array([(0,0,0),(0,0.0008,0.0005),(0,-0.0008,0)]) #geschwindigkeitsMatrix[sonne,erde,mond]
M=array([200,100,50])# massevektor[sonne,erde,mond]
P = arange(n*3).reshape(n,3) # dummy für Impulsmatrix
for i in range(0,n):# i:0,1,2 #berechne Impulsmatrix der planeten
P[i]=V[i]*M[i]
dHdQ=arange(n*3).reshape(n,3) #dummy-hamilton für ortsableitung
dHdP=arange(n*3).reshape(n,3) #dummy-hamilton für impulsableitung
#berechne abgeleitung nach impuls für planeten###############################
for i in range(0,n):
dHdP[i]=P[i]/M[i]
################################################################
#berechne ableitung nach ort für planeten
for i in range(0,n): #ite-zeile
for j in range(0,3):#jte-komponente der ite-zeile
s=0
for k in range(0,n):
if k==i:
s=s
else:
s=s+(M[i]*M[k])*(Q[k][j]-Q[i][j])/(pow(sqrt(dot(Q[k]-Q[i],Q[k]-Q[i])),3))
dHdQ[i][j]=s
################################################################
Q2=arange(n*3).reshape(n,3) #dummy Ort2
P2=arange(n*3).reshape(n,3) #dummy Impuls2
sonne=sphere(pos=Q[0]) #generiere kugel
erde=sphere(pos=Q[1]) #generiere kugel
mond=sphere(pos=Q[2])#generiere kugel
#berechne neue position nach zeitschrit dt
dt = 0.01
while 1:
rate(10)
P2=P+dt*dHdQ#euler-verfahren: neuen impuls berechnen
print(dHdQ)
#berechne abgeleitung nach impuls P2 für planeten
for i in range(0,n):
dHdP[i]=P2[i]/M[i]
#########################################
Q2=Q+dt*dHdP#euler-verfahren:neuen ort berechen
#ausgabe:neue position zuweisen
sonne.pos=Q2[0]
erde.pos=Q2[1]
mond.pos=Q2[2]
#berechne neue ableitung nach ort2 für planeten
for i in range(0,n): #ite-zeile
for j in range(0,3):#jte-komponente der ite-zeile
s=0
for k in range(0,n):
if k==i:
s=s
else:
s=s+(M[i]*M[k])*(Q2[k][j]-Q2[i][j])/(pow(sqrt(dot(Q2[k]-Q2[i],Q2[k]-Q2[i])),3))
dHdQ[i][j]=s
###########################################
P=P2
Q=Q2
############
########################################...das programm startet aber man siehr nur 3 bälle wegfliegen...das kann doch nciht sein???
gruß
nice87day