%% Figure 3. 
% Weak and Intermediate polarization 
% With Linear Fit. 
% Assumes four panels each with a different _pair_ of M and E_{K}
% Default: Ek={-25,-45} mV X M = {0.3,0.8} \muA/(cm^{2} s)
% Ramp current injection protocol
% RIR July 5, 2015

%% Dependencies
% # IniPR_db
% # NumerEquilPR_db
% # SingIntegODE23PRWithInjCurr_db

%% User Supplied Inputs
% # Array of V_{ds}^{out} to plot
%   MaxVdsOut, DeltaVdsOut, NumVdsOut 
% 


 
MaxVdsOut=10; %Default: 10 mV
DeltaVdsOut=0.5; %Default: 0.5 mV
NumVdsOut=3;   %Default: 51
Ms=0.8;
Eks=-25;
guessVsVd=[0,0];
SomaInj=true;
%VsThresh=30;
VsThresh=10;
%delay =500; 
delay=40; %prior to this model neuron should be at rest with VdsOut. Provides a good sanity check clear before and after, just remember to subtract delay from spike time to get TTFS
Tend=17500;
termaftevent=false;
Isinj=-0.5;
smp=[15,47];
f1=figure();
for i=1:size(Ms,2)
    tmpM=Ms(i);
    tmpuAmpsPermsecCm2=tmpM*1e-3;
    for j=1:size(Eks,2)
       tmpEk=Eks(j); 
       tmpPR=IniPR_db(Isinj,tmpEk);
       for k=1:NumVdsOut
           tmpVdsOut=MaxVdsOut-(k-1)*DeltaVdsOut;
           [tmpnumSSPR,diffProjFullEq,Jacob,eigJacob,nzeig] = NumerEquilPR_db(tmpPR,guessVsVd,tmpVdsOut); 
           tmpPR.SS.NumSS=tmpnumSSPR;
           tmpSingPR = SingIntegODE23PRWithInjCurr_db(tmpPR,tmpuAmpsPermsecCm2,delay,Tend,tmpVdsOut,VsThresh,SomaInj,termaftevent);
           subsample=100;
           extr=1:subsample:size(tmpSingPR.T,1);
           VsVdTrace(i,j).T(k).T=tmpSingPR.T(extr);
           VsVdTrace(i,j).Vs(k).Vs=tmpSingPR.YMultcol(extr,1);
           VsVdTrace(i,j).Vd(k).Vd=tmpSingPR.YMultcol(extr,2);
           VsVdTrace(i,j).VdsOut(k)=tmpVdsOut;
           VsVdTrace(i,j).TTFS(k)=tmpSingPR.te(1)-tmpSingPR.delay;
           VsVdTrace(i,j).NumSpks(k)=size(tmpSingPR.te,1);
           VsVdTrace(i,j).FirstToLastSpike(k)=tmpSingPR.te(size(tmpSingPR.te,1))-tmpSingPR.te(1);
           VsVdTrace(i,j).M=tmpM;
           VsVdTrace(i,j).Ek=tmpEk;
           k
       end
%        subplot(2,1,1)
%        def=linspace(1,NumVdsOut,NumVdsOut);
%  %      plot(VsVdTrace(i,j).VdsOut(nosubidx),VsVdTrace(i,j).TTFS(nosubidx),'ok','MarkerFaceColor','k','MarkerSize',2)
%        plot(VsVdTrace(i,j).VdsOut(extr(1,:)),VsVdTrace(i,j).TTFS(extr(1,:)),'ok','MarkerFaceColor','k','MarkerSize',2)
%        title({['m = ', num2str(tmpM),' \muA/(cm^{2} s)',' Ek= ',num2str(tmpEk),' mV',' VsThresh = ',num2str(VsThresh),' mV']},'FontSize',10)
%        subplot(2,1,2)
%        if mod(k,10)==0
%             plot(VsVdTrace(i,j).T(smp(1)).T,VsVdTrace(i,j).Vs(smp(1)).Vs,'ok','MarkerFaceColor','w','MarkerSize',2,'MarkerEdgeColor','w')
%             hold on;
%        end
%        ax3=gca();
%        set(ax3,'FontSize',8)
    end
end

f2=figure();

subplot(2,1,1)
plot(VsVdTrace(1,1).VdsOut(1,:),VsVdTrace(1,1).NumSpks(1,:),'sk')
title({['m = ', num2str(tmpM),' \muA/(cm^{2} s)',' Ek= ',num2str(tmpEk),' mV',' VsThresh = ',num2str(VsThresh),' mV']},'FontSize',10)
subplot(2,1,2)
plot(VsVdTrace(1,1).VdsOut(1,:),VsVdTrace(1,1).FirstToLastSpike(1,:),'-dk')

f3=figure();
title({['m = ', num2str(tmpM),' \muA/(cm^{2} s)',' Ek= ',num2str(tmpEk),' mV',' VsThresh = ',num2str(VsThresh),' mV']},'FontSize',10)
for k=1:size(VsVdTrace(1,1).VdsOut,2)
    if VsVdTrace(1,1).NumSpks(1,k)==0
        sym='ok';
        col='k';
    elseif VsVdTrace(1,1).NumSpks(1,k)==1
        sym='sk';
        col='None';
    elseif VsVdTrace(1,1).NumSpks(1,k)==2
        sym='^k';
        col='k';
    elseif VsVdTrace(1,1).NumSpks(1,k)==3
        sym='dk';
        col='None';
    elseif VsVdTrace(1,1).NumSpks(1,k) < 10 && VsVdTrace(1,1).NumSpks(1,k) > 3
        sym='+k';
        col='k';
    else
        sym='*k';
        col='None';     
    end
    plot(VsVdTrace(1,1).VdsOut(1,k),VsVdTrace(1,1).TTFS(1,k),sym,'MarkerFaceColor',col,'MarkerSize',6)
    hold on;
end