clear all
clearplot
## This script file uses the function planet10 several times and makes combined plots

vv = 1.3 ;
trange  = [ -0.3 0.0 0.3 0.6 1.0 ]  ; ##  # angle for rotating the velocity, dimensionless [1]
sun=[0 0];
j = 0 ; 
for theta = [trange] 

M = [ cos(theta) sin(theta) ;
     -sin(theta) cos(theta) ] ;  # rotation matrix, dimensionless
  j ++; 

  ###############################################
				# Duration of simulation, and time-step; and initial time 
  T = 700 ;  dt = 0.1 ; t = 0 ; 
  ################################################
				# Initial condition for position and velocity
  ################################################
  x = [ 93 , 0 ] ;
				# x has dimensions [L]
  
				# choice of initial condition
  v = [0,vv] * M ;
				# v has dimensions [L/T]
  
  ################################################
				# setting this to a value greater than zero 
				#   causes plots to occur, and defines the duration of pause while plotting (in seconds)
  plotduring=-1; ##0.015;
  velocityfactor=10; 
  ## velocityfactor has dimensions [T] - it specifies how to render the velocity
  ##                                  (with a vector of size v*velocityfactor
  
  
  [ x , v , t , logbook ] = planet10( x , v , t , T , dt , plotduring , velocityfactor ); 

  if(j==1)
    xlog(:,1) = logbook(:,1) ;
    command=sprintf("gplot sun u 1:2 w p 3 4, xlog u %d:%d w l %d", 2*j,2*j+1, j) ; 
  else
    command=sprintf("%s, xlog u %d:%d w l %d", command , 2*j,2*j+1, j) ; 
  endif
  xlog(:,2*j:2*j+1) = logbook(:,2:3) ;
  figure(20); axis('equal'); gset nokey; gset pointsize 2;
  eval(command);
  ans=input("ready");

endfor