%EA3
% 
%%---------------------------------------------------------
%%Example: 32Inductors-32Capacitors and a voltage input source
%%---------------------------------------------------------

clear all;
close all;

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

%%
%%System parameters
L       =       1;
C       =       1;
num_elements = 64;
%num_elements = 40;

for kp = 1:2:64 %if properties vary put them in individually in Property matrix
Property(kp)=L;
Property(kp+1)=C;
end

%demo to show what happens if one of the inductors is v.large
%Property(41)=1000*L;   %uncomment this line if you want to see what happens if a middle inductance is large

%State Variables:
%       X(1,:) = current in inductor 1
%       X(2,:) = voltage across capacitor 1 (equal to potential wrt ground at node 1 etc)
%       X(3,:) = current in inductor 2
% et cetera

% Initial conditions quiescent.
i      = 1;                     %time index
X(:,1) = zeros(64,1);

%time parameters
start_time = 0;
dt         = 0.5;               %time step chosen 
final_time = 75;                %end calculation after seventy-five seconds
time(i) =  start_time;          %let us keep a vector of times for plotting

                        % State matrix for use in equation Xdot = AX+F
A = zeros(64,64);       % fill with zeros and then put in the non-zero values 
A(1,2)  =  -1/Property(1);              %handle first row separately
for k=2:1:63                            %do all the interior rows
        A(k,k-1) = 1/Property(k);
        A(k,k+1) =  -1/Property(k);
end
A(64,63)        = 1/Property(64);       %handle the last row separately


% let us make sure that the A matrix looks right
disp('The first part of the first order state matrix is:'), A(1:6,1:6)
disp('The last part of the first order state matrix is:'), A(60:64,60:64)
disp('Somewhere inside A looks like:'), A(28:34,28:34)
pause
%%%
%Voltage source
F=zeros(64,1);          %fill in zeros and then put in the only non-zero term
pulse_magnitude = 1;    %magnitude of voltage pulse at input

% Integrate equations of motion using the Runge-Kutta Method

for t = start_time:dt:final_time

        %calculate voltage source
         F(1) = 0;             % no voltage source except during pulse duration

        if ( (0<t) & (t<4) )			%single pulse
                F(1) = pulse_magnitude/Property(1);             %  source value for pulses
         end
		 
		%other kinds of voltage sources
		%         if ( ((0<t) & (t<4)) | ((6<t) & (t<10)))			%double pulse
		%        if ( (t>4) )			%a constant voltage after 4seconds		
		%			end 
       
        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 voltage variation with time at all nodes as a movie 
for ktime=1:i-1
        subplot(2,1,1), plot(1:32,X(2:2:64,ktime)),Title('Voltage at nodes')            %pick off only the voltages
        axis([0,32,-2,2])
        text(20,1.6,'Time = ')
        text(24,1.6,num2str(ktime*dt))  %display time
        subplot(2,1,2), plot(1:32,X(1:2:64,ktime)),Title('Current in inductors')                        %pick off only the currents
        axis([0,32,-2,2]),xlabel('Position')
        pause(0.1)
end
pause
%do a three-d plot of only the voltages at the nodes
figure

Z=X(2:2:num_elements,1:1:i-1);  %extract voltages out of the X matrix
surf([start_time:dt:final_time],[1:num_elements/2],Z)
%% set the angle for viewing the graph
view(-45,65)
xlabel('time (seconds)')
ylabel('node #')
zlabel('volts')
title('Voltage at each node as a function of time')