function [t_EPSP,v_EPSP,nSyn,mean_EPSP,dat,PSCbase,peak,amp,rise,decay,hfw] = read_EPSCsims_mdb(fbase, gAMPA);
%
%   read_EPSCsims_mdb
% 
%   Reads binary files generated by EPSC simulations from the Amatrudo,
%   Weaver, et al. (2012) ModelDB entry.
%
%     Sample usage:
% 
%     [t_EPSP,v_EPSP,nSyn,mean_EPSP,dat,PSCbase,peak,amp,rise,decay,hfw] = ...
%               read_EPSCsims_bin('fig11_V1lowAMPA_apic',0.00017);
% 
%           and
% 
%     [t_EPSP,v_EPSP,nSyn,mean_EPSP,dat,PSCbase,peak,amp,rise,decay,hfw] = ...
%               read_EPSCsims_bin('fig11_V1sameAMPA_bas',0.00054);
%
%
%   written by Christina Weaver, christina.weaver@fandm.edu, May 2012
%

if( ~exist('fbase','var') ) fbase = 'PFC';  end;

SYN_START = 200;
SYN_DELAY = 200;
tau1 = 0.15;
tau2 = 4.0;
writeSmry = 0;

inbase = sprintf('%s_tR%.4f_tF%.2f_gAMP%.5f',fbase,tau1,tau2,gAMPA);

figttl = sprintf('PFC test: t_1 %.2f t_2 %.2f gAMP %.3f',tau1,tau2,gAMPA);
makeDistPlot = 1;
TipsOnly = 1;


[t,v]=readNRNbin_Vclamp(inbase,0);

txt_fname = sprintf('%s_dist.txt',inbase);
[dat] = dlmread(txt_fname);
nSyn = size(dat,1);

sTimes = SYN_START : SYN_DELAY : t(end);
sTimes = sTimes(1:nSyn);
idx = [min(find(t >= sTimes(1)))  max(find(t < sTimes(2)))];
tstep = [0 : idx(2)-idx(1)-1];
t_EPSP = t(1+tstep);

idxBase = find(t<sTimes(1) & t>=sTimes(1)-10);
PSCbase = mean(v(idxBase));

v_EPSP = [];
peak = [];
amp = [];
rise = [];
decay = [];
hfw = [];

for i = 1 : nSyn
    idx = min(find(t >= sTimes(i)));
    if( ~isempty(find(idx+tstep > length(v))) )
        %         We've moved past the end of the synapses.
        fprintf('No synapse for %d\n',i);
        continue;
    end;
    tmp = v(idx+tstep)';
    if( max(abs(diff(tmp))) == 0 )
        fprintf('SKIP %d:  No voltage change for proposed synapse.\n',i);
        continue;
    end;
    
    [mn,mnI]=min(tmp);
    
    % start the trace at zero
    v_EPSP(1+tstep,end+1) = tmp-tmp(1);
    peak(end+1,:) = [t_EPSP(mnI)  abs(mn-tmp(1))];
    
    [amp(end+1),rise(end+1),decay(end+1),hfw(end+1)]=analyze_EPSC(t_EPSP,tmp);
end;

nSyn = size(v_EPSP,2);
fprintf('Found %d synapses\n',nSyn);
fprintf('Mean amp\t%.2f\n',mean(amp)*1e3);
fprintf('Mean rise\t%.2f\n',mean(rise));
fprintf('Mean decay\t%.2f\n',mean(decay));
fprintf('Mean 1/2 width\t%.2f\n',mean(hfw));

mean_EPSP = mean(v_EPSP,2);
figure;
subplot(1,3,1);
plot(t_EPSP,v_EPSP,'k',t_EPSP,mean_EPSP,'r');

ttl = sprintf('%s:  superimposing all EPSCs',figttl);
title(ttl);
xlabel('time (ms)');
ylabel('Postsynaptic conductance (nA)');
fname = sprintf('%s_meanEPSCond.fig',inbase);

subplot(1,3,2);
plot(dat(:,2),peak(:,2),'.');
ttl = sprintf('%s:  peak I change');
title(ttl);
xlabel('Traversed distance from soma (microns)');
ylabel('Max EPSC amplitude');


subplot(1,3,3);
xbar=[.0005:.001:.060];
n_elements=histc(peak(:,2),xbar);
npts=max(size(peak));
c_elements = cumsum(n_elements);

plot(xbar,c_elements/npts,'o-');
xlim([0 0.06])
title(fbase);

saveas(gcf,fname);

if( writeSmry )
    fname = sprintf('%s_smry.txt',inbase);
    fout=fopen(fname,'w');
    fprintf(fout,'distance\tamp\trise\tdecay\thalfwidth\n');
    for i=1:size(dat,1)
        fprintf(fout,'%g\t%g\t%g\t%g\t%g\n',dat(i,2),1e3*amp(i),rise(i),decay(i),hfw(i));
    end;
    fclose(fout);
end;


return;