
%
%  dcab1.m
%  Steve Cox, CAAM 335
%
%  A driver for cab1
%
function dcab1

cab1f(8,'ko-')
hold on
cab1f(16,'kx-')
cab1f(32,'k+-')
cab1f(64,'k*-')
legend('N=8','N=16','N=32','N=64')
hold off
return

%  cab1f.m
%  Steve Cox, CAAM 335
%
%  Compute and plot (with linetype typ) the potentials 
%  in an N compartment cable resulting from a current 
%  stimulus at its left end
%
%  usage:    cab1f(N,typ)
%  example:  cab1f(16,'k*-')
%

function cab1f(N,typ)

ell = 1;	% fiber length (cm)
a = 0.025;	% fiber radius (cm)
rho_i = 30;	% cytoplasmic resistivity (Ohm cm)
rho_m = 1000;	% membrane resistivity (Ohm cm^2)

R_i = rho_i*(ell/N)/(pi*a*a);  	% axial resistance (Ohm)
R_m = rho_m/(2*pi*a*ell/N);	% lateral resistance (Ohm)

i_0 = .001;	% current stimulus (mA)

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

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

for j=1:N,	% build A one compartment at a time

    A(2*j-1:2*j,j:j+1) = [-1 1;0 -1];

end

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(n,1);	% assemble the stimulus vector
f(1) = i_0;

x = S \ f;	% solve for the potentials

z = 0:ell/N:ell;	% nodal positions for plotting

plot(z,x,typ)	% plot and label potential vs. distance
box off
xlabel('z  (cm)','fontsize',16)
ylabel('x  (mV)','fontsize',16)
