#!/usr/bin/octave -qf

function trajectory

 global g

 g = 9.81

%
% Enter initial position and speed (two dimensional)
%
 rx0 = input('Enter initial r_x: ');
 ry0 = input('Enter initial r_y: ');
 v0 = input('Enter initial v_0: ');
 theta0 = input('Enter initial theta_0: ');
 tf = input('Enter time of integration tf: ');
 pcolor = input('Enter the desired plot color 0-6: ');
 vx0 = v0*cos(theta0);
 vy0 = v0*sin(theta0);

%
% Create a vector of times for integration, 
% in reasonable timesteps
%
 t = linspace(0,tf,10000);

%
% initial condition.  Note the form of the four coordinates
% for position and velocity in two dimensions.
%
 coords0 = [ rx0;ry0;vx0;vy0 ];

%
% Here we use a "real" differential equation solver.
% For future assignments, we'll use a simple Euler rule
% to write out own (less accurate) one, because it will
% effectively lead us to Hamilton's equations of motion.
%
% Set ode solver tolerance
 lsode_options("relative tolerance",1e-6);

%
% Solve differential equation.  "eom" is the function
% for evaluating derivatives, a.k.a. the "equation of
% motion" and is defined below.
%
 coords = lsode("eom",coords0,t);

%
% That's it.  coords now contains four vectors, each of
% the same length as the vector t defined above.  They are
% x(t), y(t), vx(t) vy(t) BUT they are NOT INDEXED that way.
% Instead it works like:
%
%   x(t[i]) = x[i]
%
% (etc) for a common index i that runs from 0 to 9999.  So
% we get them out to make it easier to define plots.
%
 x = coords(:,1);
 y = coords(:,2);
 vx = coords(:,3);
 vy = coords(:,4);

% Now we make plots of these trajectories.  Nothin' to it.
% Note that "hold on" below ensures that you can put multiple curves
% on this plot to see the effect of e.g. varying angle.
 figure(1);
 axis([0, 100, 0, 100],"normal");
 title('Trajectories with no drag forces');
 xlabel('x axis');
 ylabel('y axis');
 plot(x,y,pcolor);

endfunction

%
% Equations of motion.  We have to write them as coupled
% FIRST order differential equations.  For each dimension,
% this is ALWAYS going to be:
%
%   dr_x/dt = v_x
%   dv_x/dt = a_x = F_x/m
%
% c is the "coords" vector passed to this function by
% the ode solver.  This function needs to return a vector of
% the same form of the derivatives of the coordinates at this
% position/time.  Note that for this problem the eom are very
% simple.
%
function dydt = eom(c,t)

 global g;

 ax = 0.0;
 ay = -g;

 dydt = [c(3);c(4);ax;ay];

endfunction

%
% This is the actual "program" (the core loop).  Note that
% all it does is run trajectory (which prompts you for
% initial conditions, solves the eom, and plots the result)
% and then ask if you're done and rerun if desired.
%

 printf("==================================================================\n");
 printf(" trajectory\n");
 printf("\n");
 printf(" trajectory plots multiple simple 2d trajectories on a single plot.\n");
 printf(" You will be prompted for initial position (rx0 and ry0) and velocity\n");
 printf(" (v0 and theta0, in polar coordinates).  In addition, you must enter a\n");
 printf(" maximum integration time and a plot color.\n");
 printf("\n");
 printf(" After completing each curve, you will be returned to a prompt and asked\n");
 printf(" if you wish to do another trajectory.  When you are done, exit.\n");

done = "y";
while(done == "y")
  trajectory;
  hold on;
  done = input('Do another graph? (y/n): ',"s");
endwhile
__gnuplot_set__ term postscript;
__gnuplot_set__ output "trajectory.ps";
replot;