% Script file cabunbact1coct.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
 % March 9,12                        Asha G                        for 1c
 % 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 ('cabunbact1coct')


% declare global values
global N l Rm Ri t tau I d Cm rinp  Ao r lambda Iapp
N = 61;
l = 400*10^-4; % in cm
r= 1.85*10^-4;% in cm
d = 2*r; % in cm
Rm = 10*10^-3;% 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
gNa0= 168;%mS/cm^2
gK0 = 12.5; % mS/cm^2
lambdaNa = -0.025*(10^4);%cm^-1
lambdaK = -0.025*(10^4);
Ao= pi*(r)^2; % crosssection area at 0 in cm^2
%niter=1;
niter = 1098900;
%niter = 274730;  
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;
alpha2 = 11;
%[ Q,M] = tridiagdoublesdoct(alpha,alpha1,f);
[Q,Qs,M,Ms]=tridiagoct(alpha,alpha1,alpha2,f);

% calling function ionact1coct
[Iion,INa,IK,IL,gK,gNa,gNabar,gtot]= ionact1coct(gK0,gNa0,lambdaNa,gL,n,m,h,V,VNa,VK,VL);
%[Iion,INa,IK,IL,gKbar,gNabar,gK,gNa] = ionactoct(gK0,lambdaK,gNa0,lambdaNa,gL,n,m,h,V,VNa,VK,VL);
%disp(gNabar)
% 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*iter*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);
elseif  ( 5< (deltaT*iter*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);
y = vvvv';

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


save -mat-binary cabunbact1coctN61octtest x vv1 vv2 z vvvv  
%save -mat-binary cabunbact1coctN21mesh x vv1 vv2 z vvvv
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)
%colormap([1 1 1])
%view(30,45)
%[xx,zz]= meshgrid(x,z);
%mesh(xx,zz,y)

xlabel('t in msec',"fontsize",40)
ylabel('V in mV',"fontsize",40)
%zlabel ('V in mV',"fontsize",40)
%ylabel('X in cm',"fontsize",40)
legend ('{\fontsize{20} voltage at N=1}','{\fontsize{20}voltage at N=N}')
  title('sealed end N = 61, V vs t ',"fontsize",40)


print ('cabunbact1coctN61octtest.eps','-depsc')
%print ('cabunbact1coctN21mesh.eps','-depsc')

hold off