function modify_segs_list(hdatafile, ldatafile, h1sciencelist, h2sciencelist, l1sciencelist, h1outfile, h2outfile, l1outfile)
%modify_segs_list creates output files that hold mean and sigma for timing
% data for LHO and LLO science segments.
%
% modify_segs_list modifies the standard science segment lists, and produces
% three files titled h1outfile, h2outfile, l1outfile each of which holds the
% matrices H1_science_timing, H2_science_timing, L1_science_timing each of
% who holds ten columns:
%
% Science segment length
% Science segment start GPS time
% Science segment end GPS time
% Science segment serial number
% Mean of EX/MX time-reading for that science segment
% Sigma of EX/MX time-reading for that science segment
% Mean of EY/MY time-reading for that science segment
% Sigma of EY/MY time-reading for that science segment
% Mean of LVEA/PEM time-reading for that science segment
% Sigma of LVEA/PEM time-reading for that science segment
%
% This function disregards invalid time-readings while calculating mean and
% sigma. If it finds no valid time-reading for a channels during a science
% segments, it returns NaN as mean and sigma for that segment.
%
% hdatafile and ldatafile are the outut files produced by HTMon_data and
% LTMon_data.
%
% h1sciencelist, h2sciencelist, and l1sciencelist should contain the two dimensional
% matrices titled H1_science_segments, H2_science_segments, and L1_science_segments each of which
% should hold four columns:
%
% Science segment length
% Science segment start GPS time
% Science segment end GPS time
% Science segment serial number
%
% Example: modify_segs_list('HTimeMon_all.mat', 'LTimeMon_all.mat', 'H1_science_segments.list', 'H2_science_segments.list', 'L1_science_segments.list', 'H1_science_timing', 'H2_science_timing', 'L1_science_timing')
%
% See also: HTMon_data, LTMon_data
%
% Rubab Khan, GECo, Columbia
% https://geco.phys.columbia.edu/~rubab/
% December 2, 2007.

% Load the data files

load(hdatafile);       % Load LHO time data
load(ldatafile);       % load LLO time data

a = H_time';                    % LHO time channel
b = L_time';                    % LLO time channel
c = [];
d = [];
e = [];

% Load science segment lists

load(h1sciencelist);
st1 = H1_science_segments(:,2); % Science segment start times
en1 = H1_science_segments(:,3); % Science segment end times
A = length(st1);                % Total number of science segments for this detector
for i = 1:A                     % For each science segment
    cc=find(a>=st1(i)&a<=en1(i));
    c = [c cc];                 % Record the time when there was valid time-reading
end

% Repeat the steps for H2 and L1

load(h2sciencelist);
st2 = H2_science_segments(:,2);
en2 = H2_science_segments(:,3);
B = length(st2);
for i = 1:B
    dd=find(a>=st2(i)&a<=en2(i));
    d = [d dd];
end

load(l1sciencelist);
st3 = L1_science_segments(:,2);
en3 = L1_science_segments(:,3);
C = length(st3);
for i = 1:C
    ee=find(b>=st3(i)&b<=en3(i));
    e = [e ee];
end

% Calculate and record means and sigmas

% For H1EX
f = [a' H1EX_mean];             % f consists GPS time and time-reading for this channel
g = f(c,1); h = f(c,2);         % g and h consists science segment GPS time and time-reading for this channel
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);             % Now g and h consists valid (corresponding GPS time and) time-reading only
k1 = [];
l1 = [];
for i = 1:A                      % For each science segment
    m=find(g>=st1(i)&g<=en1(i)); % Find valid time-reading during a science segment
    n = length(m);               % Number of valid  readings
    if n == 0
        kk = NaN; ll = NaN;          % If no valid reading, mean = 0 and sigma = 0
    else
        kk = mean(h(m)); ll = std(h(m)); % Otherwise fine mean and sigma
    end
    k1 = [k1 kk]; l1 = [l1 ll];
end
k1 = k1'; l1 = l1';              %Record means in vector K1, sigmas in vector L1


% Repeat same steps for other channels

% For H1EY
f = [a' H1EY_mean];
g = f(c,1); h = f(c,2);
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);
k2 = [];
l2 = [];
for i = 1:A
    m=find(g>=st1(i)&g<=en1(i));
    n = length(m);
    if n == 0
        kk = NaN; ll = NaN;
    else
        kk = mean(h(m)); ll = std(h(m));
    end
    k2 = [k2 kk]; l2 = [l2 ll];
end
k2 = k2'; l2 = l2';

% For H1LV
f = [a' H1LV_mean];
g = f(c,1); h = f(c,2);
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);
k3 = [];
l3 = [];
for i = 1:A
    m=find(g>=st1(i)&g<=en1(i));
    n = length(m);
    if n == 0
        kk = NaN; ll = NaN;
    else
        kk = mean(h(m)); ll = std(h(m));
    end
    k3 = [k3 kk]; l3 = [l3 ll];
end
k3 = k3'; l3 = l3';

% For H2MX
f = [a' H2MX_mean];
g = f(d,1); h = f(d,2);
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);
k4 = [];
l4 = [];
for i = 1:B
    m=find(g>=st2(i)&g<=en2(i));
    n = length(m);
    if n == 0
        kk = NaN; ll = NaN;
    else
        kk = mean(h(m)); ll = std(h(m));
    end
    k4 = [k4 kk]; l4 = [l4 ll];
end
k4 = k4'; l4 = l4';

% For H2MY
f = [a' H2MY_mean];
g = f(d,1); h = f(d,2);
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);
k5 = [];
l5 = [];
for i = 1:B
    m=find(g>=st2(i)&g<=en2(i));
    n = length(m);
    if n == 0
        kk = NaN; ll = NaN;
    else
        kk = mean(h(m)); ll = std(h(m));
    end
    k5 = [k5 kk]; l5 = [l5 ll];
end
k5 = k5'; l5 = l5';

% For H2LV
f = [a' H2LV_mean];
g = f(d,1); h = f(d,2);
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);
k6 = [];
l6 = [];
for i = 1:B
    m=find(g>=st2(i)&g<=en2(i));
    n = length(m);
    if n == 0
        kk = NaN; ll = NaN;
    else
        kk = mean(h(m)); ll = std(h(m));
    end
    k6 = [k6 kk]; l6 = [l6 ll];
end
k6 = k6'; l6 = l6';

% For L0EX
f = [b' L0EX_mean];
g = f(e,1); h = f(e,2);
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);
k7 = [];
l7 = [];
for i = 1:C
    m=find(g>=st3(i)&g<=en3(i));
    n = length(m);
    if n == 0
        kk = NaN; ll = NaN;
    else
        kk = mean(h(m)); ll = std(h(m));
    end
    k7 = [k7 kk]; l7 = [l7 ll];
end
k7 = k7'; l7 = l7';

% For L1EY
f = [b' L1EY_mean];
g = f(e,1); h = f(e,2);
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);
k8 = [];
l8 = [];
for i = 1:C
    m=find(g>=st3(i)&g<=en3(i));
    n = length(m);
    if n == 0
        kk = NaN; ll = NaN;
    else
        kk = mean(h(m)); ll = std(h(m));
    end
    k8 = [k8 kk]; l8 = [l8 ll];
end
k8 = k8'; l8 = l8';

% For L1PM
f = [b' L1PM_mean];
g = f(e,1); h = f(e,2);
j=find(h~=0&(abs(h))<1e-4);           % Disregard data when invalid time-reading was read
g = g(j); h = h(j);
k9 = [];
l9 = [];
for i = 1:C
    m=find(g>=st3(i)&g<=en3(i));
    n = length(m);
    if n == 0
        kk = NaN; ll = NaN;
    else
        kk = mean(h(m)); ll = std(h(m));
    end
    k9 = [k9 kk]; l9 = [l9 ll];
end
k9 = k9'; l9 = l9';

% Read first four columns
p = H1_science_segments;
q = H2_science_segments;
r = L1_science_segments;

% Create the ten column matrices
H1_science_timing = [p k1 l1 k2 l2 k3 l3];
H2_science_timing = [q k4 l4 k5 l5 k6 l6];
L1_science_timing = [r k7 l7 k8 l8 k9 l9];

% create output files

fid1 = fopen(h1outfile,'w');
fid2 = fopen(h2outfile,'w');
fid3 = fopen(l1outfile,'w');

for i = 1:length(H1_science_timing(:,1))
    fprintf(fid1,sprintf('%d\t%d\t%d\t%d\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\n', ... 
        H1_science_timing(i,1), H1_science_timing(i,2), ...
        H1_science_timing(i,3), H1_science_timing(i,4), ...
        H1_science_timing(i,5), H1_science_timing(i,6), ...
        H1_science_timing(i,7), H1_science_timing(i,8), ...
        H1_science_timing(i,9), H1_science_timing(i,10)));
end

for i = 1:length(H2_science_timing(:,1))
    fprintf(fid2,sprintf('%d\t%d\t%d\t%d\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\n', ... 
        H2_science_timing(i,1), H2_science_timing(i,2), ...
        H2_science_timing(i,3), H2_science_timing(i,4), ...
        H2_science_timing(i,5), H2_science_timing(i,6), ...
        H2_science_timing(i,7), H2_science_timing(i,8), ...
        H2_science_timing(i,9), H2_science_timing(i,10)));
end

for i = 1:length(L1_science_timing(:,1))
    fprintf(fid3,sprintf('%d\t%d\t%d\t%d\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\t\t%1.12f\n', ... 
        L1_science_timing(i,1), L1_science_timing(i,2), ...
        L1_science_timing(i,3), L1_science_timing(i,4), ...
        L1_science_timing(i,5), L1_science_timing(i,6), ...
        L1_science_timing(i,7), L1_science_timing(i,8), ...
        L1_science_timing(i,9), L1_science_timing(i,10)));
end

fclose(fid1);
fclose(fid2);
fclose(fid3);

h1outfile_mat = ([h1outfile '.mat']);
h2outfile_mat = ([h2outfile '.mat']);
l1outfile_mat = ([l1outfile '.mat']);

save(h1outfile_mat, 'H1_science_timing', '-MAT');
save(h2outfile_mat, 'H2_science_timing', '-MAT');
save(l1outfile_mat, 'L1_science_timing', '-MAT');