A&A Title Image Startseite | Programm: Nebensonnen | Algorithmen, Lektion, Simulationen | Unsere Vorläufer | Kontakt, Datenschutz

VdS-Journal 68 : Simulationen zu Nebensonnen

Das in Heft 68 abgedruckte Programm

Es berechnet die Strahlablenkung für 2500 zufällig einfallende Strahlen:


from turtle import *
from random import *
from math import *

def plot(x,y): # einen Punkt setzen
    penup()
    goto(x,y)
    pendown()
    dot(4)
    hideturtle()

# Hauptprogramm: Halo #5-67
n=1.310 # Brechungsindex
l=[0]*900 # Leere Liste für die Ergebnisse aller 0,1°
stichprobe=2500
w60=pi*60/180 # 60 Grad
w90=pi*90/180 # 90 Grad
for i in range(stichprobe): 
    zz=random(); 
    a1=zz*w90 # Winkel von 0-90°
    a2=asin(sin(a1)/n)
    a3=w60-a2
    h=n*sin(a3) # Totalreflexion?
    if (abs(h)<0.99999): # nein!
        a4=asin(h)
        d=a1+a4-w60
        index=int(10*180*d/pi) # eine Spalte je 1/10 Grad
        l[index]=l[index]+1
tracer(0,0) # Bildschrimaktualisierung aus
for i in range (601): # 600 Ergebnisse ab 20°
    plot(i-300, l[i])
    if (int(i/10)==i/10):
        plot(i-300, -5)
    if (int(i/100)==i/100):
        plot(i-300, -5)
        plot(i-300, -6)
        plot(i-300, -7)
    #print(i, l[i])
    #print(i/10.0, 1000*l[i]/stichprobe) # Promille
update() # Bildschirm anzeigen
name = input("Fertig?")

Zu Heft 68: Version mit modifizierten Zufallszahlen

Dieses Programm berücksichtigt, dass Eintrittswinkel um so weniger wahrscheinlich werden, je mehr sich dieser Winkel 90° nähert. Hierzu werden ungleich verteilte Zufallszahlen benutzt.

 

from turtle import *
from random import *
from math import *

def plot(x,y): # einen Punkt setzen
    penup()
    goto(x,y)
    pendown()
    dot(4)
    hideturtle()

# Hauptprogramm
n=1.310 # Brechungsindex
l=[0]*900 # Leere Liste für die Ergebnisse
stichprobe=2500
w60=pi*60/180 # 60 Grad
w90=pi*90/180 # 90 Grad
for i in range(stichprobe):
    zz=random();  # Minimum zweier ZZ ...
    zz2=random(); # ... abnehmende Häufigkeit
    if (zz2<zz):
        zz=zz2
    a1=zz*w90 # Winkel von 0-90°
    a2=asin(sin(a1)/n)
    a3=w60-a2
    h=n*sin(a3) # Totalreflexion?
    if (abs(h)<0.99999): # nein!
        a4=asin(h)
        d=a1-(w60-a4)
        index=int(10*180*d/pi) # eine Spalte je 1/10 Grad
        l[index]=l[index]+1
tracer(0,0) # Bildschirmaktualisierung aus
for i in range (601): # 600 
    plot(i-300, l[i])
    # Koordinatensystem
    if (int(i/10)==i/10):
        plot(i-300, -5)
    if (int(i/100)==i/100):
        plot(i-300, -5)
        plot(i-300, -6)
        plot(i-300, -7)
update() # Bildschirm anzeigen
name = input("Fertig?")