%EA3
% 
%%---------------------------------------------------------
%%Example: LC-oscillator with a voltage source
%%---------------------------------------------------------


clear all;
close all;

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

%%
%%System parameters
L	=	1;				%element - inductor
C	=	1;				%element - capacitor

% Voltage source
V0	=	9;				% source amplitude
%State Variables:
%	X(1,:) = voltage across capacitor element
%	X(2,:) = current through inductor element

% Initial conditions quiescent.
i	=	1;					%time index
X(1,1) = 0;
X(2,1) = 0;
%time parameters

ws=1/sqrt(L*C);				%natural frequency of oscillation of unforced system
period=2*pi/ws;				%period of oscillation


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

%State matrix
A = [0,		1/C;
	-1/L,	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,	V0/L]';				%  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;

		i = i+1;
end

%%Plot
comet(time(:),X(1,:));
pause
subplot(2,1,1), plot(time(:),X(1,:)), title('Voltage across Capacitor'), grid
subplot(2,1,2), plot(time(:),X(2,:)), title('Current in Circuit'), xlabel('Time'),grid