3

私は、使用している天体物理学モデル(具体的にはsph)からのデータの分析を実行するためにいくつかのコードに取り組んでいます。私は現在、出力の下のコード(10個)のグラフが正しいように見えても、正しい順序ではないという問題を抱えています。30度のグラフを自分自身と呼んでいるものは、50度が出力され、80が40になるというように見えますが、明らかなパターンはありません。これを修正するために提供できる助けをいただければ幸いです。

#Set inclincations of observation
id <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90)
i <- (id * pi) / 180

#Determine sign of velocity vector in direction of observer, then magnitude of said vector (using mask)
for (i in 1:length(i)) {
GD$vecdir <- (GD$Xvel * sin(i) + GD$Zvel * cos(i))
mask <- (GD$vecdir >= 0)
if (sum(mask) > 0) { 
  GD$vecmag[mask] <- sqrt((GD$Xvel^2 * sin(i)^2 ) + (GD$Zvel^2 * cos(i)^2))
}
if (sum(mask) < length(mask)) { 
  GD$vecmag[!mask] <- -(sqrt((GD$Xvel^2 * sin(i)^2 ) + (GD$Zvel^2 * cos(i)^2)))
}  
h2 <- hist(GD$vecmag, breaks=720,  plot=FALSE)
jpeg(paste("VMLOS", (i-1)*10, ".jpeg", sep = ""), width = 1280, height = 720)
print(plot(h2$mids, h2$counts, type="p", pch=4, lty=3, xlab="Velocity Vector Magnitude", ylab="counts", main=paste("VMLOS", (i-1)*10, sep="")))
print(lines(lowess(h2$mids, h2$counts, f=0.05), lwd=4, col="red"))
dev.off()
}

コンテキストでは、iは、原点にあるオブジェクトが観測されているxz平面の傾きであり、ラジアンで、コードに示されている度から変換されます。次に、vecdirセクションは、セットまたは粒子のベクトルの大きさを計算します。原点から離れる方向に移動することは、オブザーバーを基準にして正または負である必要があり、マスクセクションはそのような大きさを計算します。さらに、コードのマスク部分は警告を返します。

1: In GD$vecmag[mask] <- sqrt((GD$Xvel^2 * sin(i)^2) + (GD$Zvel^2 *  ... :
  number of items to replace is not a multiple of replacement length
2: In GD$vecmag[!mask] <- -(sqrt((GD$Xvel^2 * sin(i)^2) +  ... :
  number of items to replace is not a multiple of replacement length

使用ごとに(つまり20の警告)、コードの出力に悪影響を与えるようには見えませんが、苛立たしく、これまでのところ修正できていません。

どんな助けでもありがたいことに受け取られます。

4

1 に答える 1

2

警告については、おそらくこれが必要です (右側の式のインデックスに注意してください)。

GD$vecmag[mask] <- sqrt((GD$Xvel^2 * sin(i)^2 ) + (GD$Zvel^2 * cos(i)^2))[mask]

# ...

GD$vecmag[!mask] <- -(sqrt((GD$Xvel^2 * sin(i)^2 ) + (GD$Zvel^2 * cos(i)^2)))[!mask]

これにより、結果が変わります。

さらに、正しい値を繰り返し処理していません。1:length(i) は i の値ではありません。これを変える:

i <- (id * pi) / 180

#Determine sign of velocity vector in direction of observer, then magnitude of said vector (using mask)
for (i in 1:length(i)) {

これに:

angles <- (id * pi) / 180

#Determine sign of velocity vector in direction of observer, then magnitude of said vector (using mask)
for (i in angles) {

元のコードは、ベクトル内の値ではなく、sin(1)、sin(2)、... を使用しています。これにより、結果もシャッフルされます (間違った角度が使用されていることは明らかではありません)。

于 2012-12-09T15:13:34.570 に答える