Das deutsche Python-Forum

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 

############



########################################