%2015/06/29
%2016/10/19

\chapter{放射伝達計算}

%DCRTM では, DCPAM に実装されている放射サブルーチンを呼び出して放射伝達を計算している. 

%\begin{verbatim}
%https://www.gfd-dennou.org/library/dcpam/dcpam5/dcpam5_latest/src/
%radiation/rad_rte_two_stream_app.f90
%\end{verbatim}

%本コードは, Toon et al., 1989 による放射コードである. 太陽放射の散乱光のフラックス計算には $\delta$-Eddington 近似による二方向近似を行う. 赤外放射に対しては, 二方向近似($\delta$ 半球平均近似) によって多重散乱に寄与する放射フラックスを計算し, それを用いて天頂角余弦で離散化した放射強度を求める. 

\section{はじめに}

ここでは, 放射伝達計算について述べる. DCRTM の放射伝達計算は, Toon et al. (1989) を基本としている. 中心星放射が惑星大気により散乱されるされることによる放射は, Toon et al., (1989) の Two-Steam 法を用いて計算する. 散乱は $\delta$-Eddington 近似を用いている. 惑星大気が射出することによる放射伝達は, Toon et al., (1989) の Two-Stream Source Function 法を用いて計算する. この計算には惑星放射の散乱も考慮されている. 惑星放射の散乱は, $\delta$-半球平均近似による. この Two-Stream Source Function 法は, Two-Stream 法により計算された放射場を放射源関数として用いて放射伝達を解く方法であり, Toon et al., (1989) によれば, 惑星放射の放射伝達計算については, Two-Stream 法よりも一般に精度がよいとされている. しかし, ある条件のもとでは, Two-Stream Source Function 法による放射場で計算した加熱率が物理的に不自然な値となることがある. これは, Two-Stream 法による放射場と, それを放射源とした Two-Stream Source Function 法による放射場が整合的でない場合に起こる. 本来2つの放射場は一致すべきものであり, DCRTM では, Two-Stream Source Function 法による放射場を再度放射源として放射伝達式を解くことを繰り返し, self consistent な解を計算している. 

\section{Two-Stream 法}

DCRTM の放射伝達計算モジュールは, DCPAM の放射伝達モジュールをもとにしている. Two-Stream 法の方程式については, DCPAM のドキュメント (DCPAM5支配方程式とその離散化) を参照のこと. 

\section{Two-Stream Source Function 法 (放射源関数法)}

放射伝達方程式 (波数の添え字は省略, 方位角方向には積分済)
\begin{eqnarray} 
\mu \frac{\partial I (\tau, \mu)}{\partial \tau} = I(\tau, \mu) - S(\tau, \mu)
\end{eqnarray} 
より, 第 n 層の上端からの上向き放射 $I_n^+(0,\mu)$ は, 
\begin{eqnarray} 
I_n^+(0,\mu) = I_n^+(\tau_n,\mu) \exp [-\tau_n/\mu]+\int_0^{\tau_n} S_n^+(t,\mu) \exp[-t/\mu] \frac{dt}{\mu}. 
\end{eqnarray} 
同様に, 第 n 層の下端からの下向き放射 $I_n^-(\tau_n,-\mu)$ は, 
\begin{eqnarray} 
I_n^-(\tau_n,-\mu) = I_n^-(0,-\mu) \exp [-\tau_n/\mu]+\int_0^{\tau_n} S_n^-(t,-\mu) \exp[-(\tau_n-t)/\mu] \frac{dt}{\mu}. 
\end{eqnarray} 
$S(\tau, \mu) $ は放射源を表す関数であり, 
\begin{eqnarray} 
S(\tau, \mu) = S_e(\tau, \mu) + \frac{\varpi}{2} \int^1_{-1} P(\mu, \mu^{\prime}) I (\tau, \mu^{\prime}) d\mu^{\prime}. 
\end{eqnarray} 
右辺第1項は, 大気の射出項で, 
\begin{eqnarray} 
S_e(\tau, \mu) = 2 \pi (1-\varpi) B (T(\tau)). 
\end{eqnarray} 
右辺第2項は, 散乱項である. 

\subsection{Toon et al. (1989) による Two-Stream Source Function 法}

Toon et al. (1989) の Two-Stream Source Function 法では, 放射源関数 $S(\tau,\mu)$ 内の $I(\tau,\mu)$ に, Two-Stream 法により計算された放射場を与える. Two-Stream 法によるフラックスを $F_{TS}(\tau)$ とすると, 
\begin{eqnarray} 
I^{\pm}(\tau,\mu) = F_{TS}^{\pm}(\tau)/\mu_1,
\end{eqnarray} 
\begin{eqnarray} 
\mu_1 = \frac{1}{2}.
\end{eqnarray} 
第n層のフラックスは, 
\begin{eqnarray} 
 F_{TS, n}^{+} (\tau)&=& Y_{1,n}\left(\exp[-\lambda_n(\tau_n -\tau)] + \Gamma_n \exp[-\lambda_n \tau]\right) \\\nonumber
&&+ Y_{2,n}\left(\exp[-\lambda_n(\tau_n -\tau)] - \Gamma_n \exp[-\lambda_n \tau]\right) + C_n^+(\tau).
\end{eqnarray} 
\begin{eqnarray} 
 F_{TS, n}^{-} (\tau)&=& Y_{1,n}\left(\Gamma_n \exp[-\lambda_n(\tau_n -\tau)] + \exp[-\lambda_n \tau]\right) \\\nonumber
&&+ Y_{2,n}\left(\Gamma_n \exp[-\lambda_n(\tau_n -\tau)] -  \exp[-\lambda_n \tau]\right) + C_n^-(\tau).
\end{eqnarray} 
位相関数 $P(\mu, \mu^{\prime})$は, 
\begin{eqnarray}
P(\mu, \mu^{\prime})=\left\{ \begin{array}{ll}
 1+ g & (前方散乱) \\
 1-g & (後方散乱) 
 \end{array} \right.
\end{eqnarray}  
$\mu > 0 $ の時, 
\begin{eqnarray} 
\frac{\varpi}{2} \int^1_{-1} P(\mu, \mu^{\prime}) I (\tau, \mu^{\prime}) d\mu^{\prime} &=&
\frac{\varpi}{2} \int^0_{-1} P(\mu, \mu^{\prime}) \frac{F_{TS}^- (\tau)}{\mu_1} d\mu^{\prime} +
\frac{\varpi}{2} \int^1_0 P(\mu, \mu^{\prime}) \frac{F_{TS}^+ (\tau)}{\mu_1} d\mu^{\prime}\\
&=&
\frac{\varpi}{2} (1-g) \frac{F_{TS}^- (\tau)}{\mu_1} +
\frac{\varpi}{2} (1+g) \frac{F_{TS}^+ (\tau)}{\mu_1}. 
\end{eqnarray} 
$\mu < 0 $ の時, 
\begin{eqnarray} 
\frac{\varpi}{2} \int^1_{-1} P(\mu, \mu^{\prime}) I (\tau, \mu^{\prime}) d\mu^{\prime} &=&
\frac{\varpi}{2} \int^0_{-1} P(\mu, \mu^{\prime}) \frac{F_{TS}^- (\tau)}{\mu_1} d\mu^{\prime} +
\frac{\varpi}{2} \int^1_0 P(\mu, \mu^{\prime}) \frac{F_{TS}^+ (\tau)}{\mu_1} d\mu^{\prime}\\
&=&
\frac{\varpi}{2} (1+g) \frac{F_{TS}^- (\tau)}{\mu_1} +
\frac{\varpi}{2} (1-g) \frac{F_{TS}^+ (\tau)}{\mu_1}. 
\end{eqnarray} 
%ここで, 計算を簡単にするために $g = 0$ の条件の場合を考える. この時, 第n層の散乱項は, 
よって, 
\begin{eqnarray} 
\frac{\varpi_n}{2} \int^1_{-1} P_n(\mu, \mu^{\prime}) I_n (\tau, \mu^{\prime}) d\mu^{\prime} &=&
\frac{\varpi_n}{2\mu_1} (F_{TS,n}^+ (\tau) + F_{TS,n}^- (\tau))\\\nonumber
&&+\frac{\varpi_n g}{2\mu_1} (F_{TS,n}^+ (\tau) - F_{TS,n}^- (\tau)) \\
&=&
\varpi_n (F_{TS,n}^- (\tau) + F_{TS,n}^+ (\tau))\\\nonumber
&&+\varpi_n g(F_{TS,n}^+ (\tau) - F_{TS,n}^- (\tau))\\
&=&
\varpi_n(Y_{1,n}+Y_{2,n})(1+\Gamma_n)\exp[-\lambda_n(\tau_n-\tau)]\\\nonumber
&&+\varpi_n(Y_{1,n}-Y_{2,n})(1+\Gamma_n)\exp[-\lambda_n\tau]\\\nonumber
&&+\varpi_n 2\pi  (B_{0,n}+B_{1,n}\tau)\\\nonumber
&&\pm\varpi_n g(Y_{1,n}+Y_{2,n})(1-\Gamma_n)\exp[-\lambda_n(\tau_n-\tau)]\\\nonumber
&&\mp\varpi_n g(Y_{1,n}-Y_{2,n})(1-\Gamma_n)\exp[-\lambda_n\tau]\\\nonumber
&&\pm\varpi_n g\frac{2\pi B_{1,n}}{\gamma_{1,n}+\gamma_{2,n}}. 
\end{eqnarray} 
(ここで, $\pm, \mp$ の符号は, 上が $\mu > 0 下が, \mu < 0$.以下同様.)
さらに, 
\begin{eqnarray} 
S_{e,n}(\tau, \mu) = 2\pi(1-\varpi_n) B (T(\tau)) = 2\pi(1-\varpi_n) (B_{0,n}+B_{1,n} \tau)
\end{eqnarray} 
であるから, 
\begin{eqnarray} 
S_n^{\pm}(\tau, \pm\mu) &=& S_{e,n}(\tau, \mu) + \frac{\varpi_n}{2} \int^1_{-1} P_n(\mu, \mu^{\prime}) I_n (\tau,\mu^{\prime})d\mu^{\prime}\\ 
&=&\varpi_n(Y_{1,n}+Y_{2,n})(1+\Gamma_n)\exp[-\lambda_n(\tau_n-\tau)]\\\nonumber
&&+\varpi_n(Y_{1,n}-Y_{2,n})(1+\Gamma_n)\exp[-\lambda_n\tau]\\\nonumber
&&+2\pi (B_{0,n}+B_{1,n}\tau)\\\nonumber
&&\pm\varpi_n g(Y_{1,n}+Y_{2,n})(1-\Gamma_n)\exp[-\lambda_n(\tau_n-\tau)]\\\nonumber
&&\mp\varpi_n g(Y_{1,n}-Y_{2,n})(1-\Gamma_n)\exp[-\lambda_n\tau]\\\nonumber
&&\pm\varpi_n g\frac{2\pi B_{1,n}}{\gamma_{1,n}+\gamma_{2,n}}. 
\end{eqnarray} 
よって, Two-Stream Source Function 法による第 n 層の上端からの上向き放射 $I^+_n(0, \mu)$, 下端からの下向き放射 $I^-_n(\tau_n, -\mu)$は,  それぞれ以下のようになる. 
\begin{eqnarray} 
I^+_n(0, \mu) &=&I^+_n(\tau_n, \mu) \exp[-\tau_n/\mu] + \int^{\tau_n}_0 S_n^+(t,\mu) \exp[-t/\mu] \frac{dt}{\mu}\\
&=&I^+_n(\tau_n,\mu)\exp[-\tau_n/\mu]\\\nonumber
&&+\frac{G_n}{(\lambda_n \mu - 1)} \left(\exp[-\tau_n/\mu]-\exp[-\tau_n\lambda_n]\right)\\\nonumber
&&+\frac{H_n}{(\lambda_n \mu + 1)} \left(1-\exp[-\tau_n(\lambda_n +1/\mu)]\right)\\\nonumber
&&+\alpha_1 \left(1-\exp[-\tau_n/\mu] \right)+\alpha_2 \left(\mu - (\tau_n+\mu)\exp[-\tau_n/\mu] \right). 
\end{eqnarray} 
\begin{eqnarray} 
I^-_n(\tau_n, -\mu) &=&I^-_n(0, -\mu) \exp[-\tau_n/\mu] + \int^{\tau_n}_0 S_n^-(t,-\mu) \exp[-(\tau_n-t)/\mu] \frac{dt}{\mu}\\
&=&I^-_n(0,-\mu)\exp[-\tau_n/\mu]\\\nonumber
&&+\frac{J_n}{(\lambda_n \mu + 1)} \left(1-\exp[-\tau_n(\lambda_n +1/\mu)]\right)\\\nonumber
&&+\frac{K_n}{(\lambda_n \mu - 1)} \left(\exp[-\tau_n/\mu]-\exp[-\tau_n\lambda_n]\right)\\\nonumber
&&+\sigma_1 \left(1-\exp[-\tau_n/\mu] \right)+\sigma_2 \left(\mu \exp[-\tau_n/\mu]+\tau_n - \mu) \right). 
\end{eqnarray} 
ここで, 
 \begin{eqnarray}
G_n &=& (Y_{1,n}+Y_{2,n}) \varpi_n (1+g + \Gamma_n (1-g))\\
H_n &=& (Y_{1,n}-Y_{2,n}) \varpi_n (1-g + \Gamma_n (1+g))\\
J_n &=& (Y_{1,n}+Y_{2,n}) \varpi_n (1-g + \Gamma_n (1+g))\\
K_n &=& (Y_{1,n}-Y_{2,n}) \varpi_n (1+g + \Gamma_n (1-g))\\
\alpha_1&=&2\pi \left(B_{0,n} + \frac{\varpi_n g B_{1,n}}{\gamma_{1,n}+\gamma_{2,n}} \right)\\
\alpha_2&=&2\pi B_{1,n}\\
\sigma_1&=&2\pi \left(B_{0,n} - \frac{\varpi_n g B_{1,n}}{\gamma_{1,n}+\gamma_{2,n}} \right)\\
\sigma_2&=&2\pi B_{1,n}
 \end{eqnarray}  
さらに, 
 \begin{eqnarray}
\varpi_n (1+g) &=& 2 - \gamma_{1,n}\\
\varpi_n (1-g) &=& \gamma_{2,n}\\
\gamma_{2,n}\Gamma_n &=& \gamma_{1,n} -\lambda_n\\
\gamma_{2,n}/\Gamma_n &=& \gamma_{1,n} +\lambda_n\\
\mu_1 &=& 1/2
 \end{eqnarray}  
であることに注意すると, 
 \begin{eqnarray}
G_n &=& (Y_{1,n}+Y_{2,n}) \varpi_n (1+g + \Gamma_n (1-g))\\
&=&(Y_{1,n}+Y_{2,n}) \left( 2-\gamma_{1,n} + \gamma_{2,n}\Gamma_n \right)\\
&=&(Y_{1,n}+Y_{2,n}) \left( 2-\gamma_{1,n} + \gamma_{1,n} -\lambda_n \right)\\
&=&(Y_{1,n}+Y_{2,n}) \left( 1/\mu_1 -\lambda_n \right)
\end{eqnarray}
 \begin{eqnarray}
H_n &=& (Y_{1,n}-Y_{2,n}) \varpi_n (1-g + \Gamma_n (1+g))\\
&=&(Y_{1,n}-Y_{2,n}) \left(\gamma_{2,n} + (2-\gamma_{1,n})\Gamma_n \right)\\
&=&(Y_{1,n}-Y_{2,n}) \Gamma_n \left( \gamma_{2,n}/\Gamma_n +2- \gamma_{1,n} \right)\\
&=&(Y_{1,n}-Y_{2,n}) \Gamma_n \left( \gamma_{1,n} +\lambda_n +2- \gamma_{1,n} \right)\\
&=&(Y_{1,n}-Y_{2,n}) \left( \lambda_n -1/\mu_1\right)
\end{eqnarray}
同様に, 
 \begin{eqnarray}
J_n &=& (Y_{1,n}+Y_{2,n})\Gamma_n \left(\lambda_n+1/\mu_1\right)\\
K_n &=& (Y_{1,n}-Y_{2,n}) \left(1/\mu_1 - \lambda_n \right)
 \end{eqnarray}  
また, 
 \begin{eqnarray}
\alpha_1&=&2\pi \left(B_{0,n} + \frac{\varpi_n g B_{1,n}}{\gamma_{1,n}+\gamma_{2,n}} \right)\\
&=&2\pi \left(B_{0,n} + \frac{(1-(\gamma_{1,n}+\gamma_{2,n})/2) B_{1,n}}{\gamma_{1,n}+\gamma_{2,n}} \right)\\
&=&2\pi \left(B_{0,n} + B_{1,n}\left(\frac{1}{\gamma_{1,n}+\gamma_{2,n}}-\frac{1}{2} \right) \right)\\
&=&2\pi \left(B_{0,n} + B_{1,n}\left(\frac{1}{\gamma_{1,n}+\gamma_{2,n}}-\mu_1 \right) \right)\\
 \end{eqnarray} 
同様に, 
\begin{eqnarray}
\sigma_1= 2\pi \left(B_{0,n} - B_{1,n}\left(\frac{1}{\gamma_{1,n}+\gamma_{2,n}}-\mu_1 \right) \right)
\end{eqnarray} 

\subsection{放射源の放射場を self consistent に求める}

放射源関数内の $I(\tau. \mu)$ と放射伝達方程式を解いて得られる放射場は整合的でなければならない. ここでは, self consistent な放射場を得るために, iterative に放射場を求めることを想定して放射伝達方程式を書き下す. 

\subsubsection {上向き放射 $I^+_n(0,\mu)$}

第 n 層上端から $\mu$ の角度に射出する放射 $I^+_n(0,\mu)$ は, 以下のように書ける. 
\begin{eqnarray} 
I_n^+(0,\mu) = I_n^+(\tau_n,\mu) \exp [-\tau_n/\mu]+\int_0^{\tau_n} S_n^+(t,\mu) \exp[-t/\mu] \frac{dt}{\mu}. 
\end{eqnarray} 
\begin{eqnarray} 
S_n(\tau, \mu) = 2 \pi (1-\varpi_n) B (T(\tau_n))+ \frac{\varpi_n}{2} \int^1_{-1} P_n(\mu, \mu^{\prime}) I_n (\tau, \mu^{\prime}) d\mu^{\prime}. 
\end{eqnarray} 

今, 放射源関数内の放射場を $I_{s,n} (\tau,\mu)$ で与える. 散乱には半球平均近似を用いると, $\mu > 0 の条件で, $
\begin{eqnarray} 
\frac{\varpi_n}{2} \int^1_{-1} P_n(\mu, \mu^{\prime}) I (\tau, \mu^{\prime}) d\mu^{\prime} 
&=&
\frac{\varpi_n}{2} \int^0_{-1} P_n(\mu, \mu^{\prime}) I^-_{s,n}(\tau,-\mu^{\prime}) d\mu^{\prime} \\\nonumber
&&+
\frac{\varpi_n}{2} \int^1_0 P_n(\mu, \mu^{\prime}) I^+_{s,n} (\tau,\mu^{\prime})d\mu^{\prime}\\
&=&
\frac{\varpi_n}{2} (1-g) \int^0_{-1}I^-_{s,n}(\tau,-\mu^{\prime})d\mu^{\prime}\\\nonumber
&&+
\frac{\varpi_n}{2} (1+g)  \int^1_0I^+_{s,n}(\tau,\mu^{\prime})d\mu^{\prime}. 
\end{eqnarray} 
ここで, 天頂角積分の評価にはガウス求積法を用いる. $\omega_i$ をガウスの重み関数として, 
\begin{eqnarray} 
\frac{\varpi_n}{2} \int^1_{-1} P_n(\mu, \mu^{\prime}) I (\tau, \mu^{\prime}) d\mu^{\prime}
&\sim&
\frac{\varpi_n}{2} (1-g) \sum_i \omega_i I^-_{s,n}(\tau,-\mu^{\prime}_i)\\\nonumber
&&+\frac{\varpi_n}{2} (1+g)  \sum_i \omega_i I^+_{s,n}(\tau,\mu^{\prime}_i)
\end{eqnarray} 

さてここで, 第 n 層内の放射強度を, 層の上端と下端の放射強度の光学的厚さに対する線型な関数で近似すると, 
\begin{eqnarray} 
I^{\pm}_{s,n}(\tau, \pm\mu) &=& I^{\pm}_{0,n}(\pm\mu) + I^{\pm}_{1,n}(\pm\mu)\tau
\end{eqnarray} 
ここで, 
\begin{eqnarray} 
I^{\pm}_{0,n}(\pm\mu) &=& I^{\pm}_{s,n}(0,\pm\mu)\\
I^{\pm}_{1,n}(\pm\mu) &=& \left(I^{\pm}_{s,n}(\tau_n,\pm\mu)-I^{\pm}_{s,n}(0,\pm\mu)\right)/\tau_n. 
\end{eqnarray} 

以上から, $S^+_n(\tau, \mu) は, \mu に依存しない関数となることが分かる. 従って S^+_n(\tau, \mu) = S^+_n(\tau) として, $
\begin{eqnarray} 
S^+_n(\tau) &\sim& 2\pi(1-\varpi_n)(B_{0,n}+B_{1,n}\tau)\\\nonumber
&&+\frac{\varpi_n}{2} (1-g) \sum_i \omega_i \left(I^-_{0,n}(-\mu^{\prime}_i) + I^-_{1,n}(-\mu^{\prime}_i)\tau\right)\\\nonumber
&&+\frac{\varpi_n}{2} (1+g)  \sum_i \omega_i \left(I^+_{0,n}(\mu^{\prime}_i) + I^+_{1,n}(\mu^{\prime}_i)\tau\right)\\
&=&2\pi(1-\varpi_n)B_{0,n}+2\pi(1-\varpi_n)B_{1,n}\tau\\\nonumber
&&+\frac{\varpi_n}{2}\sum_i\left((1-g)\omega_i I^-_{0,n}(-\mu^{\prime}_i)+(1+g)\omega_i I^+_{0,n}(\mu^{\prime}_i)\right)\\\nonumber
&&+\frac{\varpi_n}{2}\sum_i\left((1-g)\omega_i I^-_{1,n}(-\mu^{\prime}_i)+(1+g)\omega_i I^+_{1,n}(\mu^{\prime}_i)\right)\tau\\
&\equiv&U_n+V_n\tau+W^+_n+X^+_n\tau
\end{eqnarray} 
ここで, 
\begin{eqnarray} 
U_n&=&2\pi(1-\varpi_n)B_{0,n}\\
V_n&=&2\pi(1-\varpi_n)B_{1,n}\\
W^{\pm}_n&=&\frac{\varpi_n}{2}\sum_i\left\{(1\pm g)\omega_i I^+_{0,n}(\mu^{\prime}_i\}+(1\mp g)\omega_i I^-_{0,n}(-\mu^{\prime}_i)\right\}\\
&=&\frac{\varpi_n}{2}\left\{(1\pm g)L^+_n(0)+(1 \mp g)L^-_n(0)\right\}\\
X^{\pm}_n&=&\frac{\varpi_n}{2}\sum_i\left\{(1\pm g)\omega_i I^+_{1,n}(\mu^{\prime}_i)+(1\mp g)\omega_i I^-_{1,n}(-\mu^{\prime}_i)\right\}\\
&=&\frac{\varpi_n}{2}\left\{(1\pm g)(L^+_n(\tau_n)-L^+_n(0))+(1\mp g)(L^-_n(\tau_n)-L^-_n(0))\right\}/\tau_n\\
L^{\pm}_n(\tau)&=&\sum_i\omega_i I^{\pm}_{s,n}(\tau,\pm\mu_i)
\end{eqnarray} 
以上より, 
\begin{eqnarray} 
I_n^+(0,\mu) &=& I_n^+(\tau_n,\mu) \exp [-\tau_n/\mu]+\int_0^{\tau_n} S_n^+(t,\mu) \exp[-t/\mu] \frac{dt}{\mu}\\
&\sim&I_n^+(\tau_n,\mu) \exp [-\tau_n/\mu]+(U_n+W^+_n)(1-\exp[-\tau/\mu])\\\nonumber
&&+(V_n+X^+_n)\left\{\mu-(\mu+\tau_n)\exp[-\tau_n/\mu]\right\}. 
\end{eqnarray} 

\subsubsection {下向き放射 $I^-_n(\tau_n,-\mu)$}

上向き放射と同様の議論により, 下向き放射 $I^-_n(\tau_n,-\mu)$ の式を導く. 

第 n 層下端から $\mu$ の角度に射出する放射 $I^-_n(\tau_n,-\mu)$ は, 以下のように書ける. 
\begin{eqnarray} 
I_n^-(\tau_n,-\mu) = I_n^-(0,-\mu) \exp [-\tau_n/\mu]+\int_0^{\tau_n} S_n^-(t,-\mu) \exp[-(\tau_n-t)/\mu] \frac{dt}{\mu}. 
\end{eqnarray} 
\begin{eqnarray} 
S^-_n(\tau, -\mu) 
&=&
2 \pi (1-\varpi_n) B (T(\tau_n))+ \frac{\varpi_n}{2} \int^1_{-1} P_n(-\mu, \mu^{\prime}) I_n (\tau, \mu^{\prime}) d\mu^{\prime}\\
&=&
2\pi(1-\varpi_n)(B_{0,n}+B_{1,n}\tau)\\\nonumber
&&+\frac{\varpi_n}{2} \int^0_{-1} P_n(-\mu, \mu^{\prime}) I^-_{s,n}(\tau,-\mu^{\prime}) d\mu^{\prime} \\\nonumber
&&+\frac{\varpi_n}{2} \int^1_0 P_n(-\mu, \mu^{\prime}) I^+_{s,n} (\tau,\mu^{\prime})d\mu^{\prime}\\
&=&
2\pi(1-\varpi_n)(B_{0,n}+B_{1,n}\tau)\\\nonumber
&&+\frac{\varpi_n}{2} (1+g) \int^0_{-1}I^-_{s,n}(\tau,-\mu^{\prime})d\mu^{\prime}\\\nonumber
&&+\frac{\varpi_n}{2} (1-g)  \int^1_0I^+_{s,n}(\tau,\mu^{\prime})d\mu^{\prime}\\
&\sim&
2\pi(1-\varpi_n)(B_{0,n}+B_{1,n}\tau)\\\nonumber
&&+\frac{\varpi_n}{2} (1+g) \sum_i \omega_i I^-_{s,n}(\tau,-\mu^{\prime}_i)\\\nonumber
&&+\frac{\varpi_n}{2} (1-g)  \sum_i \omega_i I^+_{s,n}(\tau,\mu^{\prime}_i)\\
&=&
2\pi(1-\varpi_n)(B_{0,n}+B_{1,n}\tau)\\\nonumber
&&+\frac{\varpi_n}{2} (1+g) \sum_i \omega_i \left\{I^-_{0,n}(-\mu^{\prime}_i) + I^-_{1,n}(-\mu^{\prime}_i)\tau\right\}\\\nonumber
&&+\frac{\varpi_n}{2} (1-g)  \sum_i \omega_i \left\{I^+_{0,n}(\mu^{\prime}_i) + I^+_{1,n}(\mu^{\prime}_i)\tau\right\}\\
&=&
2\pi(1-\varpi_n)B_{0,n}+2\pi(1-\varpi_n)B_{1,n}\tau\\\nonumber
&&+\frac{\varpi_n}{2}\sum_i\left\{(1+g)\omega_i I^-_{0,n}(-\mu^{\prime}_i)+(1-g)\omega_i I^+_{0,n}(\mu^{\prime}_i)\right\}\\\nonumber
&&+\frac{\varpi_n}{2}\sum_i\left\{(1+g)\omega_i I^-_{1,n}(-\mu^{\prime}_i)+(1-g)\omega_i I^+_{1,n}(\mu^{\prime}_i)\right\}\tau\\
&=&
U_n+V_n\tau+W^-_n+X^-_n\tau
\end{eqnarray} 
以上より, 
\begin{eqnarray} 
I_n^-(\tau_n,-\mu) &=& I_n^-(0,-\mu) \exp [-\tau_n/\mu]+\int_0^{\tau_n} S_n^-(t,-\mu) \exp[-(\tau_n-t)/\mu] \frac{dt}{\mu}\\
&\sim&I_n^-(0,-\mu) \exp [-\tau_n/\mu]+(U_n+W^-_n)(1-\exp[-\tau/\mu])\\\nonumber
&&+(V_n+X^-_n)\left\{\tau_n-\mu(1-\exp[-\tau_n/\mu])\right\}. 
\end{eqnarray} 

\subsubsection {iterative に放射場を求める}

放射源関数法によって求められた第 n 層からの上向き, 下向き放射 $I^+_n(0,\mu), I^-_n(\tau_n,-\mu)$ は, 本来放射源関数 $S_n(\tau, \mu) $と整合的でなければならない. しかし, 放射源関数法で求められた放射場は必ずしもそうなっていない. そこで, $I^+_n(0,\mu), I^-_n(\tau_n,-\mu)$ を放射源関数とし, 再度放射伝達を解くことを繰り返すことで, 整合的な放射場を得ることを考える. 


今, 放射場 $I_{s, n, l}(\tau,\mu)$ をもとにした放射源関数 $S_{n,l}(\tau,\mu)$ を使って計算された放射場を, $I_{n,l}(\tau,\mu)$ とする. すなわち上向き放射を例に取れば, 
\begin{eqnarray} 
I_{n,l}^+(0,\mu) = I_{n,l}^+(\tau_n,\mu) \exp [-\tau_n/\mu]+\int_0^{\tau_n} S_{n,l}^+(t,\mu) \exp[-t/\mu] \frac{dt}{\mu}. 
\end{eqnarray} 
ここで, 
\begin{eqnarray} 
I_{s,n,l+1}^+(\tau,\mu) = I_{n,l}^+(\tau,\mu) 
\end{eqnarray} 
として, 再度放射伝達計算を行う. 
そして, 
\begin{eqnarray} 
I_{n,l+1}^+(\tau,\mu) = I_{n,l}^+(\tau,\mu) 
\end{eqnarray} 
となるまで計算を繰り返すことで, 整合的な放射場を得ることができる. 下向き放射についても同様である. ちなみに, 最初の放射源関数は, Two-Stream 法による放射場を与える. すなわち, 
\begin{eqnarray} 
I^{\pm}_{s,n,1}(\tau,\pm\mu) = F^{\pm}_{TS, n}(\tau,\pm\mu) /\mu_1.
\end{eqnarray} 
