RでHR図 その1

">

(参考)
ヒッパルコス星表で太陽近傍恒星のHR図を描こう:必要項目の計算（２）恒星の表示色

上記の参考にしたサイトで、恒星の表示色の色度座標xyの計算に２種類の方法が紹介されています。この２つの方法のプログラムをRで作成して重ね合わせたグラフを作ってみました。

但し、「近似式を利用して色度図座標を計算」の近似式は、

## Planckian Locus approximate formula
x = -0.2661239*10^9/T^3-0.2343580*10^6/T^2+0.8776956*10^3/T+0.179910 （ 1667K＜T＜ 4000K）
x = -3.0258469*10^9/T^3+2.1070379*10^6/T^2+0.2226347*10^3/T+0.240390 （ 4000K＜T＜25000K）
y = -1.1063814*x^3-1.34811020*x^2+2.18555832*x-0.20219683 （ 1667K＜T＜ 2222K）
y = -0.9549476*x^3-1.37418593*x^2+2.09137015*x-0.16748867 （ 2222K＜T＜ 4000K）
y = +3.0817580*x^3-5.87338670*x^2+3.75112997*x-0.37001483 （4000K＜T＜25000K）

ですが、「T＜25000K」の条件は無視します。

比較のグラフ(1400Kから50000Kまで)

計算速度は、近似式を使った方がかなり速く1700K以上では「黒体放射の分光分布とCIEの等色関数から色度図座標を計算」した値とほぼ等しくなります。
色度座標xyの計算は、1700K以上では近似式プログラム、1700Kより低い値では「黒体放射の分光分布とCIEの等色関数から色度図座標を計算」プログラムを使う。という方法をとりたいと思います。

Rコード

自作プログラムの定義

a) 黒体放射の分光分布とCIEの等色関数から色度図座標を計算

BBR2xy=function(T){
    lambda<-c(380,385,390,395,400,405,410,415,420,425,430,435,440,445,450,455,460,465,470,475,480,
        485,490,495,500,505,510,515,520,525,530,535,540,545,550,555,560,565,570,575,580,585,590,595,
        600,605,610,615,620,625,630,635,640,645,650,655,660,665,670,675,680,685,690,695,700,705,710,
        715,720,725,730,735,740,745,750,755,760,765,770,775,780)
    x<-c(0.0014,0.0022,0.0042,0.0076,0.0143,0.0232,0.0435,0.0776,0.1344,0.2148,0.2839,0.3285,0.3483,0.3481,
        0.3362,0.3187,0.2908,0.2511,0.1954,0.1421,0.0956,0.058,0.032,0.0147,0.0049,0.0024,0.0093,0.0291,
        0.0633,0.1096,0.1655,0.2257,0.2904,0.3597,0.4334,0.5121,0.5945,0.6784,0.7621,0.8425,0.9163,0.9786,
        1.0263,1.0567,1.0622,1.0456,1.0026,0.9384,0.8544,0.7514,0.6424,0.5419,0.4479,0.3608,0.2835,0.2187,
        0.1649,0.1212,0.0874,0.0636,0.0468,0.0329,0.0227,0.0158,0.0114,0.0081,0.0058,0.0041,0.0029,0.002,
        0.0014,0.001,0.0007,0.0005,0.0003,0.0002,0.0002,0.0001,0.0001,0.0001,0)
    y<-c(0,0.0001,0.0001,0.0002,0.0004,0.0006,0.0012,0.0022,0.004,0.0073,0.0116,0.0168,0.023,0.0298,0.038,
        0.048,0.06,0.0739,0.091,0.1126,0.139,0.1693,0.208,0.2586,0.323,0.4073,0.503,0.6082,0.71,0.7932,
        0.862,0.9149,0.954,0.9803,0.995,1,0.995,0.9786,0.952,0.9154,0.87,0.8163,0.757,0.6949,0.631,0.5668,
        0.503,0.4412,0.381,0.321,0.265,0.217,0.175,0.1382,0.107,0.0816,0.061,0.0446,0.032,0.0232,0.017,
        0.0119,0.0082,0.0057,0.0041,0.0029,0.0021,0.0015,0.001,0.0007,0.0005,0.0004,0.0002,0.0002,0.0001,
        0.0001,0.0001,0,0,0,0)
    z<-c(0.0065,0.0105,0.0201,0.0362,0.0679,0.1102,0.2074,0.3713,0.6456,1.0391,1.3856,1.623,1.7471,1.7826,
        1.7721,1.7441,1.6692,1.5281,1.2876,1.0419,0.813,0.6162,0.4652,0.3533,0.272,0.2123,0.1582,0.1117,
        0.0782,0.0573,0.0422,0.0298,0.0203,0.0134,0.0087,0.0057,0.0039,0.0027,0.0021,0.0018,0.0017,0.0014,
        0.0011,0.001,0.0008,0.0006,0.0003,0.0002,0.0002,0.0001,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0)
    c2 = 1.438775 * 10^7
    M=(lambda^5) * (exp(c2/(lambda * matrix(rep(T,length(lambda)),byrow=T,ncol=length(T)) ) ) -1 )
    X=colSums(x/M)
    Y=colSums(y/M)
    Z=colSums(z/M)
    x0 = X / (X + Y + Z)
    y0 = Y / (X + Y + Z)
    return(data.frame(x=x0,y=y0))
}

(b)近似式を利用して色度図座標を計算

#(注意) T < 25000 という条件はなくした。
temp2xy=function(T){
    x=ifelse( T>1667 & T<4000 ,
        -0.2661239*10^9/T^3-0.2343580*10^6/T^2+0.8776956*10^3/T+0.179910 ,
        -3.0258469*10^9/T^3+2.1070379*10^6/T^2+0.2226347*10^3/T+0.240390 )
    y=ifelse(T>1667 & T<2222 ,
        -1.1063814*x^3-1.34811020*x^2+2.18555832*x-0.20219683 ,
        -0.9549476*x^3-1.37418593*x^2+2.09137015*x-0.16748867 )
    y=ifelse(T>=4000 ,
        +3.0817580*x^3-5.87338670*x^2+3.75112997*x-0.37001483 ,
        y )
    return(data.frame(x=x,y=y))
}

計算し、グラフ作成

# プロットする温度
T= seq(1400,50000,100)
# 黒体放射の分光分布とCIEの等色関数から色度図座標を計算
xy=BBR2xy(T)
x0 = xy$x
y0 = xy$y
# 近似式を利用して色度図座標を計算
xy=temp2xy(T)
x1 = xy$x
y1 = xy$y
# グラフ作成
# png("Ttoxy01.png",width=800,height=400)
par(mfrow=c(1,2),mar=c(4,4,3,2))
plot(x=T,y=x0,las=1,xlab="温度(K)",ylab="x",main="表面温度(K)から色度図座標xyを計算 x",bty="n")
box(bty="l",lwd=2)
lines(T,x1,col="red")
axis(1,at=1667,labels="1667",line=-0.7)
legend("topright",inset=c(0.05,0),legend="黒体放射",pch=1,bty="n",x.intersp = 1.5)
legend("topright",inset=c(0.06,0.045),legend="近似式",lty=1,col="red",bty="n")
abline(v=1667,lty=2)
#
plot(x=T,y=y0,las=1,xlab="温度(K)",ylab="y",main="表面温度(K)から色度図座標xyを計算 y",bty="n")
box(bty="l",lwd=2)
lines(T,y1,col="red")
axis(1,at=1667,labels="1667",line=-0.7)
legend("topright",inset=c(0.05,0),legend="黒体放射",pch=1,bty="n",x.intersp = 1.5)
legend("topright",inset=c(0.06,0.045),legend="近似式",lty=1,col="red",bty="n")
abline(v=1667,lty=2)
par(mfrow=c(1,1))
# dev.off()