아래 예제 simulation은 point source가 $C_n^2=1\times10^{-16}m^{-2/3}$인 거리 $\triangle z=50km$를 turbulent path를 통과하는 것이다. 단순함을 위해, Kolmogorv refractive-index Power Spectral Desnity 함수로도 충분함을 가정하였고, 관측 망원경의 지름은 $D_2=0.5m$이다.

1. Setting

1.1 Set up source and receiver geometry

parameter	value	description
$D_2$	0.5m	observation 구경 지름
$\lambda$	10^{-6}m	파장
$k$	$\cfrac{2\pi}{\lambda}$	파수
$\triangle z$	50km	전파거리
$D_{ROI}$	2m	observation 관심 영역 지름
$R$	50km	wavefront의 구면 반지름
$D_1$	$\cfrac{\lambda\ \triangle z}{D_{ROI}}$	central lobe 너비($\theta_{3dB}$)

D2 = 0.5     # diameter of the observation aperture [m]
wavelength = 1e-6 # optical wavelength [m]
k = 2*np.pi/wvl # optical wavenemumber [rad/m]
Dz = 50e3 # propagation distance [m]

DROI = 4 * D2 # diameter of observation plane region of intersets [m]
D1 = wavelength * Dz / DROI # width of central lobe [m]
R = Dz  # wavefront radius of curvature [m]

위 정보를 바탕으로 관심있는 atmospheric 변수를 계산할 수 있다. 물론, 이 변수들은 빛이 전파된 후 우리가 imaging, wavefront sensing, adaptive optics 등 뭐하는지에 따라 달려 있다.

point source의 central lobe의 너비로부터 $D_1$을 계산한다.
$D_1=\cfrac{\lambda*\triangle z}{D_{ROI}}$
이것은 source에 의해 observation 평면에서 균일하게 조명되는 관심영역의 직경을 설정하는 것에서 시작한다.

1.2. Set up turbulence related quantities

effective coherence diamter $r_0$

• $r_{0,sw}=[0.423k^2\sum\limits_{i=1}^NC_{n_i}^2 (\cfrac{z_i}{\triangle z})^{5/3} \triangle z_i]^{-3/5}=(0.423k^2C_n^2\ \triangle z \ \cfrac{1}{N}\sum\limits_{i=1}^{N}(\cfrac{ z_i}{\triangle z})^{5/3})^{-3/5}$
• $r_{0, pw}=[0.423k^2\sum\limits_{i=1}^N C_{n_i}^2\triangle z_i]^{-3/5}=[0.423k^2C_n^2\triangle z]^{-3/5}$

rytov variance $\sigma^2_{\chi}$ (log-amplitude variance)

• $\sigma^2_{\chi, sw}=0.563k^{7/6}\triangle z^{5/6}\sum\limits^N_{i=1}C^2_{n_i}(\cfrac{z_i}{\triangle z})^{5/6}(1-\cfrac{z_i}{\triangle z})^{5/6}\triangle z_i=0.563k^{7/6}C_n^2(z_2-z_1) \sum\limits_{i=1}^Nz_i^{5/6}(1-\cfrac{z_i}{\triangle z})^{5/6}$

# atmospheric properties
Cn2 = 1e-16
# SW and PW coherence diameters [m]
r0_sw = (0.423 * k**2 * Cn2 * 3/8 * Dz)**(-3/5)
r0_pw = (0.423 * k**2 * Cn2 * Dz)**(-3/5)
p = np.linspace(0, Dz, int(1e3)+1) # 0~Dz를 1001개의 평면으로 만듬

# log-amplitude variance
rytov = 0.563 * k**(7/6) * np.sum(Cn2 * (1-p/Dz)**(5/6) * p**(5/6) * (p[2]-p[1])) 

• $r_{0,sw}=[0.423k^2C^2_n\int\limits^{\triangle z}_0(\cfrac{z}{\triangle z})^{5/3}dz]^{-3/5}$
• $\sigma^2_{\chi,sw}=0.563k^{7/6}C^2_n\int\limits^{\triangle z}_0z^{5/6}(1-\cfrac{z}{\triangle z})^{5/6}dz$
  핵심 대기 매개변수 $r_{0,sw}=12.7cm$,
  $\sigma_{\chi,sw}^2=0.436$이다.

1.3 Set up screen properties

$\alpha_i=\cfrac{z_i}{\triangle z}$, $r_{0_i}=[0.423k^2C^2_{n_i}\triangle z_i]^{-3/5}$로 표현하면,

• $r_{0,sw}^{-5/3}=(\sum\limits_{i=1}^{N} r_{0_i}^{-5/3} \alpha_i^{5/3})$
• $\sigma^2_{\chi,sw}=1.33k^{-5/6}\triangle z^{5/6}\sum\limits^N_{i=1}r_{0_i}^{-5/3}\alpha_i^{5/6}(1-\alpha_i)^{5/6}$

• $\cfrac{\sigma^2_{\chi, sw}}{1.33} (\cfrac{k}{\triangle z})^{5/6}=\sum\limits^N_{i=1}r_{0_i}^{-5/3}\alpha_i^{5/6}(1-\alpha_i)^{5/6}$

여기에서 각 phase screen의 적절한 $r_0$ 값들을 계산하기 위해 행렬을 이용한다.
$r_0^{-5/3}$ 의 경우 접근법이 어려운데, 해 $r_0$ 벡터에서 음수 항목들은 비물리적이기 때문에, $r_0>0$으로 제한해야 하고, 여러 phase screen이 있는 시뮬레이션을 위해 $r_0$ 값을 계산하는데 constrained optmization을 사용해야 한다.

n_screen=11
A = np.zeros((2,n_screen))
alpha = np.arange(0, n_screen) / (nscreen-1)

A[0][:] = alpha**(5/3)
A[1][:] = (1-alpha)**(5/6) * alpha**(5/6)

b = np.array([[r0_sw**(-5/3)], [rytov/1.33*(k/Dz)**(5/6)]])

$\begin{pmatrix} r_{0,sw}^{-5/3} \\ \frac{\sigma_{\chi,sw}^2}{1.33}(\frac{k}{\triangle z})^{5/6} \end{pmatrix}=\begin{pmatrix} \alpha_1^{5/3} & ...& \alpha_N^{5/3} \\ \alpha_1^{5/6}(1-\alpha_1)^{5/6}& ...& \alpha_N^{5/6}(1-\alpha_N)^{5/6} \end{pmatrix} \begin{pmatrix} r_{0_1}^{-5/3} \\ r_{0_2}^{-5/3} \\ ... \\ r_{0_N}^{-5/3} \\ \end{pmatrix}$

최소화해야할 목적 함수

# objective function
def obj_func(X):
    # [2 x n_screen] x [n_screen x 1] - [2x1] = [2x1]
    return np.sum((A@X-b)**2)

위의 $\mathbf{b}=A\mathbf{x}$ 행렬을 통해 $r_{0_i}$값을 구하기 위해서, 초기 $\mathbf{x}_0$ 값을 다음으로 설정한다.
$\mathbf{x}_0=\begin{pmatrix} \cfrac{N}{3}r_{0,sw}^{-5/3} \\ ... \\ \cfrac{N}{3}r_{0,sw}^{-5/3} \\ \end{pmatrix}$

 # initial guess
x0 = (n_screen/3 * r0_sw * np.ones((n_screen,)))**(-5/3) # [ n_screen, ]

• $\mathbf{x}$의 상한선을 $\mathbf{x}_2$라고 지칭하면, 각 부분 전파할 때마다 rytov variance의 최대값을 0.1이라고 하면(Martin and Flatte에 의해 제안된 guideline과 관련있다.)
  $A[1]\mathbf{x}_2=\cfrac{\sigma_{\chi,sw}^2}{1.33}(\cfrac{k}{\triangle z})^{5/6}$
  $\mathbf{x}_{2}=A[1]^{-1}\cfrac{\sigma_{\chi,sw}^2}{1.33}(\cfrac{k}{\triangle z})^{5/6}$ 값을 상한선으로 잡을 수 있다.
• $\mathbf{x}_1=0$ 행렬로 하한선으로 설정, $r_0$의 무한 screen에 대응한다.

x1 = np.zeros((n_screen,))

r_max = 0.1  # maximum Rytov number per partial prop
x2 = r_max/1.33*(k/Dz)**(5/6) / A[1][:] # [ n_screen, ]


for idx, a in enumerate(A[1][:]):
    if a==0.0:
        x2[idx] = 50**(-5./3)
        
x2 = x2.reshape(-1,) # [n_screen, ]

bounds = []
for bound_per_elem in zip(x1,x2):

Sequential Least SQuares Programming 방식을 이용하여

• 목적함수, 제약조건 함수를 quadratic 형태로 이차 근사하고, 이를 반복적으로 풀어 최적해를 찾는다.
• 각 반복마다 선형 검색을 사용해, 새로운해를 찾고 해가 제약 조건(bound, constraints) 만족하는 확인

# 목적함수, 제약조건 함수를 quadratic 형태로 근사. 이를 반복적으로 풀어 최적해 찾기.
# 각 반복에서 선형 검색을 사용하여 새로운 해를 찾고, 이 해가 제약 조건을 만족하는지 확인
result = minimize(obj_func, x0, method='SLSQP', bounds=bounds)
X = result.x
#print(result.fun)

# check screen r0_i
r0_screen = X**(-3/5)
r0_screen[np.isinf(r0_screen)] = 1e6

# Checking resulting r0_sw, & rytov
bp = A@X

r0sw = bp[0]**(-3/5)
rytov = bp[1]*1.33*(Dz/k)**(5/6)

print("r0sw:", r0sw)
print("rytov:", rytov)

---

2. Analyze the sampling constraints

위에 구조와 turbulence 조건이 설정되면, 그 다음으로 grid spacings와 grid point의 수를 결정하기 위해 sampling constraint들을 분석할 수 있다.

그다음 constraint 1~3의 경계조건을 구해야 한다.
1. $\delta_n \le \cfrac{\lambda \triangle z-D'_2\delta_1}{D_1'}$
2. $N\ge\cfrac{D_1'}{2\delta_1}+\cfrac{D_2'}{2\delta_n}+\cfrac{\lambda \triangle z}{2\delta_1 \delta_n}$
3. $(1+\cfrac{\triangle z}{R})\delta_1-\cfrac{\lambda\triangle z}{D_1}\le\delta_2\le(1+\cfrac{\triangle z}{R})\delta_1+\cfrac{\lambda\triangle z}{D_1}$

이렇게 $N, \delta_1, \delta_n$이 선택되면, partial propagation에 대해

• $\triangle z_{max}=\cfrac{min(\delta_1, \delta_n)^2 N}{\lambda}$
• $n_{min}=ceil(\cfrac{\triangle z}{\triangle z_{max}})+1$

종합적인 constraint의 contour plot은 다음 그림과 같다. 그림에서 constraint 2에서 N의 하한과 constraint 1,3의 상한을 겹쳐서 보여준다.

# analysis_pt_source_atmosphere_sampling
c = 2
D1p = D1 + c * wavelength * Dz / r0_sw
D2p = D2 + c * wavelength * Dz / r0_sw

delta1 = np.linspace(0, 1.1*wavelength*Dz/D2p, 100)
deltan = np.linspace(0, 1.1*wavelength*Dz/D1p, 100)

# constraint 1
deltan_max = -D2p / D1p*delta1 + wavelength * Dz / D1p
# constraint 3
d2min3 = (1+Dz/R)*delta1 - wavelength*Dz/D1p
d2max3 = (1+Dz/R)*delta1 + wavelength*Dz/D1p

delta1, delta_n = np.meshgrid(delta1, deltan)

# constraint 2
N2 = (wavelength * Dz + D1p*deltan + D2p*delta1)/(2*delta1*deltan)

# constraint 4
d1 = 1e-3
d2 = 1e-3
N = 512
d2 = d1*d2*N/wavelength

z_max = np.min([d1, d2])**2 * N / wavelength
n_min = np.ceil(Dz/z_max)+1

print("maximum allowed propagation distance:",z_max)
print("minimum number of grid point:", n_min)

---

3. Perform a vaccum simulation

grid parameter $N, \delta_1, \delta_n$가 결정되면 다음 과정으로 진공 환경에서 시뮬레이션을 수행한다. 2가지 중요한 목적이 있는데,

1. simulation이 turbulence 없이 정확한 결과를 만들어내는지
   • 이 경우, point source로 시뮬레이션한다면, 알려진 분석된 해와 진공 시뮬레이션 결과가 일치하는지 확인해야 한다.
2. 진공 시뮬레이션과 turbulent 시뮬레이션을 비교한다. 그래야 얼마나 turbulence에 의해 degraded 되는지 성능을 알 수 있다. 예시로 Strehl ratio를 계산한다.

3.1 point source의 vaccum simulation

$U_{pt}(x,y)=e^{-i\frac{k}{2R}(x_1^2+y_1^2)}*\cfrac{1}{D_1^2}*sinc(x_1/D_1)*sinc(y_1/D_1)*e^{-\frac{x_1^2+y_1^2}{4D_1}}$

---

4. Perform the turbulent simulations

시뮬레이션이 올바르게 동작하는지 확인하기 위해, source를 구현된 수많은 turbulence를 통해 source를 전파시키고 coherence factor를 계산한후, 이론적 기대치에 대해 그래프를 그린다. 또한 phase screen들의 위치와 그들의 coherence 지름을 결정에 필요할가 있다.

마침내, phase screen의 구현들로 난류 시뮬레이션을 수행할 수 있다. 아래 코드에서는 정확한 grid spacing에서 11개의 phase screen을 생성하여 1개의 난란 경로(realization)를 만들고 전파를 시뮬레이션한다.
이 과정을 40번 반복하여 독립적이고 동일하게 분포된 대기를 통해 광학 필드의 40개 실현을 생성한다.
이러한 다양한 realization들을 수집함으로써 coherence factor, 파 구조 함수 및 로그 진폭 분산과 같은 앙상블 통계량을 추정할 수 있다.

만약 우리가 동적으로 발전하는 대기를 시뮬레이션하고 싶다면, 각 대기의 구현마다 phase screen을 시간이 경과함에 따라 횡단 차원에서 이동시켜야 한다. 이는 Taylor frozen-turbulence 가정을 명시적으로 활용하는 것.
15번 가정대로 screen의 속도는 Greenwood 주파수와 같은 시간적인 양을 통해 결정되어야 한다. 이는 시뮬레이션의 시간적 특성을 확인하고 조정 광학과 같은 adaptive optics 시스템에서 시뮬레이션을 사용할 수 있도록 한다.