# Programm Zernicke1.py
# Optische Fehler – Die Zernicke-Polynome
# rechnet lange, v.a. bei starkem Defokus
from math import *
from turtle import *
from random import *

def Plot(x,y, thickness, col): # einen Punkt setzen
    color((col, col, col))
    penup()
    goto(x,y)
    pendown()
    dot(thickness)
    hideturtle()

def plot(x,y):
    Plot(x,y,3, 'black')
    

def Line(x0,y0,x1,y1, thickness, col):
    color(col)
    width(thickness)
    penup(); goto(x0,y0)
    dx=x1-x0; dy=y1-y0
    pendown()
    if (dx<1e-10):
        setheading(90)
        forward(dy)
    else:
        setheading(180*atan(dy/dx)/pi)
        forward(sqrt(dx*dx+dy*dy))
    hideturtle()

def line(x0,y0,x1,y1):
    Line(x0,y0,x1,y1, 1, 'black')

def greenline(x0,y0,x1,y1):
    Line(x0,y0,x1,y1, 1, 'lightgreen')

# Koordinatensystem
def initKoor(x0,x1,stepX, y0,y1,stepY):
    tracer(0,0) 
    setworldcoordinates(x0-0.01*(x1-x0),y0-0.01*(y1-y0),x1+0.01*(x1-x0),y1+0.01*(y1-y0))
    line(x0,y0,x0,y1)
    line(x0,y0,x1,y0)
    x=x0;y=y0
    while (y<=y1):
        line(x,y,x+0.01*(x1-x0),y)
        write(y,font=('Arial',10))
        if y>y0+1e-10:
            greenline(x+0.01*(x1-x0),y,x1,y)
        y=y+stepY
    x=x0;y=y0
    while (x<=x1):
        line(x,y,x,y+0.01*(y1-y0))
        #line(x,y+0.01*(y1-y0),x,y1)
        write(x,font=('Arial',10))
        if x>x0+1e-10:
            greenline(x,y+0.01*(y1-y0),x,y1)
        x=x+stepX
    update()

def mycircle(r, col):
    color((col, col, col))
    penup()
    goto (0,-r)
    begin_fill()
    circle(r)
    end_fill()
    
def sq(x): # Quadrat
    return x*x


# Main

# Vorgabewerte
l=0.55e-6 # lambda 555nm
f=650 # Brennweite
A=0.105 # Apertur
O=0.0 # Obstruktion, z.B. 0.2
# die Aberrationen
D=-2 # Defokus
k=1 # Kugelgestaltsfehler
K=0 # Koma
a=0 # Astigmatismus

# Bestimmung der Rahmenbedingungen für die Rechnung - abängig vom Defokus

# die erforderliche Bildgröße
if abs(D)<1:
    w = 4.5 * 2.44 * l * f / A
else:
    w = 4.5 * D * 2.44 * l * f / A
# Anzahl der Quellpunkte
S=25
if (D > 2):
   S = 0.5*S  * D;
# Anzahl der Bildpunkte
p=50
if (D > 4):
   p = p + int((D - 4) * p / 5);

# 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 n in range(3)]

G = A * A / 4 / f; # Sagitta
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 = G + sqrt(sq(X) + sq(f) + sq(Z)) # Q ist für einen Quellpunkt unveränderlich und kann hoch
                        q = g + sqrt(sq(x - X) + sq(Y) + sq(z - Z))
                        R = r / (A / 2); # Zernicke - R
                        # alle Aberration gehen als +- in die Rechnung ein, deshlab müssen
                        # die Abweichungen halbiert werden um ptv zu erhalten
                        # Defokus
                        Q = Q + D * l * (2 * 0 * 0 - 1) / 2; 
                        q = q + D * l * (2 * R * R - 1) / 2;
                        # sph Aberration
                        Q = Q + k * l * (6 * 0 * 0 * 0 * 0 - 6 * 0 * 0 + 1) / 2; 
                        q = q + k * l * (6 * R * R * R * R - 6 * R * R + 1) / 2;
                        # Koma
                        Q = Q + K * l * (3 * 0 * 0 * 0 - 2 * 0) * sin(phi) / 2; 
                        q = q + K * l * (3 * R * R * R - 2 * R) * sin(phi) / 2;
                        # Astigmatismus
                        Q = Q + a * l * 0 * 0 * cos(2 * phi) / 2;
                        q = q + a * l * R * R * cos(2 * phi) / 2;
                        # gesamter Laufzeitunterschied
                        s = Q - 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
        Plot(i-p/2, j-p/2, 2, a**0.7) # die Potenz 0,7 verbessert den Kontrast dunkler Strukturen
        # Plot(i-p/2,j-p/2,2,a)

hideturtle();update();done()

Dateianhänge:

z.jpg [ 64.41 KiB | 3746 mal betrachtet ]