본문 바로가기
Camera

Incoherent 65D 다색광 Code 분석

by 이센 2025. 2. 20.

optical_imaging_system_using_convolution__polychromatic_incoherent.py  코드 

■ 코드 개요

  • 비코히런트(Incoherent) 광학 이미징 시스템 시뮬레이션
  • 출구 동공(Exit Pupil)을 통한 광학 전달 함수(OTF) 적용
  • 푸리에 변환(FFT) 기반 필터링
  • Rayleigh 기준을 이용한 해상도 분석
import diffractsim
##diffractsim.set_backend("CUDA")
import numpy as np
from diffractsim import PolychromaticField, ApertureFromImage, cf, mm, cm, CircularAperture

"""
MPL 2.0 License

Copyright (c) 2022, Rafael de la Fuente
All rights reserved.
"""


def get_colors_at_image_plane(F, radius, M,  zi, z0):
    from diffractsim.util.backend_functions import backend as bd
    from diffractsim.util.backend_functions import backend_name

    import numpy as np
    import time
    import progressbar
    """
    Assuming an incoherent optical system with linear response and assuming the system is only diffraction-limited by
    the exit pupil of the system, compute the field at its image plane

    
    Parameters
    ----------

    radius: exit pupil radius

    zi: distance from the image plane to the exit pupil
    z0: distance from the exit pupil to the current position

    M: magnification factor of the optical system
    (If the optical system is a single lens, magnification = - zi/z0)

    Reference:
    Introduction to Fourier Optics J. Goodman, Frequency Analysis of Optical Imaging Systems
    
    """
    pupil = CircularAperture(radius)
    F.z += zi + z0




    for j in range(len(F.optical_elements)):
        F.E = F.E * F.optical_elements[j].get_transmittance(F.xx, F.yy, 0)

    # if the magnification is negative, the image is inverted
    if M < 0:
        F.E = bd.flip(F.E)
    M_abs = bd.abs(M)

    F.E = F.E/M_abs

    Ip = F.E * bd.conjugate(F.E)
    
    fft_c = bd.fft.fft2(Ip)
    c = bd.fft.fftshift(fft_c)

    fx = bd.fft.fftshift(bd.fft.fftfreq(F.Nx, d = F.x[1]-F.x[0]))/M_abs
    fy = bd.fft.fftshift(bd.fft.fftfreq(F.Ny, d = F.y[1]-F.y[0]))/M_abs
    fx, fy = bd.meshgrid(fx, fy)
    fp = bd.sqrt(fx**2 + fy**2)

    bar = progressbar.ProgressBar()

    # We compute the pattern of each wavelength separately, and associate it to small spectrum interval dλ = (780- 380)/spectrum_divisions . We approximately the final colour
    # by summing the contribution of each small spectrum interval converting its intensity distribution to a RGB space.
    
    sRGB_linear = bd.zeros((3, F.Nx * F.Ny))

    t0 = time.time()

    for i in bar(range(F.spectrum_divisions)):
        #Definte the OTF function, representing the Fourier transform of the circular pupil function.

        fc = radius / (F.λ_list_samples[i]* nm  * zi) # coherent cutoff frequency

        H = pupil.get_optical_transfer_function(fx, fy, zi, F.λ_list_samples[i]* nm )
        #H = bd.where(fp < 2 * fc, 2/bd.pi * (bd.arccos(fp / (2*fc)) - fp / (2*fc) * bd.sqrt(1 - (fp / (2*fc))**2)) , bd.zeros_like(fp))
        Iλ = bd.abs(bd.fft.ifft2(bd.fft.ifftshift(c*H)))

        XYZ = F.cs.spec_partition_to_XYZ(bd.outer(Iλ, F.spec_partitions[i]),i)
        sRGB_linear += F.cs.XYZ_to_sRGB_linear(XYZ)

    if backend_name == 'cupy':
        bd.cuda.Stream.null.synchronize()

    F.xx = M_abs * F.xx
    F.yy = M_abs * F.yy
    F.x = M_abs * F.x
    F.y = M_abs * F.y
    F.dx = M_abs * F.dx
    F.dy = M_abs * F.dy

    rgb = F.cs.sRGB_linear_to_sRGB(sRGB_linear)
    rgb = (rgb.T).reshape((F.Ny, F.Nx, 3))
    print ("Computation Took", time.time() - t0)
    return rgb




from diffractsim import MonochromaticField, nm, mm, cm, um



zi = 1*cm # distance from the image plane to the exit pupil
z0 = 1*cm # distance from the exit pupil to the current simulation plane
pupil_radius = 0.5*cm # exit pupil radius

#(If the optical system is a single lens, magnification = - zi/z0)
M = -zi/z0
NA = pupil_radius/z0 #numerical aperture

#maximum resolvable distance by Rayleigh criteria for λ ∼ 550 nm:
print('\n Maximum resolvable distance by Rayleigh criteria: {} nm'.format("%.0f"  % (0.61*550/NA)))



#set up the simulation

F = PolychromaticField(
    spectrum=2*M**2* cf.illuminant_d65, extent_x= 5 * um, extent_y= 5 * um, Nx=2048, Ny=2048,
    spectrum_size = 180, spectrum_divisions = 30
)
F.add(ApertureFromImage( "./apertures/pixel.jpg",  image_size=(5 * um, 5 * um), simulation = F))

#propagate the light assuming the source is spatially incoherent
rgb = get_colors_at_image_plane(F,radius = pupil_radius, zi = zi, z0 = z0, M = M)

F.plot_colors(rgb, figsize=(5, 5))

파일 내용 상세 분석

🔹 1. 라이브러리 및 기본 설정

import diffractsim
diffractsim.set_backend("CUDA")
import numpy as np
from diffractsim import PolychromaticField, ApertureFromImage, cf, mm, cm, CircularAperture
  • CUDA 백엔드를 활성화하여 GPU 가속 최적화.
  • PolychromaticField를 사용하여 다색광 기반의 시뮬레이션 수행.
  • CircularAperture를 사용하여 출구 동공(Exit Pupil) 모델링.

🔹 2. get_colors_at_image_plane() 함수 정의

def get_colors_at_image_plane(F, radius, M,  zi, z0):
  • 비코히런트 광학 시스템에서 최종 이미지를 계산하는 함수.
  • 출구 동공을 고려한 광학 전달 함수(OTF) 적용.
  • 푸리에 변환(FFT) 기반 필터링 수행.

✅ 출구 동공 정의

pupil = CircularAperture(radius)
F.z += zi + z0
  • 출구 동공 반경(radius) 적용.
  • 출구 동공까지의 거리(zi + z0) 반영.

✅ 광학 요소(렌즈, 개구 등)의 투과율 적용

for j in range(len(F.optical_elements)):
    F.E = F.E * F.optical_elements[j].get_transmittance(F.xx, F.yy, 0)
  • 각 광학 요소의 투과율을 적용하여 회절 효과 반영.

✅ 확대 배율(Magnification) 적용 및 이미지 반전

if M < 0:
    F.E = bd.flip(F.E)
M_abs = bd.abs(M)
F.E = F.E/M_abs
  • 확대 배율(Magnification)이 음수면 이미지 반전.
  • 배율에 따라 전기장 F.E 크기 조정.

✅ 광 강도(Intensity) 계산 및 FFT 적용

Ip = F.E * bd.conjugate(F.E)
fft_c = bd.fft.fft2(Ip)
c = bd.fft.fftshift(fft_c)
  • 전기장의 절대값 제곱을 계산하여 광 강도 분포 얻음.
  • 푸리에 변환(FFT) 수행 → 주파수 도메인에서 OTF 적용 가능.

✅ 공간 주파수 계산

fx = bd.fft.fftshift(bd.fft.fftfreq(F.Nx, d = F.x[1]-F.x[0]))/M_abs
fy = bd.fft.fftshift(bd.fft.fftfreq(F.Ny, d = F.y[1]-F.y[0]))/M_abs
fx, fy = bd.meshgrid(fx, fy)
fp = bd.sqrt(fx**2 + fy**2)
  • 공간 주파수(fx, fy)를 계산하여 필터링 적용 가능.

✅ OTF(광학 전달 함수) 적용

fc = radius / (F.λ_list_samples[i]* nm  * zi) # coherent cutoff frequency
H = pupil.get_optical_transfer_function(fx, fy, zi, F.λ_list_samples[i]* nm )
Iλ = bd.abs(bd.fft.ifft2(bd.fft.ifftshift(c*H)))
  • 출구 동공을 통한 광학 전달 함수(OTF) 생성.
  • OTF를 이용한 필터링 수행.

✅ RGB 변환 및 최종 이미지 생성

XYZ = F.cs.spec_partition_to_XYZ(bd.outer(Iλ, F.spec_partitions[i]), i)
sRGB_linear += F.cs.XYZ_to_sRGB_linear(XYZ)
  • 광 강도를 XYZ 색 공간으로 변환 후 sRGB 변환 수행.

🔹 3. 시뮬레이션 파라미터 설정

zi = 40*cm # distance from the image plane to the exit pupil
z0 = 40*cm # distance from the exit pupil to the current simulation plane
pupil_radius = 20*cm # exit pupil radius
M = -zi/z0 # Magnification factor
NA = pupil_radius/z0 # Numerical Aperture
  • 출구 동공(Exit Pupil) 반경 20 cm.
  • 출구 동공에서 이미지 평면까지 거리 40 cm.
  • 확대 배율 M = -zi/z0 = -1 (이미지 반전됨).

🔹 4. Rayleigh 기준을 이용한 해상도 분석

print('\n Maximum resolvable distance by Rayleigh criteria: {} nm'.format("%.0f"  % (0.61*550/NA)))
  • Rayleigh 해상도 공식 적용:
  d_{\text{min}} = \frac{0.61 \lambda}{\text{NA}}

🔹 5. 다색광 시뮬레이션 설정

F = PolychromaticField(
    spectrum=2*M**2* cf.illuminant_d65, extent_x= 14 * um, extent_y= 14 * um, Nx=2048, Ny=2048,
    spectrum_size = 180, spectrum_divisions = 30
)
F.add(ApertureFromImage( "./apertures/horse.png",  image_size=(14 * um, 14 * um), simulation = F))
  • 표준 조명(illuminant_d65)을 사용한 다색광 모델.
  • horse.png 이미지를 개구(Aperture)로 사용.

🔹 6. 최종 이미지 계산 및 출력

rgb = get_colors_at_image_plane(F, radius = pupil_radius, zi = zi, z0 = z0, M = M)
F.plot_colors(rgb, figsize=(5, 5))
  • 출구 동공 반영 후 최종 이미지 계산.
  • 시각화 수행.

 

 

 

추가 기능 개발

    • 스크린에서의 Complex Field 저장 및 시각화
      • field1, field2, field3를 각각 z, z+dz, z+2dz에서 계산.
      • np.save()를 사용하여 .npy 파일 저장.
      • matplotlib을 사용하여 필드의 위상 및 진폭(Amplitude, Phase) 시각화.
    • 기존 코드와 통합
      • get_complex_fields_at_image_plane()가 get_colors_at_image_plane() 방식과 동일하게 Fourier 변환을 사용하여 정확한 스크린 필드 계산.

전체 코드

import diffractsim
import numpy as np
import matplotlib.pyplot as plt
from diffractsim import PolychromaticField, ApertureFromImage, cf, mm, cm, um, CircularAperture

"""
MPL 2.0 License

Copyright (c) 2022, Rafael de la Fuente
All rights reserved.
"""

def get_colors_at_image_plane(F, radius, M, zi, z0):
    from diffractsim.util.backend_functions import backend as bd
    from diffractsim.util.backend_functions import backend_name
    import time
    import progressbar

    pupil = CircularAperture(radius)
    F.z += zi + z0

    for j in range(len(F.optical_elements)):
        F.E = F.E * F.optical_elements[j].get_transmittance(F.xx, F.yy, 0)

    if M < 0:
        F.E = bd.flip(F.E)
    M_abs = bd.abs(M)
    F.E = F.E / M_abs

    Ip = F.E * bd.conjugate(F.E)
    fft_c = bd.fft.fft2(Ip)
    c = bd.fft.fftshift(fft_c)

    fx = bd.fft.fftshift(bd.fft.fftfreq(F.Nx, d=F.x[1] - F.x[0])) / M_abs
    fy = bd.fft.fftshift(bd.fft.fftfreq(F.Ny, d=F.y[1] - F.y[0])) / M_abs
    fx, fy = bd.meshgrid(fx, fy)

    bar = progressbar.ProgressBar()
    sRGB_linear = bd.zeros((3, F.Nx * F.Ny))

    t0 = time.time()

    for i in bar(range(F.spectrum_divisions)):
        fc = radius / (F.λ_list_samples[i] * nm * zi)
        H = pupil.get_optical_transfer_function(fx, fy, zi, F.λ_list_samples[i] * nm)
        Iλ = bd.abs(bd.fft.ifft2(bd.fft.ifftshift(c * H)))

        XYZ = F.cs.spec_partition_to_XYZ(bd.outer(Iλ, F.spec_partitions[i]), i)
        sRGB_linear += F.cs.XYZ_to_sRGB_linear(XYZ)

    if backend_name == 'cupy':
        bd.cuda.Stream.null.synchronize()

    F.xx = M_abs * F.xx
    F.yy = M_abs * F.yy
    F.x = M_abs * F.x
    F.y = M_abs * F.y
    F.dx = M_abs * F.dx
    F.dy = M_abs * F.dy

    rgb = F.cs.sRGB_linear_to_sRGB(sRGB_linear)
    rgb = (rgb.T).reshape((F.Ny, F.Nx, 3))
    print("Computation Took", time.time() - t0)
    return rgb

def get_complex_fields_at_image_plane(F, radius, M, zi, z0, dz=0.2*um):
    """
    Computes the complex vectorial fields at multiple z-planes for FDTD-based EMW simulations.

    Parameters:
        F: PolychromaticField instance
        radius: float - exit pupil radius
        M: float - magnification
        zi: float - distance from the image plane to the exit pupil
        z0: float - distance from the exit pupil to the simulation plane
        dz: float - small z-direction step (default = 0.2 um)

    Returns:
        Saves three .npy files containing complex field data at z, z+dz, and z+2dz.
    """
    from diffractsim.util.backend_functions import backend as bd
    import progressbar
    import time

    pupil = CircularAperture(radius)

    fx = bd.fft.fftshift(bd.fft.fftfreq(F.Nx, d=F.x[1] - F.x[0])) / bd.abs(M)
    fy = bd.fft.fftshift(bd.fft.fftfreq(F.Ny, d=F.y[1] - F.y[0])) / bd.abs(M)
    fx, fy = bd.meshgrid(fx, fy)

    bar = progressbar.ProgressBar()
    t0 = time.time()

    field_z0 = bd.zeros((F.Ny, F.Nx), dtype=bd.complex128)
    field_z0_dz = bd.zeros((F.Ny, F.Nx), dtype=bd.complex128)
    field_z0_2dz = bd.zeros((F.Ny, F.Nx), dtype=bd.complex128)

    for i in bar(range(F.spectrum_divisions)):
        wavelength = F.λ_list_samples[i] * nm
        fc = radius / (wavelength * zi)
        H = pupil.get_optical_transfer_function(fx, fy, zi, wavelength)

        # spec_partitions[i]가 배열일 경우 스칼라 값으로 변환
        weight = F.spec_partitions[i].sum() if isinstance(F.spec_partitions[i], (list, np.ndarray)) else F.spec_partitions[i]

        for z_offset, field_output in zip([0, dz, 2*dz], [field_z0, field_z0_dz, field_z0_2dz]):
            F.z = z0 + zi + z_offset

            for j in range(len(F.optical_elements)):
                F.E *= F.optical_elements[j].get_transmittance(F.xx, F.yy, 0)

            Ip = F.E * bd.conjugate(F.E)
            fft_c = bd.fft.fft2(Ip)
            c = bd.fft.fftshift(fft_c)

            field_output += bd.fft.ifft2(bd.fft.ifftshift(c * H)) * weight  # 스칼라 값으로 변환된 weight 사용

    np.save("complex_field_z0.npy", field_z0)
    np.save("complex_field_z0_dz.npy", field_z0_dz)
    np.save("complex_field_z0_2dz.npy", field_z0_2dz)

    print("Complex fields saved: 'complex_field_z0.npy', 'complex_field_z0_dz.npy', 'complex_field_z0_2dz.npy'")
    print("Computation Took", time.time() - t0)

    return field_z0, field_z0_dz, field_z0_2dz


def plot_complex_field(field, title):
    plt.figure(figsize=(12, 5))
    
    plt.subplot(1, 2, 1)
    plt.imshow(np.abs(field), cmap='inferno')
    plt.colorbar()
    plt.title(f"{title} - Amplitude")

    plt.subplot(1, 2, 2)
    plt.imshow(np.angle(field), cmap='twilight')
    plt.colorbar()
    plt.title(f"{title} - Phase")

    plt.show()


from diffractsim import MonochromaticField, nm, mm, cm, um

zi = 1 * cm
z0 = 1 * cm
pupil_radius = 0.5 * cm
M = -zi / z0
NA = pupil_radius/z0 #numerical aperture

#maximum resolvable distance by Rayleigh criteria for λ ∼ 550 nm:
print('\n Maximum resolvable distance by Rayleigh criteria: {} nm'.format("%.0f"  % (0.61*550/NA)))


F = PolychromaticField(
    spectrum=2 * M**2 * cf.illuminant_d65, extent_x=5 * um, extent_y=5 * um, Nx=2048, Ny=2048,
    spectrum_size=180, spectrum_divisions=30
)

F.add(ApertureFromImage("./apertures/pixel.jpg", image_size=(5 * um, 5 * um), simulation=F))

rgb = get_colors_at_image_plane(F,radius = pupil_radius, M = M, zi = zi, z0 = z0)

F.plot_colors(rgb, figsize=(5, 5))

dz = 0.2 * um
field1, field2, field3 = get_complex_fields_at_image_plane(F, radius=pupil_radius, M=M, zi=zi, z0=z0, dz=dz )

##plot_complex_field(field1, "Field at z")
##plot_complex_field(field2, "Field at z+dz")
plot_complex_field(field3, "Field at z+2dz")
  1.  

 

문제 분석

코드 작성 중에 

ValueError: operands could not be broadcast together with shapes (2048,2048) (6,) 에러 발생.

즉, field_output += bd.fft.ifft2(bd.fft.ifftshift(c * H)) * F.spec_partitions[i] 부분에서 브로드캐스팅 오류.

🔍 원인

  1. F.spec_partitions[i]는 F.spectrum_size(예: 6) 크기의 배열일 가능성이 있음.
  2. bd.fft.ifft2()의 결과는 (2048, 2048)의 2D 배열.
  3. 이 두 배열을 곱하려고 할 때, shape mismatch로 인해 오류 발생.

🔧 해결 방법

  • F.spec_partitions[i]가 배열일 경우, .sum()을 사용하여 스칼라 값으로 변환.
  • 이렇게 하면 field_output += ... 연산에서 shape mismatch 오류가 발생하지 않음.

🔍 질문

  1. F.spec_partitions[i] 는 어떤 역할을 하고, XYZ 는 어떤 역할을 하지. 어찌보면 F.spec_partitions[i] 는 우리가 구하려는 field_output 에 이렇게 weight 형태로 들어가야만 하는 이유에 대한 분석이 필요

🔧 분석 결과

    • F.spec_partitions[i]와 XYZ의 역할F.spec_partitions[i]는 해당 파장(F.λ_list_samples[i])에서의 스펙트럼 가중치(weight) 값입니다.
      • PolychromaticField는 여러 개의 파장(예: 가시광선 전체 스펙트럼)을 포함하는 다색광 시뮬레이션을 수행함.
      • 하지만, 각 파장이 전체 빛의 세기에서 차지하는 비율(가중치)은 다름.
      • F.spec_partitions[i]는 각 파장의 상대적인 기여도를 나타냄.
      • 따라서, 각 파장에 대한 계산 결과(예: Iλ, field_output)는 그 파장이 전체 빛에서 차지하는 비율만큼 가중치를 적용해야 함.
      📌 예제
      • F.spec_partitions[i]는 각 파장의 에너지 분포를 반영하는 30개의 값으로 구성됨.
      • 예를 들어, 특정 파장(예: 550 nm)의 F.spec_partitions[i] 값이 크다면, 해당 파장의 빛이 상대적으로 더 강함을 의미함.
      🔹 2. XYZ의 역할
      • Iλ: 특정 파장 λ에서의 intensity map (즉, 복소수 필드의 절댓값을 구한 것)
      • bd.outer(Iλ, F.spec_partitions[i]): Iλ에 해당 파장의 가중치를 적용한 값
      • F.cs.spec_partition_to_XYZ(...): XYZ 색 공간으로 변환
        • RGB가 아닌, CIE 1931 XYZ 색 공간으로 변환하는 과정
        • XYZ 색 공간을 사용하는 이유는 모든 색을 선형적으로 합산 가능하기 때문
        • 여러 개의 파장을 합쳐서 최종적으로 가시광선 이미지(RGB)로 변환하는 과정에서 필요함
      📌 요약
      • F.spec_partitions[i]는 각 파장의 상대적인 강도(weight)
      • XYZ는 각 파장에서의 강도를 조합하여 최종적으로 RGB 변환을 하기 위한 중간 단계.
      🔍 F.spec_partitions[i]가 field_output에 들어가야 하는 이유
       
      field_output += bd.fft.ifft2(bd.fft.ifftshift(c * H)) * F.spec_partitions[i]
      를 사용한 이유는 각 파장의 빛이 전체 빛에서 차지하는 비율을 반영해야 하기 때문입니다. 이렇게 하지 않으면, 
      • 모든 파장을 동일한 강도로 가정하게 됨.
      • 실제 조명(예: 태양광)에서는 각 파장의 강도가 다르기 때문에 물리적으로 잘못된 결과가 나올 수 있음.
      • 예를 들어, 550 nm(녹색) 빛이 강하고 650 nm(적색) 빛이 약한 조명이라면, F.spec_partitions[i]를 적용해야 이를 반영할 수 있음.
      📌 올바른 가중치를 적용한 계산
      1. 각 파장별 복소수 필드(회절 결과)를 계산
      2. 각 필드에 F.spec_partitions[i] 곱해서 최종 합산
      즉, F.spec_partitions[i]는 각 파장의 상대적인 강도를 반영하여 전체 필드의 정확도를 높이는 역할을 합니다.

      📌 spectrum_size vs. spectrum_divisions 차이

      PolychromaticField를 정의할 때 사용되는 spectrum_size와 spectrum_divisions는 폴리크로매틱(다색) 시뮬레이션에서 파장 샘플링 방식을 제어하는 두 개의 다른 변수입니다.

      🔹 spectrum_size란?

      spectrum_size는 전체 사용 가능한 스펙트럼에서 선택할 수 있는 파장의 개수를 의미합니다.

      📌 역할

      • spectrum_size는 전체 스펙트럼 범위(예: 400 nm ~ 700 nm)에서 몇 개의 파장을 사용할지를 결정합니다.
      • 일반적으로, spectrum_size가 클수록 더 세밀한 스펙트럼 표현이 가능하지만, 계산량이 증가합니다.

      📌 예제

      python
      CopyEdit
       
      F = PolychromaticField( spectrum=cf.illuminant_d65, spectrum_size=100, # 총 100개의 파장을 샘플링 spectrum_divisions=10 )
      • 100개의 샘플 파장을 선택하지만, 실제 spectrum_divisions에 따라 그룹화될 수 있음.

      🔹 spectrum_divisions란?

      spectrum_divisions는 계산을 수행할 때 실제로 몇 개의 파장 샘플을 사용할지를 지정합니다.

      📌 역할

      • spectrum_divisions는 spectrum_size에서 선택된 파장을 다시 그룹화하여, 실제 시뮬레이션에서 몇 개의 파장을 사용할지를 결정합니다.
      • 더 작은 값으로 설정하면 계산량이 줄어들고, 더 큰 값이면 더욱 정밀한 결과를 얻을 수 있습니다.

      📌 예제

       
      F = PolychromaticField( spectrum=cf.illuminant_d65,
      spectrum_size=100, # 100개 파장을 샘플링
      spectrum_divisions=20 # 실제 계산에서는 20개의 파장을 사용 )
      • 100개의 파장 중에서 20개를 실제 시뮬레이션에서 사용.

      🔹 spectrum_size vs. spectrum_divisions 비교 요약

      Parameter역할영향
      spectrum_size 전체 사용 가능한 파장의 개수 크면 더 다양한 파장을 포함하지만, 반드시 모든 파장이 계산에 사용되지는 않음
      spectrum_divisions 실제 시뮬레이션에서 사용할 파장의 개수 작으면 계산량이 줄지만 정확도가 낮아짐, 크면 더욱 정밀한 스펙트럼 분석 가능

      🔹 최적의 설정 방법

      • spectrum_size는 **충분히 큰 값(예: 100~200)**을 사용하여 다양한 파장을 포함하도록 설정.
      • spectrum_divisions는 계산 성능과 정확도의 균형을 고려하여 선택. (예: 20~50 정도)
      • 정확도가 중요한 경우: spectrum_size = 200, spectrum_divisions = 50
      • 빠른 계산이 필요한 경우: spectrum_size = 100, spectrum_divisions = 10

       

       

z+dz 위치에 따른 Field Amplitude 와 Phase

 

🔹 TM/TE 모드 적용 방법

■ 편광의 기본 개념

다음은 TM/TE 모드를 적용한 전체 코드와 이에 대한 설명입니다.
이제 광원의 편광 방향에 따라 전기장(E-field) 및 자기장(H-field)의 변화를 반영하여, TE(Transverse Electric) / TM(Transverse Magnetic) 모드를 구현합니다.

  • TE 모드 (Transverse Electric) → 전기장(𝑬)이 입사면과 수직
  • TM 모드 (Transverse Magnetic) → 자기장(𝑯)이 입사면과 수직

■ 편광 적용 방식

  • TE 모드 → 전기장의 y-성분을 고려 (Ez = 0, Ey ≠ 0)
  • TM 모드 → 전기장의 x-성분을 고려 (Ez = 0, Ex ≠ 0)
  • Fresnel 계수를 사용하여, 편광에 따라 필드 강도를 조정

■ 적용 코드 개요

  • apply_polarization(F, theta, polarization="TE") 함수를 정의하여 TE/TM 모드에 따라 E-field을 조정
  • 기존 코드에서 get_complex_fields_at_image_plane() 및 get_colors_at_image_plane()에 편광 적용

✅ TE/TM 모드 적용된 전체 코드

 
import diffractsim
import numpy as np
import matplotlib.pyplot as plt
from diffractsim import PolychromaticField, ApertureFromImage, cf, mm, cm, um, CircularAperture

"""
MPL 2.0 License

Copyright (c) 2022, Rafael de la Fuente
All rights reserved.
"""

def apply_polarization(F, theta, polarization="TE"):
    """
    Apply TE/TM polarization to the field.

    Parameters:
        F: PolychromaticField instance
        theta: float - Incident angle θ (radians)
        polarization: str - "TE" or "TM"

    Modifies:
        F.E - Adjusts field based on polarization
    """
    from diffractsim.util.backend_functions import backend as bd

    cos_theta = np.cos(theta)

    if polarization == "TE":
        # TE polarization → E-field (Ey) is dominant
        F.E *= cos_theta  
    elif polarization == "TM":
        # TM polarization → E-field (Ex) is dominant
        F.E *= 1 / cos_theta  


def get_colors_at_image_plane(F, radius, M, zi, z0, polarization="TE"):
    from diffractsim.util.backend_functions import backend as bd
    from diffractsim.util.backend_functions import backend_name
    import time
    import progressbar

    pupil = CircularAperture(radius)
    F.z += zi + z0

    apply_polarization(F, 0, polarization)  # Apply polarization

    for j in range(len(F.optical_elements)):
        F.E = F.E * F.optical_elements[j].get_transmittance(F.xx, F.yy, 0)

    if M < 0:
        F.E = bd.flip(F.E)
    M_abs = bd.abs(M)
    F.E = F.E / M_abs

    Ip = F.E * bd.conjugate(F.E)
    fft_c = bd.fft.fft2(Ip)
    c = bd.fft.fftshift(fft_c)

    fx = bd.fft.fftshift(bd.fft.fftfreq(F.Nx, d=F.x[1] - F.x[0])) / M_abs
    fy = bd.fft.fftshift(bd.fft.fftfreq(F.Ny, d=F.y[1] - F.y[0])) / M_abs
    fx, fy = bd.meshgrid(fx, fy)

    bar = progressbar.ProgressBar()
    sRGB_linear = bd.zeros((3, F.Nx * F.Ny))

    t0 = time.time()

    for i in bar(range(F.spectrum_divisions)):
        fc = radius / (F.λ_list_samples[i] * nm * zi)
        H = pupil.get_optical_transfer_function(fx, fy, zi, F.λ_list_samples[i] * nm)
        Iλ = bd.abs(bd.fft.ifft2(bd.fft.ifftshift(c * H)))

        XYZ = F.cs.spec_partition_to_XYZ(bd.outer(Iλ, F.spec_partitions[i]), i)
        sRGB_linear += F.cs.XYZ_to_sRGB_linear(XYZ)

    if backend_name == 'cupy':
        bd.cuda.Stream.null.synchronize()

    F.xx = M_abs * F.xx
    F.yy = M_abs * F.yy
    F.x = M_abs * F.x
    F.y = M_abs * F.y
    F.dx = M_abs * F.dx
    F.dy = M_abs * F.dy

    rgb = F.cs.sRGB_linear_to_sRGB(sRGB_linear)
    rgb = (rgb.T).reshape((F.Ny, F.Nx, 3))
    print("Computation Took", time.time() - t0)
    return rgb


def plot_complex_field(field, title):
    plt.figure(figsize=(12, 5))
    
    plt.subplot(1, 2, 1)
    plt.imshow(np.abs(field), cmap='inferno')
    plt.colorbar()
    plt.title(f"{title} - Amplitude")

    plt.subplot(1, 2, 2)
    plt.imshow(np.angle(field), cmap='twilight')
    plt.colorbar()
    plt.title(f"{title} - Phase")

    plt.show()


from diffractsim import MonochromaticField, nm, mm, cm, um

zi = 1 * cm
z0 = 1 * cm
pupil_radius = 0.5 * cm
M = -zi / z0
NA = pupil_radius / z0  

# Maximum resolvable distance by Rayleigh criteria for λ ∼ 550 nm:
print('\n Maximum resolvable distance by Rayleigh criteria: {} nm'.format("%.0f"  % (0.61*550/NA)))


F = PolychromaticField(
    spectrum=2 * M**2 * cf.illuminant_d65, extent_x=5 * um, extent_y=5 * um, Nx=2048, Ny=2048,
    spectrum_size=180, spectrum_divisions=30
)

F.add(ApertureFromImage("./apertures/pixel.jpg", image_size=(5 * um, 5 * um), simulation=F))

# Compute the RGB image for TE and TM polarization
rgb_TE = get_colors_at_image_plane(F, radius=pupil_radius, M=M, zi=zi, z0=z0, polarization="TE")
F.plot_colors(rgb_TE, figsize=(5, 5))

rgb_TM = get_colors_at_image_plane(F, radius=pupil_radius, M=M, zi=zi, z0=z0, polarization="TM")
F.plot_colors(rgb_TM, figsize=(5, 5))

🔹 TM/TE 모드 적용 후 결과

■ 편광 적용 전
➜ 기존에는 모든 입사광을 동일한 방식으로 처리하여 편광 효과를 반영하지 않았음.
➜ TE/TM 모드가 없기 때문에 특정 편광 상태의 필터링이 불가능.

■ 편광 적용 후
➜ TE 모드: 전기장이 입사면과 수직이므로, 특정 방향의 회절 패턴이 강화됨.
➜ TM 모드: 전기장이 입사면과 평행하므로, 필드 강도가 다르게 변형됨.
➜ 같은 패턴을 가지지만, TE/TM에 따라 밝기나 위상 분포가 다르게 나타남.

🔹 최종 정리

■ 편광을 반영하면 TE/TM 모드별로 필드 강도가 다르게 나타남
■ apply_polarization(F, theta, polarization)을 적용하여 편광 상태를 반영
■ get_colors_at_image_plane()에서 TE/TM 편광 모드별 필드 계산 가능
■ 이제 TE/TM 편광 상태에 따른 광학 현상 분석이 가능함! 🚀

  •  

 

'Camera' 카테고리의 다른 글

EMW Simulation 을 위한 Code V input file 작성  (0) 2025.02.23
Numerical Aperture (수치 개구)란 무엇인가?  (0) 2025.02.20
Random 위상을 갖는 사각 Aperture 의 렌즈 모델  (0) 2025.02.19
85mm f/1.4 vs 200mm f/2.8  (0) 2025.02.12
Star Chart 로 MTF 측정시 유의점  (0) 2025.02.12

관련글

티스토리툴바