%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Program to identify safe bands for instrumental dithers that do not              %%%
%%% interfere with aLIGO searches for known pulsars with accessible spindown limits  %%%
%%% K. Riles - October 2014                                                          %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Read in output of ATNF catalog (http://www.atnf.csiro.au/people/pulsar/psrcat/) 
% created with following data fields and conditions:
%
% Data fields:  JName F0 F1 PEpoch Dist
% Conditions:   F0>=5.0 && exist(F0) && exist(F1) && exist(PEpoch) && exist(Dist) 
%
% [also manually deleted blank lines after every five pulsars and removed entries for
%  J1748-2446ad (716.35556 Hz) and J1400-1438 (324.22969825 Hz) because F1 = 0]
%
% and plot veto windows within 0.5 Hz of the pulsars projected frequencies 
% (F and 2*F on July 1, 2015 (MJD 57204)

% Read in catalog output and parse
fid = fopen('atnf_catalog_28oct2014.txt','r');
data = textscan(fid,'%d %s %s %f %d %s %f %d %s %f %s %f');
index = data{1};
JName = data{2};
F0 = data{4};
F1 = data{7};
PEpoch = data{10};
Dist = data{12};
npulsar = length(JName);

% Compute frequency on reference day
MJDref = 57204;
F = F0 + F1.*(MJDref-PEpoch)*86400;

% Compute spindown limit
h_spindown = 2.5e-25 * sqrt( (1000.0./F) .* abs(F1)/1.e-10 ) ./ Dist; 

% Read in and parse projected aLIGO sensitivity curve
aligo_projection = load('ZERO_DET_high_P.txt');
aligo_freq = aligo_projection(:,1);
aligo_sqrtSh = aligo_projection(:,2);

% Compute h0 strain sensitivity for coherent 1-year targeted search with two detectors
aligo_h0sens_1year = 11.4*aligo_sqrtSh/sqrt(3.16e7*2.);

% Define half-width of veto band around each pulsar with a beatable spindown limit 
%vetohalfband = 0.01;
%vetohalfband = 0.10;
%vetohalfband = 0.50;
vetohalfband = 1.00;

% Define minimum instrumental frequency to consider
Fminveto = 10.;

% Compute projected strain sensitivity for each pulsar at 1*F and 2*F
aligo_hsens_interp_1F = zeros(npulsar,1);
aligo_hsens_interp_2F = zeros(npulsar,1);
for i = 1:npulsar
   if F(i)>=Fminveto
      aligo_h0sens_interp_1F(i) = interp1(aligo_freq,aligo_h0sens_1year,F(i));
   end
   if 2.*F(i)>=Fminveto
      aligo_h0sens_interp_2F(i) = interp1(aligo_freq,aligo_h0sens_1year,2.*F(i));
   end
end

% Define initial vetoed spectrum up to 2 kHz in 0.01-Hz intervals (nothing vetoed yet)
band = [Fminveto:0.01:2000.];
nband = length(band);
veto_1F = false(1,nband);
veto_2F = false(1,nband);

% Open output text file
fname = sprintf('Pulsar_Gaps_Veto_%4.2fHz.txt',vetohalfband);
fidout = fopen(fname,'w');

% Loop over pulsars to veto bands too close to 1*F or 2*F frequencies
for i = 1:npulsar
   if F(i)>= Fminveto && aligo_h0sens_interp_1F(i) < h_spindown(i)
   fprintf('Vetoing band around pulsar %d %s : 1*Freq = %f, spindown limit = %e, 1-year sens = %e\n',index(i),JName{i},F(i),h_spindown(i),aligo_h0sens_interp_1F(i));
      fprintf(fidout,'Vetoing band around pulsar %d %s : 1*Freq = %f, spindown limit = %e, 1-year sens = %e\n',index(i),JName{i},F(i),h_spindown(i),aligo_h0sens_interp_1F(i));
      test = abs(band-F(i)) < vetohalfband;
      veto_1F = or(veto_1F,test);
   end
   if 2.*F(i)>= Fminveto && aligo_h0sens_interp_2F(i) < h_spindown(i)
      fprintf('Vetoing band around pulsar %d %s : 2*Freq = %f, spindown limit = %e, 1-year sens = %e\n',index(i),JName{i},2.*F(i),h_spindown(i),aligo_h0sens_interp_2F(i));
      fprintf(fidout,'Vetoing band around pulsar %d %s : 2*Freq = %f, spindown limit = %e, 1-year sens = %e\n',index(i),JName{i},2.*F(i),h_spindown(i),aligo_h0sens_interp_2F(i));
      test = abs(band-2.*F(i)) < vetohalfband;
      veto_2F = or(veto_2F,test);
   end
end

% Define veto for either 1*F or 2*F
veto_Both = or(veto_1F,veto_2F);

% Find contiguous subbands:   
runs_1F = contiguous(veto_1F,[0 1]);
okay_intervals_1F = runs_1F{1,2};
vetoed_intervals_1F = runs_1F{2,2};

runs_2F = contiguous(veto_2F,[0 1]);
okay_intervals_2F = runs_2F{1,2};
vetoed_intervals_2F = runs_2F{2,2};

runs_Both = contiguous(veto_Both,[0 1]);
okay_intervals_Both = runs_Both{1,2};
vetoed_intervals_Both = runs_Both{2,2};

% Print results:
fprintf('Non-vetoed bands for veto half-band = %f (one times pulsar frequency)\n',vetohalfband);
fprintf(fidout,'Non-vetoed bands for veto half-band = %f (one times pulsar frequency)\n',vetohalfband);
nokay_1F = length(okay_intervals_1F);
for interval = 1:nokay_1F
   fprintf('  %7.2f-%7.2f Hz (%7.2f Hz)\n',Fminveto+(okay_intervals_1F(interval,1)-1)*0.01,Fminveto+(okay_intervals_1F(interval,2)-1)*0.01,(okay_intervals_1F(interval,2)-okay_intervals_1F(interval,1)+1)*0.01);
   fprintf(fidout,'  %7.2f-%7.2f Hz (%7.2f Hz)\n',Fminveto+(okay_intervals_1F(interval,1)-1)*0.01,Fminveto+(okay_intervals_1F(interval,2)-1)*0.01,(okay_intervals_1F(interval,2)-okay_intervals_1F(interval,1)+1)*0.01);
end
fprintf('Non-vetoed bands for veto half-band = %f (two times pulsar frequency)\n',vetohalfband);
fprintf(fidout,'Non-vetoed bands for veto half-band = %f (two times pulsar frequency)\n',vetohalfband);
nokay_2F = length(okay_intervals_2F);
for interval = 1:nokay_2F
   fprintf('  %7.2f-%7.2f Hz (%7.2f Hz)\n',Fminveto+(okay_intervals_2F(interval,1)-1)*0.01,Fminveto+(okay_intervals_2F(interval,2)-1)*0.01,(okay_intervals_2F(interval,2)-okay_intervals_2F(interval,1)+1)*0.01);
   fprintf(fidout,'  %7.2f-%7.2f Hz (%7.2f Hz)\n',Fminveto+(okay_intervals_2F(interval,1)-1)*0.01,Fminveto+(okay_intervals_2F(interval,2)-1)*0.01,(okay_intervals_2F(interval,2)-okay_intervals_2F(interval,1)+1)*0.01);
end
fprintf('Non-vetoed bands for veto half-band = %f (one or two times pulsar frequency)\n',vetohalfband);
fprintf(fidout,'Non-vetoed bands for veto half-band = %f (one or two times pulsar frequency)\n',vetohalfband);
nokay_Both = length(okay_intervals_Both);
for interval = 1:nokay_Both
   fprintf('  %7.2f-%7.2f Hz (%7.2f Hz)\n',Fminveto+(okay_intervals_Both(interval,1)-1)*0.01,Fminveto+(okay_intervals_Both(interval,2)-1)*0.01,(okay_intervals_Both(interval,2)-okay_intervals_Both(interval,1)+1)*0.01);
   fprintf(fidout,'  %7.2f-%7.2f Hz (%7.2f Hz)\n',Fminveto+(okay_intervals_Both(interval,1)-1)*0.01,Fminveto+(okay_intervals_Both(interval,2)-1)*0.01,(okay_intervals_Both(interval,2)-okay_intervals_Both(interval,1)+1)*0.01);
end

Fminplot = [10 30 100 400];
Fmaxplot = [30 100 400 2000];
hmaxplot = [4.e-24 4.e-24 4.e-24 4.e-24];

for iplot = 1:4
   subplot(2,2,iplot)
%   plot(aligo_freq,aligo_h0sens_1year,'color','blue');
   semilogy(aligo_freq,aligo_h0sens_1year,'color','blue');
   xlim([Fminplot(iplot) Fmaxplot(iplot)]);
   ylim([0 hmaxplot(iplot)]);
   xlabel('Freq (Hz)')
   ylabel('Strain sensitivity / spindown limits')
   titlestring = sprintf('Bands to avoid (veto half-band = %4.2f Hz) for %6.1f-%6.1f Hz',vetohalfband,Fminplot(iplot),Fmaxplot(iplot));
   title(titlestring)
   hold on
   for i = 1:npulsar
      if F(i) >= Fminveto 
         xmin = F(i)-vetohalfband;
         xmax = F(i)+vetohalfband;
         ymin = 1.e-29;
         ymax = h_spindown(i);
         xp = [xmin xmax xmax xmin xmin];
         yp = [ymin ymin ymax ymax ymin];
         if aligo_h0sens_interp_1F(i) < h_spindown(i)
            patch(xp,yp,'m','edgecolor','m');
         else
            patch(xp,yp,'k');
         end
      end
      if 2.*F(i) >= Fminveto 
         xmin = 2.*F(i)-vetohalfband;
         xmax = 2.*F(i)+vetohalfband;
         ymin = 1.e-29;
         ymax = h_spindown(i);
         xp = [xmin xmax xmax xmin xmin];
         yp = [ymin ymin ymax ymax ymin];
         if aligo_h0sens_interp_2F(i) < h_spindown(i)
            patch(xp,yp,'g','edgecolor','g');
         else
            patch(xp,yp,'k');
         end
      end
   end
   hold off
end

fnameplot = sprintf('Pulsar_Gaps_Veto_%4.2fHz.png',vetohalfband);
print('-dpng',fnameplot);
