EMW Simulation 을 위한 Code V input file 작성

 

■ Code V Output을 EMW에 입력 및 적용 방법

◇ Code V Interface 개요

• Code V에서 생성한 복잡한 벡터 전기장 데이터를 EMW에서 입력으로 사용할 수 있도록 함.
• TFSF (Total-Field Scattered-Field) 공식을 이용하여 Code V에서 가져온 데이터를 EMW 시뮬레이션에 적용.

◇ Code V 데이터를 EMW에 입력하는 방법

1. CodeVExcitation 섹션 설정
   • Code V의 출력 데이터를 EMW에서 사용할 수 있도록 .dat 파일을 입력.
   • 필수적인 입력 파일: 두 개의 복소 전기장 데이터 파일 (InputDataFile)

예시:

CodeVExcitation { InputDataFile = {"zf.dat", "zfPlus20.dat"}
AnchorPoint1 = {-1, -1, 0.28} # [um]
AnchorPoint2 = {-1, -1, 0.30} # [um]
Wavelength = 300 # [nm]
Signal = Harmonic
NRise = 20
Delay = 0
AmplitudeScaling = 1 PhaseShift = 0 }

1. • AnchorPoint1과 AnchorPoint2는 두 데이터 슬라이스의 위치를 지정.
2. Code V 데이터 파일 형식

.dat 파일 형식:

Wavelength: 300 nm
Grid spacing: 0.000020
0.000020 mm
Array size: 101 101

 

 

이후 각 좌표별 전기장 데이터가 나열됨:

ExRe(x1,y1) ExIm(x1,y1) EyRe(x1,y1) EyIm(x1,y1) EzRe(x1,y1) EzIm(x1,y1)
ExRe(x2,y1) ExIm(x2,y1) EyRe(x2,y1) EyIm(x2,y1) EzRe(x2,y1) EzIm(x2,y1)
...

 

 

1. Code V 데이터를 사용할 때 고려해야 할 사항
   • 지원하지 않는 기능:
     • 2D 시뮬레이션에서는 지원되지 않음.
     • z 방향 이외의 광 전파 방향은 지원되지 않음.
     • 가우시안 신호 등 비조화(Harmonic이 아닌) 신호는 사용할 수 없음.
   • 제약 사항:
     • AnchorPoint1, AnchorPoint2의 x, y 좌표는 동일해야 하며, z 좌표 차이는 EMW의 z 방향 셀 크기와 같아야 함.
     • 입력된 두 슬라이스의 데이터가 시뮬레이션 영역의 xy 경계를 초과하면 회절(Diffraction) 오류 발생 가능.
2. 지원되는 경계 조건 조합
   • 아래 표는 Code V 데이터를 사용할 때 허용되는 경계 조건 조합을 나타냄:xminxmaxyminymaxzminzmax

     CPML	CPML	CPML	CPML	CPML	CPML
     Periodic	Periodic	Periodic	Periodic	CPML	CPML

---

📌 요약

• Code V의 출력 데이터를 EMW에서 사용할 때 CodeVExcitation 섹션을 설정해야 함.
• 입력 파일은 두 개의 전기장 데이터 파일을 포함하며, 슬라이스의 위치를 AnchorPoint1, AnchorPoint2로 지정.
• Code V 데이터는 3D 시뮬레이션에서만 사용 가능하며, 광 전파 방향은 z축을 따라야 함.
• CPML 및 Periodic 경계 조건 조합이 필수이며, x, y 좌표가 동일한 두 슬라이스의 z축 차이를 EMW의 셀 크기와 맞춰야 함.

 

import diffractsim
##diffractsim.set_backend("CUDA")
import numpy as np
import matplotlib.pyplot as plt
from diffractsim import PolychromaticField, ApertureFromImage, cf, mm, cm, um, CircularAperture
import os
from diffractsim.util.backend_functions import backend as bd
import progressbar

"""
MPL 2.0 License

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


def save_field_to_file(filename, field, wavelength, grid_spacing, array_size):
    """
    Save field data in EMW format.

    Parameters:
        filename: str - Output file name
        field: numpy array - Field data with shape (Nx, Ny, 6)
        wavelength: float - Wavelength in nm
        grid_spacing: float - Grid spacing in mm
        array_size: int - Size of the field array
    """
    with open(filename, "w") as f:
        f.write(f"Wavelength: {wavelength:.1f} nm\n")
        f.write(f"Grid spacing: {grid_spacing:.6f} {grid_spacing:.6f} mm\n")
        f.write(f"Array size: {array_size} {array_size}\n\n")

        for y in range(field.shape[0]):
            for x in range(field.shape[1]):
                data_line = " ".join(f"{field[y, x, i]:.6e}" for i in range(6))
                f.write(data_line + "\n")

def save_emw_input_files(F, radius, M, zi, z0, dz=0.02*um):
    """
    Compute complex vectorial fields at two z-planes for EMW and save them in zf.dat and zfPlus20.dat.

    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 - step size in z-direction (default = 0.02 um)

    Saves:
        "zf.dat" - Field data at (z0 + zi - 1um)
        "zfPlus20.dat" - Field data at (z0 + zi - 1um + 0.02um)
    """
    pupil = CircularAperture(radius)
    
    # zf 위치 설정 (1 µm 앞쪽)
    zf = z0 + zi - 1 * um
    zf_plus = zf + dz  # zf + 0.02 um

    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()
    field_zf = bd.zeros((F.Ny, F.Nx, 6), dtype=bd.float32)
    field_zf_plus = bd.zeros((F.Ny, F.Nx, 6), dtype=bd.float32)

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

        for z_pos, field_output in zip([zf, zf_plus], [field_zf, field_zf_plus]):
            F.z = z_pos
            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)

            E_field = bd.fft.ifft2(bd.fft.ifftshift(c * H))

            field_output[:, :, 0] = bd.real(E_field)  # Ex 실수부
            field_output[:, :, 1] = bd.imag(E_field)  # Ex 허수부
            field_output[:, :, 2] = bd.real(E_field)  # Ey 실수부
            field_output[:, :, 3] = bd.imag(E_field)  # Ey 허수부
            field_output[:, :, 4] = bd.real(E_field)  # Ez 실수부
            field_output[:, :, 5] = bd.imag(E_field)  # Ez 허수부

    # 파일 저장
    wavelength_nm = np.mean(F.λ_list_samples)   # 중앙 파장 (nm)
    grid_spacing_mm = (F.x[1] - F.x[0]) * 1e3  # mm 단위 변환
    array_size = F.Nx

    save_field_to_file("zf.dat", field_zf, wavelength_nm, grid_spacing_mm, array_size)
    save_field_to_file("zfPlus20.dat", field_zf_plus, wavelength_nm, grid_spacing_mm, array_size)

    print("Files saved: 'zf.dat', 'zfPlus20.dat'")

    # 필드 데이터 RGB 변환 및 시각화
    plot_rgb_field(field_zf, title="ZF - 1 µm 앞")
    plot_rgb_field(field_zf_plus, title="ZF + 0.02 µm 앞")

    return field_zf, field_zf_plus

def plot_rgb_field(field, title="RGB Field Visualization"):
    """
    Convert the complex field to an RGB image and plot it.

    Parameters:
        field: numpy array of shape (Nx, Ny, 6) containing complex field components.
        title: str - Title for the plot.
    """
    # 진폭을 사용하여 RGB 변환
    amplitude = np.sqrt(field[:, :, 0]**2 + field[:, :, 1]**2)  # Ex 실수부와 허수부의 진폭 사용
    amplitude = amplitude / np.max(amplitude)  # 정규화

    # RGB 색 변환 (임의 매핑, 필요시 변경 가능)
    r = amplitude
    g = np.roll(amplitude, shift=1, axis=0)  # y축 기준 shift
    b = np.roll(amplitude, shift=1, axis=1)  # x축 기준 shift

    rgb_image = np.stack([r, g, b], axis=-1)

    # 플로팅
    plt.figure(figsize=(6, 6))
    plt.imshow(rgb_image, cmap='inferno')
    plt.colorbar()
    plt.title(title)
    plt.show()

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.data_550nm, extent_x= 5 * um, extent_y= 5 * um, Nx=512, Ny=512,
    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,  M = M, zi = zi, z0 = z0)

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

# EMW 입력 파일 생성 및 저장
save_emw_input_files(F, radius = pupil_radius, M = M, zi = zi, z0 = z0, dz=0.02*um)

1.