! Planetary motion with the Euler-Cromer method
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program kepler
   option nolet   ! don't need the graphics libraries for this program
   call initialize(x,v_x,y,v_y,dt)   ! set up initial conditions
   call calculate(x,v_x,y,v_y,dt)    ! do the calculation
end
! initialize variables
sub initialize(x,v_x,y,v_y,dt)
   input prompt "initial x position -> ": x
   input prompt "initial y position -> ": y
   input prompt "initial x velocity -> ": v_x
   input prompt "initial y velocity -> ": v_y
   input prompt "time step -> ": dt
                          ! now set up window for plotting
   aspect = 1.33           ! aspect ratio of screen - your's may be different
   r = 1.4 * sqr(x^2+y^2) ! pick a scale a bit larger than the planet's radius
   set window -r,r,-r/aspect,r/aspect  ! set window coordinates
   set color "black"
   clear
   plot -r,0;r,0                       ! plot axes
   plot 0,-r*aspect;0,r*aspect
end sub
! x,y = position of planet
! v_x,v_y = velocity of planet
! dt = time step
sub calculate(x,v_x,y,v_y,dt)
   do                          ! use Euler-Cromer method
      r = sqr(x^2 + y^2)
      v_x = v_x - (4 * pi^2 * x * dt) / r^3
      v_y = v_y - (4 * pi^2 * y * dt) / r^3
      set color "black"        ! keep trail of the planet black
      plot x,y
      x = x + v_x * dt
      y = y + v_y * dt
      set color "red"          ! current location of the planet is red
      plot x,y
   loop until key input        ! loop until any key is hit
end sub