프레넬 회절을 분석하고 싶습니다.

배경

근접 노광을 이용한 패턴 형성을 최근 하고 있지만, 구체적인 광 강도가 신경이 쓰이기 시작했습니다. 조사해 보면, 다음과 같은 프로그램이 이미 있는 모양.

근접 노출 계산 프로그램: htp : / / s 코모. 오. 오 7. jp / f54 / hp54_1. htm

이쪽을 그대로 사용해도 좋았습니다만, Excel 경유하고 싶지 않나, 계산 영역 더 분할하고 싶다고 생각했기 때문에, 과학 계산에 강하다고 듣는 Python으로 실장해 보기로 했습니다.
실장에 있어서는, 상기의 것을 참고로 하고 있습니다.

프레넬 회절이란?

슬릿 - 스크린 사이의 거리가 작을 때 발생하는 회절 현상입니다.
슬릿이 단일 직사각형 인 경우, 스크린상의 빛의 파동 $ U (P_n) $는 슬릿면의 빛의 파동 $ U (S_n) $를 사용하여 다음과 같이 간단하게 표현됩니다.

U(P_n)=\frac{U(S_n)}{2\mathrm{i}}\left[\int_{v_1}^{v_2}\exp\left(\frac{\pi v^2}{2}\mathrm{i}\right)dv\right]\left[\int_{w_1}^{w_2}\exp\left(\frac{\pi w^2}{2}\mathrm{i}\right)dw\right]

변수 변환은 다음과 같이 수행됩니다.
또한 화면상의 좌표를 $x$ 와 $y$ , 슬릿 개구부의 좌표를 $\xi$ 와 $\eta$ , 빛의 파장을 $\lambda$, 슬릿-스크린 사이의 거리를 $b$ 로 표시했습니다.

v = \sqrt{\frac{2}{\lambda b}}(\xi - x)

w = \sqrt{\frac{2}{\lambda b}}(\eta - y)

위의 적분을 풀려고 하면 귀찮은 녀석(프레넬 적분)이 나타납니다만, scipy.special 모듈을 사용하는 것으로 심슨 법칙을 준비하지 않고 끝난다고 합니다. 편리.

구현

※가바코드 주의보※

main.py

import numpy as np
from scipy.special import fresnel
import matplotlib.pyplot as plt

# 露光条件の定義
WL = 0.365  # 波長(μm)
gap = 5  # スリット-スクリーン間の距離(μm)
INA = 0  # 光の広がり角度(rad)
X1 = -5  # 計算領域のX方向最小値(μm)
X2 = 5  # 計算領域のX方向最大値(μm)
Nx = 100  # 計算領域のX方向分割数
Y1 = -100  # 計算領域のY方向最小値(μm)
Y2 = 100  # 計算領域のY方向最大値(μm)
Ny = 200  # 計算領域のY方向分割数
Pr = 1  # スリット面の光の実数部
Pi = 0  # スリット面の光の虚数部
Px1 = -0.5  # 矩形パターンのX方向最小値(μm)
Px2 = 0.5  # 矩形パターンのX方向最大値(μm)
Py1 = -50  # 矩形パターンのY方向最小値(μm)
Py2 = 50  # 矩形パターンのY方向最大値(μm)

# 計算領域の作成
dx = (X2 - X1) / Nx
dy = (Y2 - Y1) / Ny
Xaxis = np.arange(X1, X2, dx)
Ysaxis = np.arange(Y1, Y2, dy)
Yaxis = np.c_[Ysaxis]
Isum = np.zeros((Ny, Nx))


# フレネル積分
def cal_f(t1, t2):
    s1, c1 = fresnel(t1)
    s2, c2 = fresnel(t2)
    dif_s = s2 - s1
    dif_c = c2 - c1
    return dif_s, dif_c


# 光強度の計算
def cal_int(ox, oy):
    fx1 = lambda x: np.sqrt(2 / (WL * gap)) * (Px1 - x - ox)
    fx2 = lambda x: np.sqrt(2 / (WL * gap)) * (Px2 - x - ox)

    fy1 = lambda y: np.sqrt(2 / (WL * gap)) * (Py1 - y - oy)
    fy2 = lambda y: np.sqrt(2 / (WL * gap)) * (Py2 - y - oy)

    vfx1 = np.vectorize(fx1)
    vfx2 = np.vectorize(fx2)
    vfy1 = np.vectorize(fy1)
    vfy2 = np.vectorize(fy2)
    vf_int = np.vectorize(cal_f)

    asv, acv = vf_int(vfx1(Xaxis), vfx2(Xaxis))
    asw, acw = vf_int(vfy1(Yaxis), vfy2(Yaxis))

    ar = (acv * asw + asv * acw) / 2
    ai = (asv * asw - acv * acw) / 2

    ur = Pr * ar - Pi * ai
    ui = Pr * ai + Pi * ar
    li = (np.power(ur, 2) + np.power(ui, 2))
    return li


# 光の広がりの考慮
i = -5
k = 0
while i <= 5:
    j = -5
    while j <= 5:
        if i ** 2 + j ** 2 <= 25:
            k += 1

            Ox = i * gap * INA / 5
            Oy = j * gap * INA / 5

            Isub = cal_int(Ox, Oy)
            Isum = Isum + Isub

            j += 1
        else:
            j += 1
    else:
        i += 1

Iave = Isum / k

# グラフ描画
plt.rcParams['font.size'] = 18
plt.rcParams['font.family'] = 'sans-serif'
plt.tight_layout()

plt.imshow(Iave, extent=[X1, X2, Y1, Y2], aspect='auto', interpolation='nearest', cmap='jet')
plt.xlabel('x-position (μm)')
plt.ylabel('y-position (μm)')
plt.colorbar()
plt.clim(0, 1.0)
plt.grid(True)

plt.show()