# 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.")
 

Dateianhänge:

abb2a.jpg [ 38.8 KiB | 24485 mal betrachtet ]