%% Universal Student
% *Lab1: Circuits*
% January 2011




%%  
% In this lab we model an electric circuit related to
% approximating the behavior of nerve fibers.  The circuit is modeled
% both physically and mathematically.  Discrepancies between the two
% models are accounted for through error analysis.  The voltages at
% the nodes are observed to be linear functions of the input voltage
% to the circuit.  Matlab's publish feature will be used to document
% the write up the results.
%
% The mathematical model used in this lab is described in the course
% notes.  Code from the class, fib1.m, was modified for modeling the
% circuit.  The code from the class was designed for a current input.
% The code needed for this assignment involved a voltage input.  The
% $A$ from fib1.m had to be modified to account for the fixed voltage
% at the first node.  This led to a system of the form $A^TGAx=A^TGf$,
% where $f$ captures the voltage applied to the network.



%% Code


clear
close all
%Much of the code below is from the course notes.

N = 6;

R_i = 1e3;    % axial resistance (Ohm)
R_m = 24.9e3; % lateral resistance (Ohm)

m = 2*N;      % number of edges
n = N+1;      % number of nodes

A = zeros(m,n); % initialize adjacency matrix

% build A one compartment at a time
for j=1:N                        
    A(2*j-1:2*j,j:j+1) = [-1 1;0 -1];
end


A = A(:,2:end);  %Modification for voltage input

G = zeros(m,1);  % initialize and assemble G
G(1:2:m) = 1/R_i;
G(2:2:m) = 1/R_m;
G = diag(G);

S = A'*G*A;	 % assemble main matrix

f = zeros(m,1);  % assemble the stimulus vector
f(1) = 12;
                 % solve for the potentials for twelve volts
x12 = S \(A'*G*f) ; 
x12=[12 ; x12];
% Potentials for 6 Volts are just half of those for 12 volts
x6  = x12/2;


z = 0:1:N;      % nodal positions for plotting
data12=[11.95 10.12 8.71 7.63 6.86 6.36 6.12 ]
data6 =[6.13  5.19 4.46 3.91 3.52 3.26 3.14];



fprintf('Data for 12 Volts\n')
fprintf('Node #   Calculated Potential   Measured Potential\n')
for i=0:6
    fprintf('%3d %20.2f %20.2f\n',i,x12(i+1),data12(i+1))
end
figure(1), clf

plot(z,data12,'b.'), hold on  
plot(z,x12,'ko','markersize',10)
legend('measured','calculated')
xlabel('Node number','fontsize',16)
ylabel('Potential (V)','fontsize',16)
title('Potential for 12 Volts','fontsize',16) 



fprintf('Data for 6 Volts\n')
fprintf('Node #   Calculated Potential   Measured Potential    Ratio of Measured 6 and 12\n')
for i=0:6
    fprintf('%3d %20.2f %20.2f %24.3f\n',i,x6(i+1),data6(i+1),data6(i+1)/data12(i+1))
end
figure(2), clf
plot(z,data6,'b.'), hold on    
plot(z,x6,'ko','markersize',10)
legend('measured','calculated')
xlabel('Node number','fontsize',16)
ylabel('Potential (V)','fontsize',16)
title('Potential for 6 Volts','fontsize',16) 


%% Results and Discussion 
% The tables above shows, for the most part, that the computed data
% agrees well with the measured data.  Clearly, much of the
% discrepancy between the two results from having used solely the
% display on the power source to set the input voltage.  For the 12
% Volt scenario this led to an observed voltage of 11.95V and for the
% 6 Volt scenario an observed voltage of 6.13 V.  None the less the
% measured data does indicate that the voltages at the nodes are
% linear with respect to the input voltage.  The ratio between the two
% sets of measured voltages was roughly .513 due to the input voltages
% being 11.95V and 6.13V.


% Part 2 the Monte Carlo Simulation 
% Randomly change the
% resistor values based on .01 tolerance.  Only did simulation for 12
% Volts as 6 Volts is similar.
for j=1:1000
        
        %Build the adjacency matrix
        %First cell is special        
        R_ip(1:N)=R_i*(1+.01*randn(N,1));
        R_mp(1:N)=R_m*(1+.01*randn(N,1));
        % initialize and assemble G
        G = zeros(m,1);                 
        G(1:2:m) = 1./R_ip;
        G(2:2:m) = 1./R_mp;
        G = diag(G);
        %Solve for the voltages
        V=(A'*G*A)\(A'*G*f);
        data12m(j,:)=V.';       
end


clear datac



figure(4),clf
title('Monte Carlo Results 12 V - the Mean')
lnames={'x_1','x_2','x_3','x_4','x_5','x_6'};
boxplot(data12m-ones(1000,1)*mean(data12m),'labels',lnames),hold on
plot(1:N,data12(1,2:N+1)-mean(data12m),'*-')


%% 
% A Monte Carlo simulation was performed to account for the 1%
% tolerance of the resistors.  The box plots indicate the distribution
% of the voltages at each of the nodes.  For the 12 Volt scenario,
% with the exception of the first node, all of the measured data falls
% in the inter-quartile range, the middle fifty percent of the data.
