% Script file cabunbact2oct.m 
% Purpose : To study V vs X and t for soma- dendrite, active, sealed, unbranched.  
% Current injection at x = 0 
% Record of revisions:
% Date                             Programmer                       Description
 %June 2, 08                         Asha G                       original code
 % Define variables 
% N -- discretization by space 
% V -- output volt at points x
% l -- length of dendrite in mm
% lambda -- space constant in mm
% L -- length nondimensionalised - l/lambda
% h -- delta x = l/N
% h2 -- h*h
% f -- R.H.S of compact difference scheme
% alpha -- coeff of L.H.S compact diff scheme
% alpha1 -- coeff of R.H.S compact diff scheme
% M -- tridiag matrix
% Q -- double derivative solved by M\f'
% k -- index variable 

% calling function voltseoct_cabunpa4

pack ('cabunbact2oct')

% declare global values
global N l Rm Ri t tau I d Cm rinp gamma Ao
N = 10;
l = 400*10^-4; % in cm
r= 1.85*10^-4;% in cm
d = 2*r; % in cm
%Rm = 20* 10^5; % ohms.mm^2
Rm = 20;% Kohms.cm^2
Ri=330;%Kohms.cm
Cm = 1; % microfarad / cm^2
%rinp = 1* 10^6; % in ohm
niter = 5000;
%deltaT = 2.2*10^-2;% msec
deltaT= 0.001;% msec
t= (1:niter)*deltaT;% msec 
%t = 20; % in msec
tau = 10; % in msec
I = 450; % microamperes/cm^2
lambda = ((( Rm/Ri)* (d/4))^1/2);% cm
%rm = Rm/(pi*d);% ohm/cm
%ra = (4*Ri)*(pi*d^2);% ohm/cm^2
%Vo = - I*rinp;
Vo = -70; % mV
VNa = 55; % mV
VK = - 95; % mV
VL = -60; % mv
gL = 0.1 ; %mS/cm^2
gNao = 50 ; %mS/cm^2
lambdaNa = -0.0075*10^-4; % cm^-1
gKo = 12.5; % mS/cm^2
lambdaK = -0.0075*10^-4; % cm^-1 
gamma = 1;
Ao= pi*(r)^2; % crosssection area at 0 in cm^2

vvvv=zeros(niter,N); abmm = zeros(niter,N); abnn = zeros(niter,N); abhh = zeros(niter,N);f=zeros(niter,N); Q = zeros(niter,N); M = zeros(niter,N);Iion = zeros(niter,N);INa = zeros(niter,N);IK = zeros(niter,N); IL = zeros(niter,N); gNabar = zeros(niter,N); gKbar = zeros(niter,N); gNa = zeros(niter,N); gK = zeros(niter,N); P = zeros(niter,N); S = zeros(niter,N); V= zeros(niter,N); alphan = zeros(niter,N); alpham = zeros(niter,N); alphah = zeros(niter,N); betam = zeros(niter,N); betan = zeros(niter,N); betah= zeros(niter,N); m= zeros(niter,N); n = zeros(niter,N); h = zeros(niter,N);
% calling function voltseoct_cabunpa5
[V,h2,hd] = voltseoct_cabunpa5(lambda,Vo); 

% calling function ratcon2oct
[alphan,betan,alpham,betam,alphah,betah]= ratcon2oct(V);


% calling function gatpara2oct
[n,m,h]= gatpara2oct(alphan,betan,alpham,betam,alphah,betah);

for iter = 1:niter 

[f] = compactdiffcaboct(V,h2);

% calling function tridiagdoubleoct

alpha = 2/11;
alpha1 = 1/10;
[ Q,M] = tridiagdoubleoct(alpha,alpha1,f);

% calling function ionactoct
[Iion,INa,IK,IL,gNabar,gKbar,gNa,gK]= ionactoct(gNao,lambdaNa,gKo,lambdaK,gL,n,m,h,V,VNa,VK,VL);

% calling function timedepcabactoct

[P,S]= timedepcabactoct(Q,V,h2,Iion,deltaT);

for jji = 2:( N-1)
V(jji) = S(jji);
end
if deltaT*niter<5
I = 0;
elseif 5 == deltaT*niter <= 55
%I = 450*pi*(r)^2 ;
I = 48.4* 10^-12;% Amperes ? 
V(1) = V(2)+(Ri/Ao)*I*hd;
V(N) = (4/3)* V( N-1)-(1/3)*V(N-2);
end
disp(V)
% calling function ratcon2oct
[alphan,betan,alpham,betam,alphah,betah]=ratcon2oct(V);
abm = alpham + betam;
m = m + deltaT.*(alpham - m.*(abm));
abh = alphah +betah;
h = h +deltaT.*(alphah-h.*(abh));
abn= alphan +betan;
n= n + deltaT.*(alphan-n.*(abn));

%for kkk= 2:N-1
% calling function ratcon2oct
%[alphan,betan,alpham,betam,alphah,betah]= ratcon2oct(V);
% calling function tauinfoct
%[taum,tauh,taun,minf,hinf,ninf]= tauinfoct(alpham,betam,alphah,betah,alphan,betan);
%m = minf + (m-minf).*exp(-deltaT./(2.*taum));
%n = ninf + (n-ninf).* exp(-deltaT./(2.*taun));
%h = hinf +(h-hinf).*exp(-deltaT./(2.*tauh));
%end

vvvv(iter,:) = V;
%alphann(iter,:) = alphan;
%alphahh(iter,:) = alphah;
%alphamm(iter,:) = alpham; 
abmm(iter,:)=abm;
abnn(iter,:)=abn;
abhh(iter,:)=abh;
gNait(iter,:)= gNa;
end 

hold on 
gNait1 = gNait(:,1);
abmm5=abmm(:,5);
abhh5=abhh(:,5);
abnn5=abnn(:,5);
vv = vvvv(:,1);
%vv1 = vvvv(:,N/2);
x = linspace(0,l,N);
y=linspace(0,t,niter);
%[X,Y]= meshgrid(x,y);
%Z = vvvv;
%mesh(X,Y,Z)
subplot(2,2,1)
plot(t,abmm5)
xlabel('t')
ylabel('alpha+beta')
legend('abmm5')
subplot(2,2,2)
plot(t,abhh5)
xlabel('t')
ylabel('alpha + beta')
legend('abhh5')
subplot(2,2,3)
plot(t,abnn5)
xlabel('t')
ylabel('alpha+beta')
legend('abnn5')
%plot3(x,y,vv,'r-')
%plot(y,vv)
%grid on
print ( 'cabunbact2oct450ab5.eps', '-depsc')
hold off