% Aph 161 - PSet 5, Spring 2009
% ODE Solver Tutorial - Harmonic Oscillator
% Note: Comments for code must have a percent sign before them
% Note2: Matlab is VERY case sensitive, so BE VERY CAREFUL

% This is the main script for the program. For anyone who hasn't had a lot
% of programming experience, the "main" function is where you often put
% most of your preliminary stuff (in Matlab, that consists of setting up
% certain conditions, clearing all old variables, setting up parameters).

% First, clear all variables and close all figure windows (just good form)
 clear all;
 close all;
% Now, we will initialize "options" which will be used to setup the
% tolerances for the ODE solver program (these tolerances just let the
% program (basically) know how accurate you would like your solutions)
 options = odeset('RelTol',1e-4,'AbsTol',1e-4);

% Now, we will invoke a function that will allow our graphs to have larger
% text on the axes and for all of the axes labels (title, xlabels, ylabels,
% etc.).  This is "unecessary," but it will make the graphs easier to read.
% Note that you should leave the 0 and the 'defaultaxesfontsize' as is, but
% if you would like to change the font size (make it smaller or larger),
% just change the 20 to higher or lower values.
 set(0,'defaultaxesfontsize',20); % This is just a "bell and whistle"
                                  % and is uneccessary, but nice
 
% Now, let us setup the necessary parameters for the harmonic oscillator
 w_0  = 2; % Feel free to play around with this parameter if you'd like
           % This parameter represents -sqrt(k/m) in the HO equation
 x_0  = [1;.5]; % These are initial conditions for position and velocity
 tpts = [0;30]; % This is the initial and final time point
    % The ode solver will begin calculations for timepoint 0 and end at
    % timepoint 30
    
% Now, we will call the ode45 function, which will numerically solve our
% set of ODEs. Note that the function must be called with the parameters in
% the following order: ode45(@DE_FunctionName,Initial_FinalTimepoint,
%     InitialConditionsVector,options,ExtraParameters)
 [T,X] = ode45(@HO_ODE,tpts,x_0,options,w_0);
    % Note: This will place all of the solutions for the DEs in the vector
    % X and all of the time points used in the numerical solution in the
    % vector T

% Now, we will plot our results, labeling the graph and putting a legend
% (just for fun). Note that typing 'figure' before plot will ensure that a
% new window is opened when you are graphing (so that you do not
% accidentally record over a previous graph that you have plotted)
% The plot function must be called with the parameters in the following
% order: plot(XaxisVariableVector,YaxisVariableVector,OptionalPlotOptions)
% Note that you can plot multiple things on the same plot if you'd like:
% plot(XaxisVar1,YaxisVar1,XaxisVar2,YaxisVar2,...,PlotOptions)
 figure, plot(T,X(:,1),':r',T,X(:,2),'-.g','LineWidth',1.5)
   % X(:,1) refers to the entire first column of the X vector, which
   % contains the values for position in time from the solution to the ODE
   % ':r' will make the first curve in the graph a RED DOTTED LINE (Default
   % would be a blue solid line)
   % X(:,2) refers to the entire first column of the vector X, which
   % contains the values for velocity in time from the solution to the ODE
   % '-.g' will make a GREEN DASH-DOT LINE
   % 'LineWidth',1.5 will increase the line width of the curves so that
   % they are easier to see when printing your solution
  axis([0 30 -5 5])
   % axis(['Minimum X Val' 'Maximum X Val' 'Min Y Val' 'Max Y Val'])
  xlabel('Time (dimensionless)') % X axis label
  ylabel('Velocity and Position profiles (dimensionless)') % Y axis label
  Legend('Position','Velocity','Location','SouthEast')
   % 'Position' is the legend label for the position curve
   % 'Velocity' is the legend label for the velocity curve
   % 'Location','SouthEast' puts the legend at the bottom right corner of
   % the graph