%EA3
% 
%%---------------------------------------------------------
%%Example: Spring-mass-system with a harmonic force source
%%---------------------------------------------------------


clear all;
close all;

global A F;			%share these with function frate_fn which computes Xdot=AX+F

%%
%%System parameters
K1	=	40;				%element 1 - spring
m2	=	10;				%element 2 - mass
wn=sqrt(K1/m2);				%natural frequency of oscillation of unforced system
period=2*pi/wn;				%period of oscillation

%%Forcing function
wffrac =input('Enter forcing function frequency as multiple of wn:');
wf=wffrac*wn;
%wf	=	1.95;			%frequency of force source function F3(t)=F0*sin(wf*t)
F0	=	1;				%force source amplitude
%State Variables:
%	X(1,:) = stretch of spring element 1
%	X(2,:) = velocity of mass element 2

% Initial conditions quiescent.  Motion only due to force source
i	=	1;					%time index
X(1,1) = 0;
X(2,1) = 0;

%analytical solution
an_fac = (F0/m2)*(1/(wn^2-wf^2));		%analytical solution for spring stretch
an_X(1) = 0;

start_time = 0;
dt	=	0.05*period;			%time step chosen as a fraction of one period
final_time	=	60*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/m2,	0];  % for use in equation Xdot = AX+F

%

% Integrate equations of motion using the Runge-Kutta Method

for t = start_time:dt:final_time

		F = [0,	F0*sin(wf*t)/m2]';				% current force source value
        k1 = frate_fn(X(:,i));			% current rate
        k2 = frate_fn(X(:,i) + dt*k1/2);	% estimated mid-point rate
        k3 = frate_fn(X(:,i) + dt*k2/2);	% even better mid-point rate
        k4 = frate_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;
		%analytical solution
		an_X(i+1) = an_fac*(sin(wf*time(i+1))-(wf/wn)*sin(wn*time(i+1)));

		i = i+1;
end

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