Article Outline

# 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()
```