%Plot time series and histogram of alpha-beta during science time
%segments

t0=751651213; %10/31/03 16:00 UTC: start of S3

%list of science segments
%seg#, start, stop, duration, dqbit 
S3=load('S3L1v00_segs.txt');

%list of alpha*beta values
%Data  105573x5                  4222920  double array
%'GPS'    'Alpha*Beta'    'Alpha'    'Beta'    'CalLine/Ref'
load L1AlphaBeta.mat


%Order Data like H1, H2 so programs are identical
Data=[Data(:,1) Data(:,3) Data(:,4) Data(:,2)];


%limits to gains
lg=0.404;
hg=2.27;

%tolerance
tol=0.20;
llg=(1+tol)*lg;
hhg=(1-tol)*hg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Select the times in science mode, write out low values, plot amplitudes


%Write out science seg#, bad gps time; duration
fid=fopen('L1CalDQ.txt','w');

figure(1)  

for k=1:size(S3,1)
  i0=min(find(Data(:,1)>=S3(k,2)-60));
  ie=max(find(Data(:,1)<=S3(k,3)));
  if isempty(i0)~=1 & isempty(ie)~=1 & (i0<ie), 
    ix=[i0:ie]; 
  else
    ix=[];
    fprintf(fid,'%i\t',S3(k,1));
    fprintf(fid,'%s\n','No data');
  end
  
  if length(ix)~=0,         
    %write out too high, too low,  or zero values for alpha beta
    ibad=find((Data(ix,4)<=llg) | (Data(ix,4)>=hhg));
    for j=1:length(ibad)
       fprintf(fid,'%i\t',S3(k,1));
       fprintf(fid,'%i\t',Data(ix(ibad(j)),1));
       fprintf(fid,'%i\t',60);
       fprintf(fid,'%2.4f\n',Data(ix(ibad(j)),4));
    end
        
    %do we have good data in the beginning?
    if Data(i0,1)-S3(k,2)>0,
      z0=S3(k,2);
      dz0=Data(i0,1)-S3(k,2);
      fprintf(fid,'%i\t',S3(k,1));
      fprintf(fid,'%i\t',z0);
      fprintf(fid,'%i\n',dz0);
    else
      z0=[];
      dz0=[];
    end
   
      
    %find whether there are some gaps in time larger than 60s
    n60=find(Data(ix(2:end),1)-Data(ix(1:end-1),1)~=60);
    
    if isempty(n60)~=1, 
      
      for j=1:length(n60)
	zn(j)=Data(ix(n60(j)),1)+60;
	dzn(j)=Data(ix(n60(j)+1),1)-zn(j);
	fprintf(fid,'%i\t',S3(k,1));
	fprintf(fid,'%i\t',zn(j));
	fprintf(fid,'%i\n',dzn(j));
      end
      
    else
      zn=[];
      dzn=[];
    end
    
    
    
    %do we have good data in the end?
    if S3(k,3)-Data(ie,1)>60,
      ze=Data(ie,1)+60;
      dze=S3(k,3)-Data(ie,1)-60;
      fprintf(fid,'%i\t',S3(k,1));
      fprintf(fid,'%i\t',ze);
      fprintf(fid,'%i\n',dze);
    else
      ze=[];
      dze=[];
    end
    
    %data to plot with good and bad points
    Zgps0=[z0;zn';ze;Data(ix,1)];
    [Zgps,iz]=sort(Zgps0);
    ZData=zeros(size(Zgps));
    ZData0=zeros(size(Zgps));
    ZData0(end-length(Data(ix,1))+1:end)=Data(ix,4);
    ZData=ZData0(iz);
    
   
    [zbad,izbad]=setdiff(Zgps,Data(ix,1));
    [zbad,ixbad]=intersect(Zgps,Data(ix(ibad),1));
    izbad=union(izbad,ixbad);
    
    if 1==1
    %Plot the segment
    clf
    %add a point at the end to cover the segment
    ixp=[ix min(ie+1,size(Data,1))];
    plot(Data(ixp,1)-S3(k,2),Data(ixp,2),'b.-',...
	 Data(ixp,1)-S3(k,2),Data(ixp,3),'r.-',...
	 Data(ixp,1)-S3(k,2),Data(ixp,4),'k.-')
    legend('\alpha','\beta','\alpha\beta')
    if isempty(izbad)~=1, 
      plot(Data(ixp,1)-S3(k,2),Data(ixp,2),'b.-',...
	   Data(ixp,1)-S3(k,2),Data(ixp,3),'r.-',...
	   Zgps-S3(k,2),ZData,'k.-',...
	   Zgps(izbad)-S3(k,2),ZData(izbad),'ko',...
	   Data(ixp,1)-S3(k,2),Data(ixp,4),'k.-')
      legend('\alpha','\beta','\alpha\beta','bad \alpha\beta')
    end
    hold on
    plot([-60 S3(k,4)+60],hhg*[1 1],'k-.',...
	 [-60 S3(k,4)+60],llg*[1 1],'k-.',...
	 [-60 S3(k,4)+60],hg*[1 1],'k',...
	 [-60 S3(k,4)+60],lg*[1 1],'k',...
	 [0 0],[0 2.5],'k',...
	 S3(k,4)*[1 1],[0 2],'k','LineWidth',2)
    hold off
    
    
    title(['L1 S3 Seg#' num2str(S3(k,1)) ', ' ...
	   num2str(S3(k,2)) '-' num2str(S3(k,3))])
    xlabel(['seconds since ' num2str(S3(k,2))])

    axis([-60 S3(k,4)+60 0 2])
    shg
    eval(['print -dpng L1Plots\S3L1CalibSeg_' num2str(S3(k,1),'%03d') '.png'])
    pause(.1)
    end
    
  end
  
end
fclose(fid);