% L3_1_modeling.m - Lecture 1.3 MATLAB calculations
% RMM, 6 Oct 03

%
% Spring mass system
%

% Spring mass system parameters
m = 250; m1=m; m2=m;			% masses (all equal)
k = 50; k1=k; k2=k; k3=k;		% spring constants
b = 10;					% damping
A = 0.00315; omega = 0.75;		% forcing function

% Call ode45 routine (MATLAB 6 format; help ode45 for details)
tspan=[0 500];				% time range for simulation
y0 = [0; 0; 0; 0];			% initial conditions
[t,y] = ode45(@springmass, tspan, y0, [], k1, k2, k3, m1, m2, b, A, omega); 

% Plot the input and outputs over entire period
figure(1); plot(t, A*cos(omega*t), t, y(:,1), t, y(:,2));
print -dmeta smtransient.wmf

% Now plot the data for the final 10% (assuming this is long enough...)
endlen = round(length(t)/10);		% last 10% of data record
range = [length(t)-endlen:length(t)]';	% create vector of indices (note ')
tend = t(range);
figure(2); plot(tend, A*cos(omega*tend), tend, y(range,1), tend, y(range,2));
print -dmeta smsteady.wmf

% Compute the relative phase and amplitude of the signals
%
% We make use of the fact that we have a sinusoid in steady state,
% as well as its derivative.  This allows us to compute the magnitude
% of the sinusoid using simple trigonometry ( sin^2 + cos^2 = 1).

u = A*cos(omega*tend); udot = -A*omega*sin(omega*tend);
ampu = mean( sqrt((u .* u) + (udot/omega .* udot/omega)) );
fprintf(1, 'Amplitude = %0.5e cm', ampu*100);

%
% Predator prey system
%

% Set up the initial state
R(1) = 20; F(1) = 35;

% For simplicity, keep track of the year as well
year(1) = 1845;

% Set up parameters
br = 0.7; df = 0.5; a = 0.007;
nperiods = 208;
duration=90;

% Iterate the model
for k = 1:duration*nperiods
  b = br;					% constant food supply
% b = br*(1+0.5*sin(2*pi*k/(4*nperiods)));	% varying food supply (try it!)
  R(k+1) = R(k) + (b*R(k) - a*F(k)*R(k))/nperiods;
  F(k+1) = F(k) + (a*F(k)*R(k) - df*F(k))/nperiods;
  year(k+1) = year(k) + 1/nperiods;
end;

% Plot the populations of rabbits and foxes versus time
figure(3); plot(year, R, year, F);

