% Dec. 5, 2003
%
% Subroutines used in the example of running the Bayesian Blocks Routine.
% No version numbers, since these are only designed to be guides as to how
% one can use s-lang/isis or s-lang/sherpa to read, plot, and write 
% Chandra data. These  are meant to be starting points for custom 
% modifications.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

public define sitar_examp_read( file )
{
   variable event_times = fits_read_col(file,"TIME");
   variable tstart = fits_read_key(file, "TSTART");
   variable tstop = fits_read_key(file, "TSTOP");
   variable frame = fits_read_key(file, "TIMEDEL");
   variable dtcor = fits_read_key(file,"DTCOR");
   variable object = fits_read_key(file, "OBJECT");
%
   return event_times, tstart, tstop, frame, dtcor, object;
%
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

public define sitar_examp_bin( event_times, tstart, tstop, dtcor, plot_step )
{
   variable ev = struct{ lo_t, hi_t, rate, err };
   ev.lo_t = [ tstart : tstop : plot_step ];
   ev.hi_t = ev.lo_t + plot_step;
   ev.rate = @ev.lo_t;
   ev.err = @ev.lo_t;

   variable i=0, bincts, intvl;
   loop(length(ev.lo_t))
   {  
      bincts = length( where( event_times >= ev.lo_t[i] and 
                              event_times  < ev.hi_t[i] ) );
      intvl = ( ev.hi_t[i] - ev.lo_t[i] ) * dtcor;
      ev.rate[i] = bincts / intvl;
      ev.err[i] = ( 0.75 + sqrt( 1. + bincts ) ) / intvl;
      i++;               
   } 
   return ev;
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Two examples here.  One is a very PGPLOT specific plot, while
% the other is Chips specific (thanks to Doug Burke)

#ifnexists isis_exit

   % plot up a "dataset" as a histogram in ChIPS
   %
   %
   % Although ChIPS does allow easy plotting of histograms
   % (plot points and then change the linestyle to step)
   % it is not really designed for histograms with
   % unequal bins, hence we "do it ourselves"
   % (assuming there are no gaps)

   private define _plot_histogram_chips( data, tstart ) {

      variable n = length(data.lo_t);
      variable x = @Array_Type ( _typeof(data.lo_t), [2*n] );
      variable y = @Array_Type ( _typeof(data.rate), [2*n] );

      variable idx = [0:2*n-1:2];
      x[idx] = data.lo_t - tstart;
      y[idx] = data.rate;
      idx++;
      x[idx] = data.hi_t - tstart;
      y[idx] = data.rate;

      () = curve( x, y );

   } % _plot_histogram_chips()
#endif

#ifnexists isis_exit

   % conversion of sitar_examp_plot() to use ChIPS
   % for plotting rather than the routines in ISIS
   % Doug Burke (10/06/03)

   public define sitar_examp_plot( 
                           results, ev, object, file, plot_step, ncp_prior, 
                           tstart, tstop, plotit
                          )
   {
      variable sig = typecast( (1.-exp(-ncp_prior))*10000., Integer_Type );
      sig = string(sig/100.);
      object = strtrans(object," ","");  

      variable oflag  = chips_auto_redraw(0);
      variable ochips = @chips;

      chips.symbolstyle = _chips->none;
      chips.curvestyle  = _chips->simpleline;

      chips_clear();

      _plot_histogram_chips( ev, tstart );
      () = chips_eval( "dot" );
      _plot_histogram_chips( results, tstart );
      () = chips_eval( "red" );

      () = chips_eval( "xlabel 'Time (sec)'" );
      () = chips_eval( "ylabel 'Rate (cps)'" );
      () = chips_eval( "title '"+
      object+"/"+file+" ("+string(plot_step)+" sec. bins, "+sig+"% sig.)"+
		    "'" );
      () = chips_eval( "xlabel size 1.2" );
      () = chips_eval( "ylabel size 1.2" );
      () = chips_eval( "title size 1.2" );
 
      () = chips_set_xrange( 0.0, tstop-tstart );
      () = chips_set_yrange( 
		 min(ev.rate[ [0:length(ev.rate)-2] ]) - 0.05*max(ev.rate),
		 1.05*max( [ results.rate, ev.rate ] ) );

      if(plotit)
      {
        () = chips_eval( "print postfile " + object+"_"+file+"_"
             +string(plot_step)+"_"+sig+".ps" );
      }
      chips_redraw();
      () = chips_auto_redraw( oflag );
      set_state( "chips", ochips );

      return;
   }
#endif

#ifexists isis_exit
   public define sitar_examp_plot( 
                        results, ev, object, file, plot_step, ncp_prior, 
                        tstart, tstop, plotit
                          )
   {
      variable sig = typecast( (1.-exp(-ncp_prior))*10000., Integer_Type );
      sig = string(sig/100.);
      object = strtrans(object," ","");  

      if(plotit)
      {
         variable id = open_plot(object+"_"+file+"_"+string(plot_step)+
                                 "_"+sig+".ps/vcps");
         resize(10,0.9);
      }

      label(
            "Time (sec)", "Rate (cps)",
            "\\fr "+object+"/"+file+" ("+string(plot_step)+" sec. bins, "+
            sig+"% sig.)"
           );

      linestyle(4);
      set_line_width(3);
      charsize(1.2);

      xrange( 0., tstop - tstart );
      yrange( min(ev.rate[ [0:length(ev.rate)-2] ]) - 0.05*max(ev.rate),
              1.05*max( [ results.rate, ev.rate ] ) );

      % For plotting purposes, renorm to input start time

      hplot( ev.lo_t-tstart, ev.hi_t-tstart, ev.rate );
      linestyle(1);
      set_line_width(6);
      ohplot( results.lo_t-tstart, results.hi_t-tstart, results.rate, 2 );

      if(plotit)
      {
         close_plot(id);
      }
      return;
   }
#endif

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

public define sitar_examp_write( results, object, file, plot_step, ncp_prior )
{
   variable sig = typecast((1.-exp(-ncp_prior))*10000.,Integer_Type);
   sig = string(sig/100.);
   object = strtrans(object," ","");
   variable fp = fopen(object+"_"+file+"_"+string(plot_step)+"_"
                       +sig+".dat","w");
   variable i=0;

   % NB: Change location is specific to the output of "sitar_make_data_cells".

   loop( length(results.cpt) )
   {
      () = fprintf(fp,"%9i  %15.7e  %15.7e  %12.3e  %12.3e  \n",
           results.cpt[i], results.lo_t[i], results.hi_t[i], 
           results.rate[i], results.err[i]);
      i++;
   }
   () = fprintf(fp,"Change Location -  Lower Time - Upper Time - Rate (cps) - Error  (cps)\n");
   () = fclose(fp);
   return;
}  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%