%EA3
% 
%%---------------------------------------------------------
%%Example: Spring-mass-damper system 
%%---------------------------------------------------------


clear all;
close all;

global A;			%share these with function rate_fn which computes Xdot=AX

%%
%%System parameters
K1	=	20;				%element 1 - spring
%C2  =   2;              %element 2 - damper
C2=input('Enter damping coefficient:');
m3	=	10;				%element 3 - mass
wn=sqrt(K1/m3);				%natural frequency of oscillation of undamped system
period=2*pi/wn;				%period of oscillation

%State Variables:
%	X(1,:) = stretch of spring element 1
%	X(2,:) = velocity of mass element 3

% Initial conditions .  
i	=	1;					%time index
X(1,1) = 0.1;   %stretch of spring
X(2,1) = 0;

start_time = 0;
dt	=	0.05*period;			%time step chosen as a fraction of one period
final_time	=	20*period;		%end calculation after ten periods
time(i)	=	start_time;			%let us keep a vector of times for plotting

%State matrix
A = [0,		1;
	-K1/m3,	-C2/m3];  % for use in equation Xdot = AX+F

%

% Integrate equations of motion using the Runge-Kutta Method

for t = start_time:dt:final_time

        k1 = rate_fn(X(:,i));			% current rate
        k2 = rate_fn(X(:,i) + dt*k1/2);	% estimated mid-point rate
        k3 = rate_fn(X(:,i) + dt*k2/2);	% even better mid-point rate
        k4 = rate_fn(X(:,i) + dt*k3);	% excellent end-point rate

        X(:,i+1) = X(:,i) + dt*(k1/6 + k2/3 + k3/3 + k4/6);	% use weighted average of rates		
        time(i+1) = time(i) + dt;
		i = i+1;
end

%%Plot the stretch of spring 1
comet(time(:),X(1,:));
pause
plot(time(:),X(1,:)), title('Stretch of Spring 1'), grid
%subplot(2,1,2), plot(time(:),an_X(:)), title('Analytical'),xlabel('Time'), grid