HOME/Articles/

RでHR図 その2

Author
統計ソフトRの備忘録３
Published
2023/01/09
Updated
2023/01/09
RHRVizieR
Article Outline

RでHR図 その2

">Hits

(参考)
ヒッパルコス星表で太陽近傍恒星のHR図を描こう

(使用したデータ元)
VizieR

(使用するデータ)
hip_hr.csv

(作成したデータ)恒星の絶対等級、表面温度、太陽半径比、恒星の色(RGB)
hip_HRrgb.csv

HR図:多くの星をプロットするために条件をかなりゆるくしている(視差角:誤差100％未満,B-V色指数:hip$BV_se<0.6)

HR01.png

HR02.png

おまけ：コードはなし

(作成したデータ)hip_HRrgb.csv を使った。

hip_equatorial_T01.png

hip_galactic_T01.png

Rコード

(使用するデータ)hip_hr.csv をダウンロードし、作業フォルダに置く。

自作プログラムの定義

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
###### 最初は for文を使っていたが激遅だったので書き換えた。######
    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))
}

データの読み込み。恒星の絶対等級、表面温度、太陽半径比、恒星の色(RGB)を計算

表面温度Ｔの計算には、Ballesterosの式: T= 4600 (1/(0.92 * bv + 1.7) + 1/(0.92 * bv + 0.62))を使った。

hip<- read.csv("hip_hr.csv")
nrow(hip)
# [1] 117955
# 年周視差の値が負になっている行を削除(正のデータのみにする)
hip= hip[hip$p>=0  ,]
#
#視差角とB-V色指数について、標準誤差から精度の低いデータを削除する
# 多くの星をプロットするために条件をかなりゆるくしている。
# 視差角:誤差100％未満
hip<- hip[abs(hip$p_se/hip$p) < 1 ,]
# B-V色指数:hip$BV_se<0.6
hip<- hip[ hip$BV_se < 0.6 ,]]
#
# pがNAのデータを削除
hip<-hip[!is.na(hip$p),]
# BVがNAのデータを削除
hip<-hip[!is.na(hip$BV),]
nrow(hip)
# (1)恒星の絶対等級Ｍの計算
# d；距離（pc）
# d=1000/p    p；視差角（ミリ秒角）
# 絶対等級    
# M＝m＋5−5*log10d    m:みかけの等級
d= round(1000/hip$p,2)
M= round(hip$m+5-5*log10(d),2)
hip$d<- d
hip$M<- M
# 絶対等級Ｍが　-Inf のデータは除く
hip<- hip[!(hip$M== -Inf),]
nrow(hip)
# (2)表面温度Ｔの計算
#  Ballesterosの式: T= 4600 (1/(0.92 * bv + 1.7) + 1/(0.92 * bv + 0.62))
T <- round(4600 * ( 1/(0.92*hip$BV+1.7) + 1/(0.92*hip$BV+0.62) ),2)
hip$T<- T
# (3)恒星の半径 R/Rs（太陽半径との比）
# Ls；太陽の光度，Rs；太陽の半径，Ts；太陽の表面温度（5800K），Ms；太陽の絶対等級（4.83mag）
# L；恒星の光度，R；恒星の半径，T；恒星の表面温度，M；恒星の絶対等級
# Dm=Ms-M    Dm；絶対等級の差
# L/Ls=10^(0.4*Dm)    L/Ls；光度比  10^((Ms-hip$M)*0.4)
# R/Rs=(Ts/T)^2*(L/Ls)^(1/2)    R/Rs；半径比 (Ts/T)^2 * sqrt(L/Ls)
Ts=5800
Ms=4.83
Dm=Ms-hip$M
L_Ls=10^(0.4*Dm)
R_Rs= round((Ts/hip$T)^2*(L_Ls)^(1/2),2)
hip$R_Rs<- R_Rs
summary(hip)
# hip$Tが欠損値のデータを削除
hip<- hip[!(is.na(hip$T)) , ]
# 近似式を使うデータ
hip1<- hip[hip$T>1700,]
# 黒体放射の分光分布を使うデータ
hip2<- hip[hip$T<=1700,]
#
if (length(hip1$T) !=0 ){
    xy=temp2xy(hip1$T)
    hip1$color_x = round(xy$x,6)
    hip1$color_y = round(xy$y,6)
} else {
    hip1=NULL
}
#
if (length(hip2$T) !=0 ){
    xy=BBR2xy(hip2$T)
    hip2$color_x = round(xy$x,6)
    hip2$color_y = round(xy$y,6)
} else {
    hip2=NULL
}
# hip1とhip2をつないでhipに戻す。
hip<- rbind(hip1,hip2)
x0<- hip$color_x
y0<- hip$color_y
# 2)３刺激値XYZの計算
Y=rep(1,length(x0))
X=Y*(x0/y0)
Z=(Y/y0)*(1-x0-y0)
# (3)RGB 値の計算 ３刺激値XYZからRGB値を計算します。
# sRGB(D65)
R =　3.2410*X - 1.5374*Y - 0.4986*Z
G = -0.9692*X + 1.8760*Y + 0.0416*Z
B =　0.0556*X - 0.2040*Y + 1.0507*Z
# マイナスの値は0にする。
R= ifelse( R<0 , 0 , R )
G= ifelse( G<0 , 0 , G )
B= ifelse( B<0 , 0 , B )
# (4)γ補正(windows)
R1=R^(1/2.2)
G1=G^(1/2.2)
B1=B^(1/2.2)
# 1以上の値は1にする。
R1= ifelse( R1>1 , 1 , R )
G1= ifelse( G1>1 , 1 , G )
B1= ifelse( B1>1 , 1 , B )
hip$R = round(R1,4)
hip$G = round(G1,4)
hip$B = round(B1,4)

HR図作成

cex=0.1で統一
# 半径 R/Rs（太陽半径との比）降順に並べ替え
hip=hip[order(-hip$R_Rs),]
# png("HR01.png")
cex=0.1
par(bg="gray5", fg="white", col.main="white", col.lab="white", col.axis="white")
par(mar=c(4,4,3,1.5))
plot(NA,xlim=c(max(hip$T),min(hip$T)),ylim=c(max(hip$M),min(hip$M)*1.1),las=1,bty="n",xlab="",ylab="",log = "x")
mtext("有効放射温度(K):対数",1,line=2.2)
mtext("絶対等級",2,line=2.2)
abline(h=seq(-15,15,5),col="white",lwd=0.5)
abline(v=c(10^3*(2:9),seq(1e4,2e4+2e3,1e3)),col="white",lwd=0.5)
box(bty="l",lwd=2)
points(x=hip$T,y=hip$M,pch=21,cex=cex,
    col=rgb(hip$R,hip$G,hip$B,alpha=1),bg=rgb(hip$R,hip$G,hip$B,alpha=0.8))
# 太陽の表面温度（5800K），Ms；太陽の絶対等級（4.84）
points(x=5800,y=4.84,col="black",cex=0.3)
text(x=5800,y=4.84,labels="太陽",col="black",cex=0.8,pos=3)
title("ヒッパルコスカタログ(2007)によるH-R図")
par(bg="white", fg="black", col.main="black", col.lab="black", col.axis="black")
# dev.off()
cex=(hip$R_Rs/10^2)/2
# png("HR02.png")
cex=(hip$R_Rs/10^2)/2
par(bg="gray5", fg="white", col.main="white", col.lab="white", col.axis="white")
par(mar=c(4,4,3,1.5))
plot(NA,xlim=c(max(hip$T),min(hip$T)),ylim=c(max(hip$M),min(hip$M)*1.1),las=1,bty="n",xlab="",ylab="",log = "x")
mtext("有効放射温度(K):対数",1,line=2.2)
mtext("絶対等級",2,line=2.2)
abline(h=seq(-15,15,5),col="white",lwd=0.5)
abline(v=c(10^3*(2:9),seq(1e4,2e4+2e3,1e3)),col="white",lwd=0.5)
box(bty="l",lwd=2)
points(x=hip$T,y=hip$M,pch=21,cex=cex,
    col=rgb(hip$R,hip$G,hip$B,alpha=1),bg=rgb(hip$R,hip$G,hip$B,alpha=0.8))
# 太陽の表面温度（5800K），Ms；太陽の絶対等級（4.84）
points(x=5800,y=4.84,col="black",cex=0.3)
text(x=5800,y=4.84,labels="太陽",col="black",cex=0.8,pos=3)
# プロットする中で一番大きいガーネット・スターにラベルを付けた。
big<- head(hip[order(-hip$R_Rs),],1)
text(x=big$T,y=big$M,labels=paste0(big$NameJ,"\n太陽の",round(big$R_Rs,0),"倍"),cex=0.8,col="white")
title("ヒッパルコスカタログ(2007)によるH-R図")
# 凡例
legend("bottomright",inset=c(0.01,0),legend=c("100","300","500"),
    pch=21,pt.bg="lightgoldenrod1",col="gray80",pt.cex=c(0.5,1.5,2.5),x.intersp=1.5 ,y.intersp=1,
    box.col ="gray5",bty="n",
    text.col="white",title="R/Rs")
par(bg="white", fg="black", col.main="black", col.lab="black", col.axis="black")
# dev.off()
BACK TO TOP
© 統計ソフトRの備忘録３ 2019
Powered by GitPress.io
Sponsored by Quail.ink