標準光源Dの分光分布

標準光源D65は観測の統計から規格化された分光分布であり、色の計算ではよく使われる。
単純な黒体放射の分光分布とは異なり、太陽が大気を通過した光のことなので、
大気による波長の吸収が起きたり、大気の散乱による青空の光も含まれている。
補助光源としてD50、D55、D75等があり、印刷の色評価ではD50が使われている。
これらのDシリーズは「Illuminant series D」や「CIE昼光(D)」と呼ばれる。

色温度4000Kから25000Kに対する昼光を求める方法があり、
昼光の分光分布S0と固有ベクトルS1,S2にそれぞれ係数M1,M2を掛け合わせてS0と足し合わせることで求められる。
S(λ) = S0(λ) + S1(λ) *M1 + S2(λ) *M2

## 参考
↓ここに計算方法が書いてある
[Standard illuminant - Wikipedia](https://en.wikipedia.org/wiki/Standard_illuminant)

↓計算しているエクセルシートがあり、そこにs0,s1,s2の値も書いてある。
[RIT | Color Science | Resources | Useful Color data](http://www.rit.edu/cos/colorscience/rc_useful_data.php)
「Excel Daylight Series Calculator」をクリックすると「DaylightSeries.xlsx」がダウンロードできる

##プロットして確認
上のエクセルから作成した

```py3
import numpy as np
import matplotlib.pyplot as plt

#illuminantDの計算で使うs0,s1,s2の分光分布の入力 300-830nm 10nm step
s0 = np.array([[0.04,6,29.6,55.3,57.3,61.8,61.5,68.8,63.4,65.8,94.8,104.8,105.9,96.8,113.9,125.6,125.5,121.3,121.3,113.5,113.1,110.8,106.5,108.8,105.3,104.4,100,96,95.1,89.1,90.5,90.3,88.4,84,85.1,81.9,82.6,84.9,81.3,71.9,74.3,76.4,63.3,71.7,77,65.2,47.7,68.6,65,66,61,53.3,58.9,61.9]])
s1 = np.array([[0.02,4.5,22.4,42,40.6,41.6,38,42.4,38.5,35,43.4,46.3,43.9,37.1,36.7,35.9,32.6,27.9,24.3,20.1,16.2,13.2,8.6,6.1,4.2,1.9,0,-1.6,-3.5,-3.5,-5.8,-7.2,-8.6,-9.5,-10.9,-10.7,-12,-14,-13.6,-12,-13.3,-12.9,-10.6,-11.6,-12.2,-10.2,-7.8,-11.2,-10.4,-10.6,-9.7,-8.3,-9.3,-9.8]])
s2 = np.array([[0,2,4,8.5,7.8,6.7,5.3,6.1,2,1.2,-1.1,-0.5,-0.7,-1.2,-2.6,-2.9,-2.8,-2.6,-2.6,-1.8,-1.5,-1.3,-1.2,-1,-0.5,-0.3,0,0.2,0.5,2.1,3.2,4.1,4.7,5.1,6.7,7.3,8.6,9.8,10.2,8.3,9.6,8.5,7,7.6,8,6.7,5.2,7.4,6.8,7,6.4,5.5,6.1,6.5]])
Vaxis = np.array(range(300, 831, 10))

#色温度を指定してD光源の分光分布を合成
def calcIlluminantD(temperature):
    T = temperature
    xd = 0.0
    if (4000<= T and T<=7000):
        xd = 0.244063+0.09911*10**3/T + 2.9678*10**6/(T**2)-4.6070*10**9/(T**3)
    elif  (7000<= T and T<=25000):
        xd = 0.237040+0.24748*10**3/T + 1.9018*10**6/(T**2)-2.0064*10**9/(T**3)
    yd = -3.0*xd**2 + 2.87*xd -0.275
    m = 0.0241+0.2562*xd-0.7341*yd
    m1 = (-1.3515-1.7703*xd+5.9114*yd)/m
    m2 = (0.03-31.4424*xd+30.0717*yd)/m
    sd = s0 + s1*m1 + s2*m2
    return (sd.reshape(54))

#グラフのプロット
def plotIlluminantD(Tlist):
    for T in Tlist:
        sd = calcIlluminantD(T)
        plt.plot(Vaxis,sd,label=str(T)+'k')
    plt.title("illuminant D")
    plt.xlabel("wavelength [nm]")
    plt.ylabel("relative radiant power")
    plt.ylim(ymax=250)
    plt.xlim(xmax=830)
    plt.legend()
    plt.show()
    
Tlist = [25000,13000,9000,6500,5000,4000]
plotIlluminantD(Tlist)
```
![image.png](https://qiita-image-store.s3.amazonaws.com/0/205351/bb318576-ab4f-c2ef-9222-9676db49824c.png)

## 他参考
[Color Management addon/jp - RawPedia](http://rawpedia.rawtherapee.com/Color_Management_addon/jp) より引用
>初めの部分は“固定”の部分です。試験されたサンプル全ての平均値です。
2番目は“変化する”部分です。S1は雲があるかどうか、或いは直接的な太陽光の強さと位置で決まる“ブルー/イエロー”の変化量のことです。
3番目も“変化”する部分で、S2は蒸気や霧による湿度で決まる“ピンク/グリーン”の変化量のことです。
しかし、実際には、色温度と光源をベースにした2つの値、xDとyDを決定する単純な式で表わされます。

http://www.lrc.rpi.edu/programs/nlpip/lightinganswers/fullspectrum/comparisons.asp
[JIS Z 8720:2012 測色用の標準イルミナント（標準の光）及び標準光源](http://kikakurui.com/z8/Z8720-2012-01.html)
[Appendix 5: Relative Spectral Power Distributions of Illuminants - Measuring Colour - Hunt - Wiley Online Library](http://onlinelibrary.wiley.com/doi/10.1002/9781119975595.app5/pdf)