[PYTHON] Numerische Berechnung der Linsenpunktbildverteilungsfunktion und der MTF-Kurve (Beugungsberechnung)

Ich habe kürzlich die Punktbildverteilungsfunktion (PSF), die häufig im Zusammenhang mit dem Linsendesign auftritt, und die spezifische numerische Berechnungsmethode der MTF-Kurve untersucht und sie zusammengefasst. Ich bin kein Spezialist für Objektivdesign und komme nicht aus der Physikabteilung. Ich bin mir also sicher, dass es einige Fehler gibt, aber ich würde es begrüßen, wenn Sie auf Fehler hinweisen könnten.

Informationen zur Punktbildverteilungsfunktion (PSF)

** Punktstreufunktion --PSF ** ist eine Funktion, die die Lichtintensitätsverteilung auf der Bildebene ausdrückt, wenn eine Punktlichtquelle auf der Bildebene abgebildet wird.

psf.png

Das Bild oben ist ein Beispiel für PSF. Je näher es an Rot liegt, desto höher ist die Lichtintensität. Bei der idealen Bildgebung fokussiert eine Punktlichtquelle auf einen Punkt, sodass PSF eine Funktion ist, die nur an einem Punkt eine Spitze aufweist, wie z. B. die Delta-Funktion. In der Realität wird das Bild jedoch aufgrund der Auswirkungen von Aberration und Beugung nicht auf einen Punkt fokussiert, und es kann ein Bild mit einem bestimmten Grad an Streuung erhalten werden, wie in Fig. 1 gezeigt.

Es gibt zwei Arten von PSF: ** geometrische optische PSF **, bei denen der Beugungseffekt nicht berücksichtigt wird, und ** optische Wellen-PSF **, bei der der Beugungseffekt berücksichtigt wird. Im Allgemeinen scheint PSF wellenoptische PSF zu bedeuten. Das obige Beispiel ist eine wellenoptische PSF, die die Beugung berücksichtigt, wobei ein weiterer Peak um das Zentrum herum wellig erscheint.

Über die MTF-Kurve

Die Erklärung von hier ist leicht zu verstehen, daher werde ich sie zitieren.

MTF (Modulation Transfer Function) ist eine der Skalen zur Bewertung der Objektivleistung und drückt aus, wie genau der Kontrast des Motivs auf der Bildebene als räumliche Frequenzcharakteristik reproduziert werden kann.

MTF wird im Allgemeinen durch einen Graphen dargestellt, der wie folgt als MTF-Kurve bezeichnet wird. Dieses Diagramm stammt auch von hier.

figure_diffraction.gif

Die vertikale Achse repräsentiert den Kontrast und die horizontale Achse repräsentiert die Bildhöhe, dh den Abstand von der Mitte der Bildebene. Die roten und grünen Linien entsprechen bestimmten Raumfrequenzen, und die durchgezogenen und gestrichelten Linien entsprechen der Sagittalrichtung (S) und der Meridionalrichtung (M). Ich verstehe die Bedeutung der sagittalen und meridionalen Richtungen hier nicht wirklich, aber wahrscheinlich die dynamische Richtung $ r $ und die Azimutwinkelrichtung $ \ theta, wenn die Bildebene durch das Polarkoordinatensystem $ (r, \ theta) $ dargestellt wird. Es scheint, dass es $ entspricht.

Zusammenfassend stellt die MTF-Kurve den Kontrast zwischen der dynamischen Richtung und der Azimutwinkelrichtung bei bestimmten räumlichen Frequenzen dar, während die Bildhöhe geändert wird.

Wie PSF hat MTF ** geometrisches optisches MTF ** und ** wellenoptisches MTF **, je nachdem, ob der Einfluss der Lichtbeugung berücksichtigt wird oder nicht.

Numerische Berechnung der PSF

Es gab zwei Arten von PSF: geometrische optische PSF und wellenoptische PSF. Schauen wir uns zunächst an, wie die geometrische optische PSF berechnet wird.

Numerische Berechnung der geometrischen optischen PSF

Die geometrische und optische PSF ist eine Darstellung, wie eine Punktlichtquelle ein Bild auf der Bildebene erzeugt, ohne die Auswirkungen der Beugung zu berücksichtigen. Dies kann relativ einfach mit ** Raytracing ** berechnet werden.

Eine Punktlichtquelle wird an einer geeigneten Position platziert (im Allgemeinen im Unendlichen?), Und eine große Anzahl von Strahlen wird von dort zur Vorderseite der Linse emittiert. Jedes Mal, wenn diese Strahlen auf die Linse treffen, werden sie gebrochen, um die nächste Richtung zu berechnen, und dann wird die Berechnung der Kollisionsposition mit der nächsten Linse wiederholt. Berechnen Sie abschließend die Position, an der der Lichtstrahl die Bildebene schneidet, und addieren Sie die Strahlungshelligkeit des Lichts zu dieser Position.

Durch Wiederholen der späten Schnürung für eine große Anzahl von Strahlen auf diese Weise und geeignetes Interpolieren zwischen den Abtastpunkten kann die Intensitätsverteilung des Lichts auf der Bildebene erhalten werden. Auf diese Weise erhalten Sie eine geometrische optische PSF.

Das Diagramm, das die Position zeigt, an der der Lichtstrahl die Bildebene schneidet, wird als ** Punktdiagramm ** bezeichnet. Die folgende Abbildung ist ein Beispiel.

spot_diagram.png

Beugungsberechnung

Um die optische Wellen-PSF zu berechnen, ist es notwendig, das Licht als Welle auszudrücken und die Intensität des Lichts auf der Bildebene unter Berücksichtigung der Beugung beim Durchgang durch die Linse zu erhalten. Schauen wir uns zunächst die grundlegende Methode zur Berechnung der Beugung an. Schließlich werde ich versuchen, numerische Berechnungen mit Python durchzuführen.

Hier werden die folgenden Beugungsberechnungsmethoden betrachtet.

Wie man sich als Lichtwelle ausdrückt

Um Licht als Welle zu betrachten und in einer mathematischen Formel auszudrücken, wird der Ausdruck mit komplexen Zahlen verwendet.

Für ebene Wellen ist die komplexe Amplitude $ u (\ boldsymbol {r}, t) $ an Position $ \ boldsymbol {r} $, Zeit $ t $

u(\boldsymbol{r}, t) = A\exp[i(\boldsymbol{k}\cdot\boldsymbol{r} - \omega t + \delta)]

Wird ausgedrückt als. Wobei $ A $ die Amplitude ist, $ i $ die imaginäre Einheit ist, $ \ boldsymbol {k} $ der Wellenzahlvektor ist, $ \ omega $ die Winkelfrequenz ist und $ \ delta $ die Anfangsphase ist.

Bei sphärischen Wellen beträgt der Abstand von der Wellenquelle $ r $

u(r, t) = \frac{A}{r}\exp[i(kr - \omega t + \delta)]

Wird ausgedrückt als. Wobei $ k $ die Anzahl der Wellen ist.

In den folgenden Berechnungen haben die Begriffe $ \ omega t $ und $ \ delta $ keine Auswirkung und werden ignoriert.

Beugungsintegration

Betrachten Sie das folgende Koordinatensystem, um die Berechnung der Beugung zu erklären.

diagram-20200411.png

$ \ Sigma $ repräsentiert die Öffnung und $ P (x_p, y_p, 0) $ ist der Punkt über $ \ Sigma $. $ Q (x_q, y_q, z_q) $ ist ein Punkt auf der Bildebene. $ r $ ist die Länge der Linie $ PQ $ und $ \ theta $ ist der Winkel zwischen der Linie $ PQ $ und der Achse $ z $.

Betrachten Sie die Situation, in der eine ebene Welle, die durch die Öffnung $ \ Sigma $ geht, beim Beugen den Punkt $ Q $ erreicht. Nach dem Fresnel-Huihens-Prinzip kann die komplexe Amplitude des Punktes $ Q $ als Überlagerung von Sekundärwellen ausgedrückt werden, die auf $ \ Sigma $ erzeugt werden.

Betrachten Sie zunächst die komplexe Amplitude der quadratischen Welle, die vom Punkt $ P $ den Punkt $ Q $ erreicht. $ \frac{1}{i\lambda}\frac{e^{ikr}}{r}u(x_p, y_p, 0)\cos{\theta} $ Es wird sein. Dabei ist $ \ lambda $ die Wellenlänge des Lichts, $ k $ die Anzahl der Wellen und $ u (x_p, y_p, 0) $ die komplexe Amplitude am Punkt $ P $.

Da die komplexe Amplitude $ u (x_q, y_q, z_q) $ am Punkt $ Q $ die Überlagerung dieser Sekundärwellen über $ \ Sigma $ ist.

\begin{align}
u(x_q, y_q, z_q) &= \frac{1}{i\lambda}\int\int_{\Sigma} u(x, y, 0)\frac{e^{ikr}}{r}\cos{\theta}dxdy \tag{1} \\
&= \frac{z_q}{i\lambda}\int\int_{\Sigma} u(x, y, 0)\frac{e^{ikr}}{r^2}dxdy
\end{align}

Es wird sein. Bei der Formeltransformation haben wir $ \ cos {\ theta} = \ frac {z_q} {r} $ verwendet.

Dieses Integral heißt ** Fresnel-Kirchhoff-Beugungsintegral **. Da die komplexe Amplitude am Punkt $ Q $ durch numerische Integration berechnet werden kann, kann die Intensität auf der Bildebene durch Berechnung des Quadrats des Absolutwerts erhalten werden. Diese Berechnungsmethode ist jedoch nicht praktikabel, da der Berechnungsbetrag $ O (N ^ 4) $ beträgt, wenn die Anzahl der Teilungen zwischen der Öffnungsfläche und der Bildfläche $ N $ beträgt.

Frenel-Näherung

Stellen Sie sich eine Situation vor, in der $ \ theta $ klein und $ \ cos {\ theta} \ ungefähr 1 $ ist (Annäherung nahe der Achse).

Wenn man sich auf Gleichung 1 konzentriert, wird der Teil $ \ cos {\ theta} $ zu 1, und $ \ frac {1} {r} $ kann als $ \ frac {1} {z} $ angenähert werden. Andererseits hat der Term $ e ^ {ikr} $ einen großen Wert von $ k = \ frac {2 \ pi} {\ lambda} $. Wenn Sie also $ r $ durch $ z $ ersetzen, ist die Genauigkeit der Approximation gut. Da ist gar nichts.

Erwägen Sie daher, $ r $ zu einer Reihe zu erweitern und bis zum zweiten Term zu verwenden. Wenn $ r = \ sqrt {(x_q - x) ^ 2 + (y_q - y) ^ 2 + z_q ^ 2} $ auf den zweiten Term erweitert wird

r = z_q + \frac{(x_q - x)^2 + (y_q - y)^2}{2z_q} + \cdots

Durch Einsetzen in die Gleichung für die Beugungsintegration

u(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q} \int\int_{\Sigma} u(x, y, 0)e^{ik\frac{(x_q - x)^2 + (y_q - y)^2}{2z_q}}dxdy \tag{2}

Es wird sein. Diese Näherung wird als ** Frenel-Näherung ** bezeichnet, und Gleichung 2 wird als ** Frenel-Beugungsintegration ** bezeichnet.

Erweitern des Inhalts der $ e $ -Schulter in Gleichung 2

u(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x_p, y_p, 0)e^{-ik\frac{xx_q + yy_q}{z_q}}e^{ik\frac{x^2 + y^2}{2z_q}}dx dy \tag{3}

Es wird sein. In Gleichung 3 ist $ f_x = \ frac {x_q} {\ lambda z_q} $, $ f_y = \ frac {y_q} {\ lambda z_q} $

\begin{align} 
u(x_q, y_q, z_q) &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x, y, 0)e^{-ik\frac{xx_q + yy_q}{z_q}}e^{ik\frac{x^2 + y^2}{2z_q}}dx dy \\ 
&= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x, y, 0)e^{ik\frac{x^2 + y^2}{2z_q}}e^{-2\pi i(f_x x + f_y y)}dx dy \\ 
&= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \mathcal{F}\{u(x, y, 0)e^{ik\frac{x^2 + y^2}{2z_q}}\}_{f_x = \frac{x_q}{\lambda z_q}, f_y = \frac{y_q}{\lambda z_q}} \end{align}

Daher kann es durch die zweidimensionale Fourier-Transformation von $ u (x, y, 0) e ^ {ik \ frac {x ^ 2 + y ^ 2} {2z_q}} $ berechnet werden.

Damit die Fresnel-Näherung gilt, wird $ r $ mit dem dritten Term multipliziert, der eine Reihenerweiterung darstellt, und $ k $ ist erforderlich.

\frac{k}{8z_q^3}[(x_q - x)^2 + (y_q - y)^2]^2 \ll 1

Muss für $ z_q $ erfüllt und neu geschrieben werden

z_q^3 \gg \frac{\pi}{4\lambda}[(x_q - x)^2 + (y_q - y)^2]^2

Es wird sein.

Fraunhofer-Näherung

In Gleichung 3 ist $ \ frac {x ^ 2 + y ^ 2} {2z_q} $ klein genug und liegt nahe an $ e ^ {ik \ frac {x ^ 2 + y ^ 2} {2z_q}} \ append Dann

u(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x, y, 0)e^{-ik\frac{xx_q + yy_q}{z_q}}dxdy \tag{4}

Es wird sein. Diese Näherung wird als ** Fraunhofer-Näherung ** bezeichnet, und Gleichung 4 wird als ** Fraunhofer-Beugungsintegral ** bezeichnet.

In Gleichung 4 ist $ f_x = \ frac {x_q} {\ lambda z_q} $, $ f_y = \ frac {y_q} {\ lambda z_q} $

\begin{align} u(x_q, y_q, z_q) &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x, y, 0)e^{-2\pi i(f_x x + f_y y)}dx dy \\ &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \mathcal{F}\{u(x, y, 0)\}_{f_x = \frac{x_q}{\lambda z_q}, f_y = \frac{y_q}{\lambda z_q}}
\end{align}

Daher können wir sehen, dass es durch die zweidimensionale Fourier-Transformation von $ u (x_p, y_p, 0) $ berechnet werden kann. Damit die Fraunhofer-Näherung gilt, muss $ \ frac {k} {2z_q} (x ^ 2 + y ^ 2) \ ll 1 $ erfüllt sein. Dies kann umgeschrieben werden, damit $ z_q $ $ z_q \ ist Es wird gg \ frac {\ pi} {\ lambda} (x ^ 2 + y ^ 2) $.

Die Fraunhofer-Näherung ist eine gültige Näherung, wenn die Apertur klein genug ist und das Bild sehr weit von der Apertur entfernt beobachtet wird. Andererseits ist die Fresnel-Näherung eine gültige Näherung, die noch näher als die Fraunhofer-Näherung liegt.

Winkelspektrum-Methode

Die inverse Fourier-Transformation von $ U (f_x, f_y, 0) $, die die Fourier-Transformation der komplexen Amplitude $ u (x, y, 0) $ der offenen Oberfläche ist, ist

u(x, y, 0) = \int\int U(f_x, f_y, 0)\exp[2\pi i(x f_x + y f_y)]df_x df_y

ist. In dieser Gleichung ist die komplexe Amplitude $ u (x, y, 0) $ der Öffnungsfläche die Summe der ebenen Welle $ U (f_x, f_y, 0) \ exp [2 \ pi i (x f_x + y f_y)] $. Es zeigt, dass es ausgedrückt werden kann als.

Finden jeder Komponente des Wellenzahlvektors $ \ boldsymbol {k} $ jeder ebenen Welle

\begin{align}
k_x &= 2\pi f_x \\
k_y &= 2\pi f_y \\
k_z &= \sqrt{k^2 - k_x^2 - k_y^2}
\end{align}

Es wird sein.k_zIst\|\boldsymbol{k}\| = kIch werde von gehorchen.

Wenn sich jede ebene Welle $ z_q $ in Richtung $ z $ bewegt, ist ihre komplexe Amplitude

U(f_x, f_y, 0)\exp[2\pi i(x f_x + y f_y)] \exp[i k_z z_q]

Es wird sein. Dann ist die komplexe Amplitude $ u (x ', y', z_q) $ auf der Bildebene die Summe dieser ebenen Wellen.

\begin{align}
u(x', y', z_q) &= \int\int U(f_x, f_y, 0)\exp[2\pi i(x f_x + y f_y)] \exp[i k_z z_q] df_x df_y \\
&= \mathcal{F}^{-1}\{ U(f_x, f_y, 0) \exp[i k_z z_q] \}
\end{align}

Es ist ersichtlich, dass die komplexe Amplitude $ u (x ', y, z_q) $ auf der Bildebene durch die inverse Fourier-Transformation von $ U (f_x, f_y, 0) \ exp [i k_z z_q] $ berechnet werden kann. Diese Berechnungsmethode heißt ** Angular Spectrum Method **.

Die Winkelspektrummethode kann die Beugung in dem Bereich berechnen, der näher an der Öffnungsfläche liegt als die Frenel-Näherung oder die Fraunhofer-Näherung. Da die Berechnung für die Fourier-Transformation und die inverse Fourier-Transformation nur einmal durchgeführt wird, kann sie außerdem mit hoher Geschwindigkeit berechnet werden.

Numerisches Berechnungsbeispiel für die Beugung

Öffnung $ \ Sigma $ mit Radius $ r = 1 ~ \ mathrm {[mm]} $ kreisförmige Öffnung, $ \ lambda = 550 ~ \ mathrm {[nm]} $, $ z_q = 5 ~ \ mathrm {[m] Berechnen wir numerisch als} $. Das für die Berechnung verwendete Jupyter-Notizbuch finden Sie unter Google Colab.

Numerische Berechnung der Fraunhofer-Beugung

Berechnen wir zunächst numerisch mit der Fraunhofer-Näherung. Wenn Sie die Näherungsbedingungen überprüfen, sind $ x $ und $ y $ maximal $ r $.

z_q \gg \frac{\pi}{\lambda}r^2 \approx 5.71

Es wird sein. Es entspricht nicht den Annahmen, aber ich werde mit der Berechnung so fortfahren, wie sie ist.

Die Fraunhofer-Beugung kann berechnet werden, indem eine zweidimensionale Fourier-Transformation an der komplexen Amplitude $ u (x, y, 0) $ an der Apertur $ \ Sigma $ durchgeführt wird. Um numerische Berechnungen programmgesteuert durchzuführen, teilen Sie den Berechnungsbereich mit einer Seitenlänge von $ D $ in ein Gitter von $ N \ mal N $, um die Öffnungsfläche abzudecken, und speichern Sie die komplexen Amplituden auf jedem Gitter. Kann durch Ausführen einer zweidimensionalen diskreten Fourier-Transformation berechnet werden.

Berechnen wir numerisch mit Python. Bereiten Sie zunächst ein Array vor, das komplexe Amplituden im Raster speichert.

import numpy as np
import matplotlib.pyplot as plt

l = 0.55e-6 #Wellenlänge[m]
k = 2 * np.pi / l #Anzahl der Wellen
zq = 5 #Abstand von der Öffnungsfläche zur Bildfläche[m]
N = 200 #Anzahl der Gitterunterteilungen
r = 0.001 #Öffnungsradius[m]
D = 0.05 #Rechenbereich auf der offenen Fläche[m]

fscale = l*zq

#Gibt die komplexe Amplitude auf der offenen Oberfläche zurück
def P(x, y):
    if x**2 + y**2 < r**2:
        return 1
    else:
        return 0

#Netzerzeugung
X, Y = np.meshgrid(np.linspace(-D/2, D/2, num=N), np.linspace(-D/2, D/2, num=N))
Pv = np.vectorize(P)
Z = Pv(X, Y)

Lassen Sie es uns planen.

plt.imshow(Z, extent=[-D/2, D/2, -D/2, D/2], cmap='gray')
plt.xlabel('$x_p$ [m]')
plt.ylabel('$y_p$ [m]')

download.png

Der weiße Teil entspricht der Öffnungsfläche.

Als nächstes wird dies einer diskreten Fourier-Transformation unterzogen, aber vorher wird der Bereich auf der Bildebene berechnet, der dem Berechnungsbereich entspricht.

v = np.fft.fftfreq(N, d=D/N)
v = fscale * np.fft.fftshift(v)
#Einheiten[mm]
extent= [1e3 * v[0], 1e3 * v[-1], 1e3 * v[0], 1e3 * v[-1]]

Das Fraunhofer-Beugungsbild wird durch diskrete Fourier-Transformation der komplexen Amplitude auf der offenen Oberfläche erhalten.

#Berechnen Sie die komplexe Amplitude in der Bildebene
U_fraun = np.fft.fft2(Z)
U_fraun = 1/(1.0j * l * zq) * np.exp(1.0j * k * zq) * np.exp(1.0j * k * (v**2 + v**2)/(2*zq)) * np.fft.fftshift(U_fraun)

#Stärke berechnen
I_fraun = np.abs(U_fraun)**2
#Normalisierung
I_fraun = I_fraun / np.max(I_fraun)

Zeichnen wir die Intensität.

plt.imshow(I_fraun, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fraunhofer Diffraction')

download (1).png

Es gibt auch einen Intensitätspeak um das Zentrum, der zeigt, dass der Beugungseffekt berechnet werden kann.

Numerische Berechnung der Fresnel-Beugung

Überprüfen Sie zunächst die Annäherungsbedingungen. Der bedingte Ausdruck für die Approximation der Fresnel-Beugung ist

z_q^3 \gg \frac{\pi}{4\lambda}[(x_q - x)^2 + (y_q - y)^2]^2

war. Dies ist etwas kompliziert, da im Gegensatz zur Fraunhofer-Beugung die Größe der Beobachtungsfläche berücksichtigt werden muss. Vorerst die Größe einer Seite der BeobachtungsflächelIch werde es als belassen. Dann|x_q - x| \le r + \frac{l}{2}、yIst das gleiche, also wenn Sie dies ersetzen

z_q \gg \sqrt[3]{\frac{\pi}{\lambda}(r + \frac{l}{2})^4} \approx 0.21

Es wird sein.

Die Frennell-Beugung kann sowohl mit der zweidimensionalen diskreten Fourier-Transformation als auch mit der Fraunhofer-Beugung berechnet werden.

#Berechnen Sie die komplexe Amplitude in der Bildebene
U_fresnel = np.fft.fft2(Z * np.exp(1.0j * k * (X**2 + Y**2)/(2*zq)))
U_fresnel = 1/(1.0j * l * zq) * np.exp(1.0j * k * zq) * np.exp(1.0j * k * (v**2 + v**2)/(2*zq)) * np.fft.fftshift(U_fresnel)

#Stärke berechnen
I_fresnel = np.abs(U_fresnel)**2
#Normalisierung
I_fresnel = I_fresnel / np.max(I_fresnel)

Zeichnen wir die Intensität.

plt.imshow(I_fresnel, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fresnel Diffraction')

download (2).png

Es scheint, dass der Peak um das Zentrum stärker ist als in der Fraunhofer-Näherung.

Direkte Berechnung der Fresnel-Kirchhoff-Beugungsformel

Fresnel Kirchhoffs Beugungsformel

\frac{z_q}{i\lambda}\int\int_{\Sigma} u(x, y, 0)\frac{e^{ikr}}{r^2}dx dy

Berechnen wir durch direkte numerische Integration. Diese Methode ist relativ streng, da keine Näherungswerte verwendet werden, sie ist jedoch aufgrund des überwältigenden Rechenaufwands nicht praktikabel. Ich wollte es mit der Approximationsmethode vergleichen, also habe ich es berechnet.

#Fresnel-Kirchhoffs diffraktive numerische Integration
def Uf(x_q, y_q):
    r = np.sqrt((x_q - X)**2 + (y_q - Y)**2 + zq**2)
    return zq/(1.0j * l) * (Z * np.exp(1.0j * k * r) / r**2).sum()

#Berechnungsbereich auf der Bildebene
X_Q, Y_Q = np.meshgrid(np.linspace(extent[0]*1e-3, extent[1]*1e-3, num=N), np.linspace(extent[2]*1e-3, extent[3]*1e-3, num=N))
#Numerische Integration(Eine Skalierung ist angemessen)
U_strict = np.vectorize(Uf)(X_Q, Y_Q)

#Stärke berechnen
I_strict = np.abs(U_strict)**2
#Normalisierung
I_strict = I_strict / np.max(I_strict)

Zeichnen wir die Intensität.

plt.imshow(I_strict, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fresnel Kirchhoff Diffraction')

download (3).png

Numerische Berechnung nach der Winkelspektrummethode

Lassen Sie uns abschließend eine numerische Berechnung mit der Winkelspektrummethode durchführen. Die Prozedur besteht darin, zuerst die Fourier-Transformation $ U (f_x, f_y, 0) $ der komplexen Amplitude $ u (x, y, 0) $ der offenen Oberfläche und dann $ U (f_x, f_y, 0) $ und $ \ zu finden. Inverses Fourier transformiert das Produkt von exp [ik_z z_q] $.

#Fourier-Transformation der komplexen Amplitude der offenen Oberfläche
U = np.fft.fft2(Z)

#Berechnung des Wellenzahlvektors
k_x = np.fft.fftfreq(N, d=D/N) * 2 * np.pi
k_y = np.fft.fftfreq(N, d=D/N) * 2 * np.pi
K_X, K_Y = np.meshgrid(k_x, k_y)
k_z = np.sqrt(k**2 - K_X**2 - K_Y**2)

#Berechnen Sie die komplexe Amplitude der Bildebene
U_angular = np.fft.ifft2(U * np.exp(1.0j * k_z * zq))

#Stärke berechnen
I_angular = np.abs(U_angular)**2
#Normalisierung
I_angular = I_angular / np.max(I_angular)

Lassen Sie es uns planen.

plt.imshow(I_angular, cmap='jet', extent=[-D/2, D/2, -D/2, D/2])
plt.xlabel('$x_q$ [m]')
plt.ylabel('$y_q$ [m]')
plt.title('Angular Spectrum Diffraction')

download (4).png

Im Gegensatz zur Fresnel-Näherung und zur Fraunhofer-Näherung entspricht die Bildebenenfläche der Berechnungsfläche, sodass sie klein aussieht.

Es ist lange her, also schreibe ich den Rest auf eine andere Seite.

Referenzen etc.

Recommended Posts

Numerische Berechnung der Linsenpunktbildverteilungsfunktion und der MTF-Kurve (Beugungsberechnung)
[Einführung in Scipy] Berechnung der Lorenzkurve und des Gini-Koeffizienten ♬