輻射輸送を考慮した空の色の試算

空の色を計算する試みの続き。内容的には、ほぼ独立している。

前回の計算では、輻射の方向性など全く考慮してない雑な計算だった。しかし、レイリーの時代の計算は、あれより極端に高度ということもなかっただろうと思う。

 

太陽輻射には明らかに方向性があり、レイリー散乱にも異方性があるので、方向も込みで、真面目に散乱を考慮したら、どうなるかという計算。具体的には、輻射輸送方程式を解いて空の色を計算できるか検討していく。

輻射強度は、位置$${\mathbf{r}}$$、方向$${w}$$、周波数$${\nu}$$、時間$${t}$$に依存し、それを$${I(\mathbf{r} , w , \nu, t)}$$と書くことにする。

大気圏外では、方向wについて積分すると、黒体輻射スペクトルで近似できるようになる。

 

輻射輸送方程式は$${c}$$を光速として、

$${ \dfrac{1}{c} \dfrac{\partial I}{\partial t} + (w \cdot \nabla) I = A_{放射} - A_{吸収} I - A_{散乱} \displaystyle \int S(w,w') I(\mathbf{r} ,w , \nu,t) dw' + A_{散乱} \displaystyle \int S(w,w') I(\mathbf{r} ,w' , \nu,t) dw' }$$

と書ける。$${A_{放射},A_{吸収},A_{散乱}}$$は、放射、吸収、散乱に関わる係数。

放射は、位置、方向、周波数、時間全てに依存する可能性がある。最も単純には、位置の温度に応じた等方的な黒体輻射を仮定するものだろう。

吸収はいいだろう。水は光を吸収するし、空気も少しは可視光を吸収する。

散乱は、例えば、レイリー散乱などを想定している。$${S(w,w')}$$は散乱確率密度で、方向$${w}$$から$${w'}$$へ散乱される確率を表す。散乱が完全に等方的なら、単に$${1/4 \pi}$$に等しい。等方的でない場合は、微分散乱断面積で表示される量となる。レイリー散乱は、方向依存性を持つ。

$${A_{散乱}}$$は、位置、時間、周波数に依存しうる。例えば、レイリー散乱では、散乱断面積は、波長の4乗に逆比例するので、周波数の4乗に比例する。波長が凄く短い時には、レイリー散乱の前提条件(散乱体は波長に比べて十分に小さい)が妥当性を失う。

 

そのままでは、変数が多すぎて辛い。まずは時間微分項は落とすことにする。太陽の地球に対する相対位置は、時間によって変動しているが、多分、それが空の青さに影響してるわけではないだろう。

放射項と吸収項も無視しよう。放射は、大気中で強い可視光が放射されてたりはせんだろうという勘に基づく。吸収項は無視してはいけない気もするが、可視領域では、酸素や窒素など大気主成分の影響は比較的小さいように思う。

太陽光は鉛直真上から入射してきて、入射面は無限に広い平面とする。従って、水平方向の依存性はないとする。実際は、鉛直方向のスケールに比べて、水平方向のスケールが十分大きいかは明らかではないと思う。

朝焼けや夕焼けを解析したい場合は、この仮定は不適切だが、真上から入射してくる時でも、空は青い。

 

これで大分変数が減って、残った項は、

$${ w_{z} \dfrac{\partial I}{\partial z} = -A_{散乱}\displaystyle \int S(w,w') I(z,w) dw' + A_{散乱} \displaystyle \int S(w,w') I(z,w') dw'}$$

方向成分について、球面座標で

$${w = (w_1,w_2,w_3) = (\cos \psi \sin \theta , \sin \psi \sin \theta , \cos \theta)}$$

と変数を導入すると、水平方向の依存性は無視できるとしたので、$${I}$$は$${\theta}$$のみに依存する。

 

球面調和関数で軸対称性を仮定したので、$${\cos(n \theta)}$$で展開すればいいだろう。$${\theta}$$の周期関数と見て、フーリエ展開すると考えてもいい。

$${I(z,\theta,\nu) = \displaystyle \sum_{n=0}^{N} c_{n}(z,\nu) \cos(n \theta)}$$

という形。$${N}$$はフーリエ級数を打ち切る項数(+1)。

左辺で$${ \cos \theta }$$を掛けているので、フーリエ級数を打ち切っても、高次の項が出てくるが、そのような項は全部落とす。

 

これで、$${\cos(n\theta)}$$の係数ごとに比較していけば、最終的には、$${\mathbf{c} = (c_{0} , \cdots , c_{N})}$$に対して、

$${ B \dfrac{d \mathbf{c}}{dz} = A_{散乱} D \mathbf{c} }$$

という常備分方程式の形になるはず。$${B}$$が逆行列持ち、$${A_{散乱}}$$が$${z}$$に依存しなければ、行列のexponentialで計算できる。$${A_{散乱}}$$が$${z}$$に依存しても、積分は、難しくないだろう。

常微分方程式ではあるが、初期値はどうすればいいだろうか。つまり、大気圏外で$${I}$$の方向依存性に冠するデータがないという話。

太陽から地球に向かってきてる光だし、$${\theta=0}$$付近に集中しているだろうとは思う。デルタ関数的にしようと思って、$${c_{n}=1/N}$$を初期値にすると、輻射強度が負になる場所がでてきて嬉しくない。

今回は単純に$${c_{n}=1/N}$$で、負の輻射強度が出ないように、$${c_{0}}$$に余分な定数を足して、最後に積分が全輻射強度に等しくなるように規格化した。

 

レイリー散乱の微分断面積によって

$${S(\theta , \theta') = \dfrac{1}{3 \pi^2} (1 + \cos(\theta - \theta')^2)}$$

定数倍は、$${S}$$が確率密度なので、$${\theta,\psi}$$に冠する積分をして、丁度1になるように規格化するため。

あとは、

$${dw = \sin \theta d \theta d \psi}$$

で、積分すればいい。

計算すると分かるけど、$${n}$$が1と3の時は

$${\displaystyle \int_{0}^{\pi} S(\theta,\theta') \cos(n \theta') \sin(\theta') d \theta'}$$

は$${\sin(2 \theta)}$$に比例する。従って、$${c_{1}=c_{3}=0}$$でないといけない。

しかし、たとえ、初期値を0にしても、他の項の影響で、$${c_1=c_3=0}$$は維持できなくなる。とりあえず、$${c_{1}=c_{3}=0}$$となるように射影行列を両辺から掛けてしまうか、単に、ここの項はなかったことにしてしまうかの二択。後者の方が手間が少ないのは明らか。

 

$${n \neq 1,3}$$の時

$${ \displaystyle \int_{0}^{\pi} S(\theta,\theta') \cos(n \theta') \sin(\theta') d\theta' = \dfrac{1}{3\pi^2} \left( -\dfrac{1+\cos(n\pi)}{n^2-1} - \dfrac{1}{2}(P_{n}-Q_{n}) + \dfrac{1}{2}(P_{n}-Q_{n}) \cos (2 \theta) \right)}$$

$${P_{n} = \displaystyle \int_{0}^{\pi} \cos(x)^2\sin(x) \cos(n x) dx}$$

$${Q_{n} = \displaystyle \int_{0}^{\pi} \sin(x)^3 \cos(n x) dx}$$

となる。$${P_{n},Q_{n}}$$は解析的に出るだろうが数値積分してしまえばいいだろう。

 

最初の積分の方は、

$${\displaystyle \int_{0}^{\pi} S(\theta,\theta') \sin \theta' \cos(n \theta) d \theta' = \dfrac{1}{9 \pi^2} \left( \cos(|n-2| \theta) + \cos( (n+2)\theta) \right) }$$

なので、ここからも、より高次の項が出てしまう。勿論、これも無視する。

 

最後に、$${A_{散乱}}$$であるが、波長$${\lambda = \dfrac{c}{\nu}}$$に対して

$${A_{散乱} = \dfrac{K}{\lambda^4} }$$

となる。$${K}$$は定数。$${K}$$は本来$${z}$$に依存する(厳密には波長にも依存するが、可視光領域では波長依存性は小さい)が、ここでは定数として扱う。

決め方としては、波長$${\lambda = 600 \mathrm{nm}}$$で、高度$${H=100 \mathrm{km}}$$の時、$${KH = \lambda^4}$$になるとする。数字にそれほど意味はないが、オーダーとして、これくらいということ。

以上の準備を行って、光の進行距離を変えて鉛直真下方向に入射する放射照度の計算をしたのが以下のグラフ

空の色の近似スペクトルかもしれないもの

なんだか不思議な挙動をしているが、2km地点では、波長の短い成分が強くて、青そうな感じはしている。方向を考慮せず減衰だけ見てた時とは違いそう。

5km,10kmと進むにつれて、長波長成分が出てくる。つまり、赤くなる。

夕焼けとか朝焼けで、単純に距離が長くなるので赤くなるという説明は本当に正しいのかもしれない。

5kmとかエヴェレストより低いが、そもそも、$${A_{散乱}}$$は高度に依存して小さくなるので、それは考慮すべきだと思う。

物凄く簡略化した計算だが、きちんと輻射輸送方程式を解けば、散乱項はレイリー散乱で、空の色が計算できる可能性は出てきた。

しかし、一般には、もっと面倒だし、レイリーがこんな計算をできたわけもない。なにしろ、輻射輸送方程式が書かれたのは20世紀初頭になってからなので。

一つ言えることは、空の色は難しい。