Liebe Sternfreunde, mithilfe der Wellenoptik berechnet dieses Programm die Airy-Scheibe. So wie abgedruckt, ist das Instrument obstruktionsfrei. Durch eine kleine Modifikation kann man die Wirkung der zentralen Obstruktion studieren.
Code:
# Programm Airy1.py
# Stern-Abbildung eines Newton-Teleskops
# Nur Obstruktion, keine Aberrationen eingebaut
from math import *
from turtle import *
from random import *
from plot import *
def sq(x): # Quadrat
return x*x
# Vorgabewerte
l=0.55e-6 # lambda 555nm
f=0.650 # Brennweite
A=0.105 # Apertur
O=0.0 # Obstruktion, z.B. 0.2
# Bestimmung der Rahmenbedingungen für die Rechnung - Vorläufig eine reine Festlegung
w = 4.5 * 2.44 * l * f / A; # die erforderliche Bildgröße
S=25 # Anzahl der Quellpunkte in jede Koordinatenrichtung
p=50 # Anzahl der Bildpunkte in jede Koordinatenrichtung
# das Ergebnisfeld für die Punkte: zwei Phasenlagen und das Ergebnis in einem 3D-Feld
I = [[[0 for i in range(p)] for j in range(p)] for k in range(3)]
G = A * A / 4 / f; # Pfeiltiefe
print(p,":")
for i in range(p): # jeder Bildpunkt X
print(i, end=" ", flush=True)
X = -w / 2 + i * w / (p-1)
for j in range (p): # jeder Bildpunkt Z
Z = -w / 2 + j * w / (p-1)
for m in range(S): # Quellkoordinate x
x = -A / 2.0 + m * A / S
for n in range(S): # Quellkoordinate z
z = -A / 2.0 + n * A / S
r = sqrt(x * x + z * z)
phi = atan2(z, x)
if (r <= A / 2): # nur der Kreis
if (r >= O * A / 2): # zentrale Obstruktion berücksichtigen
y = sq(r) / 4 / f
g = G - y
Y = f - y
Q = sqrt(sq(X) + sq(f) + sq(Z))
q = sqrt(sq(x - X) + sq(Y) + sq(z - Z))
s = (G + Q) - (g + q)
for o in range(2): # beide Phasenwinkel
phasenwinkel = o * pi / 2
a = cos(2 * pi * s / l + phasenwinkel) / sq(q)
I[o][i][j] = I[o][i][j] + a
# beide Phasen summieren nach I[2]
M = 0
for i in range(p):
for j in range(p):
for o in range(2):
I[o][i][j] = sq(I[o][i][j])
a = I[o][i][j]
I[2][i][j] += a
if (I[2][i][j] > M):
M = I[2][i][j]
# Normalisierung und Ausgabe
tracer(0,0)
bgcolor("black")
for i in range(p):
for j in range(p):
a = I[2][i][j] / M
# die Potenz 0,7 verbessert den Kontrast dunkler Strukturen:
PLOT(i-p/2, j-p/2, 2, a**0.7)
hideturtle();update()
ende = input("Fertig.")