% Script file : cabunpa2sevoltNoct2.m
% Purpose : To study V/Vo vs X for soma- dendrite,passive, sealed end.Correction for current inj cond. For various values of N
% Current Injection at x = 0;
% Record of revisions :
% Date                                         Programmer                         Description
% June 1, 2007                                  Asha G                              original code
% June 14, 2007                                 Asha G                              removed II  
% August 31, 2008                              Asha G                              correction for current inj condn
 % October 26, 2010                             Asha G                              For various values of N 
% October 29,2010                               Asha G                              tau changed and I changed to nanoamp
% Nov 5,2010                                   Asha G                             added linewidth and markersize and printpdf
% Nov 12,2010                                  Asha G                             Multiply I by taufactor 
% 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 

% calling function voltseoct_cabunpa4

clear all 
pack('cabunpa2sevoltNoct2')

% declare global values
global N l Rm Ri  tau I d Cm ra rinp Cfactor d Iinj ri
N = 10;
l = 400*(10^-4);%in cm
%d = 10*(10^-4);% in cm
d= 2*1.85*(10^-4);% in cm
%r= 1.85*(10^-4);% in cm
Rm = 20000;% ohms.cm^2
%Rm = 2500*(10^-3);%kohms.cm^2
Ri = 330;% ohms.cm
%Ri = 70*(10^-3);%kohms.cm
Cm = 1;% microfarad/cm^2
%Cm = 2;% microfarad/cm^2
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) % ohms 
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
Iinj = I* lambda; % microamperes/cm

niter = 8000;
%niter = 200;

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


for n = 1: niter

[f]= compactdiffcaboct(V,h2);
%[f]=compactdiffcab4oct(V,h2);

% calling function tridiagdoubleoct

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

%[Q,M]= tridiagdoubleoct(alpha,alpha1,f);
% calling function tridiagoct
[ Q,Qs,M,Ms] = tridiagoct( alpha,alpha1,f);
%disp(M)
%disp(Qs)
%disp(Ms)
% calling function timedepcabvoltoct

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

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

endfor 

  %Icomp = I/pi*d*h;
%V(1) = V(2)+ (lambda*Ri*(10^3)*4*I*h*(10^-6))/(pi*(d)^2);
V(1) = V(2)+(4*Ri*lambda*I*h)/(pi*(d^2));
%V(1) = V(2)+(ri*I*h*lambda*(10^3));
%V(1) = V(2)+ tau* I *h/(Cm*pi*d);
%disp(V(1))
V(N) = (4/3)* V(N-1)-(1/3)* V(N-2);
%V(N) = V(N-1);
%disp(V)

nnnn(n,:)= V;

endfor 
%disp(nnnn)
z = V; 
q = Q;
disp(z)
Vo = nnnn(niter,1);
%disp(Vo)
y = z/Vo;
%disp(y)
%y = z;

% calling function volt_cabunpaseoct2
[Va,Vo1]= volt_cabunpaseoct2(lambda);
disp(Vo1)

y1 = Va/Vo1;
 
hold on 
x = linspace(0,l,N);
X= x/lambda;

save -mat-binary cabunpa2sevoltNoct2N10 X y y1 
set (findall (gcf, "-property", "linewidth"), "linewidth",2)
set (findall (gcf, "-property", "markersize"), "markersize",6)
set(gca(),"fontsize",20)

plot (X,y,'rs',"markersize",6)
plot (X,y1,'k',"linewidth",2)
%plot(X,full(q),'m*',"markersize",6)
%plot(X,z,'gs',"markersize",6)

xlabel('X(lambda)',"fontsize",20)
ylabel ('V/Vo',"fontsize",20)
%ylabel ('V')
legend ("on")
%gset nokey;
title('sealed end N = 10, V/Vo vs X ')
legend ('{\fontsize{14} compact}', '{\fontsize{14} analytical}')
print ('cabunpa2sevoltNoct2N10.eps','-depsc')
%print cabunpa2sevoltNoct2pubN30.pdf
hold off