>"Quite precise data for thousands of NEOs is available, but its needs converted to position and movement vectors to use in this simulation." https://github.com/strainer/fancy

Thanks for the good discussion. I'm still unconvinced that asking for climate models to report on conservation of energy is "undue speculation", but I really don't know how big a deal it is. Anyway, I looked at my old code and here is some R to get the JPL data and convert to state vectors, hopefully you can use/translate it to easier test your sim:

  dlHorizons <- function(ObjName = "Apophis", 
                         sd = "2000-Jan-01", st = "00:00:00",
                         ed = "2000-Jan-02", et = "00:00:00", 
                         stepSize = "1", stepUnit = "d"){
  
    root     = "http://ssd.jpl.nasa.gov/horizons_batch.cgi?batch=1&"
    start    = paste0("START_TIME%20=%20%27", sd, "%20", st, "%27&")
    end      = paste0("STOP_TIME%20=%20%27", ed, "%20", et, "%27&")
    tble     = "TABLE_TYPE%20=%20%27Vector%27&"
    RefPlne  = "REF_PLANE%20=%20%27Ecliptic%27&"
    RefCntr  = "CENTER%20=%20%27@000%27&"
    VecTab   = "VECT_TABLE=%20%272%27&"
    step     = paste0("STEP_SIZE= %27", stepSize, "%20", stepUnit, "%27&")
    Cmd      = paste0("COMMAND=%20%27",ObjName, "%27")
  
    res = scan(paste0(root, start, end, tble, RefPlne, RefCntr, VecTab, step, Cmd), what="")
    return(res)
  }
  
  
  extractEphemeris <- function(input){
    ends  = grep("(\\$\\$SOE)|(\\$\\$EOE)", input)
    input = input[ends[1]:ends[2]]
  
    dates = input[grep("(TDB)", input) - 2]
    times = input[grep("(TDB)", input) - 1]
  
    start = grep("(TDB)", input) + 1
    end   = start + 5
    ret   = as.data.frame(matrix(nrow = length(start), ncol = 8), stringsAsFactors=F)
    colnames(ret) = c("Date","Time", "px", "py", "pz", "vx", "vy", "vz")

    for(i in 1:length(start)){
      ret[i,3:8] = as.numeric(input[start[i]:end[i]])*1000
    }
  
    ret$Date = dates
    ret$Time = times
    return(ret)
  }