function data =  plot_bs(bs_scanfile)
%data =  plot_bs(bs_scanfile)
% 
% bs_scafile is an output file which results from running scanbs.py 
%   on a scan file: the default of the python program is to just append
%  ".bs" to the original scanfile name
%
%
% "data" is a structure with


dd = read_asc(bs_scanfile);

% 1     2     3    4      5     6      7        8       9          10
%[y startdd enddd beam    N   MEAN  STD_DEV     SEM   minimum    maximum]
num_ens = length(dd(:,1))/4;

year       = dd(:,1);
startdd    = dd(:,2);
enddd      = dd(:,3);

start_date = to_date(year, startdd, 'n');
end_date   = to_date(year, enddd ,'n');


%reshape:

data.startdd = reshape(to_day(year(1), start_date), 4, num_ens)';
data.enddd   = reshape(to_day(year(1), end_date), 4, num_ens)';
data.beam       = reshape(dd(:,4), 4, num_ens)';
data.num        = reshape(dd(:,5), 4, num_ens)';
data.meanamp    = reshape(dd(:,6), 4, num_ens)';
data.stddev     = reshape(dd(:,7), 4, num_ens)';
data.stderrm    = reshape(dd(:,8), 4, num_ens)';
data.minamp     = reshape(dd(:,9), 4, num_ens)';
data.maxamp     = reshape(dd(:,10),4, num_ens)';

figure
dt = [data.enddd(:,1)-data.startdd(:,1)];
t1 = data.startdd(:,1) + 0*dt;
t2 = data.startdd(:,1) + .1*dt;
t3 = data.startdd(:,1) + .2*dt;
t4 = data.startdd(:,1) + .3*dt;

plot([t1 t1]',  [data.minamp(:,1) data.maxamp(:,1)]','r')
hold on
h=plot(t1, data.meanamp(:,1), 'r*');set(h,'markersize', 10)

plot([t2 t2]',  [data.minamp(:,2) data.maxamp(:,2)]','g')
h=plot(t2, data.meanamp(:,2), 'g*');set(h,'markersize', 10)

plot([t3 t3]',  [data.minamp(:,3) data.maxamp(:,3)]','b')
h=plot(t3, data.meanamp(:,3), 'b*');set(h,'markersize', 10)

plot([t4 t4]',  [data.minamp(:,4) data.maxamp(:,4)]','k')
h=plot(t4, data.meanamp(:,4), 'k*');set(h,'markersize', 10)


title('beam statistics: [min, mean, max] rgbk = 1,2,3,4')
xlabel('dday')

orient portrait
print -dpsc beamstats.ps
