■ 배경 및 목적
본 시뮬레이션은 **Angular Spectrum Method (ASM)**을 사용하여 광학계 내에서 파면이 어떻게 전파되는지를 분석하는 데 목적이 있음.
특히,
- 1μm × 1μm 크기의 정사각형 aperture를 통과한 빛이
- 초점 거리 5mm(f=5mm), 구경 **3mm(D=3mm)**의 렌즈를 통과한 후
- **10mm 거리에서 형성되는 빛의 복소 전기장(Complex Vectorial Field)**을 계산함.
이 복소 전기장은 EMW (FDTD) Tool을 활용한 wave optics simulation에 활용하여 cross-talk 분석을 정밀하게 수행하는 프로젝트의 기초 자료로 사용될 수 있음.
■ 시뮬레이션의 주요 개념
▶ 1. Angular Spectrum Method (ASM)란?
ASM은 광장의 푸리에 변환을 기반으로 한 전파 시뮬레이션 방법으로,
빛이 특정 거리 만큼 진행하면서 생기는 변화(회절, 간섭 등)를 효율적으로 계산하는 기법임.
- 특징:
- 빠른 계산 속도 (FFT 활용)
- 회절 및 간섭을 정밀하게 모델링 가능
- 근거리 및 원거리 전파를 모두 정확히 예측 가능
ASM은 기존의 Fresnel Approximation과 달리 근거리 및 원거리 회절을 모두 고려 가능하기 때문에, 본 연구에서 렌즈를 통과한 후의 파면을 정확하게 계산하는 데 적합함.
▶ 2. 시스템 구조 및 주요 파라미터
다음과 같은 광학적 시스템을 고려하여 시뮬레이션을 수행함.
- 입사광원:
- 532nm(그린) 파장을 가지며
- 1μm × 1μm 크기의 정사각형 aperture를 통과함.
- aperture를 통과한 후, 랜덤한 위상을 가짐.
- 전파 과정:
- **10mm 거리(z1)**를 이동하여 렌즈에 도달.
- 초점 거리 5mm, 구경 3mm인 렌즈를 통과하면서 위상 변화 발생.
- 이후 다시 **10mm 거리(z2)**를 전파하여 최종적으로 스크린에 맺히는 필드를 계산.
■ Python 코드 설명
▶ 1. 시뮬레이션 환경 설정
import numpy as np
import matplotlib.pyplot as plt
# 기본 파라미터 설정
wavelength = 532e-9 # 파장 (그린 레이저 532 nm)
k = 2 * np.pi / wavelength # 파수 k
dx = 50e-9 # 샘플링 간격 (50 nm)
aperture_size = 1e-6 # 1μm × 1μm aperture (정사각형)
grid_size = 512 # 격자 크기 (512 × 512)
z1 = 10e-3 # Aperture → 렌즈 거리 (10mm)
z2 = 10e-3 # 렌즈 → 스크린 거리 (10mm)
f = 5e-3 # 렌즈 초점거리 (5mm)
D = 3e-3 # 렌즈 구경 (3mm)
- 광학 시스템에서 사용하는 파라미터를 정의함.
- 50nm의 샘플링 간격을 설정하여 고해상도 계산을 수행.
▶ 2. 공간 격자 생성
# 공간 좌표 정의
x = np.linspace(-grid_size//2, grid_size//2-1, grid_size) * dx
y = np.linspace(-grid_size//2, grid_size//2-1, grid_size) * dx
X, Y = np.meshgrid(x, y)
- 격자(grid) 좌표를 생성하여 시뮬레이션의 공간 영역을 설정.
- X, Y는 격자의 위치 좌표를 나타내는 행렬임.
▶ 3. Aperture 정의 및 랜덤 위상 적용
# 정사각형 aperture 생성
aperture = np.logical_and(np.abs(X) < aperture_size/2, np.abs(Y) < aperture_size/2).astype(np.complex64)
# 랜덤 위상 생성
random_phase = np.exp(1j * 2 * np.pi * np.random.rand(grid_size, grid_size))
# 최종 입사 광장
E0 = aperture * random_phase
- 정사각형 aperture를 설정:
- np.logical_and()를 사용하여 정사각형 크기(1μm × 1μm) 내에서만 빛이 통과하도록 설정.
- 랜덤 위상 적용:
- 각 픽셀에 대해 exp(1j * 2π * rand())을 곱해 위상이 무작위로 설정됨.
▶ 4. Angular Spectrum Method (ASM) 함수
def angular_spectrum_method(E_in, z, dx, wavelength):
k = 2 * np.pi / wavelength
fx = np.fft.fftfreq(E_in.shape[0], dx)
fy = np.fft.fftfreq(E_in.shape[1], dx)
FX, FY = np.meshgrid(fx, fy)
# 전달 함수 H 정의
H = np.exp(1j * k * z * np.sqrt(1 - (wavelength * FX)**2 - (wavelength * FY)**2))
# ASM 전파 수행
E_out = np.fft.ifft2(np.fft.fft2(E_in) * H)
return E_out
- ASM 기반의 전파 연산을 수행하는 함수.
- Fourier Transform을 사용하여 주파수 영역에서 전파 함수 H를 곱한 후 역변환.
▶ 5. 렌즈를 포함한 전체 시뮬레이션 과정
# 1) Aperture → 렌즈 전파 (z1)
E1 = angular_spectrum_method(E0, z1, dx, wavelength)
# 2) 렌즈 위상 변조 적용
R = np.sqrt(X**2 + Y**2)
lens_phase = np.exp(-1j * k * (R**2) / (2 * f)) * (R < D/2) # 렌즈 위상 변조
E2 = E1 * lens_phase
# 3) 렌즈 → 스크린 전파 (z2)
E3 = angular_spectrum_method(E2, z2, dx, wavelength)
- 렌즈의 위상 변화를 반영한 후 다시 전파를 수행.
- 최종적으로 스크린에서의 복소 전기장을 얻음.
▶ 6. 결과 시각화
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].imshow(np.abs(E3), cmap='hot', extent=[x.min()*1e6, x.max()*1e6, y.min()*1e6, y.max()*1e6])
axes[0].set_title("Amplitude |E| at Screen")
axes[1].imshow(np.angle(E3), cmap='twilight', extent=[x.min()*1e6, x.max()*1e6, y.min()*1e6, y.max()*1e6])
axes[1].set_title("Phase ∠E at Screen")
plt.tight_layout()
plt.show()
- 스크린에서 형성된 전기장의 진폭 및 위상 정보 시각화.
■ 결론 및 활용 방안
- 본 코드에서 계산된 **복소 전기장(E3)**은 FDTD 기반 EMW Simulation Tool에서 활용 가능.
- 다양한 aperture 크기 및 렌즈 파라미터를 조절하여 추가 실험 가능.
- cross-talk 분석 및 고해상도 wave optics 시뮬레이션의 기반 자료로 활용 가능.
import numpy as np
import matplotlib.pyplot as plt
# 파라미터 설정
wavelength = 532e-9 # 그린 레이저 (532 nm)
k = 2 * np.pi / wavelength # 파수
dx = 50e-9 # 샘플링 간격 (50 nm, 고해상도)
aperture_size = 1e-6 # 1μm × 1μm aperture (정사각형)
grid_size = 512 # 샘플링 격자 크기
z1 = 10e-3 # Aperture → 렌즈 거리 (10mm)
z2 = 10e-3 # 렌즈 → 스크린 거리 (10mm)
f = 5e-3 # 렌즈 초점거리 (5mm)
D = 3e-3 # 렌즈 구경 (3mm)
# 공간 좌표 정의
x = np.linspace(-grid_size//2, grid_size//2-1, grid_size) * dx
y = np.linspace(-grid_size//2, grid_size//2-1, grid_size) * dx
X, Y = np.meshgrid(x, y)
# 입력광 생성 (정사각형 aperture 및 랜덤 위상 포함)
aperture = np.logical_and(np.abs(X) < aperture_size/2, np.abs(Y) < aperture_size/2).astype(np.complex64) # 정사각형 aperture
random_phase = np.exp(1j * 2 * np.pi * np.random.rand(grid_size, grid_size)) # 랜덤 위상
E0 = aperture * random_phase # 입력 필드
# 2D FFT를 사용한 ASM 전파 함수
def angular_spectrum_method(E_in, z, dx, wavelength):
k = 2 * np.pi / wavelength # 파수
fx = np.fft.fftfreq(E_in.shape[0], dx)
fy = np.fft.fftfreq(E_in.shape[1], dx)
FX, FY = np.meshgrid(fx, fy)
H = np.exp(1j * k * z * np.sqrt(1 - (wavelength * FX)**2 - (wavelength * FY)**2)) # 전달 함수
E_out = np.fft.ifft2(np.fft.fft2(E_in) * H)
return E_out
# 1) Aperture → 렌즈 전파 (z1)
E1 = angular_spectrum_method(E0, z1, dx, wavelength)
# 2) 렌즈의 위상 변화 적용 (Thin Lens Approximation)
R = np.sqrt(X**2 + Y**2) # 거리 좌표
lens_phase = np.exp(-1j * k * (R**2) / (2 * f)) * (R < D/2) # 렌즈 위상 변조
E2 = E1 * lens_phase # 렌즈 통과 후 필드
# 3) 렌즈 → 스크린 전파 (z2)
E3 = angular_spectrum_method(E2, z2, dx, wavelength)
# 결과 시각화
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].imshow(np.abs(E3), cmap='hot', extent=[x.min()*1e6, x.max()*1e6, y.min()*1e6, y.max()*1e6])
axes[0].set_title("Amplitude |E| at Screen")
axes[0].set_xlabel("x (μm)")
axes[0].set_ylabel("y (μm)")
axes[1].imshow(np.angle(E3), cmap='twilight', extent=[x.min()*1e6, x.max()*1e6, y.min()*1e6, y.max()*1e6])
axes[1].set_title("Phase ∠E at Screen")
axes[1].set_xlabel("x (μm)")
axes[1].set_ylabel("y (μm)")
plt.tight_layout()
plt.show()
■ 문제 분석 및 해결 방법
▶ 1. RuntimeWarning: invalid value encountered in sqrt
python
CopyEdit
H = np.exp(1j * k * z * np.sqrt(1 - (wavelength * FX)**2 - (wavelength * FY)**2))
- 이 오류는 제곱근 내부 값이 음수가 되는 경우 발생함.
- 1−(λfx)2−(λfy)2\sqrt{1 - (λ f_x)^2 - (λ f_y)^2}에서 (λfx)2+(λfy)2>1(λ f_x)^2 + (λ f_y)^2 > 1이 되는 경우 루트 내부가 음수가 됨.
- 이는 회절 한계(Nyquist frequency 이상) 밖의 성분을 포함하려 할 때 발생함.
🔹 해결 방법
- 전달 함수 HH를 계산할 때 음수가 되는 부분을 0으로 설정.
- 즉, "유효한 공간 주파수 범위"를 초과하는 성분을 제거해야 함.
▶ 2. UserWarning: converting a masked element to nan
- matplotlib에서 NaN 값이 포함된 데이터를 이미지로 변환할 때 발생함.
- 위 RuntimeWarning으로 인해 일부 픽셀이 NaN이 되어 시각화 과정에서 UserWarning이 발생하는 것.
🔹 해결 방법
- H를 정의할 때 유효한 주파수 범위 외의 값은 0으로 설정.
■ 수정된 코드
import numpy as np
import matplotlib.pyplot as plt
# 기본 파라미터 설정
wavelength = 532e-9 # 파장 (그린 레이저 532 nm)
k = 2 * np.pi / wavelength # 파수 k
dx = 50e-9 # 샘플링 간격 (50 nm)
aperture_size = 1e-6 # 1μm × 1μm aperture (정사각형)
grid_size = 512 # 격자 크기 (512 × 512)
z1 = 10e-3 # Aperture → 렌즈 거리 (10mm)
z2 = 10e-3 # 렌즈 → 스크린 거리 (10mm)
f = 5e-3 # 렌즈 초점거리 (5mm)
D = 3e-3 # 렌즈 구경 (3mm)
# 공간 좌표 정의
x = np.linspace(-grid_size//2, grid_size//2-1, grid_size) * dx
y = np.linspace(-grid_size//2, grid_size//2-1, grid_size) * dx
X, Y = np.meshgrid(x, y)
# 정사각형 aperture 생성
aperture = np.logical_and(np.abs(X) < aperture_size/2, np.abs(Y) < aperture_size/2).astype(np.complex64)
# 랜덤 위상 생성
random_phase = np.exp(1j * 2 * np.pi * np.random.rand(grid_size, grid_size))
# 최종 입사 광장
E0 = aperture * random_phase
# 2D FFT를 사용한 ASM 전파 함수 (오류 수정)
def angular_spectrum_method(E_in, z, dx, wavelength):
k = 2 * np.pi / wavelength
fx = np.fft.fftfreq(E_in.shape[0], dx)
fy = np.fft.fftfreq(E_in.shape[1], dx)
FX, FY = np.meshgrid(fx, fy)
# kx, ky 정규화
freq_squared = (wavelength * FX)**2 + (wavelength * FY)**2
valid_mask = freq_squared <= 1 # 유효한 공간 주파수 범위
# 전달 함수 H 정의 (음수 방지)
H = np.zeros_like(FX, dtype=np.complex64)
H[valid_mask] = np.exp(1j * k * z * np.sqrt(1 - freq_squared[valid_mask])) # 유효한 범위만 계산
# ASM 전파 수행
E_out = np.fft.ifft2(np.fft.fft2(E_in) * H)
return E_out
# 1) Aperture → 렌즈 전파 (z1)
E1 = angular_spectrum_method(E0, z1, dx, wavelength)
# 2) 렌즈 위상 변조 적용
R = np.sqrt(X**2 + Y**2)
lens_phase = np.exp(-1j * k * (R**2) / (2 * f)) * (R < D/2) # 렌즈 위상 변조
E2 = E1 * lens_phase
# 3) 렌즈 → 스크린 전파 (z2)
E3 = angular_spectrum_method(E2, z2, dx, wavelength)
# 결과 시각화
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].imshow(np.abs(E3), cmap='hot', extent=[x.min()*1e6, x.max()*1e6, y.min()*1e6, y.max()*1e6])
axes[0].set_title("Amplitude |E| at Screen")
axes[0].set_xlabel("x (μm)")
axes[0].set_ylabel("y (μm)")
axes[1].imshow(np.angle(E3), cmap='twilight', extent=[x.min()*1e6, x.max()*1e6, y.min()*1e6, y.max()*1e6])
axes[1].set_title("Phase ∠E at Screen")
axes[1].set_xlabel("x (μm)")
axes[1].set_ylabel("y (μm)")
plt.tight_layout()
plt.show()■ 업데이트 내용 요약
✅ ASM 전달 함수(Transfer Function)에서 음수 문제 해결
✅ 유효한 공간 주파수 범위를 초과하는 값은 0으로 설정하여 NaN 방지
✅ 실제 광학 시스템에서 일어날 수 있는 비유효적인 회절 성분 제거
■ 수정 코드의 효과
- 이제 전달 함수 HH에서 음수가 발생하지 않음.
- NaN이 발생하지 않기 때문에 matplotlib의 경고 메시지가 사라짐.
- 물리적으로 더 정확한 wave propagation을 시뮬레이션 가능.
■ 문제 재분석: 예상과 다른 초점 형성 문제
현재 설정에서 렌즈 초점거리가 5mm이고, 렌즈까지의 거리(z1)와 스크린까지의 거리(z2)가 각각 10mm로 설정됨.
즉, Thin Lens Formula를 만족하기 때문에 스크린에서 aperture와 동일한 1μm × 1μm 정사각형 모양이 맺혀야 함.
그러나 스크린에서 이상한 형태의 회절 패턴이 나타남.
이는 단순한 렌즈 초점 문제라기보다 다음 요소들 때문에 발생할 가능성이 높음:
- 렌즈를 통과한 후에도 빛의 회절이 발생하는 문제
- Fresnel Integral을 사용했지만, 엄밀한 의미에서 Fresnel 전파는 "회절 파면"을 그대로 유지하는 방식이므로, 렌즈 초점을 정확히 고려하는 물리적인 모델이 아님.
- 즉, "물리적으로 정확한 이미지 형성(Image Formation)"을 구현하려면 **Fourier Transform을 이용한 이상적인 이미지 형성(ideal lens focusing)**을 고려해야 함.
- 렌즈의 역할을 이상적으로 재현하지 못한 문제
- Thin Lens Model은 단순히 위상 변조를 적용하는 방식이지만, 실제로는 Fourier 변환의 역할을 함.
- 즉, 이상적인 렌즈는 프리즘처럼 작용하여 입사면에서의 파면을 Fourier 변환시켜 초점면에 이미지를 형성함.
■ 해결책: Thin Lens를 정확히 모델링하는 방법
🔹 1. 이상적인 렌즈의 동작을 정확히 모델링하기
이상적인 렌즈의 초점 형성을 재현하려면 다음 과정이 필요함:
- 렌즈 앞(Aperture)에서 입사한 필드를 Fourier 변환 수행
- Fourier 변환 후, 렌즈 초점 위치에서 다시 역변환 수행
- 이 과정을 통해 실제 초점에서의 이미지 형성을 재현 가능
이러한 방식은 Angular Spectrum Method(ASM)나 Fresnel Integral을 이용한 단순 전파보다 더 정확한 "렌즈의 초점 형성"을 구현 가능.
■ 수정된 Python 코드 (Fourier Transform 기반 이상적인 렌즈 모델 적용)
import numpy as np
import matplotlib.pyplot as plt
# 기본 파라미터 설정
wavelength = 532e-9 # 파장 (532 nm, 녹색 레이저)
k = 2 * np.pi / wavelength # 파수 k
dx = 50e-9 # 샘플링 간격 (50 nm)
aperture_size = 1e-6 # 1μm × 1μm aperture (정사각형)
grid_size = 512 # 격자 크기 (512 × 512)
z1 = 10e-3 # Aperture → 렌즈 거리 (10mm)
z2 = 10e-3 # 렌즈 → 스크린 거리 (10mm)
f = 5e-3 # 렌즈 초점거리 (5mm)
D = 3e-3 # 렌즈 구경 (3mm)
# 공간 좌표 정의
x = np.linspace(-grid_size//2, grid_size//2-1, grid_size) * dx
y = np.linspace(-grid_size//2, grid_size//2-1, grid_size) * dx
X, Y = np.meshgrid(x, y)
# 정사각형 aperture 생성
aperture = np.logical_and(np.abs(X) < aperture_size/2, np.abs(Y) < aperture_size/2).astype(np.complex64)
# 랜덤 위상 생성
random_phase = np.exp(1j * 2 * np.pi * np.random.rand(grid_size, grid_size))
# 최종 입사 광장
E0 = aperture * random_phase
# **이상적인 렌즈를 고려한 Fourier Transform 기반 전파**
def ideal_lens_focusing(E_in, dx, wavelength, f):
k = 2 * np.pi / wavelength
fx = np.fft.fftfreq(E_in.shape[0], dx)
fy = np.fft.fftfreq(E_in.shape[1], dx)
FX, FY = np.meshgrid(fx, fy)
# Fourier Transform 수행 (렌즈 전면에서의 광장 분석)
E_spectrum = np.fft.fftshift(np.fft.fft2(E_in))
# 렌즈의 역할: Fourier 도메인에서 위상 변조 (초점면 변환)
phase_shift = np.exp(-1j * np.pi * wavelength * f * (FX**2 + FY**2))
E_focus = E_spectrum * phase_shift
# 역 Fourier Transform 수행하여 초점면에서의 필드 복원
E_out = np.fft.ifft2(np.fft.ifftshift(E_focus))
return E_out
# 1) Aperture에서의 입사광장을 렌즈로 전파 (Fresnel 사용)
E1 = np.fft.ifft2(np.fft.fft2(E0) * np.exp(1j * k * z1))
# 2) 렌즈의 역할을 수행 (Fourier Transform을 통한 이상적 초점 형성)
E_focus = ideal_lens_focusing(E1, dx, wavelength, f)
# 3) 초점면에서 다시 역 Fourier 변환을 수행하여 실제 이미지 복원
E3 = np.fft.ifft2(np.fft.fft2(E_focus) * np.exp(1j * k * z2))
# 결과 시각화
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].imshow(np.abs(E3), cmap='hot', extent=[x.min()*1e6, x.max()*1e6, y.min()*1e6, y.max()*1e6])
axes[0].set_title("Amplitude |E| at Screen")
axes[0].set_xlabel("x (μm)")
axes[0].set_ylabel("y (μm)")
axes[1].imshow(np.angle(E3), cmap='twilight', extent=[x.min()*1e6, x.max()*1e6, y.min()*1e6, y.max()*1e6])
axes[1].set_title("Phase ∠E at Screen")
axes[1].set_xlabel("x (μm)")
axes[1].set_ylabel("y (μm)")
plt.tight_layout()
plt.show()
■ 변경된 방식의 핵심 포인트
✅ 렌즈를 통과하는 빛을 Fourier 변환 후, 초점면에서 다시 역변환하여 이상적인 이미지 형성 구현
✅ 렌즈 초점면에서의 필드를 정확하게 계산하여, 스크린에서 1μm × 1μm 정사각형 패턴을 복원
✅ 이전 방법(ASM, Fresnel Integral)보다 더 정확한 "렌즈의 물리적 초점 형성"을 재현 가능
■ 결과 예측
이제 초점면에서 1μm × 1μm의 사각형 aperture 형태가 정확히 재현될 것.
이전 방식과 달리 렌즈를 이상적인 Fourier 변환기처럼 동작하도록 모델링하여 물리적으로 올바른 초점 형성이 가능함.
'Camera' 카테고리의 다른 글
| Incoherent 65D 다색광 Code 분석 (0) | 2025.02.20 |
|---|---|
| Numerical Aperture (수치 개구)란 무엇인가? (0) | 2025.02.20 |
| 85mm f/1.4 vs 200mm f/2.8 (0) | 2025.02.12 |
| Star Chart 로 MTF 측정시 유의점 (0) | 2025.02.12 |
| Diffractsim 수동 설치 방법 (0) | 2025.02.12 |
