function [t_EPSP,v_EPSP,nSyn,mean_EPSP,dat,PSCbase,peak,amp,rise,decay,hfw] = read_IPSCsims_mdb(fbase,gGABA);
%
%   read_IPSCsims_mdb
% 
%   Reads binary files generated by IPSC 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_IPSCsims_bin('fig12_V1lowAMPA_apic',0.00067);
% 
%           and
% 
%     [t_EPSP,v_EPSP,nSyn,mean_EPSP,dat,PSCbase,peak,amp,rise,decay,hfw] = ...
%               read_IPSCsims_bin('fig12_V1sameAMPA_bas',0.00115);
%
%
%   written by Christina Weaver, christina.weaver@fandm.edu, May 2012
%

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


SYN_START = 200;
SYN_DELAY = 200;
tau1 = 2.5;
tau2 = 7.5;
writeSmry = 0;

inbase = sprintf('%s_INtR%.4f_tF%.2f_gGAB%.5f',fbase,tau1,tau2,gGABA);

figttl = sprintf('PFC test: t_1 %.2f t_2 %.2f gAMP %.3f',tau1,tau2,gGABA);
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;

    tmp = tmp*-1;
    [mn,mnI]=min(tmp);
    
    % corrected 4/10/12:  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,-1*v_EPSP,'k',t_EPSP,-1*mean_EPSP,'r');
% plot(t_EPSP,v_EPSP-PSCbase,'k',t_EPSP,mean_EPSP-PSCbase,'r');
xlim([0 50]);
ttl = sprintf('%s:  superimposing all IPSCs',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 IPSC amplitude');

subplot(1,3,3);
xbar=[.0005:.001:.080];
n_elements=histc(peak(:,2),xbar);
npts=max(size(peak));
c_elements = cumsum(n_elements);
% c_elements = histc(peak(:,2),xbar);
plot(xbar,c_elements/npts,'o-');
xlim([0 0.08])
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);
    % fname = sprintf('%s_amplitudes.txt',inbase);
    % fout=fopen(fname,'w');
    % for i=1:size(dat,1)
    %     fprintf(fout,'%g\t%g\n',dat(i,2),peak(i,2));
    % end;
    % fclose(fout);
end;


return;