3D-Animation mit VPython

Antworten
asterix
Beiträge: 228
Registriert: 23.05.2018, 08:24

31.07.2019, 11:41

Wenn man vektororientierte Simulationen schreibt, dann kommt man schnell auf die Idee, den Ablauf der Simulation selbst graphisch darzustellen. Wenn die Simulation fertig ist, dann sollte es doch ziemlich einfach sein, die ganzen Vektoren mit dem originalen Code auch als Grafik auf den Bildschirm zu bringen. Fehler in der Vektorrechnung müssten dann sofort sichtbar werden. Um das einmal auszutesten, war der erste Versuch den Ablauf der Schleifensimulation über javafx darzustellen. In javafx ist eine recht einfache 3D-Darstellung realisierbar. Der Haken kam dann beim Versuch zwei ineinander geschachtelte Schleifen zu animieren ans Licht. Die einfachste Möglichkeit die ich gefunden habe, war dazu beide Schleifen in separaten Animationsthreads laufen zu lassen. Annimationsthreads werden in javafx über die Bildwiederholung getaktet. Die laufen aber parallel ab und müssen synchronisiert werden, damit die Reihenfolge in der Berechnung stimmt. Das geht bei normalen Threads recht einfach. Bei Annimationsthreads leider nicht. Es funktioniert zwar, aber da kommt ein für eine Fehlerüberprüfung untauglich kompliziertes Programm heraus, das selbst noch einmal ziemlich fehleranfällig ist. Also habe ich einmal nachgeschaut, was es denn sonst noch so gibt. Bei VPython bin ich dann fündig geworden. Installiert habe ich die open-source Anaconda Distribution Anaconda 2019.07 for Windows Installer mit der Python 3.7 Version. Da gab es dann dummerweise auch einen Haken, der mich einen Tag unnötige Rumprobiererei gekostet hat. Python-Programme haben prima funktioniert und alle VPython-Programme sind mit der Fehlermeldung " KeyError: 'default/interpreter/dedicated' " abgebrochen Für einen Python-Neuling, ziemlich nichtssagend. Nach einem Tag rumprobieren bin ich dann drauf gekommen, dass VPython nur richtig läuft, wenn Anaconda Navigator mit Adminrechten gestartet wird. Es hat sich aber gelohnt. Animationen mit VPhyton zu realisieren ist wirklich sehr einfach.

Die Installationsvorgehensweise ist hier recht anschaulich beschrieben: https://www.datacamp.com/community/tuto ... da-windows

Wichtig ist, dass die PATH-Variable nach der Installation richtig gesetzt ist. Wenn noch kein Python auf dem Rechner installiert ist, dann geht das am Einfachsten über das obere Häkchen im Installationsmenü:

Abb_1.gif
Abb_1.gif (45.21 KiB) 274 mal betrachtet

Das ist u.A. wichtig, damit der Rechner weiß, wo die Programme conda, python und pip zu finden sind.

Ob alles gut gelaufen ist kann man dann testen:

pip --version
conda --version
python --version

Abb_2.gif
Abb_2.gif (16.73 KiB) 274 mal betrachtet

Als nächstes wird Anaconda Navigator mit Adminrechten (rechte Maustaste, als Administrator ausführen) gestartet:

Abb_3.gif
Abb_3.gif (60.78 KiB) 274 mal betrachtet

Beim Erststart wird noch einiges vom Programm erledigt.

Jetzt fehlt noch das Paket VPython.
Pakete werden mit pip installiert.

Die Anweisung "pip install VPhyton" in einer Eingabeaufforderung erledigt das.

Abb_4.gif
Abb_4.gif (15.08 KiB) 274 mal betrachtet

Aus dem Navigator heraus kann man dann Spyder starten.

Abb_5.gif
Abb_5.gif (59.86 KiB) 274 mal betrachtet

Jetzt kann man endlich eine erste Animation testen.
Einfach den Code in das linke Fenster kopieren und "Ausführen" im Menü auswählen.

Code: Alles auswählen

from vpython import *
scene = canvas() #in Jupyter-Notebooks und -Labors erforderlich, damit Programme problemlos erneut ausgeführt werden können


scene.title = "\nSimulationstest\nRechte Maustaste -> drehen\nMittlere Maustaste -> zoomen\n"
scene.width = 640
scene.height = 400
scene.background = color.gray(0.7)


scene.range = 5                                                             #Kamara in der Entfernung aufstellen, in der ein Bereich von +-5 abgedeckt wird
scene.center = vector(0,0,0)                                                #Koordinatenursprung in die Mitte
scene.forward = vector(1,2,0.1)                                             #Blickrichtung
scene.camera.rotate(angle=pi/1.6, axis=vector(0,1,0), origin=vector(0,0,0)) #Kamera so drehen, dass X-Achse nach rechts zeigt
beleuchtung=distant_light(direction=vector(0,-10,0),color=color.gray(0.8))

#Koordinatenkreuz
x_achse = arrow(pos=vector(0,0,0), axis=vector(5,0,0), color=color.yellow,shaftwidth=0.1)
y_achse = arrow(pos=vector(0,0,0), axis=vector(0,5,0), color=color.green,shaftwidth=0.1)
z_achse = arrow(pos=vector(0,0,0), axis=vector(0,0,5), color=color.blue,shaftwidth=0.1)

#Verschachtelte Rotation
rotor = arrow(pos=vector(0,0,0),axis=vector(3,0,0),color=color.gray(0.8),radius=0.2)
schleife=ring(pos=rotor.axis, axis=vector(1,0,0), radius=0.5, thickness=0.1)
subrotor = arrow(pos=rotor.axis,axis=vector(0,1,0),color=color.orange,radius=0.1)

angle = 0
da=0.1
da1=0.1
while True:
  rate(50) #maximal 50 Durchläufe/sec für die Schleife 
  rotor.rotate(angle=da, axis=vector(0,1,0), origin=rotor.pos)    #Rotor f. Rotation um die y-Achse
  schleife.rotate(angle=da, axis=vector(0,1,0), origin=rotor.pos) #Schleife f. Rotation um die y-Achse
  subrotor.rotate(angle=da, axis=vector(0,1,0), origin=rotor.pos) #Rotor f. Rotation um die rotor-Achse
  angle1 = 0
  while angle1 < 2*pi:
    rate(50) #maximal 50 Durchläufe/sec für die Schleife 
    subrotor.rotate(angle=da1, axis=rotor.axis, origin=rotor.pos)
    angle1 += da1
  angle += da

Die kleine Animation startet dann in einem Browser:

Abb_6.gif
Abb_6.gif (32.75 KiB) 274 mal betrachtet

Eine knappe Übersicht zu den 3d-Basisobjekten findet man hier: https://www.glowscript.org/docs/VPython ... tives.html

Beispielprogramme die mit der Version laufen findet man hier: https://mybinder.org/v2/gh/BruceSherwoo ... ndex.ipynb


Dass man VPhyton nur mit Admin-Rechten benutzen kann, ist nicht gerade das gelbe vom Ei. Aber vielleicht hat ja noch einer eine Idee dazu.
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
asterix
Beiträge: 228
Registriert: 23.05.2018, 08:24

02.08.2019, 21:02

Noch ein Nachtrag.

Wenn es nach der Erstinstallation auch mit Admin-Rechten nicht funktioniert, dann ist noch ein Upgrate erforderlich:       pip install vpython --upgrade

Es funktioniert auch aus normalen Kennungen heraus, wenn "Spider (Anakonda3)" im Windows-Menü, anstatt
Spider aus dem Anaconda-Navigator, gestartet wird.

 
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
asterix
Beiträge: 228
Registriert: 23.05.2018, 08:24

06.08.2019, 15:51

Bevor die Vektoren zum Algorythmus animiert werden können, muss ersteinmal der Algorithmus selbst nach Python konvertiert werden.  Die konvertierung der Vektorrechnung selbst kann man sich sparen, weil die in VPython bereis implementiert ist. Beachten muss man nur, dass Vektoren dort im Urzeigersinn gdreht werden. Um keine zur Java-Version Y-achsenmäßig gespiegelte Darstellung zu bekommen muss man mit dem negativen Winkel drehen. Den Progammteil zur Plot-Ausgabe kann man sich auch sparen. Da stellt Python bereits etwas zur Verfügung. Beim Recherchieren bin ich da auch bei Java auf ein ziemilich gutes Paket gestossen. Aber dazu später in einem anderen Thread. Der Haken bei Python ist die Geschwindigkeit. Python ist eine Interpretersprache. Rechenschleifenintensive Algorithmen laufen da etwa 20 mal langsamer ab als in Java. Ein Simulationszenario das in Java z.b. 1.1 Sekunden braucht, ist in Python erst nach 19.1 Sekunden fertig. Bei der Animation gibt es dagegen praktisch keinerlei Unterschied. Aber zuest die Simulation, zu der ein Ablauf animiert werden soll.

Der Code ist recht kurz geraten.

Code: Alles auswählen

#S_simu_2s_rot_vpy.py
#Simulation zu https://www.power2world.net/viewtopic.php?p=1019#p1019
from vpython import *
import matplotlib.pyplot as plt
import textwrap as tw
import matplotlib as mpl
import time

muo = 4.0*pi*1.0e-7
steps = 256         #   integration steps Rotor
stepsR = 32         #   integration steps Shleife
r1=0.021            #   outer radius rotor #1
r2=0.030            #   outer radius rotor #2
Spalt=0.004         #      gap between rotors
phi0=-5             #   angle shift between rotors
I1=880              #   current in loop #1
I2=880              #   current in loop #2             (3)
R1=0.005            #   radius loop #1
R2=0.005            #    radius loop #2 muo = 4.0*pi*1.0e-7

#Arrays deklarieren
Mphi1=[ 0 for x in range(steps) ]
Mphi2=[ 0 for x in range(steps) ]
Mphi=[ 0 for x in range(steps) ]
IMphi=[ 0 for x in range(steps) ]

def Skalarmultiplikation(v,f):
  return vector(v.x*f,v.y*f,v.z*f) #

def Flussdichte_B1(Pb1,phi2): #Biot Savart Flussdichte B1 verursacht durch S2
  H=vector(0,0,0)
  dphidl2=2.0*pi/stepsR
  Pdl2=vector(r2,R2,0) #(1) S2 linke Schleife
  Pdl2=rotate(Pdl2,-phi2,vector(0,1,0))
  n=1
  while n <= stepsR:#bis einschl Pdlp_nplus1 == 2*pi
    phidl2=n*dphidl2    # links herum
    Pdl2_nplus1=vector(r2,R2*cos(phidl2),R2*sin(phidl2)) #(1) S2 linke Schleife
    Pdl2_nplus1=rotate(Pdl2_nplus1,-phi2,vector(0,1,0))
    dl2=Pdl2_nplus1-Pdl2 #(2)
    rd1=Pdl2-Pb1  #(12)
    abst=rd1.mag
    dH=Skalarmultiplikation(cross(dl2,rd1),1.0/(abst*abst*abst))#()
    H=H+dH
    Pdl2=Pdl2_nplus1
    n+=1
  return Skalarmultiplikation(H,I2*muo/(4.0*pi))  #(5)

def dphi_rotieren_S1(phi2,phi1,n_): # Rotor-Teildrehung dphi
  Mphi1[n_]=0 #(14)
  Fs=vector(0,0,0)
  dphidl1=2.0*pi/stepsR
  Pdl1=vector(-r1,R1,0) # (3)   S1 rechte Schleife
  Pdl1=rotate(Pdl1,-phi1,vector(0,1,0)) # ()
  n=1
  while n <= stepsR:  #bis einschl Pdls_nplus1 == 2*pi
    phidl1=n*dphidl1      # links herum
    Pdl1_nplus1=vector(-r1,R1*cos(phidl1),R1*sin(phidl1)) #(3)   S1 rechte Schleife
    Pdl1_nplus1=rotate(Pdl1_nplus1,-phi1,vector(0,1,0))  #um Phi rotieren
    dl1=Pdl1_nplus1-Pdl1 #()
    Pb1=vector(r1+r2+Spalt,0,0)+Pdl1 # (12)
    B1=Flussdichte_B1(Pb1,phi2) # (11)
    dFs1=Skalarmultiplikation(cross(dl1,B1),I1) #(13)
    Fs=Fs+dFs1
    Mphi1[n_]+=cross(Pdl1,dFs1).y #(14)
    Pdl1=Pdl1_nplus1
    n+=1
  return Fs

def Flussdichte_B2(Pb2,phi1): #Biot Savart Flussdichte B2 verursacht durch S1
  H=vector(0,0,0)
  dphidl1=2.0*pi/stepsR
  Pdl1=vector(-r1,R1,0) # (3)   S1 rechte Schleife
  Pdl1=rotate(Pdl1,-phi1,vector(0,1,0))
  n=1
  while n <= stepsR:#bis einschl Pdlp_nplus1 == 2*pi
    phidl1=n*dphidl1    # links herum
    Pdl1_nplus1=vector(-r1,R1*cos(phidl1),R1*sin(phidl1)) #(3)   S1 rechte Schleife
    Pdl1_nplus1=rotate(Pdl1_nplus1,-phi1,vector(0,1,0))
    dl1=Pdl1_nplus1-Pdl1 #(4)
    rd2=Pb2+Pdl1   #(6)
    abst=rd2.mag
    dH=Skalarmultiplikation(cross(dl1,rd2),1.0/(abst*abst*abst))#(5)
    H=H+dH
    Pdl1=Pdl1=Pdl1_nplus1
    n+=1
  return Skalarmultiplikation(H,I1*muo/(4.0*pi))  #(5)

def dphi_rotieren_S2(phi2,phi1,n_): # Rotor-Teildrehung dphi
  Mphi2[n_]=0    #(14)
  Fs=vector(0,0,0)
  dphidl2=2.0*pi/stepsR
  Pdl2=vector(r2,R2,0)      # (1) S2 linke Schleife
  Pdl2=rotate(Pdl2,-phi2,vector(0,1,0))
  n=1
  while n <= stepsR: #bis einschl Pdlp_nplus1 == 2*pi
    phidl2=n*dphidl2   # links herum
    Pdl2_nplus1=vector(r2,R2*cos(phidl2),R2*sin(phidl2))  #(1) S2 linke Schleife
    Pdl2_nplus1=rotate(Pdl2_nplus1,-phi2,vector(0,1,0))   #um Phi rotieren
    dl2=Pdl2_nplus1-Pdl2
    Pb2=vector(r1+r2+Spalt,0,0)-Pdl2   # (6)
    B2=Flussdichte_B2(Pb2,phi1)   # (5)
    dFs2=Skalarmultiplikation(cross(dl2,B2),I2)#(9)
    Fs=Fs+dFs2
    Mphi2[n_]+=cross(Pdl2,dFs2).y#(10)
    Pdl2=Pdl2_nplus1
    n+=1
  return Fs

dphi=2.0*pi/steps
n=0
winkel= [0.0 for x in range(steps)]

t0 = time.time()
delta=2
tf = time.time()+delta

while n < steps:
  phi2=n*dphi
  phi1=-n*dphi+phi0/360.0*2.0*pi
  phi1+=pi
  phi2+=pi

  dphi_rotieren_S2(phi2,phi1,n)
  dphi_rotieren_S1(phi2,phi1,n)
  winkel[n]=phi2/(2*pi)*360-360
  Mphi[n]=Mphi1[n]-Mphi2[n]
  if n==0:
    IMphi[n]=Mphi[n]
  else:
    IMphi[n]=IMphi[n-1]+Mphi[n]
  n+=1
  t1 = time.time()
  t  = time.time()
  if t > tf or n==steps:
    txt = 'n={:d} '.format(n)
    print(txt,end = '')
    tf=t+delta;
print(' ')
print('Laufzeit =',time.time()-t0,' sec')



fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 9
fig_size[1] = 5

plt.rcParams["figure.figsize"] = fig_size
plt.ylabel('Drehmoment/Nm')
plt.legend(['Drehmoment'])
plt.xlabel('Winkel/°')
plt.grid(True)
plt.plot(winkel,Mphi,color="blue", linewidth=2.5, linestyle="-")
comment2_txt = 'steps={:d} stepsR={:d} R2={: .1f}mm R1={:.1f}mm r2={:.1f}mm r1={:.1f}mm Spalt={:.1f}mm phi0={:.1f}mm I1={:.1f}A I2={:.1f}A'.format(steps,stepsR,R2*1000.0,R1*1000.0,r2*1000.0,r1*1000.0,Spalt*1000.0,phi0,I1,I2)
fig_txt = tw.fill(tw.dedent(comment2_txt.rstrip() ), width=60)
plt.figtext(0.5, -0.11, fig_txt, horizontalalignment='center',fontsize=12, multialignment='left',bbox=dict(boxstyle="round", facecolor='#D8D8D8',ec="0.5", pad=0.5, alpha=1), fontweight='bold')
plt.show()

plt.ylabel('Summe Drehmomente/Nm°')
plt.xlabel('Winkel/°')
plt.legend(['Summe Drehmomente'])
plt.grid(True)
plt.plot(winkel,IMphi,color="green", linewidth=2.5, linestyle="-")
comment2_txt = 'steps={:d} stepsR={:d} R2={: .1f}mm R1={:.1f}mm r2={:.1f}mm r1={:.1f}mm Spalt={:.1f}mm phi0={:.1f}mm I1={:.1f}A I2={:.1f}A'.format(steps,stepsR,R2*1000.0,R1*1000.0,r2*1000.0,r1*1000.0,Spalt*1000.0,phi0,I1,I2)
fig_txt = tw.fill(tw.dedent(comment2_txt.rstrip() ), width=60)
plt.figtext(0.5, -0.11, fig_txt, horizontalalignment='center',fontsize=12, multialignment='left',bbox=dict(boxstyle="round", facecolor='#D8D8D8',ec="0.5", pad=0.5, alpha=1), fontweight='bold')
plt.show()


Die Ausgabe macht auch einen guten Eindruck. Skaliert wird automatisch.

Abb1.gif
Abb1.gif (13.88 KiB) 242 mal betrachtet

Abb2.gif
Abb2.gif (13.9 KiB) 242 mal betrachtet



S_simu_2s_rot_vpy.zip
(1.91 KiB) 4-mal heruntergeladen
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
Uatu
Beiträge: 360
Registriert: 23.05.2018, 15:53

06.08.2019, 19:54

@asterix:

asterix hat geschrieben:
06.08.2019, 15:51
Der Haken bei Python ist die Geschwindigkeit. Python ist eine Interpretersprache. Rechenschleifenintensive Algorithmen laufen da etwa 20 mal langsamer ab als in Java. Ein Simulationszenario das in Java z.b. 1.1 Sekunden braucht, ist in Python erst nach 19.1 Sekunden fertig.

Das ist natürlich unschön. Ich hätte nicht erwartet, dass der Unterschied so gross ist. Java ist ja auch schon kein Renner.
asterix
Beiträge: 228
Registriert: 23.05.2018, 08:24

07.08.2019, 11:15

@uatu
Vielleich gibt es ja irgendwann einmal für Python einen Compiler der Maschinencode erzeugt. Dann wäre das Problem weg.

Zur Visualisierung der beteiligten Vektoren müssen eigentlich nur noch die Grafikausgaben in die Iterationsschleifen eingefügt werden. In Python besonders einfach, weil die Grafikelemente von Haus aus schon über Verktoren beschrieben werden. Allerdings stellt VPython nur Pfeile mit eckigem Schaft zur Verfügung. Eine Klasse für Ausgabe von Pfeilen mit rundem Schaft lässt sich aber leicht realiieren.

Code: Alles auswählen

#g_2s_rot_vpy.py
from vpython import *
import time


muo = 4.0*pi*1.0e-7
steps = 40          #   integration steps Rotor
stepsR = 8          #   integration steps Shleife
r1=0.021            #   outer radius rotor #1
r2=0.030            #   outer radius rotor #2
Spalt=0.004         #  	gap between rotors
phi0=-5             #   angle shift between rotors
I1=880              #   current in loop #1
I2=880              #   current in loop #2 			(3)
R1=0.005            #   radius loop #1
R2=0.005            #	radius loop #2 muo = 4.0*pi*1.0e-7

scene = canvas()
scene.width = 1100
scene.height = 500
scene.background = color.gray(0.7)

kposx=0.0194786
kposy=-0.07*4
kposz=0.0265079
scene.autoscale = False
scene.camera.pos=vector(kposx,kposy,kposz);
scene.camera.axis=vector(kposx/100,-kposy,-kposz);
scene.range=r1+r2
beleuchtung=distant_light(direction=scene.camera.pos,color=color.gray(0.8))


Mphi1=[ 0 for x in range(steps) ]
Mphi2=[ 0 for x in range(steps) ]
Mphi=[ 0 for x in range(steps) ]
IMphi=[ 0 for x in range(steps) ]


def DummyFun(c):
  return

reset=False
def ResetFun(c):
  global reset
  reset=True
  return

stop=True
def StopFun(c):
  global stop
  if stop:
    stop=False
    c.text='Stop'
  else:
    stop=True
    c.text='Start'
  return

grate=5
def setspeed(s):
    global grate
    grate = s.value

#scene.caption('\n')
scene.append_to_caption('\n')
scene.append_to_caption('\n')
sp=button(bind=ResetFun, text=' Reset ')
scene.append_to_caption('      ')
st=button(bind=StopFun, text='Start')
scene.append_to_caption('\n')
scene.append_to_caption('\n')
scene.append_to_caption("Speed:")
slider(min=1, max=50, value=grate, length=350, bind=setspeed)
scene.append_to_caption('\n')
bt_Pdl2 = radio(bind=DummyFun, checked=True, text='Pdl2 ')
scene.append_to_caption('\t\t')
bt_Pdl1 = radio(bind=DummyFun, checked=True, text='Pdl1 ')
scene.append_to_caption('\n')
bt_dl2  = radio(bind=DummyFun, checked=True, text=' dl2 ')
scene.append_to_caption('\t\t')
bt_dl1  = radio(bind=DummyFun, checked=True, text=' dl1 ')
scene.append_to_caption('\n')
bt_rd2  = radio(bind=DummyFun, checked=True, text=' rd2 ')
scene.append_to_caption('\t\t')
bt_rd1  = radio(bind=DummyFun, checked=True, text=' rd1 ')
scene.append_to_caption('\n')




vradius=0.0002
"""
def rotateY(v, a):
  return rotate(v,angle=-a,axis=vector(0,1,0))
"""
def Skalarmultiplikation(v,f):
  return vector(v.x*f,v.y*f,v.z*f)

def Flussdichte_B1(Pb1,phi2): #Biot Savart Flussdichte B1 verursacht durch S2
  gPdl2=None
  gdl2=None
  grd1=None
  gdl2_list = []
  H=vector(0,0,0)
  dphidl2=2.0*pi/stepsR
  Pdl2=vector(r2,R2,0) #(1) S2 linke Schleife
  Pdl2=rotate(Pdl2,-phi2,vector(0,1,0))
  n=1
  while n <= stepsR and not reset:#bis einschl Pdlp_nplus1 == 2*pi
    while stop and not reset: continue
    gPdl2=gdraw(gPdl2,vector(0,0,0),Pdl2,color.red,vradius,bt_Pdl2.checked)
    phidl2=n*dphidl2    # links herum
    Pdl2_nplus1=vector(r2,R2*cos(phidl2),R2*sin(phidl2)) #(1) S2 linke Schleife
    Pdl2_nplus1=rotate(Pdl2_nplus1,-phi2,vector(0,1,0))
    dl2=Pdl2_nplus1-Pdl2 #(2)
    gdl2 =gadd(gdl2,gdl2_list,Pdl2,dl2,color.red,vradius,bt_dl2.checked)
    rd1=Pdl2-Pb1  #(12)
    grd1 = gdraw(grd1,Pb1,rd1,color.blue,vradius,bt_rd1.checked)
    abst=rd1.mag
    dH=Skalarmultiplikation(cross(dl2,rd1),1.0/(abst*abst*abst))#()
    H=H+dH
    Pdl2=Pdl2_nplus1
    rate(grate)
    n+=1
  gdel(gPdl2)
  glist_del(gdl2_list)
  gdel(grd1)
  return Skalarmultiplikation(H,I2*muo/(4.0*pi))  #(5)

def dphi_rotieren_S1(phi2,phi1,n_): # Rotor-Teildrehung dphi
  gPdl1=None
  gdl1=None
  grd2=None
  gdl1_list = []
  Mphi1[n_]=0 #(14)
  Fs=vector(0,0,0)
  dphidl1=2.0*pi/stepsR
  Pdl1=vector(-r1,R1,0) # (3)   S1 rechte Schleife
  Pdl1=rotate(Pdl1,-phi1,vector(0,1,0)) # ()
  n=1
  while n <= stepsR and not reset:  #bis einschl Pdls_nplus1 == 2*pi
    while stop and not reset: continue
    gPdl1=gdraw(gPdl1,vector(r1+r2+Spalt,0,0),Pdl1,color.red,vradius,bt_Pdl1.checked)
    phidl1=n*dphidl1      # links herum
    Pdl1_nplus1=vector(-r1,R1*cos(phidl1),R1*sin(phidl1)) #(3)   S1 rechte Schleife
    Pdl1_nplus1=rotate(Pdl1_nplus1,-phi1,vector(0,1,0))  #um Phi rotieren
    dl1=Pdl1_nplus1-Pdl1 #()
    gdl1 =gadd(gdl1,gdl1_list,vector(r1+r2+Spalt,0,0)+Pdl1,dl1,color.red,vradius,bt_dl1.checked)
    Pb1=vector(r1+r2+Spalt,0,0)+Pdl1 # (12)
    B1=Flussdichte_B1(Pb1,phi2) # (11)
    dFs1=Skalarmultiplikation(cross(dl1,B1),I1) #(13)
    Fs=Fs+dFs1
    Mphi1[n_]+=cross(Pdl1,dFs1).y #(14)
    Pdl1=Pdl1_nplus1
    rate(grate)
    n+=1
  gdel(gPdl1)
  glist_del(gdl1_list)
  return Fs

def Flussdichte_B2(Pb2,phi1): #Biot Savart Flussdichte B2 verursacht durch S1
  gPdl1=None
  gdl1=None
  grd2=None
  gdl1_list = []
  H=vector(0,0,0)
  dphidl1=2.0*pi/stepsR
  Pdl1=vector(-r1,R1,0) # (3)   S1 rechte Schleife
  Pdl1=rotate(Pdl1,-phi1,vector(0,1,0))
  n=1
  while n <= stepsR and not reset:#bis einschl Pdlp_nplus1 == 2*pi
    while stop and not reset: continue
    gPdl1=gdraw(gPdl1,vector(r1+r2+Spalt,0,0),Pdl1,color.red,vradius,bt_Pdl1.checked)
    phidl1=n*dphidl1    # links herum
    Pdl1_nplus1=vector(-r1,R1*cos(phidl1),R1*sin(phidl1)) #(3)   S1 rechte Schleife
    Pdl1_nplus1=rotate(Pdl1_nplus1,-phi1,vector(0,1,0))
    dl1=Pdl1_nplus1-Pdl1 #(4)
    gdl1 =gadd(gdl1,gdl1_list,vector(r1+r2+Spalt,0,0)+Pdl1,dl1,color.red,vradius,bt_dl1.checked)
    rd2=Pb2+Pdl1   #(6)
    grd2 = gdraw(grd2,vector(r1+r2+Spalt,0,0)-Pb2,rd2,color.blue,vradius,bt_rd2.checked)
    abst=rd2.mag
    dH=Skalarmultiplikation(cross(dl1,rd2),1.0/(abst*abst*abst))#(5)
    H=H+dH
    Pdl1=Pdl1=Pdl1_nplus1
    rate(grate)
    n+=1
  gdel(gPdl1)
  glist_del(gdl1_list)
  gdel(grd2)
  return Skalarmultiplikation(H,I1*muo/(4.0*pi))  #(5)

def dphi_rotieren_S2(phi2,phi1,n_): # Rotor-Teildrehung dphi
  gFs=None
  gPdl2=None
  gdl2=None
  gdl2_list = []
  Mphi2[n_]=0    #(14)
  Fs=vector(0,0,0)
  dphidl2=2.0*pi/stepsR
  Pdl2=vector(r2,R2,0)      # (1) S2 linke Schleife
  Pdl2=rotate(Pdl2,-phi2,vector(0,1,0))
  n=1
  while n <= stepsR and not reset: #bis einschl Pdlp_nplus1 == 2*pi
    while stop and not reset: continue
    gPdl2=gdraw(gPdl2,vector(0,0,0),Pdl2,color.red,vradius,bt_Pdl2.checked)
    phidl2=n*dphidl2   # links herum
    Pdl2_nplus1=vector(r2,R2*cos(phidl2),R2*sin(phidl2))  #(1) S2 linke Schleife
    Pdl2_nplus1=rotate(Pdl2_nplus1,-phi2,vector(0,1,0))   #um Phi rotieren
    dl2=Pdl2_nplus1-Pdl2
    gdl2 =gadd(gdl2,gdl2_list,Pdl2,dl2,color.red,vradius,bt_dl2.checked)
    Pb2=vector(r1+r2+Spalt,0,0)-Pdl2   # (6)
    B2=Flussdichte_B2(Pb2,phi1)   # (5)
    dFs2=Skalarmultiplikation(cross(dl2,B2),I2)#(9)
    Fs=Fs+dFs2
    Mphi2[n_]+=cross(Pdl2,dFs2).y#(10)
    Pdl2=Pdl2_nplus1
    rate(grate)
    n+=1
  gdel(gPdl2)
  glist_del(gdl2_list)
  return Fs

def glist_del(lst):
  for element in lst:
    gdel(element)


class r_arrow:
  co=None
  cy=None
  slen=1;
  radius=1
  def __init__(self,p,ax,c,r):
    slen=10*r;
    l=ax.mag-slen;
    if l > 0:
      self.cy = cylinder(pos=p,axis=ax,color=c,radius=r)
      self.cy.length=self.cy.length-slen
      self.co=cone(pos=p+self.cy.axis,axis=ax,color=c,radius=2*r)
      self.co.length=slen
    else:
      self.co=None
      self.co=cone(pos=p,axis=ax,color=c,radius=2*r)
      self.co.length=slen
    self.radius=r

  def set_rc(self,r,c):
    self.cy.radius=r
    self.co.radius=2*r
    self.cy.color=c
    self.co.color=c

  def __del__(self):
    self.cy.visible = False
    del self.cy
    self.co.visible = False
    del self.co

def gadd(old,lst,p,ax,c,sw,on):
  if on:
    if old !=None:
      old.set_rc(old.radius/2,color.gray(0.4))
      old.color=color.gray(0.5)
    gobj = r_arrow(p,ax,c,sw)
    lst.append(gobj)
    return gobj

def gdraw(gobj,p,ax,c,sw,on):
  if gobj !=None:
    gdel(gobj)
  if on:
    gobj = r_arrow(p,ax,c,sw)
  return gobj

def gdel(gobj):
  if gobj !=None:
    gobj.visible = False
    del gobj
  return None


xyz=r1/2
rko=0.0005
x_achse = r_arrow(vector(0,0,0),vector(xyz,0,0),color.yellow,rko)
y_achse = r_arrow(vector(0,0,0),vector(0,xyz,0),color.green,rko)
z_achse = r_arrow(vector(0,0,0),vector(0,0,xyz),color.blue,rko)


camax=None;
campos=None;
dphi=2.0*pi/steps
n=0
while True:
  rate(grate)
  if reset:
    n=0
    reset=False
  phi2=n*dphi
  phi1=-n*dphi+phi0/360.0*2.0*pi
  if phi2 < pi/2:
    scene.camera.axis=vector(kposx/100,-kposy,-kposz);
    scene.camera.pos=vector(kposx+0.01,kposy,kposz+0.013)
    scene.range=r1*0.9
  elif phi2 < 2*pi-pi/2:
    scene.camera.axis=vector(kposx/100,-kposy,-kposz);
    scene.camera.pos=vector(kposx,kposy,kposz)
    scene.range=r1*1.7
  else:
    scene.camera.axis=vector(kposx/100,-kposy,-kposz);
    scene.camera.pos=vector(kposx,kposy,kposz-0.013)
    scene.range=r1*0.9
  if not reset:
    dphi_rotieren_S2(phi2,phi1,n)
  if not reset:
    dphi_rotieren_S1(phi2,phi1,n)
  scene.autoscale = False
  n+=1
  if n >= steps:  n=0


Bei der Gelegenheit habe ich mich noch umgeschaut wo man seine Videos noch so ablegen kann. Bei YouTube hatte ich schon zweimal Datentotalverlust, weil die mir ohne sinnvolle Begründung und ohne jegliche nachvollziehare Gründe die Kennung gesperrt haben. Nur meine Daten haben sie behalten. Einschließlich der Telefonnummer . Bei dailymotion kann man auch über eine e-mail Adresse eine uploadfähige Kennung bekommen. Ohne Telefonnummer. Da ist eigentlich ziemlich optimal für Videos die in ein Forum eingebettet werden sollen. Das habe ich dann auch gleich ausprobiert.





g_2s_rot_vpy.zip
(2.46 KiB) 4-mal heruntergeladen
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
Uatu
Beiträge: 360
Registriert: 23.05.2018, 15:53

08.08.2019, 00:10

@asterix:

asterix hat geschrieben:
07.08.2019, 11:15
Vielleich gibt es ja irgendwann einmal für Python einen Compiler der Maschinencode erzeugt. Dann wäre das Problem weg.

Ich weiss nicht, ob Du Dich mit der Thematik schon mal auseinandergesetzt hast, aber das ist aufgrund einiger konzeptioneller Eigenschaften von Python (insbesondere der dynamischen Typisierung) problematisch. Ich habe einige Erfahrung im Compilerbau (auch wenn das schon lange zurückliegt) und kenne das Problem von der Implementierungsseite. Bei statisch typisierten Programmiersprachen (u.a. den meisten C- und Pascal-basierten Sprachen) reicht i.d.R. eine einfache lokale Analyse des Sourcecodes, um effizienten Maschinencode zu erzeugen, während das bei dynamisch typisierten Sprachen i.d.R. erheblich schwieriger ist, und meist trotzdem zu wesentlich ineffizienterem Maschinencode führt. Es gibt einige Compiler für Python, aber die Berichte darüber, die ich bisher gelesen habe, waren wenig überzeugend. PyPy könnte evtl. einen Versuch wert sein, wobei das ein Just-in-Time-Compiler, also ein Zwischending ist, und es ein paar Einschränkungen gegenüber normalem Python gibt.

Auf eine andere Möglichkeit der Geschwindigkeitssteigerung weist die folgende Aussage im offiziellen VPython-Architektur-Dokument hin:

Computationally intensive GlowScript VPython programs run about an order of magnitude faster than VPython 7 programs, because they are compiled to (fast) JavaScript (but there is no access to Python modules).

(VPython Architecture (pdf))

Bei einer normalen VPython-Installation ist VPython eine Python-Bibliothek, die für die Grafikausgabe eine als GlowScript bezeichnete JavaScript-Bibliothek im Webbrowser ansteuert. Man kann das Ganze aber auch vollständig im Webbrowser ablaufen lassen, wofür es eine Reihe verschiedener Möglichkeiten gibt. Das hat den o.g. Geschwindigkeitsvorteil, und man braucht auch keine lokale Python-Installation. Nachteil ist, dass man keine Python-Bibliotheken, sondern nur den vorgegebenen Sprachumfang plus JavaScript-Bibliotheken verwenden kann.

Für die direkte Nutzung von GlowScript gibt es drei Möglichkeiten:

1. Man kann von JavaScript aus die GlowScript-Bibliothek ansprechen. Man programmiert also im Prinzip in JavaScript. Das ist die einfachste Lösung, die Syntax ist allerdings weniger elegant als bei VPython.

2. Man kann in RapydScript programmieren, was im wesentlichen der Python-Syntax entspricht, und in JavaScript umgesetzt wird. Von dieser Variante halte ich nichts, weil ich bezweifle, dass RapydScript langfristig verfügbar sein wird.

3. Man kann in VPython programmieren, das mit einem Verfahren ähnlich wie bei Möglichkeit 2 in JavaScript umgesetzt wird. Obwohl dieses Verfahren aktuell auf RapydScript beruht, gelten meine Bedenken zu Möglichkeit 2 in diesem Fall nicht. Zum einen ist das Umsetzungsverfahren in diesem Fall eine Black-Box, d.h. es könnte irgendwann auch völlig anders als mit RapydScript implementiert werden, ohne dass der Sourcecode angepasst werden muss. Zum anderen sollte der Sourcecode (zumindest sehr weitgehend) auch unter einer normalen (Python-basierten) VPython-Installation lauffähig sein, d.h. man hat ggf. eine Ausweichmöglichkeit.

Zum Vergleich der Syntax hier als Beispiel (von der offiziellen GlowScript-Website) der Sourcecode für eine einfache astronomische Animation in JavaScript, RapydScript und VPython.

Alle drei Varianten kann man sowohl cloudbasiert als auch offline nutzen. Das wird auf der offiziellen VPython-Website unter Using VPython without installing any software erläutert.

Aus theoretischer Sicht erscheinen mir die Lösungen 1 und 3 etwa gleichwertig (persönliche Programmiersprachen-Präferenzen mal aussen vor gelassen). Beide sind mit Vor- und Nachteilen verbunden, die sich aber m.E. ungefähr gegeneinander aufheben. In der Praxis kann das natürlich anders aussehen.
asterix
Beiträge: 228
Registriert: 23.05.2018, 08:24

09.08.2019, 10:32

Ich habe noch eine Möglichkeit gefunden mit der man die Geschwindigkeit steigern kann. Es gibt noch die Möglichkeit Extensions in C zu schreiben und die dann als Modul zu importieren. Ein funktionierendes Beispiel ist hier beschrieben. Darauf aufbauend habe ich die Methode einmal getestet.

Zuerst wird das c/c++ Programm geschrieben. Im Beispiel wird das Modul "winkelfunktionen" erzeugt. Zu einem mit Winkeln gefüllten Array berechnet das Modul cosinus/sinus zu jedem Winkel und gibt diese in einem anderen Array zurück.

Zuerst wird das c-Programm (winkelfunkt.c) geschrieben:

winkelfunkt.c

Code: Alles auswählen

#include <Python.h>
#include <numpy/arrayobject.h>
#include <math.h>

static PyObject* sin_func_np(PyObject* self, PyObject* args)
  {
  PyArrayObject *in_array;
  PyObject      *out_array;
  NpyIter *n_in;
  NpyIter *n_out;
  if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &in_array))
      return NULL;
  out_array = PyArray_NewLikeArray(in_array, NPY_ANYORDER, NULL, 0);
  if (out_array == NULL)
      return NULL;
  n_in = NpyIter_New(in_array, NPY_ITER_READONLY, NPY_KEEPORDER,NPY_NO_CASTING, NULL);
  if(n_in == NULL)
    {
    Py_XDECREF(out_array);
    return NULL;
    }
  n_out = NpyIter_New((PyArrayObject *)out_array, NPY_ITER_READWRITE,NPY_KEEPORDER, NPY_NO_CASTING, NULL);
  if (n_out == NULL)
    {
    NpyIter_Deallocate(n_in);
      {
      Py_XDECREF(out_array);
      return NULL;
      }
    }
  NpyIter_IterNextFunc *in_next = NpyIter_GetIterNext(n_in, NULL);
  NpyIter_IterNextFunc *out_next = NpyIter_GetIterNext(n_out, NULL);
  if (in_next == NULL || out_next == NULL)
    {
    NpyIter_Deallocate(n_in);
    NpyIter_Deallocate(n_out);
      {
      Py_XDECREF(out_array);
      return NULL;
      }
    }
  double ** in_dataptr = (double **) NpyIter_GetDataPtrArray(n_in);
  double ** out_dataptr = (double **) NpyIter_GetDataPtrArray(n_out);
  do
    {
    **out_dataptr = sin(**in_dataptr);
    } while(in_next(n_in) && out_next(n_out));

  NpyIter_Deallocate(n_in);
  NpyIter_Deallocate(n_out);
  Py_INCREF(out_array);
  return out_array;
  }

static PyObject* cos_func_np(PyObject* self, PyObject* args)
  {
  PyArrayObject *in_array;
  PyObject      *out_array;
  NpyIter *n_in;
  NpyIter *n_out;
  if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &in_array))
      return NULL;
  out_array = PyArray_NewLikeArray(in_array, NPY_ANYORDER, NULL, 0);
  if (out_array == NULL)
      return NULL;
  n_in = NpyIter_New(in_array, NPY_ITER_READONLY, NPY_KEEPORDER,NPY_NO_CASTING, NULL);
  if(n_in == NULL)
    {
    Py_XDECREF(out_array);
    return NULL;
    }
  n_out = NpyIter_New((PyArrayObject *)out_array, NPY_ITER_READWRITE,NPY_KEEPORDER, NPY_NO_CASTING, NULL);
  if (n_out == NULL)
    {
    NpyIter_Deallocate(n_in);
      {
      Py_XDECREF(out_array);
      return NULL;
      }
    }
  NpyIter_IterNextFunc *in_next = NpyIter_GetIterNext(n_in, NULL);
  NpyIter_IterNextFunc *out_next = NpyIter_GetIterNext(n_out, NULL);
  if (in_next == NULL || out_next == NULL)
    {
    NpyIter_Deallocate(n_in);
    NpyIter_Deallocate(n_out);
      {
      Py_XDECREF(out_array);
      return NULL;
      }
    }
  double ** in_dataptr = (double **) NpyIter_GetDataPtrArray(n_in);
  double ** out_dataptr = (double **) NpyIter_GetDataPtrArray(n_out);
  do
    {
    **out_dataptr = cos(**in_dataptr);
    } while(in_next(n_in) && out_next(n_out));

  NpyIter_Deallocate(n_in);
  NpyIter_Deallocate(n_out);
  Py_INCREF(out_array);
  return out_array;
  }

/*  Funktionen die exportiert werden in Array schreiben */
static PyMethodDef winkelfunktionen[] =
  {
     {"cosinus", cos_func_np, METH_VARARGS,"cos in numpy array"},
     {"sinus", sin_func_np, METH_VARARGS,"sin in numpy array"},
     {NULL, NULL, 0, NULL}
  };


#if PY_MAJOR_VERSION >= 3  /* Python version 3*/
static struct PyModuleDef cModPyDem =
  {
  PyModuleDef_HEAD_INIT,"winkelfunktionen", "Doku",
  -1,
  winkelfunktionen
  };
PyMODINIT_FUNC PyInit_winkelfunktionen(void)  //PyInit_[modulname]
  {
  PyObject *m;
  m = PyModule_Create(&cModPyDem);
  if(m==NULL) return NULL;
  import_array();
  if (PyErr_Occurred()) return NULL;
    return m;
  }

#else /* Python version 2 */
PyMODINIT_FUNC initcos_module_np(void)
  {
  PyObject *m;
  module = Py_InitModule("winkelfunktionen", winkelfunktionen);
  if(m==NULL)
    return;
  import_array();
  return;
  }

#endif

Mit Python lassen sich c-Progamme in ein Modulformat kompilieren. Dazu muss eine setup-Datei angefertigt werden:

setup.py

Code: Alles auswählen

from distutils.core import setup, Extension
import numpy

winkelfunktionen = Extension('winkelfunktionen', sources=['winkelfunkt.c'],include_dirs=[numpy.get_include()])
setup(ext_modules=[winkelfunktionen])

#Ausführen mit:
#python setup.py build_ext --inplace


Neben dem Modulnamen kann dort noch angegeben werden wo eventuell im c-Programm aufgerufene .h-Dateien zu finden sind. numpy.get_include() gibt hier den Pfad zu den .h-Dateien von numpy zurück. Die Datei wird anschließend, vom Ordner mit der c-Ouelle aus, folgendermassen aufgerufen:

python setup.py build_ext --inplace


Python ruft dann einen c-Compiler auf und erzeugt das Modul winkelfunktionen.cp37-win_amd64.pyd. Das kann man dann in winkelfunktionen.pyd umbenennen
Abb1.gif
Abb1.gif (35.36 KiB) 150 mal betrachtet

Mit einem kleinen Python-Programm kann man das neue Modul testen:

c_in_py_test.py

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt
import os, sys
from os.path import dirname
sys.path.append(dirname(__file__))  #Pfad vom aktuellen .py wo auch winkelfunktionen.pyp liegt
import winkelfunktionen

x = np.arange(0, 2 * np.pi, 0.1)
y = winkelfunktionen.cosinus(x)
plt.plot(x, y)
plt.show()

x = np.arange(0, 2 * np.pi, 0.1)
y = winkelfunktionen.sinus(x)
plt.plot(x, y)
plt.show()



Die Ausgabe sollte so aussehen:
 
Abb2.gif
Abb2.gif (10.45 KiB) 150 mal betrachtet



ProgQuellen.zip
(7.43 KiB) 1-mal heruntergeladen

 
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
asterix
Beiträge: 228
Registriert: 23.05.2018, 08:24

10.08.2019, 12:27

Ich habe den Algorithmus jetzt einmal als cpp-Extension implementiert und mit der original Python-Version und der Java-Version verglichen.


Test_py.gif
Test_py.gif (4.56 KiB) 122 mal betrachtet

Dazu ergaben sich die folgenden Laufzeiten:
Python ................... 341.9 sec
Python mit cpp-Extension ... 3.1 sec
Java ...................... 12.4 sec


Was rechenaufwendige Chart-Anwendungen angeht, da ist es wahrscheinlich am Einfachsten, wenn man die zunächst in Java realisiert. NetBeans macht das ziemlich bequem. Es gibt dazu auch recht einfach handhabbare Chart-Erweiterungen . Eine ggf. anschließende Umsetzung in eine cpp-Extension ist auch ziemlich einfach. Man kann da den größten Teil des Java-codes 1:1 übernehmen. Und Python unterstützt die Einbindung eigentlich auch recht gut.


S_simu_c_py.zip
(23 KiB) 2-mal heruntergeladen
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
Antworten
  • Information