Kepler. Interaction between two planets

The Kepler class

library(rODE)

# This code can also be found in the `examples` folder under this name:
#
# Kepler.R
#

setClass("Kepler", slots = c(
    GM = "numeric"
    ),
    contains = c("ODE")
)

setMethod("initialize", "Kepler", function(.Object, ...) {
    .Object@GM <- 1.0                 # gravitation constant times combined mass
    .Object@state <- vector("numeric", 5)  # x, vx, y, vy, t
    return(.Object)
})


setMethod("getState", "Kepler", function(object, ...) {
    # Gets the state variables.
    return(object@state)
})

setMethod("getRate", "Kepler", function(object, state, ...) {
    # Computes the rate using the given state.
    r2 <- state[1] * state[1] + state[3] * state[3]  # distance squared
    r3 <- r2 * sqrt(r2)   # distance cubed
    object@rate[1] <- state[2]
    object@rate[2] <- (- object@GM * state[1]) / r3
    object@rate[3] <- state[4]
    object@rate[4] <- (- object@GM * state[3]) / r3
    object@rate[5] <- 1   # time derivative

    object@rate
})

# constructor
Kepler <- function(r, v) {
    kepler <- new("Kepler")
    kepler@state[1] = r[1]
    kepler@state[2] = v[1]
    kepler@state[3] = r[2]
    kepler@state[4] = v[2]
    kepler@state[5] = 0
    return(kepler)
}

Run the application KeplerApp

# This code can also be found in the `examples` folder under this name:
# KeplerApp.R
#
KeplerApp <- function(verbose = FALSE) {
    
    # set the orbit into a predefined state.
    r <- c(2, 0)
    v <- c(0, 0.25)
    dt <- 0.1
    
    planet <- Kepler(r, v)
    solver <- RK45(planet)
    
    # solver <- step(solver)
    
    while (planet@state[5] <= 10) {
        solver <- step(solver)
        planet <- solver@ode
        if (verbose)
            cat(sprintf("state[1]=%10f, state[2]= %10f,  
                        state[3]=%10f, state[5]=%10f\n",
                    planet@state[1],
                    planet@state[2], planet@state[3], planet@state[5]))
    }
    
# at t=100, dt=0.1,  c(2.131958,     1.105316,   100.000000)
# Java: state[0] = 0.444912, state[1]= -1.436203, state[2]= 0.459081, state[4]= 10.033245
#       currentStep=    0.061646
# R:    state[1] = 0.444912, state[2]= -1.436203, state[3]= 0.459081, state[5]= 10.033245
#       currentStep= 0.06164632
}

KeplerApp(verbose = TRUE)

## state[1]=  1.999987, state[2]=  -0.002500,  
##                         state[3]=  0.002500, state[5]=  0.010000
## state[1]=  1.998487, state[2]=  -0.027511,  
##                         state[3]=  0.027493, state[5]=  0.110000
## state[1]=  1.948511, state[2]=  -0.162156,  
##                         state[3]=  0.158503, state[5]=  0.639561
## state[1]=  1.842313, state[2]=  -0.290237,  
##                         state[3]=  0.270214, state[5]=  1.111057
## state[1]=  1.699885, state[2]=  -0.413368,  
##                         state[3]=  0.359093, state[5]=  1.517417
## state[1]=  1.538370, state[2]=  -0.532916,  
##                         state[3]=  0.425287, state[5]=  1.859899
## state[1]=  1.374619, state[2]=  -0.647215,  
##                         state[3]=  0.470134, state[5]=  2.138174
## state[1]=  1.218724, state[2]=  -0.755967,  
##                         state[3]=  0.497571, state[5]=  2.360880
## state[1]=  1.072153, state[2]=  -0.861865,  
##                         state[3]=  0.512006, state[5]=  2.542429
## state[1]=  0.937302, state[2]=  -0.965097,  
##                         state[3]=  0.516394, state[5]=  2.690303
## state[1]=  0.813892, state[2]=  -1.066667,  
##                         state[3]=  0.513150, state[5]=  2.811963
## state[1]=  0.701810, state[2]=  -1.166754,  
##                         state[3]=  0.504087, state[5]=  2.912462
## state[1]=  0.600201, state[2]=  -1.265693,  
##                         state[3]=  0.490567, state[5]=  2.996103
## state[1]=  0.508216, state[2]=  -1.363525,  
##                         state[3]=  0.473613, state[5]=  3.066145
## state[1]=  0.424811, state[2]=  -1.460284,  
##                         state[3]=  0.453937, state[5]=  3.125270
## state[1]=  0.348893, state[2]=  -1.555914,  
##                         state[3]=  0.431985, state[5]=  3.175648
## state[1]=  0.280902, state[2]=  -1.648019,  
##                         state[3]=  0.408535, state[5]=  3.218115
## state[1]=  0.223580, state[2]=  -1.730067,  
##                         state[3]=  0.385487, state[5]=  3.252065
## state[1]=  0.173098, state[2]=  -1.804540,  
##                         state[3]=  0.362226, state[5]=  3.280635
## state[1]=  0.128954, state[2]=  -1.869423,  
##                         state[3]=  0.339146, state[5]=  3.304667
## state[1]=  0.089757, state[2]=  -1.923913,  
##                         state[3]=  0.316039, state[5]=  3.325330
## state[1]=  0.054703, state[2]=  -1.965983,  
##                         state[3]=  0.292799, state[5]=  3.343347
## state[1]=  0.022867, state[2]=  -1.992815,  
##                         state[3]=  0.269038, state[5]=  3.359421
## state[1]= -0.002734, state[2]=  -1.999878,  
##                         state[3]=  0.247592, state[5]=  3.372236
## state[1]= -0.025046, state[2]=  -1.987905,  
##                         state[3]=  0.226705, state[5]=  3.383416
## state[1]= -0.044253, state[2]=  -1.955638,  
##                         state[3]=  0.206593, state[5]=  3.393146
## state[1]= -0.061107, state[2]=  -1.900866,  
##                         state[3]=  0.186790, state[5]=  3.401876
## state[1]= -0.077340, state[2]=  -1.811154,  
##                         state[3]=  0.165112, state[5]=  3.410606
## state[1]= -0.090121, state[2]=  -1.700256,  
##                         state[3]=  0.145495, state[5]=  3.417872
## state[1]= -0.101946, state[2]=  -1.546659,  
##                         state[3]=  0.124349, state[5]=  3.425139
## state[1]= -0.110487, state[2]=  -1.386681,  
##                         state[3]=  0.106306, state[5]=  3.430950
## state[1]= -0.117991, state[2]=  -1.189309,  
##                         state[3]=  0.087270, state[5]=  3.436762
## state[1]= -0.124234, state[2]=  -0.952724,  
##                         state[3]=  0.067308, state[5]=  3.442573
## state[1]= -0.128992, state[2]=  -0.678836,  
##                         state[3]=  0.046545, state[5]=  3.448384
## state[1]= -0.132065, state[2]=  -0.374557,  
##                         state[3]=  0.025178, state[5]=  3.454196
## state[1]= -0.133309, state[2]=  -0.051968,  
##                         state[3]=  0.003465, state[5]=  3.460007
## state[1]= -0.132663, state[2]=   0.273254,  
##                         state[3]= -0.018297, state[5]=  3.465819
## state[1]= -0.130160, state[2]=   0.584886,  
##                         state[3]= -0.039804, state[5]=  3.471630
## state[1]= -0.125918, state[2]=   0.869447,  
##                         state[3]= -0.060784, state[5]=  3.477442
## state[1]= -0.120124, state[2]=   1.118322,  
##                         state[3]= -0.081018, state[5]=  3.483253
## state[1]= -0.112996, state[2]=   1.328113,  
##                         state[3]= -0.100358, state[5]=  3.489065
## state[1]= -0.104762, state[2]=   1.499625,  
##                         state[3]= -0.118720, state[5]=  3.494876
## state[1]= -0.095634, state[2]=   1.636314,  
##                         state[3]= -0.136076, state[5]=  3.500687
## state[1]= -0.085801, state[2]=   1.742877,  
##                         state[3]= -0.152435, state[5]=  3.506499
## state[1]= -0.075425, state[2]=   1.824251,  
##                         state[3]= -0.167834, state[5]=  3.512310
## state[1]= -0.064638, state[2]=   1.885043,  
##                         state[3]= -0.182323, state[5]=  3.518122
## state[1]= -0.053547, state[2]=   1.929270,  
##                         state[3]= -0.195961, state[5]=  3.523933
## state[1]= -0.042240, state[2]=   1.960296,  
##                         state[3]= -0.208811, state[5]=  3.529745
## state[1]= -0.024762, state[2]=   1.988206,  
##                         state[3]= -0.226986, state[5]=  3.538590
## state[1]= -0.007117, state[2]=   1.999149,  
##                         state[3]= -0.243668, state[5]=  3.547435
## state[1]=  0.010569, state[2]=   1.998339,  
##                         state[3]= -0.259032, state[5]=  3.556280
## state[1]=  0.028211, state[2]=   1.989426,  
##                         state[3]= -0.273232, state[5]=  3.565126
## state[1]=  0.055438, state[2]=   1.965209,  
##                         state[3]= -0.293315, state[5]=  3.578890
## state[1]=  0.082275, state[2]=   1.933609,  
##                         state[3]= -0.311301, state[5]=  3.592653
## state[1]=  0.108648, state[2]=   1.898277,  
##                         state[3]= -0.327516, state[5]=  3.606417
## state[1]=  0.134523, state[2]=   1.861355,  
##                         state[3]= -0.342217, state[5]=  3.620181
## state[1]=  0.174201, state[2]=   1.802907,  
##                         state[3]= -0.362768, state[5]=  3.641839
## state[1]=  0.212629, state[2]=   1.746108,  
##                         state[3]= -0.380695, state[5]=  3.663497
## state[1]=  0.249854, state[2]=   1.692020,  
##                         state[3]= -0.396459, state[5]=  3.685154
## state[1]=  0.302892, state[2]=   1.617570,  
##                         state[3]= -0.416548, state[5]=  3.717216
## state[1]=  0.353649, state[2]=   1.549701,  
##                         state[3]= -0.433484, state[5]=  3.749278
## state[1]=  0.402328, state[2]=   1.487825,  
##                         state[3]= -0.447862, state[5]=  3.781339
## state[1]=  0.475441, state[2]=   1.400570,  
##                         state[3]= -0.466397, state[5]=  3.831999
## state[1]=  0.544414, state[2]=   1.323993,  
##                         state[3]= -0.480848, state[5]=  3.882659
## state[1]=  0.609733, state[2]=   1.256032,  
##                         state[3]= -0.492060, state[5]=  3.933319
## state[1]=  0.705551, state[2]=   1.163272,  
##                         state[3]= -0.504486, state[5]=  4.012611
## state[1]=  0.794548, state[2]=   1.083352,  
##                         state[3]= -0.512007, state[5]=  4.091904
## state[1]=  0.877613, state[2]=   1.013237,  
##                         state[3]= -0.515692, state[5]=  4.171196
## state[1]=  0.997346, state[2]=   0.918281,  
##                         state[3]= -0.515466, state[5]=  4.295343
## state[1]=  1.106172, state[2]=   0.836806,  
##                         state[3]= -0.509573, state[5]=  4.419491
## state[1]=  1.205532, state[2]=   0.765306,  
##                         state[3]= -0.499300, state[5]=  4.543638
## state[1]=  1.349394, state[2]=   0.664692,  
##                         state[3]= -0.475493, state[5]=  4.745231
## state[1]=  1.474449, state[2]=   0.577891,  
##                         state[3]= -0.445016, state[5]=  4.946823
## state[1]=  1.583049, state[2]=   0.500925,  
##                         state[3]= -0.409547, state[5]=  5.148416
## state[1]=  1.722552, state[2]=   0.395386,  
##                         state[3]= -0.347392, state[5]=  5.460505
## state[1]=  1.831013, state[2]=   0.301195,  
##                         state[3]= -0.278926, state[5]=  5.772593
## state[1]=  1.911323, state[2]=   0.214463,  
##                         state[3]= -0.206142, state[5]=  6.084682
## state[1]=  1.981504, state[2]=   0.096526,  
##                         state[3]= -0.095744, state[5]=  6.537692
## state[1]=  1.999406, state[2]=  -0.017214,  
##                         state[3]=  0.017210, state[5]=  6.990702
## state[1]=  1.965799, state[2]=  -0.131685,  
##                         state[3]=  0.129715, state[5]=  7.443712
## state[1]=  1.879222, state[2]=  -0.251998,  
##                         state[3]=  0.238683, state[5]=  7.896721
## state[1]=  1.735613, state[2]=  -0.384820,  
##                         state[3]=  0.340309, state[5]=  8.349731
## state[1]=  1.581971, state[2]=  -0.501704,  
##                         state[3]=  0.409948, state[5]=  8.697452
## state[1]=  1.413517, state[2]=  -0.620264,  
##                         state[3]=  0.461113, state[5]=  8.998642
## state[1]=  1.257756, state[2]=  -0.728500,  
##                         state[3]=  0.491933, state[5]=  9.230141
## state[1]=  1.107293, state[2]=  -0.835986,  
##                         state[3]=  0.509484, state[5]=  9.422887
## state[1]=  0.970161, state[2]=  -0.939291,  
##                         state[3]=  0.516090, state[5]=  9.577646
## state[1]=  0.843554, state[2]=  -1.041514,  
##                         state[3]=  0.514565, state[5]=  9.705677
## state[1]=  0.728870, state[2]=  -1.141798,  
##                         state[3]=  0.506822, state[5]=  9.810873
## state[1]=  0.624641, state[2]=  -1.241086,  
##                         state[3]=  0.494300, state[5]=  9.898458
## state[1]=  0.530379, state[2]=  -1.339153,  
##                         state[3]=  0.478132, state[5]=  9.971599
## state[1]=  0.444912, state[2]=  -1.436203,  
##                         state[3]=  0.459081, state[5]= 10.033245