## This is a function file.  The script RUNplanet10.m
##                           shows an example of how to use this function
##
function [ x , v , t , logbook ] = planet10( x , v , t , T , dt , plotduring , velocityfactor ) 
  ## 
  ## x has dimensions [L]
  ## v has dimensions [L/T]
  ## T and dt and t have dimensions [T]
  ## if plotduring > 0,  show the state as we run the simulation, 
  ##                   and use a pause of duration plotduring seconds, every loop
  ## velocityfactor has dimensions [T] - it specifies how to render the velocity
  ##                                  (with a vector of size v*velocityfactor
  ## This function returns the final 
  ## position and velocity and time (x,v,t), 
  ## and also returns a big logbook matrix giving all the [t,x,v]

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

					    ## Prepare the logbook matrix.
					    ## Setting a matrix like this to zero first is recommended, for speed.
					    ## The number of columns in the logbook will 1(for t)+2(for x)+2(for v) = 5.
  logbook=zeros( (T/dt) , 1+2+2 ) ;
  j = 0 ; ## index for log

  if(plotduring>0)
    ## prepare the plot
    ##    guess an appropriate range
    range = sqrt( sum(x.^2) ) * 1.1 ; 
    figure(1); axis([-range range -range range]); gset pointsize 3 ; gset nokey ;
    axis("equal");
    sun = [ 0 0 ] ;
  endif
  
  while( t <= T ) 
    if(plotduring>0)
      ## every loop, plot the state:  
      mytitle=sprintf("t=%5.1f",t);
      title(mytitle);
      if(velocityfactor>0) ## create a little vector to plot
	vel_arrow = [ x ; x + v*velocityfactor ] ; ## [L]+[LT^-1 * T]
      endif
      if(velocityfactor>0) 
	gplot sun with points 3 3, x with points 5 4, vel_arrow with line 2 
      else
	gplot sun with points 3 3, x with points 5 6
      endif
      pause(plotduring);
    endif
    ## keep log  
    j ++ ;
    logbook( j , : ) = [ t , x , v ] ;

    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
endfunction