空の色を計算する試み
今日の空は、雲ひとつない快晴だったので、帰ったら、空の色を計算しようと思った。
空が青いのは、レイリー散乱のせいというのは、常識みたいに語られている。そう言われても、何も分かった気がしないし、「エライ人が言うならそういうもんか」という以上の感想しかない。
レイリー散乱で散乱された光が目に入っているのだという説明が正しいんだとしたら、レイリー散乱で消散した分の放射照度を色に変換してやれば、空の色になるはずじゃないだろうか。
黒体輻射の色
まずは、大気圏外の太陽光のスペクトル分布を計算する。
太陽光のスペクトル分布は、約6000Kの黒体輻射でよく近似できるとされている。そのスペクトル分布の式は、1900年にプランクが得た。Wikipediaの式を持ってくると、放射輝度は
$${ I(\lambda) = \dfrac{2 h c^2}{\lambda^5} \dfrac{1}{e^{hc/\lambda k_B T} - 1} }$$
となる。$${\lambda,h,c,k_B,T}$$は、それぞれ、光の波長、プランク定数、光速、ボルツマン定数、温度。全体として、単位は、$${\mathrm{W/m^3 \cdot sr}}$$となる。単位時間、単位面積、単位立体角あたり単位波長の光が通過するエネルギー量という解釈。
これは太陽表面での計算式なので、地球に到達する時には、(太陽半径)/(軌道半径)の2乗を掛けた分まで減っている。また、方向についても(ランバート則を考慮して)積分しておく必要もある。
全部考慮して計算したグラフが以下。
ASTM G173-03というのが、よく使われる観測データらしいので、一緒にプロットしておいた。いうほど黒体輻射か？という気もするが。
折角なので、色に変換する。スペクトル分布から色に変換するには、XYZ色空間なるものに変換してから、RGBに変換することが推奨されているらしい。放射輝度からXYZ色空間に変換するための等色関数(color matching function)というものが定義されている。詳しくは、Wikipedia。
誰かが1931年に頑張って決めたらしい。適当に検索して、CIE 1931 COLOUR-MATCHING FUNCTIONSから数表を持ってきて、数値積分する。
そして、それをRGBに変換する手順もWikipedia。
波長ごとの等色関数のXYZ値をRGBに変換してやった(つまり、各波長を平均とするデルタ関数的なスペクトル分布の色を計算したことに相当)ものが以下。期待通り、虹ができた。
で、黒体輻射近似した太陽輻射の色を計算する。生のスペクトル強度を使うと、真っ白になってしまう。人間も太陽を直視できないことを思えば、仕方ない。というわけで、強度を一定率弱める必要がある。ここでは、50%にした。
計算してできた太陽輻射の色。
なんか肌色。輻射温度を5770Kとして計算した太陽輻射のRGB値は、(255,234,182)。検索して出てきたカラーサイト.comによると、肌色のRGBは(254,220,189)。
人間の肌の色は、主にメラニン色素による(ヘモグロビンとかの影響もあるかもしれないが)とはいうが、メラニンは光を吸収し、人間が見るのは、吸収されず反射された光。なので、メラニン量が少ない人の肌の色は、太陽輻射の色に近付くのかもしれない。従って、肌色が太陽輻射の色なのは、妥当な気もする。
太陽輻射自体は、RGBの内、赤が最も強い。リファレンスデータであるAM1.5の色も計算した。
黒体輻射より少し暗いが、同じ系統の色。この色が直接、空の色になってるというわけではない。
特に本題とは関係ないけど、輻射温度を変えてみた。放射照度を全波長で合計すると、温度の4乗に比例するので、強度を弱める割合は、それに合わせて変えている。
レイリー散乱による消散の計算
大気には、散乱とは別に吸収による減衰もある。散乱と吸収は同時に進むので、吸収が大きい場合は、無視すると、散乱による消散分を過大に見積もることになってしまう。が、レイリー散乱のみで、どこまで説明できるか見るために、ここでは考えない。
レイリーは1871年に"On the Light from the Sky, Its Polarization and Color"という論文を書いたそうだ。
論文の冒頭に
IT is now, I believe, generally admitted that the light which we receive from the clear sky is due in one way or another to which divert the light from its regular course. On this point the experiments of Tyndall with precipitated clouds seem quite decisive.
と書いてあって、当時、空の色が青いのは、小さな粒子(small suspended particles)で散乱されるからという認識はあったらしい。
1899年にも"On the transmission of light through an atmosphere containing small particles in suspension, and on the origin of the blue of the sky"という長いタイトルの論文を書いていて、ここでは、この粒子は分子だと主張した。
レイリー散乱は、散乱球体が波長に比べて十分小さい時成立する。光の波長は500nm前後な一方で、ボーア半径は0.05nm。また、標準状態で22.4リットルの空気中にアボガドロ数個の分子があるとして、一つの分子が占める領域の平均体積は、半径2nmの球と等しい。気体は、粒子間距離が広いことも踏まえれば、条件を満たしていると考えていいだろう。
例によって、Wikipediaをカンニングすると、微小球一つあたりの全散乱断面積$${\sigma_R}$$は
$${ \sigma_{R} = \dfrac{8 \pi}{3} \left( \dfrac{2\pi}{\lambda} \right)^4 d^6 \left( \dfrac{n^2-1}{n^2+2} \right)^2 }$$
で書ける。$${\lambda,d,n}$$は、光の波長、球半径、屈折率。"断面積"という名前の通り、単位は、長さの2乗。単位体積あたりの粒子数密度$${N}$$を掛けて
$${\alpha = N \sigma_{R}}$$
がレイリー散乱による消散係数になる。長さの逆数の次元なので、$${1/\alpha}$$進んだ時に、光の強度が$${1/e}$$になる。
数値例を見たいところだが、散乱球半径の6乗に比例するので、この項が2倍違うだけで、散乱断面積は64倍も違ってくる。
レイリーは1899年の論文の中で、波長600nmでの試算を試みているが、当時はアボガドロ数が分かっていなかった(まだ原子や分子の存在に批判があった時代)。それでも、$${1 / \alpha = 83 \mathrm{km}}$$という数字を得ている。大体、$${1.2 \times 10^{-5} \mathrm{m}^{-1}}$$に相当する。
現在はアボガドロ数が分かっているが、光の波長600nm、大気屈折率1.0003、散乱球半径1nmとした時、
$${\sigma_{R} \approx 4.03 \times 10^{-33} \mathrm{m}^2}$$
で、アボガドロ数を$${N_{A}}$$とすると、地表付近の粒子数密度は
$${N = N_{A} \times \dfrac{1000}{22.4} \approx 2.69 \times 10^{25} \mathrm{m}^{-3}}$$
と計算できて、消散係数は
$${\alpha = N \sigma_{R} \approx 1.08 \times 10^{-7} \mathrm{m}^{-1}}$$
となる。標準状態で気体分子が占める体積は、半径2nm程度の球で、これを散乱体とした場合、
$${\alpha \approx 6.93 \times 10^{-6} \mathrm{m}^{-1}}$$
で、レイリーの数値の半分くらい。144kmで1/eになる程度。
現代から考えると、マイクロ波が粉塵で散乱されるみたいな状況と、可視光が分子に散乱されるみたいな状況を同じように扱っていいかは怪しい。前者は、レイリーの時代の物理学で扱えて、レイリーの散乱公式を適用して問題ないはず。
後者については、レイリーの時代の物理学で扱える範囲を超えてるので、結果の妥当性は自明でない。少なくとも、分子の半径とかいうパラメータは出てくるべきでないだろう。
とりあえず今は、レイリー散乱公式の導出に煩わされたくないので、波長や屈折率への依存性は可視光でもOKとされてるのを信じて、
$${\sigma_R = \alpha_{0} \dfrac{1}{\lambda^4} \left( \dfrac{n^2-1}{n^2+2} \right)^2 }$$
という形で、$${\alpha_0}$$というパラメータを定義して、このパラメータを勘で与えることにする。
レイリーの古典的結果では、
$${\alpha_{0} = \dfrac{128 \pi^5 d^6}{3}}$$
なので、おおまかな目安は決められる。
レイリー散乱のみがあるとして地表でのスペクトル分布がどうなるか計算する。大気圏外から地表まで波長ごとの減衰を計算してやればいいだろう。しかし、高度によって、空気は薄くなり、分子数が減るので、それを考慮して積分してやる必要があるように思う。
また、空気屈折率には波長依存性がある。屈折率も、温度や圧力に依存するので、高度依存性がある。めんどくせー
湿度依存性なんかもあるだろうが、経験的に、湿度のせいで、空の色が全然違うということはないので、無視していいだろう。
傾向を見るために、そういうのを全部無視して計算してみる。消散係数は波長の4乗に逆比例し、定数$${K}$$を使って
$${\alpha = \dfrac{K}{\lambda^4}}$$
とする。距離$${H}$$だけ進んだ時の放射強度は
$${ I(\lambda) e^{-\dfrac{K}{\lambda^4} H} }$$
で計算できる。
$${K = 1296000 (\mathrm{nm^3/m)}}$$
$${H = 100 \mathrm{km}}$$
とすると、$${\lambda=600 \mathrm{nm}}$$の時に$${ \dfrac{K}{\lambda ^4} H = 1}$$になる。
散乱によって失われた放射強度
$${ S(\lambda) = I(\lambda) (1 - e^{-\dfrac{K}{\lambda^4} H}) }$$
をプロットしよう。
距離$${H}$$を変えた時に、失われた放射強度のRGB値を見ると
のようになる。距離が大きくなるにつれて、青色は飽和して、赤や緑に抜かれる。色にすると
左から右にいくほど、距離が大きくなる。どこを見ても、あんまり空っぽい色は出ない。
屈折率の波長依存性や高度依存性などを考慮していくと、空の色が出たりするだろうか。
屈折率の波長依存性は、Sellmeierの式という経験式(多分)があるそうだ。常温常圧空気については、The Refractive Index of Airとかいう論文の式は、その変種だろう。
この式を使って計算すると、$${\dfrac{n^2-1}{n^2+2}}$$の値は380nmと780nmで、0.3%ほどしか違わない。屈折率の波長依存性は可視光に対しては無視して構わないだろう。
温度$${T}$$と圧力$${p}$$依存性は、ローレンツ・ローレンツの式をTaylor展開した
$${n \approx \sqrt{1+ \dfrac{3 A p}{RT}} \approx 1 + \dfrac{3 A p}{2 R T}}$$
でいいだろう。$${R}$$は気体定数。$${A}$$はモル屈折率で、
$${A \approx \dfrac{n^2-1}{3} \dfrac{RT}{p}}$$
なので、常温常圧の屈折率から決定できる。
次に、温度、圧力、分子数の高度依存性を知りたい。温度は100mにつき、0.65度上昇とかいうのは対流圏の話で、成層圏、中間圏、熱圏というのがあって、そこの温度はあがったり下がったりするらしい。その上を外気圏というそうだ。
熱圏は、途中(100kmあたりより上空)から宇宙空間扱いで、気体分子も希薄らしいので、その下で、減衰を考えればいいだろう。空の色を出すのに、どこまでのの高度の散乱が効いてくるかは直感的に分からない。成層圏はよく聞くし、なんか考慮した方がよさそう。
というか、比較的高緯度限定で夜間とはいえ、オーロラとかが見えるんだから、電離層からの光も割と届くのかもしれない。電離層は熱圏の上部。てことで、高度100kmからの散乱を全部入れて良さそうに思える。
超高層大気についてという1958年の論文の表5に、高度10km刻みで、分子数、密度、気圧、温度のデータがあるので、この表の数値を線形補間して使うことにする。古いデータだが、1950年代にも、空は青かっただろう。
消散係数$${\alpha}$$が定数の時は、進行距離$${x}$$に対して、
$${I(x) = I(0) e^{-\alpha x}}$$
だが、$${\alpha}$$が場所に依存する場合は、どうすればいいのか。仮に、進行区間を$${x_{0}=0 \lt \cdots \lt x_{n}}$$と分割すると
$${I(x_{n}) \approx I(x_0) \exp \left( -\displaystyle \sum_{i=0}^{n-1} (x_{i+1}-x_{i})\alpha(x_{i}) \right)}$$
とかだろう。分割を細かくして極限を取れば
$${I(x) = I(0) \exp\left( -\displaystyle \int_{0}^{x} \alpha(u) du \right)}$$
でよさそう。
光が斜めから入射してくる場合、光路長が長くなり、消散が大きいと考えられるが、垂直に入射してくるとして、波長毎に消散率を計算した。リファレンスデータの減衰率と比較したのが、以下のグラフ。
減衰要因は他にも沢山あるだろうから、特に矛盾はない。それで、レイリー散乱で失われた分の放射照度を色に変換した。そのままでは暗すぎるので、4倍にスケールして明るくしている。
色々頑張って計算したけど、高度依存性を完全無視した時と様子は変わらない。RGB値は、(127, 128, 151)と、青が一番強いが、かなり均衡している。空色のRGBを検索すると、(152,187,219)とか出てくるので、赤味が強すぎ、青味が足りない。
元々の太陽光は赤成分が強いので、減衰率が低くても、絶対的な減衰量は多いってことだと思う。
言い訳とか
結局、空の色は出なかった。レイリー散乱だけで、空の色を説明していいのか、私には分からない。
2012年の論文The Color of the Skyには、レイリーが$${\dfrac{1}{\lambda^4}e^{-a/\lambda^4}}$$に比例する減衰則を導いたと書いてある(オゾンによる吸収と、朝焼け、夕焼けはエアロゾルの影響もあるみたいなことも書いてあるが)。
1871年のレイリーの論文を見ると
$${ I = I_{0} e^{-k \lambda^{-4} x} }$$
$${ I \propto \lambda^{-4} e^{-k \lambda^{-4} x} }$$
の両方の式がある。前者は普通の減衰なので理解できるが、後者の式の導入は、ad hocで言葉による謎の説明のみがある。
散乱された光がそのまま眼球に到達すると仮定して(前者の式を導いて)きたが、何らかのmodificationを受けるみたいなことが書いてあるが、何を想定した文章か意味が分からない。1899年の論文には、後者の式は出てない。
指数関数の外の項によって、短い波長と長い波長の両方で、減衰が出るようになる。
何らかの理由で、眼球に届く光には、そういう因子が追加されるなら、比例係数を計算できるべきだと思う。私が何か勘違いしてるか見落としてるのかもしれないけど。
大体、レイリー散乱自体は、火星や金星の大気でも起きるはず。火星と金星の大気圧は、地球の1/100と100倍らしいから、分子数密度も違うだろうし、大気組成も違うので、屈折率も違うだろうが、波長依存性自体は同じ。
レイリーの時代には知り得なかった事柄だろうけど、火星や金星の空が青く見えたことは多分ない。私も火星や金星に行ったことないので、あれが砂塵とかのせいで、それがなければ、青く見えたりするのかもしれないけど。
太陽輻射の色計算プログラムたち
グラフを作ったりするのに使用したpythonプログラムたち。
import math
def xyz2rgb(xyz):
def gamma_correct(c):
if c<=0.0031308:
return 12.92*c
else:
return 1.055*math.pow(c , 1/2.4) - 0.055
MATRIX_SRGB_D65 = np.array( [[3.2404542, -1.5371385, -0.4985314],[-0.9692660, 1.8760108, 0.0415560],[0.0556434, -0.2040259, 1.0572252]] )
rgbL = np.dot(MATRIX_SRGB_D65 , xyz)
srgb = np.vectorize(gamma_correct)(rgbL)
crgb = np.clip(srgb,0,1) * 255
return crgb
def I(L,T):
"""
L: 光の波長 in nm
T: 輻射温度 in K
"""
A = 2 * 6.62607015e-34 * ( 299792458 ** 2) #-- h c^2
B = 6.62607015e-34 * 299792458 / 1.380649e-23 #-- h c / k_B
R = 6.960e8 #-- 太陽半径(m)
D = 1.496e11 #-- 軌道半径(m)
val = math.pi * 1.0e-9*A/(L*1.0e-9)**5 * 1/(math.exp(B/(L*1.0e-9*T)) - 1) #-- W/m^2/nm
return (val * (R/D)**2)
n_tbl = [2.55e25, 8.60e24, 1.85e24, 3.71e23, 8.32e22, 2.25e22, 7.26e21, 2.09e21, 4.50e20, 8.31e19, 1.65e19] #-- per m^3
p_tbl = [1.013e5 , 2.649e4 , 5.528e3 , 1.185e3 , 2.997e2 , 8.784e1 , 2.581e1 , 6.320e0 , 1.223e0 , 2.257e-1 , 4.629e-2] #-- Pa
t_tbl = [288.16,223.26,216.66,231.24,260.91,282.66,257.55,219.33,196.86,196.86,203.06] #-- K
def loc_att(L , h,alpha_0):
"""
h : height / m
L: wavelength of light / nm
alpha_0: "standard attenuation coefficient"
"""
H = 10000
idx = int(h/H)
n = (n_tbl[idx] * (H*idx+H-h) + n_tbl[idx+1]*(h-H*idx)) / H #-- numer of molecules at given height
p = (p_tbl[idx] * (H*idx+H-h) + p_tbl[idx+1]*(h-H*idx)) / H #-- pressure at given height
t = (t_tbl[idx] * (H*idx+H-h) + t_tbl[idx+1]*(h-H*idx)) / H #-- temperature at given height
A = 4.85306070022225e-07 #--standard molar refractive index of air / gas constant
mu = 1 + 1.5 * A * p / t #-- refractive index of air at given height
ret = alpha_0 * (n/(L*1.0e-9)**4)*((mu**2-1)/(mu**2+2))**2
return ret
def att(L,alpha_0=1.0e-48):
max_height = 100e3 #-- 100km
a_tot = np.sum(np.vectorize(lambda x:loc_att(L , x , alpha_0))(np.arange(0,max_height)))
return math.exp(-a_tot)
"""
converting blackbody radiation spectrum to rgb color
"""
class BBColor:
def __init__(self , s , ratio=0.5):
Ls,xf,yf,zf = [],[],[],[]
for line in s.split("\n"):
dat=line.split(",")
if len(dat)!=4:continue
Ls.append(float(dat[0]))
xf.append(float(dat[1]))
yf.append(float(dat[2]))
zf.append(float(dat[3]))
self.Ls = Ls
self.xf = xf
self.yf = yf
self.zf = zf
self.x_sum = sum(xf)
self.y_sum = sum(yf)
self.z_sum = sum(zf)
self.ratio = ratio
def __call__(self,T,residual=False):
X,Y,Z = 0 , 0, 0
Ls,xf,yf,zf = self.Ls, self.xf, self.yf, self.zf
for n,(L0,L1) in enumerate(zip(Ls[:-1],Ls[1:])):
assert(L1 > L0)
dL = L1-L0
L_mid = (L1+L0)/2
Ival = I( L_mid , T)
if residual:
Ival *= (1 - att(L_mid))
X += 0.5*(xf[n]+xf[n+1])*Ival * dL * self.ratio
Y += 0.5*(yf[n]+yf[n+1])*Ival * dL * self.ratio
Z += 0.5*(zf[n]*zf[n+1])*Ival * dL * self.ratio
xyz = np.array([X/self.x_sum , Y /self.y_sum , Z/self.z_sum])
return xyz2rgb(xyz).astype(np.uint8)
import urllib
import cv2
import numpy as np
if __name__=="__main__":
dataurl = "https://files.cie.co.at/CIE_xyz_1931_2deg.csv"
T = 5770 #-- radiation temperature in K
s = urllib.request.urlopen(dataurl).read() #-- to get color matching function table
s = s.decode('utf-8')
rgb = BBColor(s)
#-- color of black body radiation
r,g,b = rgb(T)
img = np.ones((100,400,3),dtype=np.uint8)
img[:,:,0] *= np.uint8(r)
img[:,:,1] *= np.uint8(g)
img[:,:,2] *= np.uint8(b)
img = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)
cv2.imwrite("solar_color_5770K.png",img)
#-- temperature variation
for temp in [5000,9000,10000]:
rgb.ratio = 0.5* (5770/temp)**4
r,g,b = rgb(temp)
img = np.ones((100,400,3),dtype=np.uint8)
img[:,:,0] *= np.uint8(r)
img[:,:,1] *= np.uint8(g)
img[:,:,2] *= np.uint8(b)
img = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)
cv2.imwrite("solar_color_{}K.png".format(temp),img)
おまけ: 太陽輻射温度の計算
輻射温度$${T}$$の計算は、どこにでも書いてあるが、太陽定数1366$${\mathrm{W/m^2}}$$から出す。
光の周波数を$${\nu}$$として、
$${ I(\lambda) d \lambda = - \tilde{I}(\nu) d \nu}$$
となる$${\tilde{I}(\nu)}$$を定義しておく。マイナスが付いてるのは、$${\lambda = \dfrac{c}{\nu}}$$で積分区間の大小が逆転するから事前に考慮している。物理的には$${ \tilde{I}(\nu) }$$が先に導けて、そこから$${I(\lambda)}$$が出る。
$${ \displaystyle \int_{0}^{\infty} \tilde{I}(\nu) d\nu \displaystyle \int_{半球} \cos(\theta) d\Omega }$$
を計算すると、単位時間、単位面積あたりに太陽から放射されるエネルギー量が分かる。$${d\Omega}$$は立体角。黒体表面がランベルト面だとして、放射面から垂線方向に放射されるエネルギーを数えている。で、
$${ \displaystyle \int_{半球} \cos(\theta) d\Omega = \displaystyle \int_{0}^{2\pi} d\phi \displaystyle \int_{0}^{\pi/2} \cos(\theta) \sin(\theta) d \theta = \pi}$$
$${ \displaystyle \int_{0}^{\infty} I(\dfrac{c}{\nu}) d\nu = \dfrac{2 k_{B}^4 T^4}{c^2 h^3} \displaystyle \int_{0}^{\infty} \dfrac{x^3}{e^{x}-1} dx }$$
多少初等的でない積分
$${ \displaystyle \int_{0}^{\infty} \dfrac{x^3}{e^{x}-1} dx =\dfrac{\pi^4}{15} }$$
を飲み込めば、
$${ \displaystyle \int_{0}^{\infty} \tilde{I}(\nu) d\nu \displaystyle \int_{半球} \cos(\theta) d\Omega = \dfrac{2 \pi^5 k_{B}^4}{15 c^2 h^3} T^4}$$
で、ステファン・ボルツマン則を得る。太陽半径、太陽定数、地球と太陽の距離を$${R_{s} , E_{s} , D}$$とすると
$${4 \pi D^2 E_{s} = 4 \pi R_{s}^2 \dfrac{2 \pi^5 k_{B}^4}{15 c^2 h^3} T^4}$$
なので、$${E_{s} = 1366 \mathrm{W/m^2}}$$として、天文定数と物理定数から、輻射温度$${T \approx 5770 \mathrm{K}}$$が決まる。