简洁版本
正如评论中所建议的，我将在此处留下一个总结的简短版本的问题。原始的完整解释可以在下面找到。
总之，我目前正在使用 deSolve 包，并且有兴趣实现响应根函数的事件。我知道我可以使用此类事件来引入模型状态变量的突然变化，但我也想修改模型的参数值以响应此类根函数。
这可能吗？
长版
我已经在 R 中实现了一个轨道数值传播器（一个在给定初始位置和速度状态的情况下计算空间卫星位置的函数）。这个问题被表述为一组 6 个 ODE（位置和速度的 X、Y 和 Z 分量）。我实现的模型计算任何给定时间的加速度，然后我使用 deSolve 包执行积分并计算轨迹。
进行此类计算时必须确定的一个关键参数是参考系的中心，该参考系的中心通常放置在对卫星施加最显着引力影响的天体的质心处。这是因为，尽管原则上可以使用任意参考系进行积分和计算轨迹，但实际上我们只有将坐标中心放在施加主要引力影响的天体上才能获得合理的结果（即，地球用于地球轨道卫星，月球用于月球轨道卫星，等等），如本 SE 问题中所述。
最初，我的实现使用了一个恒定的坐标中心，或者由用户提供，或者根据不同主要天体的影响范围自动确定。
然而，这不适用于模拟行星际空间任务，因为施加主要引力影响的天体在轨迹期间会发生变化。一个很好的例子是阿波罗任务，卫星从地球轨道开始，然后移动到月球轨道。
我已经设法检测到中央天体的这种变化何时发生，并将其作为积分器结果的一部分返回。然而，为了实际执行正确的建模，当检测到这些变化时，需要更改集成期间使用的中心体。这个“改变中心体”的过程涉及两个任务（注意只是坐标中心的移动，不涉及旋转）：
- 从卫星坐标中减去要用作新坐标中心的天体坐标（通过这样做，卫星坐标现在称为新天体）。
- 修改指定传递给计算加速度的函数的中央天体的参数的值，这是提供给定义 ODE 模型的函数的参数列表的元素之一。
我相信任务 1 可以通过使用根激活事件轻松解决。为此，我在专门为此目的创建的环境中定义了一个变量，该变量存储自动计算的天体的值，该天体在积分器的每次迭代中施加主要的引力影响。在新的迭代中，计算一个新值，并与之前的值进行比较。如果相同，则不会发生任何事情。但如果不同，根函数将返回 0，从而触发事件函数。然后，事件函数将返回减去新中央天体坐标的位置。
但是，我不确定如何执行任务 2，因为这将涉及更改提供给 ODE 模型的初始参数之一。对此的任何想法将不胜感激！（要么是我的方法的延续，要么是完全不同的方法）。
我留下了所涉及代码的简化版本。
我的 main 函数被称为hpop，并且是传递初始状态向量和其他参数的用户级函数。它看起来像这样：
hpop <- function(position, velocity, dateTime, times, satelliteMass, dragArea,
radiationArea, dragCoefficient, radiationCoefficient,
earthSphericalHarmonicsDegree = 130, solidEarthTides=TRUE,
oceanTides=TRUE, moonSphericalHarmonicsDegree = 150,
centralBody="Earth", ...) {
extraArgs <- list(...)
## This is the environment used to hold the variable that keeps track of what
## the central body should be at each iteration
propagationGlobalVars <- new.env()
## This is the initial value of such variable, which is the user-provided central body, by default Earth
propagationGlobalVars$latestCentralBody <- centralBody
## This is the initial state, composed by 6 variables (X, Y and Z components of position and velocity)
initial_state <- c(position, velocity)
## This is the list of parameters required for trajectory calculation
## There is quite a few, but the 2 most relevant for this question are the last 2
## centralBody is the body that will be used as the center of coordinates, and
## globalVarsEnv is the environment that will be containing the variable that keeps track of what the central body should be
parameters = list(
dateTime = dateTime,
solarArea = radiationArea,
satelliteMass = satelliteMass,
satelliteArea = dragArea,
Cr = radiationCoefficient,
Cd = dragCoefficient,
earthSPHDegree = earthSphericalHarmonicsDegree,
moonSPHDegree = moonSphericalHarmonicsDegree,
SETcorrections = solidEarthTides,
OTcorrections = oceanTides,
centralBody = centralBody,
globalVarsEnv = propagationGlobalVars)
## This calles function ode from the deSolve package, passing the previously defined initial state,
## integration times and the function defining the ode model (code below)
integration_results <- ode(y=initial_state, times=times, func=odeModel,
parms=parameters, method="radau", rtol=1e-13,
atol=1e-16, hini=0.01, ...)
numeric_results <- integration_results[, 1:7]
central_bodies <- names(centralBodiesNum[integration_results[, 8]])
output <- cbind(as.data.frame(numeric_results), central_bodies)
colnames(output) <- c("time", "X", "Y", "Z", "dX", "dY", "dZ", "Central body")
return(output)
}
odeModelODE模型func的代码的简化版本ode是：
odeModel <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
state_vector <- state
## Function accel calculates the acceleration and velocity at time t. It returns a list with
## two elements. The first is a numeric vector with the X, Y and Z components of velocity,
## and the X, Y and Z components of acceleration, in this order.
## The second element is the celestial body that, given the position at that iteration,
## exerts the main gravitational influence on the satellite
results <- accel(t, state_vector, dateTime, solarArea, satelliteMass,
satelliteArea, Cr, Cd, earthSPHDegree, SETcorrections,
OTcorrections, moonSPHDegree, centralBody)
centralBody <- results[[2]]
## Now we can compare the central body with that of the previous iteration,
## and assign the new value to the tracking variable
if(centralBody != globalVarsEnv$latestCentralBody) {
message(strwrap(paste("A transition from the sphere of influence of ",
globalVarsEnv$latestCentralBody, " to that of ",
centralBody, " has been detected.", sep=""), initial="", prefix="\n"))
assign("latestCentralBody", centralBody, envir = globalVarsEnv)
## here I would also trigger the root function, to perform Task 1 of the two
## tasks required to change the central body as described above. And also,
## I would need a way to change the value of centralBody in the parameters argument, which is the main issue of this question
}
acceleration <- results[[1]]
dx <- acceleration[1, 1]
dy <- acceleration[1, 2]
dz <- acceleration[1, 3]
d2x <- acceleration[2, 1]
d2y <- acceleration[2, 2]
d2z <- acceleration[2, 3]
## This is the return value of the odeModel function. As specified in the
## documentation for the func argument of the ode function, the first element
## of the list are the values of the derivatives of the model, and the
## second one can be any other variable. Note that since such other variables
## must be numeric, I actually access a named vector of numbers, and then convert
## back to the proper names when outputting final results to the users.
## It is important to provide the central body of each output step so that
## we know what the center of coordinates are at each step
list(c(dx, dy, dz, d2x, d2y, d2z),
centralBodiesNum[centralBody])
})
}