
%% Load all data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bands =  [10, 20; 20, 30; 30, 40; 43, 58; 70, 100; 140, 160; 200, 220; 240, 260; 305, 370; ...
         340, 360; 400, 450; 700, 900];
files = dir('LHO/lho_blrms*.txt');
nfiles = length(files);

b = [];

for i=1:nfiles
    % get GPS time from file name
    x = strsplit(files(i).name, {'_', '.'});
    gps0 = str2double(x(4));
    % load data
    bb = load(['LHO/', files(i).name]);
    if ~isempty(bb)
        % first column is time from start, add initial GPS time
        bb(:,1) = bb(:,1) + gps0;
        b = [b; bb];
    end
    fprintf(2, 'Read file %s (%d)\n', files(i).name, gps0);
end

% time in GPS and day from 2015-09-19 00:00 UTC
t = b(:,1);
T = (t - 1126569617)/86400 + datenum('09-19-2015');
b = b(:,2:end);
Tticks = floor(min(T)):ceil(max(T));


%% plot all trends %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ymax = [5e-21, 5e-22, 1e-19,   1e-20,   3e-19,   1e-18,   2e-18,   2e-18,   2e-16,   2e-17, 1e-17,   1e-16];
ymin = [0e-21, 0e-22, 0.5e-19, 0.1e-20, 0.4e-19, 0.2e-18, 0.5e-18, 0.8e-18, 0.8e-16, 0e-17, 0.5e-17, 0.5e-16];

for i=1:size(b,2)
    bsmooth = running_average(t, b(:,i), 60, 60);
    figure('Position', [100 100 1500 800])
    plot(T, b(:,i), '.')
    hold all
    plot(T, bsmooth, '.')
    xlabel('Date [at 00:00 UTC]')
    datetick('x')
    set(gca, 'xtick', Tticks)
    ax = gca;
    ax.XTickLabelRotation = -45;
    ylabel('BLRMS [a.u.]')
    title(sprintf('BLRMS  %d-%d Hz', bands(i,1), bands(i,2)))
    legend('60s average', '1 hour average')
    xlim([min(T), max(T)])
    ylim([ymin(i) ymax(i)])
    hgexport(gcf, sprintf('lho_trend_%d_%d.png', bands(i,1), bands(i,2)),  hgexport('factorystyle'), 'Format', 'png');
end

