%
% file: make_streak_plot.sl
% 
% Usage:
%  isis> .source make_streak_plot
%
% reads in focus_trends_v2.dat
% and outputs streak_plots.ps
% 
% 10/18/08 - dd


% read the information from focus_trends_v2.dat:
variable obsids, yrs, mons, days, fwhms, errs;
(obsids,yrs,mons,days,fwhms,errs)=readcol("focus_trends_v2.dat",[2,4,5,6,13,14]);

% valid ones have fhhms > 0:
variable sel = where(fwhms > 0.0);

obsids=obsids[sel];
yrs=yrs[sel];
mons=int(mons[sel]);
days=days[sel];
fwhms=fwhms[sel];
errs=errs[sel];

% Convert the 90% confidence errors into 1-sigma:
errs = errs/1.65;
% limit error to 1 um rms:
sel = where(errs < 1.0);
errs[sel]=1.0;

% Create decimal year dates:
variable daysinmons,daysbeforemon, decidate, iloop, ilbl;
daysinmons=31.0+Float_Type[13];
daysinmons[0]=0;
daysinmons[9]=30.0;
daysinmons[4]=30.0;
daysinmons[6]=30.0;
daysinmons[11]=30.0;
daysinmons[2]=28.25;
daysbeforemon=cumsum([0.0,daysinmons]);
%
decidate = Float_Type[length(obsids)];
_for iloop (0,length(obsids)-1,1)
{
     decidate[iloop]= yrs[iloop] + (daysbeforemon[mons[iloop]]+days[iloop]-0.5)/365.5;
}

% Open plot file for hard copy:
PSwin=open_plot("streak_plot.ps/CPS");

% Set the plot range:
xrange(1999.5,2011.1);

xlin;
% have the yrange be 24 um = one pixel
yrange(20.0-0.5,44.0-0.5);
ylin;
label("Date (year)", "FWHM (um)", "FWHM of HETG Streak Core vs Time  (TGCat processed)");
plot(decidate,fwhms,15);
charsize(0.6);
color(14);
xylabel(1999.0,17.3,"created: "+time);
charsize(1.0);
% individual points etc
charsize(0.6);
ilbl=0;
_for iloop (0,length(obsids)-1,1)
{
   oplot(decidate[iloop]+[0.0,0.0],fwhms[iloop]+errs[iloop]*[-1.0,1.0],5);
   set_line_width(3);
   oplot(decidate[iloop]+0.07*[-1.0,0.0,1.0,0.0,-1.0],
     fwhms[iloop]+0.26*[0.0,1.0,0.0,-1.0,0.0],1);
   set_line_width(1);
   color(3);
   % add the obsid:
   xylabel(decidate[iloop],20.3+1.7*ilbl,string(obsids[iloop]),90.0);
   ilbl=ilbl+1;
   if (ilbl > 4) ilbl=0;
}
charsize(1.0);

% Calc and print the weighted mean value
% weights are 1/sigma^2
variable weights = 1.0/(errs)^2;
variable mean = sum(weights*fwhms)/sum(weights);
variable mean_err = sqrt(1.0/sum(weights));
message(" FWHM mean and error = "+string(mean)+",  "+string(mean_err));
% Calc the Chi-square
variable chi2n = sum( ((fwhms-mean)/errs)^2 ) / (length(fwhms)-1.0);
message(" chi-squared / N-1  = "+string(chi2n));

% show average and error limits on plot
variable date_range = [min(decidate)-0.5,max(decidate)+0.5];
% launch:
date_range[0] = 1999.0 + (daysbeforemon[7]+23.0-0.5)/365.5;
% could use today's date:
%% date_range[1] = 2010.0 + (daysbeforemon[1]+5.0-0.5)/365.5;
oplot(date_range,mean+[0.0,0.0],3);
oplot(date_range,mean+[0.0,0.0]+mean_err,5);
oplot(date_range,mean+[0.0,0.0]-mean_err,5);

color(1);
xylabel(2000.0, 41.0, "Ave FWHM = "+string((int(100.0*mean))/100.0)+
  " +/- "+string(int(100.0*mean_err)/100.0)+" um");
charsize(0.7);
xylabel(2000.0, 40.0, "Chi^2 / (N-1) = "+string(int(100.0*chi2n)/100.0)+
  "  [1um min err]");
charsize(1.0);

% all done
close_plot(PSwin);