% Script file cabunpanumcent2oct.m 
% Purpose : To study V/Vo vs X for soma - dendrite, passive, sealed end. 
% Current injection at x = 0 - comparison between numerical, analytical and central difference solutions
% Record of revisions:
% Date                             Programmer                       Description
% June 5, 07                         Asha G                       original code
% June 14,07                         Asha G                       remove II
% October 20,08                     Asha G                        correct for current inj at N= 1
% 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 
pack ('cabunpanumcent2oct')

% declare global values
global N l Rm Ri t tau I d Cm rinp ra
N = 20;
%l = 400; % in micron
l = 400*( 10^-4); % in cm 
%d = 2*1.85; % in micron
d = 2*(1.85*(10^-4));%cm
r = 1.85*(10^-4); % in cm
Rm = 20000; % ohms.cm^2
Ri = 330; % ohms.cm
Cm = 1; % microfarad / cm^2
%rinp = 1* 10^6; % in ohm
%rinp = (d^(-3/2))*((4*Rm*Ri)^(1/2))*(1/pi);
lambda = ((Rm/Ri)* (d/4))^(1/2); % cm
L= l/lambda;
%rinp = 1* 10^6; % in ohms
%rinp = (d^(-3/2))*((4*Rm*Ri)^(1/2))*(1/pi);
rinf = ((Rm*Ri)^(1/2))*(2/(pi*d^(3/2)));
%rinp = (Ri/pi*(d/2)^2)*lambda*coth(L) % kohms 
ginp = (pi*((d)^(3/2))*tanh(L))/((2*Ri*Rm*10^6)^(1/2));% milliSiemen
rinp = (1/ginp); % kohms
ri = 8.9*(10^7) ; % ohmscm^-1
%t = (1:niter)*deltaT; % in msec
tau = (Rm)*(Cm*(10^-6))*(10^3); %msec  
%tau = 5;% msec
taufactor = (tau)/(Cm*pi*d);
I= 0.1*(10^-9); % amperes

niter = 1; 

 V=zeros(niter,1);

%Calling voltseoct_cabunpa4
[V,h,h2] = voltseoct_cabunpa4(lambda); 
%disp(V)

for iter = 1:niter


[f]=compactdiffcab4oct(V,h2);

% calling function tridiagdoubleoct
% calling function tridiagoct


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

% calling function timedepcabvoltoct
% calling function timedepcabspvoltoct

gamconst = 1;
%[P,S] = timedepcabvoltoct ( Q,V,h2,gamconst);
[P,S]=timedepcabspvoltoct(Q,Qs,V,h2,gamconst);

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

endfor


%V(1)= V(2)+ rinp*I*h;
V(1) = V(2)+(4*Ri*lambda*I*h)/(pi*(d^2));
%V(1) = V(2)+(ri*I*h*lambda*(10^3));
V(N) = (4/3)* V( N-1)-(1/3)*V(N-2);
nnnn(iter,:) = V; 
endfor  
Vo = nnnn(iter,1);
y = V/Vo ;


% calling function voltoct_cabunpaseoct2
pack ('voltoct_cabunpase')
global N l rinp I ra Rm Ri t tau  d Cm 

N = 20;
%l = 400; % in micron
l = 400*( 10^-4); % in cm 
%d = 2*1.85; % in micron
d = 2*(1.85*(10^-4));%cm
r = 1.85*(10^-4); % in cm
Rm = 20000; % ohms.cm^2
Ri = 330; % ohms.cm
Cm = 1; % microfarad / cm^2
%rinp = 1* 10^6; % in ohm
%rinp = (d^(-3/2))*((4*Rm*Ri)^(1/2))*(1/pi);
lambda = ((Rm/Ri)* (d/4))^(1/2); % cm
L= l/lambda;
%rinp = 1* 10^6; % in ohms
%rinp = (d^(-3/2))*((4*Rm*Ri)^(1/2))*(1/pi);
rinf = ((Rm*Ri)^(1/2))*(2/(pi*d^(3/2)));
%rinp = (Ri/pi*(d/2)^2)*lambda*coth(L) % kohms 
ginp = (pi*((d)^(3/2))*tanh(L))/((2*Ri*Rm*10^6)^(1/2));% milliSiemen
rinp = (1/ginp); % kohms
ri = 8.9*(10^7) ; % ohmscm^-1
%t = (1:niter)*deltaT; % in msec
tau = (Rm)*(Cm*(10^-6))*(10^3); %msec  
%tau = 5;% msec
taufactor = (tau)/(Cm*pi*d);
I= 0.1*(10^-9); % amperes

 Va = zeros(niter,1);


[Va,Vo1] = volt_cabunpaseoct2( lambda);

y1 = Va/Vo1;
%disp(y1(2))

% Beginning the section on central difference 
% Define variables
% N -- discretization by space
% V -- output voltage at points x 
% l -- length of dendrite in mm
% lambda -- space constant in mm
% L -- length nondimensionalised - l/lambda
% h -- delta x = l ( N-1)
%h2 -- h*h
% f -- R.H.s of compact diff 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 

%pack ('cabunpanumcent2oct')

% declare global values
global N l Rm Ri t tau I d Cm rinp ra

N = 20;
%l = 400; % in micron
l = 400*( 10^-4); % in cm 
%d = 2*1.85; % in micron
d = 2*(1.85*(10^-4));%cm
r = 1.85*(10^-4); % in cm
Rm = 20000; % ohms.cm^2
Ri = 330; % ohms.cm
Cm = 1; % microfarad / cm^2
%rinp = 1* 10^6; % in ohm
%rinp = (d^(-3/2))*((4*Rm*Ri)^(1/2))*(1/pi);
lambda = ((Rm/Ri)* (d/4))^(1/2); % cm
L= l/lambda;
%rinp = 1* 10^6; % in ohms
%rinp = (d^(-3/2))*((4*Rm*Ri)^(1/2))*(1/pi);
rinf = ((Rm*Ri)^(1/2))*(2/(pi*d^(3/2)));
%rinp = (Ri/pi*(d/2)^2)*lambda*coth(L) % kohms 
ginp = (pi*((d)^(3/2))*tanh(L))/((2*Ri*Rm*10^6)^(1/2));% milliSiemen
rinp = (1/ginp); % kohms
ri = 8.9*(10^7) ; % ohmscm^-1
%t = (1:niter)*deltaT; % in msec
tau = (Rm*(Cm*(10^-6)))*(10^3); %msec  
%tau = 5;% msec
taufactor = (tau)/(Cm*pi*d);
I= 0.1*(10^-9); % amperes

niter = 1; 
Vc = zeros(niter,1);

% calling voltsecentoct_cabunpa4
[Vc,h,h2] = voltsecentoct_cabunpa4(lambda); 

     for iter = 1:niter
% calling function centraldiffcaboct
%[Q,h2] = centraldiffcaboct(V,h,lambda);
% calling function centraldiffcab4oct
[Qc,h2]= centraldiffcab4oct(Vc,h,lambda);
% calling function timedepcentvoltoct
gamconst = 1;
[PP,SS] = timedepcentvoltoct(Qc,Vc,h2,gamconst);
for jj = 2: N-1
Vc(jj) = SS(jj);
endfor 

%V(1)= V(2) + rinp*I*h;
Vc(1) = Vc(2)+(4*Ri*lambda*I*h)/(pi*(d^2));
%V(1) = V(2)+(ri*I*h*lambda*(10^3));
Vc(N) = (4/3)*Vc(N-1)-(1/3)*Vc(N-2);
 nnnn(iter,:) = Vc;
endfor
 Vco2 = nnnn(iter,1); 
y2 = Vc/Vco2;
%disp(Vc(2))
disp(y)
disp(y1)
disp(y2)
e1 = y1-y;
e2 =y1-y2;
e1a = y1(2)-y(2);
er1 = (e1a) / (y1(2));
ep1 = 100*er1;
%disp(e1)
%disp(e2)
%disp(e1a)
%disp(er1)
%disp(ep1)
e2a = y1(2)-y2(2);
er2 = (e2a)/(y1(2));
ep2 = 100*er2;
%disp (e2a)
%disp(er2)
%disp(ep2)

hold on 

x = linspace (0,l,N);
X = x/lambda;
save -mat-binary cabunpanumcent2oct4N20 X y y1 y2 
%save -mat-binary cabunpanumcent2oct4er10 X e1 e2 e1a er1 ep1 e2a er2 ep2

set (findall (gcf, "property", "linewidth"),"linewidth",2)
set (findall (gcf, "property", "markersize"),"markersize",20)
set(gca(),"fontsize",20)

plot (X,y,'k-*',"markersize",6) 
plot (X,y1,'k-s',"linewidth",2,"markersize",6)
plot (X,y2,'k-d',"linewidth",2)
%plot (X,e1,'k-*',"markersize",20)
%plot(X,e2,'k-s',"markersize",20)
xlabel('X(lambda)',"fontsize",20)
ylabel('V/Vo',"fontsize",20)
%ylabel('error',"fontsize",20)
title ('b')
legend ("on")
legend ('{\fontsize{20}compact diff}', '{\fontsize{20} analytical}','{\fontsize{20}central diff}') 
%legend('{\fontsize{20}err compact}','{\fontsize{20}err central}')
print ( 'cabunpanumcent2oct4pubN20.eps', '-depsc')
%print ('cabunpanumcent2octer4pubN10.eps','-depsc')
%print cabunpanumcent2octpubN10.pdf
hold off