Physics starter - Kepler

MODEL "Physics starter - Kepler"
DESCRIPTION "One eccentric orbit with GM=1 and semi-major axis=1; check energy, angular momentum and return position."
VERSION "1.0.0"
AUTHOR "Grid Team"
TAGS "physics", "starter", "kepler"

# GM=1 and semi major axis=1 give period 2*pi and specific energy -1/2.
input Eccentricity = 0.6
  VALIDATE BETWEEN 0 AND 0.8 MESSAGE "Use eccentricity from 0 to 0.8 for this finite-time orbit example."
Period = 2*PI()
InitialState = [1-Eccentricity,0,0,SQRT((1+Eccentricity)/(1-Eccentricity))]
output Orbit = ODE_SOLVE((time,state) => [INDEX(state,1,3),INDEX(state,1,4),-INDEX(state,1,1)/(INDEX(state,1,1)^2+INDEX(state,1,2)^2)^1.5,-INDEX(state,1,2)/(INDEX(state,1,1)^2+INDEX(state,1,2)^2)^1.5],InitialState,0,Period,{samples:401,rtol:1e-10,atol:1e-12})
output ReturnError = SQRT((INDEX(Orbit,401,2)-(1-Eccentricity))^2+INDEX(Orbit,401,3)^2)
output EnergyError = MAX(MAP(SEQUENCE(401),(row) => ABS((INDEX(Orbit,row,4)^2+INDEX(Orbit,row,5)^2)/2-1/SQRT(INDEX(Orbit,row,2)^2+INDEX(Orbit,row,3)^2)+0.5)))
output AngularMomentumError = MAX(MAP(SEQUENCE(401),(row) => ABS(INDEX(Orbit,row,2)*INDEX(Orbit,row,5)-INDEX(Orbit,row,3)*INDEX(Orbit,row,4)-SQRT(1-Eccentricity^2))))
A1 = "Orbit closure error"
B1 = ReturnError
A2 = "Maximum absolute specific energy error"
B2 = EnergyError
A3 = "Maximum absolute angular momentum error"
B3 = AngularMomentumError

# Scientific checks and interpretation

output Assumptions = "One eccentric orbit with GM=1 and semi-major axis=1; check energy, angular momentum and return position. Solver completion alone does not establish global accuracy or physical validity."
output ScenarioGuide = "Save a Baseline input scenario, change an admitted input, inspect numerical checks and the unsmoothed plot, then restore Baseline. Save and reopen a .gridoc archive. Plot selection also accepts the Orbit numeric array; choose its X and Y columns in the chart dialog."

Charts_1!type = "charts"
Charts_1!config = <toml>
  version = 1
  kind = "charts"
  title = "Kepler"
  [layout]
  x_axis = "x / semi-major axis"
  y_axis = "y / semi-major axis"
  [[charts]]
  type = "line"
  smooth = false
  bind = "Orbit"
  x = "1"
  y = "2"
  label = "Orbit"
  color = "#2563eb"
</toml>

END MODEL