광 전파에서 대기 turbulence 효과를 분석적 결과로 도출하는 것은 간단한 통계적 모델을 가정할 때 가능하다. 하지만, 적응 광학 시스템을 사용하여 더 복잡한 시나리오를 고려하는 것을 원하는 경우, 보통 보정된 optical field의 통계는 closed form에서 계산될 수 없다.
수식 단순화를 위해, 흔한 테크닉은 turbulence를 유한한 이산적인 layer로 다루는 것이다. 이 분석적 계산, simulation, turbulence를 실험실에서 모방할 때 흔한 방법이다.

layered model은 해당 extended 매질의 굴절률 spectrum과 scintillation 특성에 맞다면 유용하다. 각 layer는 더욱 큰 두께의 turbulent 부피를 표현하는 얇은 phase screen 단일 진폭이다. 만약 phase screen의 두께가 screen 향하는 전파 거리보다 훨씬 작다면, phase screen은 얇게 고려된다. phase screen은 대기 phase perturbation(작은 변화)의 하나 구현이고 $\mathcal{T}[z_i, z_{i+1}]=\exp[-\phi(\mathbf{r}_{i+1})]$ 굴절 연산자를 구현할 수 있다. 이것은 대기 전파를 시뮬레이션 하기 위해 대기 phase screen을 split-step beam 전파 방법을 어떻게 통합하는지 보여준다.

1. split-step beam propagation

진공이 아닌 매질에 대한 전파를 시뮬레이션은 split-step beam propagation method를 통해 달성할 수 있는데, 많은 material 종류들에서 유용하다. (불균등하고, 이방성, 비선형적)

여기서는 대기에 국한하여 선형적이고, 등방적인 물질인데 불균등한 굴절률(n)을 $n=n(x,y,z)$에 대해서 논한다. $\delta n=n-1$이 작을 때,
i+1번째 평면에서 보이즌 field는
$U(\mathbf{r}_{i+1})\simeq\mathcal{R}[\cfrac{\triangle z_i}{2}, \mathbf{r}_i,\tilde{\mathbf{r}}_{i+1}]\ \mathcal{T}[z_i, z_{i+1}]\ \mathcal{R}[\cfrac{\triangle z_i}{2}, \mathbf{r}_i,\tilde{\mathbf{r}}_{i+1}]\ \{U(\mathbf{r}_i)\}$

$\mathcal{T}[z_i, z_{i+1}]$은 phase의 축적을 표현하는 연산자이며, $\tilde{\mathbf{r}}_{i+1}$은 i, i+1 평면 상이 중간 평면의 좌표이다.

$\mathcal{T}[z_i, z_{i+1}]=\exp[-i\phi(\mathbf{r}_{i+1})]$
축적된 phase는 $\phi(\mathbf{r}_i)=k\int^{z_{i+1}}_{z_i}\delta n(\mathbf{r}_i) dz$, 즉 굴절률 차이만큼의 누적의 파수 k배이다.

매질을 통해 통한 전파를 2가지 효과로 분리할 수 있다.

1. 회절(diffraction)
   자유공간의 회절은 연산자 $\mathcal{R}$
2. 굴절(refraction)
   굴절은 연산자 $\mathcal{T}$

이 방법은 atmospheric turbulence를 통한 전파를 시뮬레이션하는데 자주 사용된다. 실제로, 광학 실험실에서 turbulence를 통한 propgation을 emulate(모사)할 때 사용하기도 한다.
빛과 매질 사이 상호작용이 있는 것에 대해 partial vaccum propagation을 대체 방법이다.

vaccum propgation 알고리즘을 살짝 고쳐 아래처럼 표현한다. interaction이 있는 평면 n개와 n-1개의 전파가 있을 때, vaccum 상황은 각 step에서 $\mathcal{T}=1$일 때고, turbulence가 있는 경우 어떻게 $\mathcal{T}$ 구현을 생성할지가 중요하다.

$U(\mathbf{r}_n)=\mathcal{Q}[\cfrac{m_{n-1}-1}{m_{n-1}\ \triangle z_{n-1}}, \mathbf{r}_n]$
$\times\prod\limits^{n-1}_{i=1}\{\mathcal{T}[z_i,z_{i+1}] \mathcal{F}^{-1}[\mathbf{f}_i,\cfrac{\mathbf{r}_{i+1}}{m_i}]\ \mathcal{Q}_2[-\cfrac{\triangle z_i}{m_i}, \mathbf{f}_i]\ \mathcal{F}[\mathbf{r}_i,\mathbf{f}_i]\cfrac{1}{m_i}\}$
$\times\{\mathcal{Q}[\cfrac{1-m_1}{\triangle z_1}, \mathbf{r}_1]\ \mathcal{T}[z_1,z_2] \ U(\mathbf{r}_1)\}$

• $\mathcal{Q}[c,\mathbf{r}]\{U(\mathbf{r})\}\equiv e^{i\frac{k}{2}c|\mathbf{r}|^2} U(\mathbf{r})$
• $\mathcal{Q}_2[d,\mathbf{r}] \{U(\mathbf{r})\} \equiv e^{i\pi^2\frac{2d}{k}|\mathbf{r}|^2}U(\mathbf{r})$
• $\mathcal{F}[\mathbf{r},\mathbf{f}]\{U(\mathbf{r})\}\equiv \int\limits^{\infty}_{-\infty}U(\mathbf{r})e^{-i2\pi f\cdot\mathbf{r}}d\mathbf{r}$

def ang_spec_multi_prop(Uin, wavelength, delta1, deltan, z, t):
    N = Uin.shape[0]

    n_1d = np.arange(-N/2,N/2)
    nx, ny = np.meshgrid(n_1d, n_1d)
    k = 2*np.pi/wavelength

    # super-Gaussian absorbing bounday
    nsq = nx**2+ny**2
    w = 0.47*N
    sg = np.exp(-nsq**8/w**16)
    #################################
    
    z = np.insert(z, 0, 0)    
    n = len(z) # 11

   
    # propgation distacne
    Delta_z = z[1:n]-z[0:n-1]

    # grid spacings
    alpha = z / z[n-1]
    
    delta = (1-alpha) * delta1 + alpha * deltan
    
    m = delta[1:n] / delta[0:n-1]
    print(m)
    x1 = nx * delta[0]
    y1 = ny * delta[0]
    rlsq = x1**2+y1**2

1. $Q_1=e^{i\frac{k}{2}\frac{1-m_1}{\triangle z_1}*(x_1^2+y_1^2)}*U_{in}(x_1,y_1)$

2. $\mathcal{T}(x)=e^{-(x/\sigma)^{2n}}$
   - $n=8$
   - $\sigma = 0.47*N*d1$
   sg = np.exp(-(x1/(0.47*N*d1))**16) * np.exp(-(y1/(0.47*N*d1))**16)
   

	Q1 = np.exp(1j*k/2*(1-m[0])/Delta_z[0]*rlsq)    
    Uin = Uin * Q1 * t[0,:,:]

2. $Q_2=e^{i\pi^2\frac{-2\triangle z_i/m_i}{k}|f_x^2+f_y^2|}$
3. $1/m_i → \mathcal{F} → \mathcal{Q}_2 → \mathcal{F}^{-1} → \mathcal{T}$

	for idx in range(0, n-1):
        # spatial frequencies (of ith plane)
        delta_f = 1 / (N*delta[idx])
        fX = nx * delta_f
        fY = ny * delta_f
        fsq = fX**2+fY**2
        dZ = Delta_z[idx] # propagation distance
        # quadratic pphase factor
        Q2 = np.exp(-1j*np.pi**2 *2*dZ/(m[idx]*k)*fsq)

        # compute th e propagated field
        Uin = sg * t[idx+1,:,:] * np.fft.ifft2(np.fft.fftshift(Q2)*np.fft.fft2(Uin/m[idx],  norm='ortho'), norm='ortho')

4. $U_{out}=e^{i\frac{k}{2}\frac{m_{n-1}-1}{m_{n-1}\triangle z_{n-1}}*(x_n^2+y_n^2)}*U_{in}(x_n,y_n)$
   $\mathcal{Q}[\cfrac{m_{n-1}-1}{m_{n-1}\ \triangle z_{n-1}}, \mathbf{r}_n]$

    xn = nx*delta[n-1]
    yn = ny*delta[n-1]
    rnsq = xn**2+yn**2
    Q3 = np.exp(1j*k/2*(m[n-2]-1)/(m[n-2]*dZ) *rnsq) # 0~9
    Uout = Q3 * Uin
    return xn, yn, Uout
    

---

2. Theory

이론적으로 대기 phase screen을 나타내기 위해, 유효 구조 상수 $C_{n_i}^2$의 항으로 turbulence profile을 나타낸다.

• $z_i:$전파 거리
• $\triangle z_i:$ i번째 phase screen에 표현되는 확장된 turbulence의 판 두께

$C_{n_i}^2$ 값은 연속 모델의 모델의 중요한 통계적 특성(평균, 분산, 갯수, low order moments)를 layer model에 맞추기 위해 선택된다.

$\int\limits_0^{\triangle z_i}C_n^2(z') (z')^m dz' = \sum\limits_{i=1}^nC_{n_i}^2z_i^m\triangle z_i$

• $n:$ 사용되는 phase screen 갯수 $0\le m\le7$

이러한 방법으로 layered 모델의 $r_0, \theta_0, \sigma_\chi^2$이 원래의 turbulence의 매개변수와 맞춘다.

• layered turbulence model을 위한 대기 parameter들은 discrete-sum 버전을 사용해 계산된다.

1. $r_{0,pw}=[0.423k^2\int\limits^{\triangle z}_0 C_n^2(z) dz]^{-3/5}$
   $r_{0,pw} = (0.423k^2\sum\limits_i C^2_{n_i}\triangle z_i)^{-3/5}$
2. $r_{0,sw}=[0.423k^2\int\limits^{\triangle z}_0 C_n^2(z)(\cfrac{z}{\triangle z})^{5/3} dz]^{-3/5}$
   $r_{0,sw} = [0.423k^2\sum\limits_{i=1}^n C^2_{n_i}(\cfrac{z_i}{\triangle z})^{5/3}\triangle z_i]^{-3/5}$
3. $\sigma^2_{\chi, pw}=0.563k^{7/6}\triangle z^{5/6}\int\limits^{\triangle z}_0C_n^2(z)(1-\cfrac{z}{\triangle z})^{5/6}dz$
   $\sigma^2_{\chi, pw}=0.563k^{7/6}\triangle z^{5/6}\sum\limits_{i=1}^nC_{n_i}^2(1-\cfrac{z_i}{\triangle z})^{5/6}\triangle z_i$
4. $\sigma^2_{\chi, sw}=0.563k^{7/6}\int\limits^{\triangle z}_0C_n^2(z)z^{5/6}(1-\cfrac{z}{\triangle z})^{5/6}dz$
   $\sigma^2_{\chi,sw}=0.563k^{7/6}\triangle z^{5/6} \sum\limits_{i=1}^nC_{n_i}^2 (\cfrac{z_i}{\triangle z})^{5/6}(1-\cfrac{z_i}{\triangle z})^{5/6}\triangle z_i$

$r_{0,pw} = (0.423k^2\sum\limits_i C^2_{n_i}\triangle z_i)^{-3/5}$을 이용하여, i번째 layer는 유효 coherence diameter $r_{0_i}$를 다음과 같이 표현하여 더 간단히 표현하면

$r_{0_i}=[0.423k^2C_{n_i}^2\triangle z_i]^{-3/5}$

1. $r_{0,pw}=(\sum\limits_{i=1}^nr_{0_i}^{-5/3})^{-3/5}$
2. $r_{0,sw}=(\sum\limits_{i=1}^nr_{0_i}^{-5/3}(\cfrac{z_i}{\triangle z})^{5/3})^{-3/5}$
3. $\sigma^2_{\chi, pw}=1.33k^{-5/6}\triangle z^{5/6}\sum\limits_{i=1}^nr_{0_i}^{-5/3}(1-\cfrac{z_i}{\triangle z})^{5/6}$
4. $\sigma^2_{\chi, sw}=1.33k^{-5/6}\triangle z^{5/6}\sum\limits_{i=1}^nr_{0_i}^{-5/3}(\cfrac{z_i}{\triangle z})^{5/6}(1-\cfrac{z_i}{\triangle z})^{5/6}$

요구된 atmospheric 조건 set($r_{0,sw}, \sigma_{\chi, sw}^2$)이 주어지면, 예로 이 방정식들은 경로를 따른 위치와 phase screen 특성을 결정하는데 사용할 수 있다.
5~10 개의 보통 phase screen 수를 사용하면 10~20개의 모르는(각 screen 마다 $r_0,z_i$) 변수가 있다. 그래서 2개의 방정식으로 시스템을 결정할 수 없다.
partial-propagation planes의 균일한 spacing은 일관성을 유지할 수 있어 이러한 방법으로 phase screen 위치를 고정하여 향상시킬 수 해결할 수 있다. 각 partial-rpopagation plane에서 phase screen을 놓는 것을 선택한 다음 방정식을 단순화하는 $\alpha_i=\cfrac{z_i}{\triangle z}$을 이용.
예로 5개의 screen이면,
$\left[ \begin{matrix} r_{0,sw}^{-5/3} \\ \frac{\sigma_{\chi,sw}^2}{1.33}(\frac{k}{\triangle z})^{5/6} \\ \end{matrix} \right] = \left[ \begin{matrix} 0 & \alpha_2^{5/3} & \alpha_3^{5/3} & \alpha_4^{5/3} & \alpha_5^{5/3} \\ 0 & \alpha_2^{5/6}(1-\alpha_2)^{5/6} & \alpha_3^{5/6}(1-\alpha_3)^{5/6} & \alpha_4^{5/6}(1-\alpha_4)^{5/6} & \alpha_5^{5/6}(1-\alpha_5)^{5/6} \\ \end{matrix} \right] \left[ \begin{matrix} r_{0_1}^{-5/3} \\ r_{0_2}^{-5/3} \\ r_{0_3}^{-5/3} \\ r_{0_4}^{-5/3} \\ r_{0_5}^{-5/3} \\ \end{matrix} \right]$

위 방법으로 왼쪽 항이 $r_{0,sw}, \sigma_{\chi, sw}^2$가 결정된다.

1. $\lambda, \triangle z, C_n^2(z)$ 모델이 주어지면, 우리는 simulation을 위해 요구되는 atmospheric 매개변수를 계산한다.
2. 위 행렬식과 같이, 적절한 시스템의 방정식을 풀어서, phase screen의 $r_0$ 값을 구한다.
3. 이 접근 방식의 어려움은 $r_0$ 벡터에서 -5/3 승수가 있다는 점이다. $r_0$ 벡터에서 음수 항목은 물리적으로 불가능하므로, 해결된 $r_0$ 벡터의 값은 양수로 제한되어야 한다.

몇 개의 phase screen을 시뮬레이션하기 위한 $r_0$ 계산하기 위한 constrained optmization의 사용을 볼 수 있다.