% Script file cabunbact1aoct.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
 %July 4, 08                         Asha G                       original code 
 % Feb 22,12                         Asha G                       rewritten again
 % 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 
clear all
%close all
%format long 

pack ('cabunbact1aoct')


% declare global values
global N l Rm Ri t tau I d Cm rinp  Ao r lambda Iapp
N = 10;
l = 400*10^-4; % in cm
r= 1.85*10^-4;% in cm
d = 2*r; % in cm
Rm = 20000;% ohms.cm^2
Ri=333.33;% ohms.cm
Cm = 1*(10^-6); % farad / cm^2
lambdaDC = ((Rm/Ri)* (d/4))^(1/2);%cm
rinf = ((Rm*Ri)^(1/2))*(2/(pi*d^(3/2)));
ri = 8.9*(10^7) ; % ohmscm^-1
%t = (1:niter)*deltaT; % in msec
tau = ((Rm*10^-3)*Cm*10^6);% msec
f= 1/(2*pi*tau*(10^-3)); % Hz
lambda = (1/2)*(d/(pi*f*Ri*Cm))^(1/2);%cm
L= l/lambda;
ginp = (pi*((d)^(3/2))*tanh(L))/((2*Ri*Rm)^(1/2));% Siemen
rinp = (1/ginp); % ohms
%I= 0.1*(10^-9); % amperes
tmax = 55;% msec
%I = 450*(10^-6); % amperes/cm^2
Iapp = (0.1*10^-9)*ones(1,N);
rm = Rm/(pi*d);% ohm/cm
ra = (4*Ri)*(pi*d^2);% ohm/cm^2
VNa = 55; % mV
VK = - 95; % mV
VL = -60; % mV
gL = 0.1 ; %mS/cm^2
gNa0 = 50 ; %mS/cm^2
gK0 = 12.5; % mS/cm^2
Ao= pi*(r)^2; % crosssection area at 0 in cm^2

%niter = 6182;
niter = 6182;  
vvvv=zeros(niter,N); abmm = zeros(niter,N); abnn = zeros(niter,N); abhh = zeros(niter,N); V= zeros(1,N);alphan = zeros(1,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);%U = zeros(1,N);abm= zeros(1,N); abn = zeros(1,N); abh = zeros(1,N);vv = zeros(niter,1); 


% calling function voltseoct_cabunpa4
[V,hd,h2] = voltseoct_cabunpa4(lambda); 

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

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

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

 
for iter = 1:niter 
% calling function compactdiffcaboct
[f] = compactdiffcaboct(V,h2);

% calling function tridiagoct

alpha = 2/11;
alpha1 = 1/10;
%[ Q,M] = tridiagdoublesdoct(alpha,alpha1,f);
[Q,Qs,M,Ms]=tridiagoct(alpha,alpha1,f);

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


% calling timedepcabspvoltoct
gamconst = 1;

% calling function timedepcabspvoltoct
[P,S,deltaT]= timedepcabspvoltoct(Q,Qs,V,h2,Iion,gamconst);
%disp(deltaT)

for jji = 2:( N-1)
V(jji) = S(jji);
endfor 

if (deltaT*niter*tau)< 5
I = 0;
V(1)= (4/3)*V(2)-(1/3)*V(3)+(2/3)*(4*Ri*lambda*I*hd)*(10^3)/(pi*(d^2));% 2nd order formula
V(N) = (4/3)*V(N-1)-(1/3)*V(N-2);

%**************************************next line changed by Doug Stewart
elseif (   5== (deltaT*niter*tau) <= 55)
I = 4.8384*(10^-11);% amperes 
V(1)= (4/3)*V(2)-(1/3)*V(3)+(2/3)*(4*Ri*lambda*I*hd)*(10^3)/(pi*(d^2));% 2nd order formula
V(N)= (4/3)*V(N-1)-(1/3)*V(N-2);
endif 



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


abm = alpham + betam;
%disp(abm)
abh = alphah + betah;

abn = alphan + betan;
% calling function tauinfoct
[taum,tauh,taun,minf,hinf,ninf]= tauinfoct(alpham,betam,alphah,betah,alphan,betan);

% calling function timedepmnhoct
gamconst = 1;

[n,m,h] = timedepmnhoct(ninf,n,taun,minf,m,taum,hinf,h,tauh,h2,gamconst);


vvvv(iter,:) = V;
abmm(iter,:)=abm;
abnn(iter,:)=abn;
abhh(iter,:)=abh;


endfor 

hold on 
x= linspace(0,tmax,niter);
z = linspace(0,l,N);


vv1 = vvvv(:,1);
vv2 = vvvv(:,N);

y = vv1;


 
save -mat-binary cabunbact1aoctmesh x z vvvv vv1
set (findall (gcf, "-property", "linewidth"), "linewidth",2)
set (findall (gcf, "-property", "markersize"), "markersize",6)
set(gca(),"fontsize",20)

%plot (x,vv1,'r',"linewidth",2)
%plot(x,vv2,'b--',"linewidth",2)
[xx,zz]= meshgrid(x,z);
mesh(xx,zz,vvvv')
%mesh(vvvv)
%surf(vvvv)
xlabel('t in msec',"fontsize",20)
ylabel ('V in mV',"fontsize",20)
zlabel('X in cm',"fontsize",20)
%legend ('{\fontsize{14} N=1}','{\fontsize{14}N=N}')
title('sealed end N = 10, V vs t & X ')

print ('cabunbact1aoctmesh.eps','-depsc')

hold off