function [st,brk,dl]=GuildTwoIfoStats(filename1,filename2,outname)
%[st,brk,dl]=GuildTwoIfoStats(filename1,filename2)
%LockStatistics infered from minute trends from guild queries of two ifos
%stored in "filename1" and "filename2"
%The file should be an ascii file and the name should have no extensions.
%
%st and brk are start and break times in minutes from the starting time.
%dl is the length of each segment. 

eval(['load -ascii ' filename1])
eval(['Data1=' filename1 ';']);
p11=Data1(:,1);
p21=Data1(:,2);

%rows that indicate a breakpoint
dp=find(p11(2:end)-p11(1:end-1)~=200);

%start and end GPS times of locked segments
st1=[p11(1);p11(dp+1)];
brk1=[p21(dp);p21(end)];

eval(['load -ascii ' filename2])
eval(['Data2=' filename2 ';']);
p12=Data2(:,1);
p22=Data2(:,2);

%rows that indicate a breakpoint
dp=find(p12(2:end)-p12(1:end-1)~=200);

%start and end GPS times of locked segments
st2=[p12(1);p12(dp+1)];
brk2=[p22(dp);p22(end)];

t0=693597613; %12/28/01 18:00 UTC: start of E7

%define a time vector, and a logical channel
t=[0:max(p21(end)-t0,p22(end)-t0)];
z1=zeros(size(t));
z2=zeros(size(t));


for i=1:length(st1)
   z1(st1(i)-t0+1:brk1(i)-t0+1)=ones(size([st1(i)-t0+1:brk1(i)-t0+1]));
end

for i=1:length(st2)
   z2(st2(i)-t0+1:brk2(i)-t0+1)=ones(size([st2(i)-t0+1:brk2(i)-t0+1]));
end

z=z1.*z2;

%%Find points with  two ifos locked

lck0=find(z~=0);

pts=find(lck0(2:end)-lck0(1:end-1)~=1); %break points: start or ends?

if z(1)==1, %it started locked
   brk=[lck0(pts) lck0(end)]; 
   st=[1 lck0(pts+1)]; 
else
   st=[lck0(1) lck0(pts+1)]; 
   brk=[lck0(pts) lck0(end)];   
end



N=length(brk); 	%number of locked segments: 
dl=brk-st+1;		%length of segments  
dnl=st(2:end)-brk(1:end-1);	%length of unlocked segments

st=st+t0;
brk=brk+t0;

disp(['Total locked time: ' num2str(sum(dl)/3600,3) ' hrs'])
disp(['Duty cycle: ' num2str(sum(dl)/t(end)*100,3) ' %'])
disp([])
disp(['Total time locked with locks longer than 15min: ' num2str(sum(dl(find(dl>15*60)))/3600,3) ' hrs'])
disp(['Duty cycle for long locks: ' num2str(sum(dl(find(dl>15*60)))/t(end)*100,3) ' %'])


%Plot length of locked segments vs time

%day labels
days={'12/28','12/29','12/30','12/31',...
      '01/01','01/02','01/03','01/04',...
      '01/05','01/06','01/07','01/08',...
      '01/09','01/10','01/11','01/12',...
      '01/13','01/14'};

figure(1)
Nmx=ceil(max(dl/600))*10;
bar((st-t0)/(3600*24),dl/60),
axis([-.5 17.5 0 Nmx])
hold on
plot((st-t0)/(3600*24),cumsum(dl/60)/60,'r',...
(st(find(dl>15*60))-t0)/(3600*24),cumsum(dl(find(dl>15*60))/60)/60,'k')
hold off
title([filename1 ' and ' filename2])
set(gca,'XTick',[.5:1:17])
set(gca,'XTickLabel',[])
for i=1:18
text(i-1.4,-Nmx/20,days{i},'FontSize',8)
end
legend(['cum sum of locked segments(hrs), total=' num2str(sum(dl/3600),3)],...
       ['cumsum of locked segments > 15min (hrs),total=' num2str(sum(dl(find(dl>15*60))/3600),3)],...
       'locked segments (min)')
grid on
shg

figure(2)
xx=[2.5:5:ceil(max(dl)/600)*10];
hist(dl/60,xx)
ylabel('# of locked segments')
xlabel('minutes')
title([filename1 ' and ' filename2])
shg

figure(3)
xx=[2.5:5:ceil(max(dnl)/600)*10];
hist(dnl/60,xx);
ylabel('# of unlocked segments')
xlabel('minutes')
title([filename1 ' and ' filename2])
shg


%mean lock length depending on the time of the day
for i=1:2:24
   z=find(mod((st-t0)/3600+18,24)>=(i-1) & mod((st-t0)/3600+18,24)<i+1);
   Mx((i-1)/2+1)=mean(dl(z)/60);
end

figure(4)
bar([1:12],Mx)
set(gca,'XTickLabel',{'0-2';'2-4';'4-6';'6-8';'8-10';...
      '10-12';'12-14';'14-16';'16-18';'18-20';'20-22';'22-0'},'FontSize',8)
xlabel('UTC hour start')
ylabel('mean lock length (min)')
title([filename1 ' and ' filename2])
shg


%%%%%%%%
%Print out gps times in a file, if asked


if nargin>2,
   
   E7days={'Dec 28','Dec 29','Dec 30','Dec 31',...
      'Jan 01','Jan 02','Jan 03','Jan 04',...
      'Jan 05','Jan 06','Jan 07','Jan 08',...
      'Jan 09','Jan 10','Jan 11','Jan 12',...
      'Jan 13','Jan 14'};

stday=floor(((st-t0)/3600+18)/24);
sthr=floor((st-t0)/3600+18-24*stday);
stmin=round((st-t0)/60+18*60-sthr*60-24*60*stday);

brkday=floor(((brk-t0)/3600+18)/24);
brkhr=floor((brk-t0)/3600+18-24*brkday);
brkmin=round((brk-t0)/60+18*60-brkhr*60-24*60*brkday);


eval(['fid=fopen(''' outname ''',''w'');'])

fprintf(fid,'%s\t','Seg#');
fprintf(fid,'%s ','Start Time UTC');
fprintf(fid,'\t %s','Stop Time UTC');
fprintf(fid,'\t %s \t \t %s \t %s \n','DT(min)','Start GPS','Stop GPS');
fprintf(fid,'\n');

for i=1:length(st)
   
   if dl(i)>5*60 %only write segments longer than 5 minutes
   
   fprintf(fid,'%d \t%',i);
   fprintf(fid,'%s ',E7days{stday(i)+1});
   fprintf(fid,'%02.0f',sthr(i));
   fprintf(fid,'%s',':');
   fprintf(fid,'%02.0f',stmin(i));
   
   fprintf(fid,'\t %s ',E7days{brkday(i)+1});
   fprintf(fid,'%02.0f',brkhr(i));
   fprintf(fid,'%s',':');
   fprintf(fid,'%02.0f',brkmin(i));

	fprintf(fid,'\t %3.0f \t \t %d \t %d \n',dl(i)/60,st(i),brk(i));
   
	end

if i<length(st) & stday(i+1)~=stday(i), fprintf(fid,'\n'); end


end

fclose(fid);

end