% Script file cabranbact13oct.m 
% Purpose : To study V vs X and t for soma- dendrite, active, sealed, branched, problem Toth 1999, figure 3.   
% 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 1 b 
    % Dec 4,12                      Asha G                        change Rm value from 20,000 to 10*10^-3;% ohms.cm^2
    % Jan 28,2012                   Asha G                         rewritten for Toth 1999,fig.3. 
 % 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 ('cabranbact13oct')

%Part 1: L1
% declare global values
global N l Rm Ri t tau I d Cm rinp  Ao rp lambda Iapp lb1 lp rb1 db1 lb2 rb3 lb4
N = 11;
lp = 10*10^-4; % in cm
rp= 5.0*10^-4;% in cm
dp = 2*rp; % in cm
Rm = 10*10^-3;% ohms.cm^2
Ri=333.33;% ohms.cm
Cm = 1*(10^-6); % farad / cm^2
lambdaDC = ((Rm/Ri)* (dp/4))^(1/2);%cm
%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)*(dp/(pi*f*Ri*Cm))^(1/2);%cm
L= lp/lambda;% check this 
ginp = (pi*((dp)^(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*dp);% ohm/cm
ra = (4*Ri)*(pi*dp^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
lambdaNa = 0; % cm^-1
lambdaK = 0;% cm^-1
Ao= pi*(rp)^2; % crosssection area at 0 in cm^2

%niter = 488400;
%niter = 1953600;
niter = 2;
  
%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 voltoct_cabranpap
[VVp,hp,hp2] = voltoct_cabranpap(lp,lambda);

% calling function volt_cabranpab1 ( L3)
lb1 = 100*(10^-4); % in cm
rb1 = 1.8*(10^-4); % in cm
db1 = rb1*2; % in cm 
lambdab1DC = ((Rm/Ri)*(db1/4))^(1/2);% cm
tau = ((Rm*10^-3)*Cm*10^6);%msec
f= 1/(2*pi*tau*(10^-3)); % Hz
lambdab1 = (1/2)*(db1/(pi*f*Ri*Cm))^(1/2);%cm
%Lb1 = lb1/lambdab1; % check 
rab1 = (4*Ri)/(pi*db1^2);% ohm/cm
rmb1 = Rm/(pi*db1);% ohm.cm
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
lambdaNab1 = -2; % cm^-1
lambdaKb1 = -2;% cm^-1
Aob1= pi*(rb1)^2; % crosssection area at 0 in cm^2

[VVb1,hb1,hb12]= volt_cabranpab1(VVp,lp,lb1,lambda);

% calling function voltoct_cabranpab2 L2
lb2 = 100*(10^-4); % in cm
rb2o = 2.5*(10^-4); % in cm
db2o= rb2o*2;
rb2l = 2.0*(10^-4);% in cm
db2l = rb2l*2;
rhob2 = (rb2l-rb2o)/lb2;
lambdab2DC = ((Rm/Ri)* (db2o/4))^(1/2);% cm
%rab2 = (4*Ri)/(pi*db2^2);
%rmb2 = Rm/(pi *db2);
x = linspace(0,lb2,N);% cm
rxb2= rb2o+rhob2.*x ;%cm
lambdab2DCtap = lambdab2DC*((rxb2/rb2o)*(1/2))*(1+(rhob2)^2)^(-1/4);% cm
tau = ((Rm*10^-3)*Cm*10^6);%msec
f= 1/(2*pi*tau*(10^-3)); % Hz
lambdab2AC = (1/2)*(db2o./(pi*f*Ri*Cm))^(1/2);%cm 
lambdab2ACtap = lambdab2AC.*(1+(2*rhob2.*x)./db2o).^(1/2);%cm
lambdab2 = (lambdab2ACtap(1)+lambdab2ACtap(N))/2;% cm
X= x/lambda ; 
%Lb2 = lb2/lambdab2
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
lambdaNab1 = -2; % cm^-1
lambdaKb1 = -2;% cm^-1
Aob2o= pi*(rb2o)^2; % crosssection area at 0 in cm^2


[VVb2,hb2,hb22]= volt_cabranpab2(VVp,lp,lb2,lambda);


% calling function volt_cabranpab3 L5
lb3 = 50*(10^-4); % in cm
rb3 = 1.8*(10^-4); % in cm
db3 = rb3*2; % cm
lambdab3DC = ((Rm/Ri)*(db3/4))^(1/2);% cm
tau = ((Rm*10^-3)*Cm*10^6);%msec
f= 1/(2*pi*tau*(10^-3)); % Hz
lambdab3 = (1/2)*(db3/(pi*f*Ri*Cm))^(1/2);%cm
%Lb3 = 
rab1 = (4*Ri)/(pi*db1^2);% ohm/cm
rmb1 = Rm/(pi*db1);% ohm.cm
VNa = 55; % mV
VK = - 95; % mV
VL = -60; % mV
gL = 0.1 ; %mS/cm^2
gNa0b3 = 0 ; %mS/cm^2
gK0b3 = 2.1; % mS/cm^2
lambdaNab3 = 0; % cm^-1
lambdaKb3 = 0;% cm^-1
Aob3= pi*(rb3)^2; % crosssection area at 0 in cm^2

[VVb3,hb3,hb32]= volt_cabranpab3(VVp,lp,lb3,lambda);

% calling function volt_cabranpab4 L4
lb4 = 300*(10^-4); % in cm
rb4o = 1.0*(10^-4); % in cm
db4o= rb4o*2;
rb4l = 0.3*(10^-4);% in cm
db4l = rb4l*2;
rhob4 = log(rb4o/rb4l)/lb4;%cm^-1
xb4 = linspace(0,lb4,N);% cm
rxb4 = rb4o.*(exp(-rhob4.*xb4));% in c
lambdab4DC = ((Rm/Ri)* (db4o/4))^(1/2);% cm
%rab2 = (4*Ri)/(pi*db2^2);
%rmb2 = Rm/(pi *db2);
lambdab4DCtap = lambdab4DC*((rxb4/rb4o)*(1/2))*(1+(rhob4)^2)^(-1/4);% cm
tau = ((Rm*10^-3)*Cm*10^6);%msec
f= 1/(2*pi*tau*(10^-3)); % Hz
lambdab4AC = (1/2)*(db4o./(pi*f*Ri*Cm))^(1/2);%cm 
lambdab4ACtap = lambdab4AC.*(1+(2*rhob4.*x)./db4o).^(1/2);%cm
lambdab4 = (lambdab4ACtap(1)+lambdab4ACtap(N))/2;% cm
X= x/lambda ;  
%Lb4 = 
VNa = 55; % mV
VK = - 95; % mV
VL = -60; % mV
gL = 0.1 ; %mS/cm^2
gNa0b4 = 2.1 ; %mS/cm^2
gK0b4 = 8.4; % mS/cm^2
lambdaNab4 = -2; % cm^-1
lambdaKb4 = -2;% cm^-1
Aob4o= pi*(rb4o)^2; % crosssection area at 0 in cm^2

[VVb4,hb4,hb42]= volt_cabranpab4(VVp,lp,lb4,lambda);


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


% calling function gatpara2oct13
[np,mp,hp,nb1,mb1,hb1,nb2,mb2,hb2,nb3,mb3,hb3,nb4,mb4,hb4]= gatpara2oct13(alphan,betan,alpham,betam,alphah,betah);

for iter = 1:niter 

% Calling function compactdiffcaboctp

[fp] = compactdiffcaboctp(VVp,hp2);

% Calling function tridiagoctp

alpha = 2/11;
alpha1 = 1/10;
alpha2 = 11;

[Qp,Qps,M,Ms]=tridiagoctp(alpha,alpha1,alpha2,fp);

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


% Calling function timedepcabranvoltoctp
gamconst = 1;

[Pp,Sp,deltaT] = timedepcabranvoltoctp(Qp,VVp,hp2,Iion,gamconst );

%disp(deltaT)
%for jj = 2: N-1
%VVp(jj) = Sp(jj);
%end 

jj = 2: N-1;
VVp(jj) = Sp(jj);

%VVp(1) = VVp(2)+rinp*I*taufactor*hp;
VVp(1)= (4/3)*VVp(2)-(1/3)*VVp(3); %( assuming sealed end - 2nd order )

% calling function ratcon2octp
[alphanp,betanp,alphamp,betamp,alphahp,betahp]= ratcon2octp(VVp);


% calling function tauinfoctp
[taump,tauhp,taunp,minfp,hinfp,ninfp]= tauinfoctp(alphamp,betamp,alphahp,betahp,alphanp,betanp);

% calling function timedepmnhoctp
gamconst = 1;

[np,mp,hp] = timedepmnhoctp(ninfp,np,taunp,minfp,mp,taump,hinfp,hp,tauhp,hp2,gamconst);

vvvvp(iter,:) = VVp;

%endfor

% section  branch L3 :(b1)

% calling function compactdiffcaboctb1

[fb1] = compactdiffcaboctb1(VVb1,hb12);

% calling function tridiagoctb1

alpha = 2/11;
alpha1 = 1/10;
alpha2 = 11;

%[Qb1,M] = tridiagdoubleoctb1(alpha,alpha1,fb1);
[ Qb1,Qb1s,M,Ms] = tridiagoctb1(alpha,alpha1,alpha2,fb1);

gK0b1 = 12.5;
gNa0b1 = 50;
lambdaKb1 = -2;
lambdaNab1 = -2;

 
% calling function ionactoctb1
[Iionb1,INab1,IKb1,ILb1,gKbarb1,gNabarb1,gKb1,gNab1]= ionactoctb1(gK0b1,lambdaKb1,gNa0b1,lambdaNab1,gL,nb1,mb1,hb1,VVb1,VNa,VK,VL);
%disp( Iionb1)
% calling function timedepcabranvoltoctb1
gamconst = 1;


[Pb1,Sb1,deltaT] = timedepcabranvoltoctb1 ( Qb1,VVb1, hp2,Iionb1,gamconst);
%disp(deltaT)

%for jj = 2:N-1
%VVb1(jj) = Sb1(jj);
%endfor 
jj = 2:N-1;
VVb1(jj) = Sb1(jj);

if (deltaT*iter*tau)< 5
I = 0;
VVb1(1)= VVp(N);
VVb1(N)= (4/3)*VVb1(N-1)-(1/3)*VVb1(N-2);%+(2/3)*(4*Ri*lambdab1*I*hb1)*(10^3)/(pi*(db1^2));

elseif  ( 5== (deltaT*iter*tau) <= 45)
I = 1.9849*(10^-11);% amperes 
VVb1(1)= VVp(N);
VVb1(N)= (4/3)*VVb1(N-1)-(1/3)*VVb1(N-2)+(2/3)*(4*Ri*lambdab1*I*hb1)*(10^3)/(pi*(db1^2));
endif 

% calling function ratcon2octb1
[alphanb1,betanb1,alphamb1,betamb1,alphahb1,betahb1]= ratcon2octb1(VVb1);

% calling function tauinfoctb1
[taumb1,tauhb1,taunb1,minfb1,hinfb1,ninfb1]= tauinfoctb1(alphamb1,betamb1,alphahb1,betahb1,alphanb1,betanb1);

% calling function timedepmnhoctb1
gamconst = 1;

[nb1,mb1,hb1] = timedepmnhoctb1(ninfb1,nb1,taunb1,minfb1,mb1,taumb1,hinfb1,hb1,tauhb1,hb12,hp2,gamconst); 

vvvvb1(iter,:) = VVb1;
%endfor 

% section on L2 (b2)

% calling function compactdiffcaboctb2

[fb2] = compactdiffcaboctb2(VVb2,hb22);

% calling function tridiagoctb2

alpha = 2/11;
alpha1 = 1/10;
alpha2 = 11;

[ Qb2,Qb2s,M,Ms] = tridiagoctb2(alpha,alpha1,alpha2,fb2);


% calling function compactdiffcabfdoctb2
[ffdb2] = compactdiffcabfdoctb2(VVb2,hb2);

alphafd = 1/3;
%alphafd = 1/2;
alpha1fd = 1/4;
alphafdb = 2;
%alphafd2 = 5;

% calling function tridiagoctaperedb2
[ Qfdb2,Qfdb2s,Mfd,Mfds] = tridiagoctaperedb2(alphafd,alpha1fd,alphafdb,ffdb2);
 
%disp(Qfd)

gK0b2 = 12.5;
gNa0b2 = 50;
lambdaKb2 = -2;
lambdaNab2 = -2;



% calling function ionactoctb2
[Iionb2,INab2,IKb2,ILb2,gKbarb2,gNabarb2,gKb2,gNab2]= ionactoctb2(gK0,lambdaKb2,gNa0,lambdaNab2,gL,nb2,mb2,hb2,VVb2,VNa,VK,VL);
 

% calling function timedepcabranvoltoctb2
gamconst = 1;
[Pb2,Sb2,deltaT] = timedepcabranvoltoctb2( Qb2,Qfdb2,VVb2, hb22,Iionb2,gamconst);


%for jji = 2:( N-1)
%VVb2(jji) = Sb2(jji);
%endfor 
jji = 2:( N-1);
VVb2(jji) = Sb2(jji);

VVb2(1) = VVp(N);
VVb2(N) % branchpoint 

% calling function ratcon2octb2
[alphanb2,betanb2,alphamb2,betamb2,alphahb2,betahb2] = ratcon2octb2(VVb2);


% calling function tauinfoctb2
[taumb2,tauhb2,taunb2,minfb2,hinfb2,ninfb2]= tauinfoctb2(alphamb2,betamb2,alphahb2,betahb2,alphanb2,betanb2);
% calling function timedepmnhoctb2
gamconst = 1;

 
 [nb2,mb2,hb2] = timedepmnhoctb2(ninfb2,nb2,taunb2,minfb2,mb2,taumb2,hinfb2,hb2,tauhb2,hb22,hp2,gamconst);

vvvvb2(iter,:) = VVb2;


%endfor 

% section on L5 ( b3) 

% calling function compactdiffcaboct

[fb3] = compactdiffcaboctb3(VVb3,hb32);

% calling function tridiagoctb3

alpha = 2/11;
alpha1 = 1/10;
alpha2 = 11;

[ Qb3,Qb3s,M,Ms] = tridiagoctb3(alpha,alpha1,alpha2,fb3);


gK0b3 = 2.1;
gNa0b3 = 0;
lambdaKb3 = 0;
lambdaNab3 = 0;



% calling function ionact4octb3

[Iionb3,INab3,IKb3,IL,gKb3,gNab3,gtotb3]= ionact4octb3(gNa0b3,gK0b3,gL,nb3,mb3,hb3,VVb3,VNa,VK,VL);


% calling function timedepcabranvoltoctb3
gamconst = 1;
[Pb3,Sb3,deltaT] = timedepcabranvoltoctb3(Qb3,VVb3,hb32,Iionb3,gamconst);

%for jji = 2:( N-1)
%VVb3(jji) = Sb3(jji);
%endfor 

jji = 2:( N-1)
VVb3(jji) = Sb3(jji);

VVb3(1) = VVb2(N);
VVb3(N) = (4/3)*VVb3(N-1)-(1/3)*VVb3(N-2);

% calling function ratcon2octb3
[alphanb3,betanb3,alphamb3,betamb3,alphahb3,betahb3]= ratcon2octb3(VVb3);

% calling function tauinfoctb3
[taumb3,tauhb3,taunb3,minfb3,hinfb3,ninfb3]= tauinfoctb3(alphamb3,betamb3,alphahb3,betahb3,alphanb3,betanb3);

% calling function timedepmnhoctb3
gamconst = 1;

[nb3,mb3,hb3] = timedepmnhoctb3(ninfb3,nb3,taunb3,minfb3,mb3,taumb3,hinfb3,hb3,tauhb3,hb32,hp2,gamconst); 

vvvvb3(iter,:) = VVb3;


%endfor 

% section L4 ( b4)

% calling function compactdiffcaboctb4

%[fb4] = compactdiffcaboctb4(VVb4,hb42);

% calling function tridiagoct

alpha = 2/11;
alpha1 = 1/10;
alpha2 = 11;

%[ Qb4,Qb4s,M,Ms] = tridiagoctb4(alpha,alpha1,alpha2,fb4);


gK0b4 = 2.1;
gNa0b2 = 8.4;
lambdaKb2 = -2;
lambdaNab2 = -2;

% calling function compactdiffcabfdoctb4
%[ffdb4] = compactdiffcabfdoctb4(VVb4,hb4);

alphafd = 1/3;
%alphafd = 1/2;
alphafd1 = 1/4;
alphafd2 = 2;
%alphafd2 = 5;

% calling function tridiagoctaperedb4
%[ Qfdb4,Qfdb4s,Mfd,Mfds] = tridiagoctaperedb4(alphafd,alphafd1,alphafd2,ffdb4);
 
%disp(Qfd)


% calling function ionactoctb4
[Iionb4,INab4,IKb4,ILb4,gKbarb4,gNabarb4,gKb4,gNab4]= ionactoctb4(gK0b4,lambdaKb4,gNa0b4,lambdaNab4,gL,nb4,mb4,hb4,VVb4,VNa,VK,VL);

% calling function timedepcabranvoltoctb4
gamconst = 1;

%[Pb4,Sb4,deltaT]= timedepcabranvoltoctb4(Qb4,Qfdb4,VVb4,hb42,Iionb4,gamconst);

%for jji = 2:( N-1)
%VVb4(jji) = Sb4(jji);
%endfor 
jji = 2:( N-1);
VVb4(jji) = Sb4(jji);

if (deltaT*iter*tau)< 10
I = 0;
VVb4(1)= VVb2(N);
VVb4(N)= (4/3)*VVb4(N-1)-(1/3)*VVb4(N-2);%+(2/3)*(4*Ri*lambdab4*I*hb4)*(10^3)/(pi*(db4l^2));

elseif  ( 10 == (deltaT*iter*tau) <= 50)
I = 8.4823*(10^-11);% amperes 
VVb4(1)= VVb2(N);
VVb4(N)= (4/3)*VVb4(N-1)-(1/3)*VVb4(N-2)+(2/3)*(4*Ri*lambdab4*I*hb4)*(10^3)/(pi*(db4l^2));
endif 


% calling function ratcon2oct
[alphanb4,betanb4,alphamb4,betamb4,alphahb4,betahb4]= ratcon2octb4(VVb4);

% calling function tauinfoctb4
[taumb4,tauhb4,taunb4,minfb4,hinfb4,ninfb4]= tauinfoctb4(alphamb4,betamb4,alphahb4,betahb4,alphanb4,betanb4);

% calling function timedepmnhoctb4

gamconst = 1;

[nb4,mb4,hb4] = timedepmnhoctb4(ninfb4,nb4,taunb4,minfb4,mb4,taumb4,hinfb4,hb4,tauhb4,hb42,hp2,gamconst); 

vvvvb4(iter,:) = VVb4;

[Sbp1] = branchpointoctbr1(VVp,VVb1,VVb2,hp,hb1,hb2,dp,db1,db2o);

 
[Sbp2]= branchpointoctbr2(VVb2,VVb3,VVb4,hb2,hb3,hb4,db2l,db3,db4o);


%VVp(N) = Sbp1(N);
VVb1(1) = VVp(N); 
VVb2(1) = VVp(N);

%VVb2(N)= Sbp2(N);
VVb3(1)= VVb2(N);
VVb4(1)= VVb2(N);


nnnn(iter,:) = VVp;
oooo(iter,:)= VVb1;
pppp(iter,:) = VVb2;
qqqq(iter,:) = VVb3;
ssss(iter,:)= VVb4;

endfor 

tmax = 50; % msec 
hold on 
x= linspace(0,tmax,niter);
z = linspace(0,lp +lb2 +lb4,N);


vv1p = vvvvp(:,1);
vv2p = vvvvp(:,N);
vv1b1 = vvvvb1(:,1);
vv2b1 = vvvvb1(:,N);
vv1b2 = vvvvb2(:,1);
vv2b2 = vvvvb2(:,N);
vv1b3 = vvvvb3(:,1);
vv2b3 = vvvvb3(:,N);
vv1b4 = vvvvb4(:,1);
vv2b4 = vvvvb4(:,N);
yp = vvvvp';
yb1 = vvvvb1';
yb2 = vvvvb2';
yb3 = vvvvb3';
yb4 = vvvvb4';

%q = [X(:) y(:)]; q1 =[ X1(:) y1(:)]; q2= [X2(:) y2(:)]; 
save -mat-binary cabranbact13N11oct x vv1p vv2p vv1b1 vv2b1 vv1b2 vv2b2 vv1b3 vv2b3 vv1b4 vv2b4 z vvvvp vvvvb1 vvvvb2 vvvvb3 vvvvb4 
set (findall (gcf, "-property", "linewidth"), "linewidth",2)
set (findall (gcf, "-property", "markersize"), "markersize",6)
set(gca(),"fontsize",20)

plot(z,vv1p,'r',"linewidth",2)
plot(z,vv1b1,'r-.',"linewidth",2)
plot(z,vv1p,'b',"linewidth",2)
plot(z,vv1b2,'b-.',"linewidth",2)
plot(z,vv1b4,'b*',"linewidth",2, "markersize",6)
plot(z,vv1p,'g',"linewidth",2)
plot(z,vv1b2,'g-.', "linewidth",2)
plot(z,vv1b3,'g*',"linewidth",2,"markersize",6)
%%colormap([1 1 1])
%view(45,45)
%[xx,zz]= meshgrid(x,z);
%mesh(xx,zz,y)xlabel('t in msec',"fontsize",40)

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

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

hold off