我一直在编写一些代码,以对我一直在使用的天体物理模型(sph specally)的数据进行一些分析。我目前遇到了一个问题,下面代码输出的图形(其中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个警告),这似乎不会对代码的输出产生负面影响,但令人恼火,到目前为止我还无法修复它。
任何帮助我们都将不胜感激。
对于警告,您可能需要以下内容(注意右侧表达式上的索引):
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:长度(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)。。。,而不是矢量中的值。这也会打乱结果(使用错误的角度可能并不明显)。