3

我一直在编写一些代码来对我一直在使用的天体物理学模型（特别是 sph）的数据进行一些分析。我目前遇到一个问题，即下面的代码输出图（其中十个）虽然它们看起来是正确的，但顺序不正确，即。一个自称为 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 回答