%
%  Gabbiani & Cox, Mathematics for Neuroscientists
%
%   cabextrainv.m
%
%  recover cable and stimulus parameters from knowledge of extracelluar
%  potentials from      trapcabextraphi.mat
%  generated by running trapcabextra.m
%
%   usage:  cabextrainv
%
function cabextrainv

close all
load trapcabextraphi

for j=1:9,
   Phi0(j,1) = trapz(tim',squeeze(phi(1,j,:)));       % zeroth moment of data
   Phi1(j,1) = trapz(tim',tim'.*squeeze(phi(1,j,:))); % first moment of data
end
    
[~,loc] = max(abs(Phi0));
Phi0rat = Phi0/Phi0(loc);

for j=1:6
    for k=1:99,
        
        y0 = k*0.001;
        lam = 0 + j*0.015;
        val = PhiMom([y0 lam],xhat(1),yhat,ell,Phi0rat,loc);
        E(k) = norm(val)^2;
        
    end
    
    plot(0.001*[1:99],E,'k')
    hold on
    
end

ylim([0 100])
xlabel('y_0  (cm)','fontsize',14)
ylabel('E(y_0,0.015m), m=1,...,6','fontsize',14)
text(0.09,90,'(A)','fontsize',20)
hold off

% now fix y0 and graph E wrt lam

for k=1:99,
    
    val = PhiMom([0.06 k*0.001],xhat(1),yhat,ell,Phi0rat,loc);
    E(k) = norm(val)^2;
    
end

figure(2)
plot(0.001*[1:99],E,'k')

ylim([0 0.1])
xlabel('\lambda  (cm)','fontsize',14)
ylabel('E(0.06,\lambda)','fontsize',14)
text(0.09,0.09,'(B)','fontsize',20)

% now solve for x = [tau; I1], 

lam = trueX(2);
J0 = trueX(1);
y0 = trueX(3);
X = trueX;

L = ell/lam;
L0 = y0/lam;
a1 = cosh(L-L0); a2 = cosh(L0);
b1 = (J0/2)*(L0*sinh(L-L0) + L*cosh(L0)/sinh(L) - cosh(L-L0));
b2 = (J0/2)*(L*cosh(L)*cosh(L0)/sinh(L) - cosh(L0) - L0*sinh(L0));
g1 = -(J0/2)*cosh(L-L0)/lam;
g2 = -(J0/2)*cosh(L0)/lam;

for i=1:9
    S1 = integral(@(y) ig1(y,X,xhat(1),yhat(i)),0,y0);
    S2 = integral(@(y) ig2(y,X,ell,xhat(1),yhat(i)),y0,ell);
    S3 = integral(@(y) ig3(y,X,xhat(1),yhat(i)),0,y0);
    S4 = integral(@(y) ig4(y,X,ell,xhat(1),yhat(i)),y0,ell);
    A(i,1) = a1*S1 + a2*S2 - lam*sinh(L)/sqrt(xhat(1)^2+(y0-yhat(i))^2);
    A(i,2) = b1*S1 + g1*S3 + b2*S2 + g2*S4;
end

x = A\Phi1;

J1 = x(1);
I1 = 4*pi*0.003*lam*sinh(L)*J1;
tau = x(2);

mc = A(:,2)./A(:,1);
rhs = Phi1./A(:,1);
figure(3)
plot(mc,rhs,'kx')   % raw data
hold on
LSA = ones(9,2);   % find best line fit through raw data
LSA(:,2) = mc;
x = LSA\rhs;
eta0 = -0.015;
eta1 = 0.02;
plot([eta0 eta1],x(1)+x(2)*[eta0 eta1],'k')
hold off
xlabel('\eta','fontsize',14)
ylabel('\zeta','fontsize',14)

return

function val = PhiMom(X,xhat,yhat,ell,Phi0rat,loc)
val = zeros(9,1);
for i=1:9,
    int1 = integral(@(y) ig1(y,X,xhat,yhat(i)),0,X(1));
    int2 = integral(@(y) ig2(y,X,ell,xhat,yhat(i)),X(1),ell);
    val(i) = (cosh((ell-X(1))/X(2))*int1 + cosh(X(1)/X(2))*int2 - ...
           X(2)*sinh(ell/X(2))/sqrt(xhat^2+(X(1)-yhat(i))^2));
end
val = val/val(loc);
val = val - Phi0rat;
return

function val = ig1(y,X,xhat,yhat)
val = cosh(y/X(2))./sqrt(xhat^2+(yhat-y).^2);
return

function val = ig2(y,X,ell,xhat,yhat)
val = cosh((ell-y)/X(2))./sqrt(xhat^2+(yhat-y).^2);
return

function val = ig3(y,X,xhat,yhat)
val = y.*sinh(y/X(2))./sqrt(xhat^2+(yhat-y).^2);
return

function val = ig4(y,X,ell,xhat,yhat)
val = (ell-y).*sinh((ell-y)/X(2))./sqrt(xhat^2+(yhat-y).^2);
return


    
