%
% compute the potential in a cell during
% somatic current injection
%
% usage
%
%  [indexes,len] = pots(neurolucida_file,param_file,out_file)
%

function [indexes,len,t,v0] = pots(neurolucida_file,param_file,out_file)

close all

eval(param_file)

gamma = 2*Ra*(Gm-Gs);

[As,indexes,list,len,lpoint] = convert(neurolucida_file);

As = (4/3)*(1e-8)*As;	% cm

[M,K] = makeMK2(h,list,indexes,len,lpoint);

%return

M = (1e-8)*M;
K = (1e-4)*K;

%save MaK M

M(1,1) = M(1,1) + As;
K(1,1) = K(1,1) - As*gamma;

n = size(M,2);

B = 2*Ra*Cm*M;

A = B + dt*K + 2*Ra*Gm*dt*M;
disp('A built')
[L,U] = lu(A);
disp('A factored')

v = zeros(n,1);

fid = fopen(out_file,'w');
fwrite(fid,dt,'double');
fwrite(fid,Nt,'double');
fwrite(fid,n,'double');
fwrite(fid,v,'double');

for j=1:Nt,

        rhs = B*v;
        t(j) = j*dt;
        rhs(1) = rhs(1) + Ra*dt*stim(t(j))/pi;
        v = U \ (L \rhs);
        v0(j) = v(1);
	fwrite(fid,v,'double');
        fprintf(1,'%d',j);

end

fclose(fid);