%
% recover cable properties
%
% compute the potential in a cell during
% somatic current injection
%

clear all
close all

dt = .1;   	     % time step, ms
Nt = 200;	     % number of time steps
Ra = .3;  	     % k Ohm micron
Cm = 1;       % micro F / micron^2
Gm = 0.3;    % mS / micron^2
q = 2*Ra*Gm;
Gs = 2*Gm;
gamma = 2*Ra*(Gm-Gs);
drec = 180;		% site of dendritic recording

[As,indexes,list,len,lpoint] = convert('antena.asc');

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

h = 1;			% space step, microns

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

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

%save MaK M K

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;

v = zeros(n,1);

disp(' ')
disp('Solving the cable equations')
disp('on iteration ')

for j=1:Nt,

        rhs = B*v;
        t(j) = j*dt;
        stim(j) = t(j)*exp(-t(j))/2000;
        rhs(1) = rhs(1) + Ra*dt*stim(j)/pi;
        v = A\rhs;
        v1(j) = v(1);
        v3(j) = v(drec);
        fprintf(1,'%d',j);

end

fprintf(1,'\n');

figure(3)
plot(t,v1)
hold on
plot(t,v3,'--')
xlabel('t (ms)','fontsize',16)
ylabel('v (mV)','fontsize',16)
hold off


mv1 = sum(v1)*dt;
mv3 = sum(v3)*dt;
Tlevel = mv3/mv1;

disp(['Tlevel = ' num2str(Tlevel)])

K(1,1) = K(1,1) + As*gamma;	% return K to native
opts.disp = 0;
[Q,D] = eigs(K,M,80,-1,opts);

e = diag(Q'*M*Q);
Qe = fliplr(Q*diag(1./sqrt(e)));
sn2 = Qe(1,:);
pn2 = sn2.^2;
pN2 = Qe(drec,:);
mun = fliplr(diag(D)');

disp(['mun(1:4) = ' num2str(mun(1:4))])
figure(4)
plot(Qe(:,1:4))

q = 0.12:0.01:0.24;
lq = length(q);

for j=1:lq,

    Tnum = sum( sn2.*pN2./(mun+q(j)) );
    Tden = sum( pn2./(mun+q(j)) );
    T(j) = Tnum/Tden;

end

figure(5)
plot(q,T)
hold on
plot(q,Tlevel*ones(1,lq))
xlabel('q (cm)^{-1}','fontsize',16)
ylabel('H','fontsize',16)

% now for dT and such

qq = 0.18;   % exact value

dTnum = sum(pn2./(mun+qq) )*sum( sn2.*pN2./((mun+qq).^2) ) -...
        sum(pn2./((mun+qq).^2))*sum( sn2.*pN2./(mun+qq) );
dTden = sum( pn2./(mun+qq) )^2;
dT = dTnum/dTden;

m1v1 = sum(t.*v1)*dt;
m1v3 = sum(t.*v3)*dt;
dTlevel = (mv1*m1v3 - m1v1*mv3)/mv1^2;
twiceRaCm = dTlevel/dT;

% now Phi to get gamma

io = stim;

K11 = K(1,1);
gc = 0;
Philevel = sum(t.*io)/sum(io) - sum(t.*v1)/sum(v1);

for gamma = -0.2:0.01:-0.16,

    K(1,1) = K11 - gamma*As;
	[Q,D] = eigs(K,M,80,'sm',opts);
	e = diag(Q'*M*Q);
	Qe = Q*diag(1./sqrt(e));
	Qe = fliplr(Qe);
	sn2 = Qe(1,:);
	pn2 = sn2.^2;
	mun = fliplr((diag(D))');
	gc = gc + 1;
	Phi(gc) = -twiceRaCm*sum(pn2./((mun+qq).^2))/sum(pn2./(mun+qq));
	gg(gc) = gamma;
end

figure(3)
plot(gg,Phi)
hold on
plot(gg,Philevel*ones(size(gg)))
%xlabel('\gamma    (cm)^{-1}','fontsize',16)
%ylabel('P','fontsize',16)

gamma = -.18;

K(1,1) = K11 - gamma*As;
[Q,D] = eigs(K,M,80,'sm',opts);
e = diag(Q'*M*Q);
Qe = Q*diag(1./sqrt(e));
Qe = fliplr(Qe);
sn2 = Qe(1,:);
pn2 = sn2.^2;
mun = (diag(D))';
mun = fliplr(mun)+qq;
    
Ra = pi*(sum(v1)/sum(io))/sum(pn2./mun)
Gm = qq/2/Ra
Cm = twiceRaCm/2/Ra
Gs = Gm - gamma/2/Ra

