## To run this script, type 
##      planet0
##
## This simple planet demonstration
## simulates the motion of the planet and plots 
## it as it goes, WITHIN THE LOOP.
## Every equation has its dimensions shown. This script 
## has a fixed initial condition. See planet10 for a more 
## general approach.
##
## In planet0a, the "plot" command is used.
## In planet0,  the "gplot" command is used.
##    For this example, "gplot" works much faster, 
##    and the command is simpler too.
## In planet0C, gplot is used, with a "pause" command to slow it down

# Gravitational parameter of sun 
A = 238 ; #  [L^3 T^{-2}] 

################################################
# Duration of simulation, and time-step; and initial time 
T = 1000 ;  dt = 1 ; t = 0 ; 
# T and dt and t have dimensions [T]
################################################
# Initial condition for position and velocity
################################################
 x = [ 93 , 0 ] ;
# x has dimensions [L]
 v = [0,1.1] ;
# v has dimensions [L/T]
################################################

more off
## prepare the plot
figure(1); axis([-100 100 -100 100]); gset pointsize 3 ; gset nokey ;
axis("equal"); ## make lengths in both directions look the same
sun = [ 0 0 ] ;
frametime = 0.02 ; ## number of seconds between plots. (1/50) seconds is
## a good value to give 50 frames per second
oldtime = time ;
timefornextplot = time + frametime ;
while( t < T ) 
  ## every loop, plot the state: (but wait such that the time between
  ## plots is frametime
  pause( timefornextplot - time )
  gplot sun with points 3 3, x with points 5 4         
  timefornextplot = timefornextplot  + frametime ;
  
  t = t + dt ;
  ##  [T] +=  [T] 
  x  = x + v * dt ;
  ##  [L] += [L/T] * [T] 
  
  v = v - A * x / (x*x')**(3.0/2.0) * dt ;
  ## [L/T] -= ( [L^3 T^{-2}] * L / L^3 )  			      
  
endwhile