%
%  fib1.m
%  Steve Cox, CAAM 335, August 1, 2000
%
%  Compute and plot the potentials in an N compartment
%  fiber resulting from a current stimulus.
%

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)

N = 64;		% number of compartments

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

y = -G*A*x;	% compute the associated currents

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

figure(1)
plot(z,x)	% plot and label potential vs. distance

xlabel('z (cm)','fontsize',16)
ylabel('x (mV)','fontsize',16)
tlab = ['Potential along ' num2str(N) ' compartment fiber'];
title(tlab,'fontsize',16) 

figure(2)
z = z(2:N+1);
y = 1000*y;		% scale for plotting purposes
plot(z,y(1:2:m)) 	% plot axial current vs. distance

xlabel('z (cm)','fontsize',16)
ylabel('y (\mu A)','fontsize',16)
tlab = ['Axial current along ' num2str(N) ' compartment fiber'];
title(tlab,'fontsize',16) 

figure(3)
plot(z,y(2:2:m)) 	% plot membrane current vs. distance

xlabel('z (cm)','fontsize',16)
ylabel('y (\mu A)','fontsize',16)
tlab = ['Membrane current along ' num2str(N) ' compartment fiber'];
title(tlab,'fontsize',16)