% SITAR - S-lang/ISIS Timing Analysis Routines 
% ***  Version 0.8.0 *** Oct. 15, 2007  ***
%    ***  Sherpa Compatible Version ***

% (Mostly) Alphabetical List of Routines in this File

% epfold_rate_use (v. 2.0.00)
% fstat           (v. 1.0.00)
% import_mutils
% glob_opt_use    (v. 2.4.00)
% log_gamma       (v. 1.0.00)
% log_post        (v. 0.4.02)
% mdc_use         (v. 2.1.00)
% msg
% ns_pop
% null_it         (v. 1.0.00)
% pfold_rate_use  (v. 2.0.00)
% readasm_use     (v. 2.0.00)
% reb_rat_use     (v. 2.1.00)

% Sherpa, in Ciao 3.3.01, imports ISIS v1.3.0, so you should be OK.

% sitar_avg_cpd               (v. 0.2.00)
% sitar_avg_psd               (v. 0.2.00)
% sitar_define_psd            (v. 0.1.01)  => ISIS Version
% sitar_define_psd            (v. 0.1.01s) => Sherpa Version
% sitar_epfold_rate           (v. 0.2.00) !!Requires ISIS > v1.1.3 !!
% sitar_global_optimum        (v. 0.4.01) !!Requires ISIS > v1.0.48!!
% sitar_lags                  (v. 0.1.00)
% sitar_lbin_cpd              (v. 0.2.00)
% sitar_lbin_psd              (v. 0.2.00)
% sitar_make_data_cells       (v. 0.4.00)
% sitar_pfold_rate            (v. 0.2.00) !!Requires ISIS > v1.1.3 !!
%  _sitar_pfold_profile       (v. 1.0.00)
%  _sitar_pfold_rate          (v. 1.0.00)
%  _sitar_pfold_rate_deriv    (v. 1.1.00)
%  _sitar_pfold_rate_noderiv  (v. 1.0.00)
% sitar_readasm               (v. 0.2.00)
% sitar_rebin_rate            (v. 0.4.00) !!Requires ISIS > v1.1.3 !!

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

% If you've got the s-lang GNU Scientific Module loaded,
% uncomment the following *before* running.  Necessary for statistics
% of epoch fold (dummy function used if left uncommented), and for
% additional functionality of the Bayesian Blocks routine (i.e.,
% `event mode' with alpha & beta > 0, **which won't crash, but won't
% give real results, either, if used without GSL**).

%require("gsl");


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

% The following is adapted from David Huenemoeder, and is designed 
% to get those utilities from ISIS that are necessary to run the
% routines.  If you're running ISIS, then you've got them already, so
% no worries.  If you're running Sherpa, this should provide you with
% everything you need.

% nspop will take a namespace and a function name and re-define it in
% the current namespace.

define nspop( ns, name )
{
   if ( is_defined("Global->" + name ) )
   {
      return; 
   }
   else
   {
      eval(sprintf("public define %s() {variable args = __pop_args(_NARGS); return %s->%s(__push_args(args));}", 
           name, ns, name ) );
   }
}       

define import_mutils()
{

% isis_exit is a function name that should probably only exist in
% ISIS, and nowhere else.  So, we use that for our check.  If it
% doesn't exist, presume that we're in Sherpa, and import ISIS.
% Then start importing functions, but make sure they don't clobber
% anything in Sherpa.

#ifnexists isis_exit
   import("isis","_ns_ISIS");
#endif

#ifnexists cumsum
   nspop("_ns_ISIS", "cumsum"      ) ;
#endif

#ifnexists fft
   nspop("_ns_ISIS", "fft"         ) ;
#endif

#ifnexists histogram
   nspop("_ns_ISIS", "histogram"   ) ;
#endif

#ifnexists linear_grid
   nspop("_ns_ISIS", "linear_grid" ) ;
#endif

#ifnexists make_hi_grid
   nspop("_ns_ISIS", "make_hi_grid"         ) ;
#endif

#ifnexists max
   nspop("_ns_ISIS", "max"         ) ;
#endif

#ifnexists mean
   nspop("_ns_ISIS", "mean"        ) ;
#endif

#ifnexists min
   nspop("_ns_ISIS", "min"         ) ;
#endif

#ifnexists moment
   nspop("_ns_ISIS", "moment"      ) ;
#endif

#ifnexists ones
   nspop("_ns_ISIS", "ones"       ) ;
#endif

#ifnexists reverse
   nspop("_ns_ISIS", "reverse"     ) ;
#endif

#ifnexists shift
   nspop("_ns_ISIS", "shift"       ) ;
#endif

#ifnexists sum
   nspop("_ns_ISIS", "sum"         ) ;
#endif

}

% Grab stuff from ISIS

import_mutils;

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

provide("sitar");

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

private define msg(str_array)
{
   () = printf("\n");
   foreach(str_array)
   {
      variable str = ();
      () = printf(" %s\n", str);
   }
   () = printf("\n");
   return;
}

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

private define null_it(a,i)
{

%       Version 1.0.00: 2003/05/26 MAN 

   variable j = Char_Type[ length(a) ];
   j[i] = 1;
   return a[ where(j==0) ];
}

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

private define mdc_use()
{

%       Version 2.2.00: 2004/03/30 MAN 

   message(" ");
   message("Usage:");
   message(" ");
   message(" cell = sitar_make_data_cells(tt,type,max_delt,frame,tstart,tstop);");
   message(" ");
   message(" Inputs: ");
   message("   tt: The times of the detected events");
   message("   type: Describe cells as 'events' (type=1 or 2), with the");
   message("         size (i.e., normalized duration) being the distance");
   message("         from halfway from the previous event to halfway to");
   message("         the next event (type=1), or from the current event");
   message("         to right before the subsequent event (type=2).");
   message("         Alternatively, events can be assigned to bins of");
   message("         uniform size (type=3).");
   message("   max_delt: Events closer together than max_delt are grouped");
   message("           together in a single cell.");
   message("   frame: The output cell sizes are in units of frame, or this");
   message("          is the bin size for type=3.");
   message("   tstart/tstop: start/stop times for making the output cells.");
   message(" ");
   message(" Outputs: ");
   message("   cell.pops: Array with the number of events per cell");  
   message("   cell.size: Event mode: size (i.e. duration) of a cell in units of");
   message("              frame.  Deadtime/efficiency should be folded into this");
   message("              definition.  E.g., size= total \"good time\" / total \"good");
   message("              time\" for an individual frame.  (Here, we use the equivalent");
   message("              for uniform deadtime, size = total time / total frame time)");
   message("              Binned mode: size is mean (over whole observation) counts");
   message("              per bin.  (Again, here is where you would put in efficiency"); 
   message("              factors, *if* efficiency/deadtime is not uniform bin to bin)");
   message("   cell.lo_t,.hi_t:  cell start and stop times");
   message("   cell.dtcor: Array, currently set to unity, for storing `dead"); 
   message("               time' or 'efficiency' corrections for the cells.");
   message("               Must be set, after the fact, by the user.");
   message("               Currently only used in sitar_global_optimum to");
   message("               give output block rates corrected for dead time");    
   message(" ");

   return;
}
  

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

public define sitar_make_data_cells( )
{ 

% NAME: 
%     sitar_make_data_cells
%
% PURPOSE:
%     Taking time-tagged event data, create structure with cell 
%     populations and sizes, suitable to send to global_optimum(),
%     and bin time boundaries, useful for plotting later on.
%
% CALLING SEQUENCE:
%     cell = sitar_make_data_cells(tt,type,max_delt,frame,tstart,tstop);
%
% INPUTS:
%     tt         : Time tags of *detected* events.
%     type       : Event data (type = 1 or 2), or binned data (type = 3)
%                  type = 1: cell is half-way from previous event to
%                            half-way to next event. (Voronoi)
%                  type = 2: cell is from first event to right before next
%                            event. (Delaunay)
%	           type = 3: Turn the event list into bins, with spacing
%                            of frame (max_delt then becomes a dummy 
%                            variable).  NOTE:  If you already have binned
%	                     counts, then set cell.pops = counts, and
%		             cell.size = an array of 1's, with length of
%                            cell.pops.
%                  type 1 vs. 2 makes the biggest difference for sparse
%	           lightcurves.  Note: I handle the endpoints slightly
%	           differently than the Scargle routines.
%     max_delt   : Maximum separation of times between nearest neighbors
%	           in any given cell.
%	           dt <= max_delt in the same cell (cell top to bottom
%                                                   can be > max_delt)
%	           dt > max_delt in different cells
%		   (dt = 0 corresponds to duplicate time tags) **OR**
%		   dummy variable for type 3
%     frame      : Length of a clock "tick", to determine how many have
%	           occured in a given interval **OR** the bin size for
%                  binned mode (type=3).
%     tstart     : Start of the lightcurve (not necessarily first detected
%                  event).
%     tstop      : Stop of the lightcurve (not necessarily last detected
%                  event).
%
% OUTPUTS:
%     cell       : Structure with cell.pops,cell.size,cell.lo_t,cell.hi_t,
%                  cell.dtcor
%     cell.pops  : Array of cell populations
%     cell.size  : Array of cell sizes (in units of 'frame').
%     cell.lo_t  : Times associated with bottom edge (>=) of cells
%     cell.hi_t  : Times associated with upper edge  (<) of cells
%     cell.dtcor : For future use as "dead time correction"
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     shift      : ISIS function to shift the elements of an array
%     mean       : ISIS function to calculate the mean of an array
%     ones       : ISIS function to create a 1-D array of ones of
%                  a given length.
%     linear_grid: ISIS function to create a linear grid, suitable 
%                  for making a histogram
%     histogram  : ISIS function to create a histogram (to make 
%                  type=3 data out of event list)
%     null_it    : Function, to remove indices from an array, 
%                  included above.
% 
% MODIFICATION HISTORY:
%     Based upon Jeff Scargle's MATLAB routines to perform a 
%     Bayesian blocks decomposition of an ordered array.
%
%     Paper: "Studies in Astronomical Time Series Analysis. VI.
%             Optimum Partition of the Interval: Bayesian Blocks,
%	      Histograms, and Triggers", J.D. Scargle et al., 2003
%             (to be submitted).
%		
%     Web Site:  http://astrophysics.arc.nasa.gov/~jeffrey
%
%     Version 0.1.00: 2003/02/25 Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.1.01: 2003/02/27 MAN - Changed to the cell structure output
%     Version 0.1.02: 2003/04/10 MAN - Added tstart and tstop, fixed bug with
%                                      Type 2, and changed to cell.lo/hi_t 
%                                      output
%     Version 0.1.03: 2003/04/15 MAN - Changed type=3 slightly
%     Version 0.1.04: 2003/04/17 MAN - Further bug in Type=2 repaired
%     Version 0.1.05: 2003/04/22 MAN - Changed the way endpoints are handled
%     Version 0.1.06: 2003/04/24 MAN - More diddling with endpoints.
%     Version 0.2.00: 2003/05/20 MAN - Final diddles with endpoints
%     Version 0.3.00: 2003/05/26 MAN - Beginning to think about GTI, 
%                                      define binned rate to be one per
%                                      cell size, and some style changes
%     Version 0.3.01: 2003/09/16 MAN - Use message added
%     Version 0.4.00: 2003/12/12 MAN - Changed the definition of
%                                      cell.size for binned mode!


   switch ( _NARGS )
   {
      case 0:
      mdc_use( );
      return;
   }
   {
      variable tt,type,max_delt,frame,tstart,tstop;
      (tt,type,max_delt,frame,tstart,tstop) = ( );
   }

   %  Initially, one datum per cell.

   variable cell = struct{ pops, size, lo_t, hi_t, dtcor };
   tt = tt[ where( tt >= tstart and tt <= tstop ) ];
   cell.pops = ones( length(tt) );

   variable mintt = min(tt);
   variable maxtt = max(tt);
   variable difftt = shift(tt,1) - tt;
   difftt = difftt[ [0:length(difftt)-2] ];

   variable ii_close = where( difftt < max_delt );
   variable ii_start, ii_beyond, ii_end, ii_clump;
   variable  clump_pop, clump_tt, ind_null;

   if ( type !=3 ) while( length(ii_close) )
   {
      ii_start = ii_close[0];
      ii_beyond = where( shift(ii_close,1) - ii_close > 1 );

      if ( length(ii_beyond) == 0 )
      {
         ii_end = ii_start + length(ii_close);
      }
      else
      {
         ii_end = ii_start + ii_beyond[0] + 1;
      }

      ii_clump = [ii_start:ii_end];

      clump_pop = sum( cell.pops[ii_clump] );
      clump_pop = typecast(clump_pop,Integer_Type);
      clump_tt = mean( tt[ii_clump] );

      % Put the members of a clump into one cell

      cell.pops[ii_start] = clump_pop;
      tt[ii_start] = clump_tt;

      % Null the cells evacuated by this operation

      if ( ii_end > ii_start+1 )
      {
         ind_null = [ii_start+1:ii_end];
      }
      else
      {
         ind_null = [ii_end:ii_end];
      }
      cell.pops = null_it(cell.pops,ind_null);
      tt = null_it(tt,ind_null);

      difftt = shift(tt,1) - tt;
      difftt = difftt[ [0:length(difftt)-2] ];
      ii_close = where( difftt < max_delt );
   }

   % Now define the cell sizes and locations.

   variable dt = shift(tt,1) - tt;
   variable ndt = length(dt)-1;
   dt = dt[[0:ndt-1]];
   variable cstart = 1, cstop = ndt-1;

   switch (type)
   { 
    case 1:

    % type = 1: Set to midpoints

      cell.size = 0.5 * (dt[[0:ndt-2]] + dt[[1:ndt-1]]);
      cell.lo_t = tt[[1:ndt-1]]-dt[[0:ndt-2]]/2.;
      cell.hi_t = tt[[1:ndt-1]]+dt[[1:ndt-1]]/2.;

      % If tstart is < first event, add in the first event.  Design feature
      % for sparse lightcurves, where a long wait time until the first event
      % might mean something.

      if ( tstart < mintt )
      {
         cell.size = [ tt[0]-tstart+dt[0]/2., cell.size ];
         cell.lo_t = [ tstart, cell.lo_t ];
         cell.hi_t = [ tt[0]+dt[0]/2., cell.hi_t ];
         cstart = 0;
      }

      % If tstop is > last event, add in the last event.  Design feature
      % for sparse lightcurves, where a long drought at the end might
      % mean something.

      if ( tstop > maxtt )
      {
         cell.size = [ cell.size, tstop-tt[ndt]+dt[ndt-1] ];
         cell.lo_t = [ cell.lo_t, tt[ndt-1]+dt[ndt-1]/2. ];
         cell.hi_t = [ cell.hi_t, tstop ];
         cstop = ndt;
      }

      % Now make cell.pops match the above.

      cell.pops = cell.pops[ [cstart:cstop] ];
   }

   {
    case 2:

    % type =2: Intervals instead of midpoints. 

      cell.size = dt;
      cell.lo_t = tt[ [0:ndt-1] ];
      cell.hi_t = tt[ [1:ndt] ];
      cstart = 0;
      cstop = ndt-1;

      % If tstop is > last event, add in the last event.  Design feature
      % for sparse lightcurves, where a long drought at the end might
      % mean something.

      if ( tstop > maxtt )
      {
         cell.size = [ cell.size, tstop-tt[ndt] ];
         cell.lo_t = [ cell.lo_t, tt[ndt] ];
         cell.hi_t = [ cell.hi_t, tstop ];
         cstop = ndt;
      }

      % I'm not entirely happy with any choices for the endpoints for
      % sparse lightcurves.  But one has to choose *something*.  For
      % tstart < mintt, append that time to the first bin.  Thought for
      % alternative and type=2: add a bin with 0 photons going from 
      % tstart to right before the first photon?  (Should be OK with
      % the statistic?)

      if ( tstart < mintt )
      {
         cell.size[0] = cell.size[0] + ( tt[0]-tstart );
         cell.lo_t[0] = tstart;
      }
      cell.pops = cell.pops[ [cstart:cstop] ];
   }

   {
    case 3:

    % type = 3: Bin using ISIS functions

      variable nbins = typecast( (tstop-tstart)/frame, Integer_Type );
      ( cell.lo_t, cell.hi_t ) = linear_grid( tstart, tstop, nbins);

      cell.pops = histogram( tt, cell.lo_t, cell.hi_t );
      cell.dtcor = ones( length(cell.pops) );
      cell.size = @cell.dtcor;

      % Based upon my thoughts of applying this to gratings data,
      % I've decided upon the following modification.  This forces
      % the "rate per cell size" to be typically less than one.  Or, 
      % in terms of the marginalized likelihood, this puts a sum
      % of "model counts per block" in the denominator (presuming that
      % your zeroth order model is `constant rate over the whole
      % observation').

      cell.size = cell.size * ( sum(cell.pops)/length(cell.size) );
      return cell;
   }

   {

    % If it's not type 1, 2, or 3, something is wrong! 

      error("Type not 1, 2, or 3 in 'sitar_make_data_cells'");
   }

   % The return for Type 1 or 2

   cell.size = cell.size/frame;
   cell.dtcor = ones( length(cell.pops) );	
   return cell;
}

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

private define log_gamma(z)
{

%       Version 1.0.00: 2003/05/26 MAN 

   variable c = [76.18009173,-86.50532033,24.01409822,-1.231739516,
                 1.20858003e-3,-5.36382e-6];
   variable csum=1.;
   foreach(c)
   {  
      csum+=()/z;
      z+=1.;
   }
   csum = log( sqrt(2.*PI)*csum ) - (z-1.5) + (z-6.5)*log(z-1.5);
   return csum;
}

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

#ifnexists beta_inc
   private define beta_inc(a,b,x)
   {
      % This is in **no way meant to be an approximation**.
      % Merely a device to keep the program from crashing.
      %
      return x;
   }
#endif

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
private define log_post( data_pops, data_size, ncp_prior, type, alpha, beta, lo )
{

% NAME: 
%     log_post
%
% PURPOSE:
%     Used in Bayesian decomposition of *time tagged event* data,
%     or *binned* data into blocks with statistically uniform "rate".
%     Note: this does not have to be used on a lightcurve, but instead
%     could be used, e.g., on a spectrum to identify potential lines.  
%     The only requirement is that the data be sequential, with an
%     independent variable (e.g., time, wavelength, etc.), and a dependent
%     variable (e.g., counts, rate).
%
%     This function calculates the log of the probability function (or 
%     fitness function) for a given decomposition.  Called as a subroutine
%     from "global_optimum".
%
% CALLING SEQUENCE:
%     log_prob = log_post(data_pops, data_size, ncp_prior, type, alpha, beta, lo);
%
% INPUTS:
%     data_pops  : Array with cell populations (i.e., number of "events" 
%                  per cell)
%     data_size  : Array with cell sizes (i.e., number of "ticks" 
%                  per cell; can contain window function/efficiency)
%     ncp_prior  : log parameter for number of change points.  It is
%                  assumed that the prior probability of # of change
%                  points goes as gamma^N_cp, then:
%	              ncp_prior = -log(gamma)
%                  MAN suggests: ncp_prior = log(# of bins in lightcurve/
%                                             maximum # of desired cells)
%                  ncp_prior=7 is ~ 10^-3 significance for each new block
%     type       : Event data (type = 1 or 2), or binned data (type = 3)
%                  type = 1: cell is half-way from previous event to
%                            half-way to next event
%                  type = 2: cell is from first event to right before next
%                            event.
%                  (1 or 2 makes no difference in this subroutine) 
%     alpha, beta: For binned mode, the prior on the bin rate goes as
%                  (Lambda)^(alpha-1) Exp(-beta*Lambda), where Lambda = rate
%                  No error checking to make sure that they're sensible
%                  choices!
%
% OUTPUTS:
%     log_prob   : The logarithm of the probability function for a given
%                  decomposition of the lightcurve.
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     log_gamma  : Logarithm of the gamma function, included above.
%	           Eventually to be replaced by a GSL routine
% 
% MODIFICATION HISTORY:
%     Based upon Jeff Scargle's MATLAB routines to perform a 
%     Bayesian blocks decomposition of an ordered array.
%
%     Paper: "Studies in Astronomical Time Series Analysis. VI.
%             Optimum Partition of the Interval: Bayesian Blocks,
%	      Histograms, and Triggers", J.D. Scargle et al., 2003
%             (to be submitted).
%		
%     Web Site:  http://astrophysics.arc.nasa.gov/~jeffrey
%
%     Version 0.1.00: 2003/02/25 Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.1.01: 2003/02/26 MAN - Changed 'jdl' lgamma to log_gamma
%     Version 0.1.02: 2003/02/27 MAN - D'oh! bug in log_gamma fixed.
%     Version 0.2.00: 2003/05/26 MAN - Added alpha, beta to input.
%     Version 0.3.00: 2003/12/12 MAN - Changed prior for event mode to
%                                      (alpha+1) (1-p_1)^alpha
%     Version 0.4.00  2004/03/30 MAN - Added 'lo' variable, and toggles
%                                      on alpha, beta variables which allow
%                                      various choices for prior in event mode
%     Version 0.4.01  2004/04/06 MAN - Fixes to probabilities
%     Version 0.4.02  2004/05/10 MAN - Erroneous print statement removed

   % Note: no error checking, since this subroutine shouldn't be called 
   % directly, only through "sitar_global_optimum();"

   variable log_prob;

   % If time tagged event data ...

   if ( type == 1 or type == 2 )
   {
      variable arg, ii;

      arg = data_size - data_pops;
      ii = where( arg > 0 );

      % Slightly different tack than Scargle here.  His MATLAB code
      % sets this to ~0.  I set this M=N in Eq. (23) of Scargle (1998,
      % ApJ, 504, pp. 405-418).  Neither is strictly correct, but
      % I don't know if it makes a substantial difference.  (If it does,
      % you should be using the binned version instead.)  Furthermore,
      % instead of marginalizing over a uniform prior for p_1, 
      % choose P(p_1) = (alpha+1) (1-p_1)^alpha

      if(alpha == 0.)
      {
         % lo will be set to 1 / not 1 based upon beta = 0 / not 0
         % (done in sitar_global_optimum)

         log_prob = log_gamma( data_size + 1./lo ) 
                  - log_gamma( data_size + 1./lo + 1 ) - log(lo);

         log_prob[ii] = log_gamma( data_pops[ii] + 1 ) 
                      + log_gamma( arg[ii] + 1./lo  )
                      - log_gamma( data_size[ii] + 1./lo + 1 )
                      - log(lo);
      }
      else if(beta == 0.)
      {
         log_prob = log_gamma( data_size + 1 ) + log_gamma( alpha + 1 )
                  - log_gamma( data_size + alpha + 2 ) + log(alpha+1);

         log_prob[ii] = log_gamma( data_pops[ii] + 1 ) 
                      + log_gamma( arg[ii] + alpha + 1 )
                      - log_gamma( data_size[ii] + 2 + alpha )
                      + log( alpha + 1. );
      }
      else
      {

         % The following represent my own thoughts for this case.

         variable ealpha = @data_size, ebeta = @ealpha;
         ealpha[*] = exp(-alpha);
         ebeta[*] = exp(-beta);
      
         % This is only an approximate expression, to make up for
         % my lack of integration skills
         
         log_prob = data_size * log(1.-ebeta) + log(beta)
                    - log( beta - alpha );

         % This, however, I believe, is correct.  It *will not*
         % evaluate correctly if the GSL module is not initialized

         log_prob[ii] = log_gamma( data_size[ii] - data_pops[ii] ) +
                        log_gamma( data_pops[ii] + 1 ) -
                        log_gamma( data_size[ii] + 1 ) +
                        log( 
                            beta_inc(arg[ii],data_pops[ii]+1,ealpha[ii]) -
                            beta_inc(arg[ii],data_pops[ii]+1,ebeta[ii])  
                            )
                        - log( beta - alpha );
      }
   }

   % If it is instead binned data ...

   else
   {
      % Note that one can modify this to take account of window functions, 
      % exposures, etc.

      log_prob = log_gamma( data_pops + alpha ) 
                    - ( data_pops + alpha ) * log( data_size + beta )
                    + alpha*log( beta ) - log_gamma( alpha );
   }

   % This goes in whether it's event or binned data.

   log_prob = log_prob - ncp_prior;

   return log_prob;   
}

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

private define glob_opt_use()
{

%       Version 2.4.00: 2006/01/21 MAN 

   message(" ");
   message("Usage: ");
   message(" ");
   message(" ans = sitar_global_optimum( cell, ncp_prior, type [, first, alpha, beta] );");
   message(" ");
   message(" Inputs:");
   message("   cell: A structure, with cell.pops, cell.size, cell.lo_t");
   message("         cell.hi_t,cell.dtcor (see sitar_make_data_cells)");
   message("   ncp_prior: Parameter for prior on number of 'blocks'"); 
   message("   type: Identical to type from sitar_make_data_cells");
   message("   first: (Optional) first > 0, code runs in 'trigger' mode,");
   message("          and returns upon the first sign of a change, so long");
   message("          as it is greater than 'first' cells from the");
   message("          beginning of the lightcurve");
   message("   alpha,: EVENT MODE: ");
   message("   beta     alpha=0, beta=0 (default): implies prior on p_1 (probability");
   message("              of one or more photons in a frame) is uniform from 0->1.");
   message("            alpha!=0, beta=0: prior is (1+alpha) (1-p_1)^alpha");
   message("            alpha=0, beta!=0: prior on *rate* is exp(-Lambda/Lambda_0), ");
   message("              where \Lambda_0 is the mean rate from the entire observation");
   message("            alpha!=0,beta!=0: prior on *rate* is uniform over bounds");
   message("              given by alpha and beta."); 
   message("            alpha = beta > 0: then auto choose bounds to be 1/3 & 3 times"); 
   message("              the *mean* rate of the input lightcurve");
   message("           BINNED MODE: ");
   message("            The prior on the bin rate (sort of) goes as");
   message("            (Lambda)^(alpha-1) Exp(-beta*Lambda), where Lambda=rate");
   message(" ");
   message(" Outputs:");
   message("   results.cpt: Array of change point locations for the maximum"); 
   message("                likelihood solution (indices are specific to the");
   message("                cell input);");
   message("   results.last, .best: (Diagnostic purposes only) Arrays of the");
   message("                        location of the last change point and the");
   message("                        associated maximum log probability");
   message("   results.cts: Counts in each block");
   message("   results.rate: Rate in each block");
   message("   results.err: (Gherels) poisson error for the block rate");
   message("   results.lo_t, .hi_t: Times of lower (>=) and upper (<) block");
   message("                        boundaries.");
   message(" ");
   return;
}

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

public define sitar_global_optimum()
{

% NAME: 
%     global_optimum
%
% PURPOSE:
%     Used in Bayesian decomposition of *time tagged event* data,
%     or *binned* data into blocks with statistically uniform "rate".
%     Note: this does not have to be used on a lightcurve, but instead
%     could be used, e.g., on a spectrum to identify potential lines.  
%     The only requirement is that the data be sequential, with an
%     independent variable (e.g., time, wavelength, etc.), and the 
%     dependent variable (i.e., counts)
%
%     This is the main routine to determine the optimal partitioning of 
%     the sequential data.
%
% CALLING SEQUENCE:
%     ans = sitar_global_optimum( cell, ncp_prior, ev_type, first );
%
% INPUTS:
%     cell       : Structure with cell.pops (i.e., number of "events" 
%                  per cell), cell.size (i.e., number of "ticks" per cell),
%                  cell.lo_t, cell.hi_t (i.e., time boundaries of cell),
%                  cell.dtcor (i.e., "dead time correction" info for cells)
%     ncp_prior  : log parameter for number of change points.  It is
%                  assumed that the prior probability of # of change
%                  points goes as (N_cp)^-gamma, then:
%		      ncp_prior = log(gamma)
%                  MAN suggests: ncp_prior = log(# of bins in lightcurve/
%                                             maximum # of desired cells)
%                  ncp_prior=7 is ~ 10^-3 significance for each new block
%                  NOTE: The latter is a potentially *very bad* frequentist
%                  interpretation, and should be verified via Monte Carlo
%                  simulations.
%     ev_type    : Event (ev_type = 1 or 2), or binned data (ev_type = 3)
%                  ev_type = 1: cell is half-way from previous event to
%                            half-way to next event
%                  ev_type = 2: cell is from first event to right before next
%                            event.
%                  (1 or 2 treated the same in this subroutine) 
%
% OPTIONAL INPUTS:
%
%     first      : Toggle, =0 for retrospective mode (i.e., analyze whole 
%                  lightcurve), >0 for trigger mode (i.e., return at first
%                  sign of a change, so long as it's greater than 'first'
%                  cells from the beginning.)
%     alpha, beta: For binned mode, the prior on the bin rate goes as
%                  (Lambda)^(alpha-1) Exp(-beta*Lambda), where Lambda=rate
%                   Event mode- alpha=0, beta=0 (default) implies prior on p_1
%                  (probability of one or more photons in a frame) is uniform
%                  from 0->1. alpha!=0, beta=0 prior is (1+alpha) (1-p_1)^alpha
%                  alpha=0, beta!=0 prior on *rate* is exp(-Lambda/Lambda_0), 
%                  where \Lambda_0 is the mean rate from the entire observation
%                  alpha!=0,beta!=0 prior on *rate* is uniform over bounds
%                  given by alpha and beta
%                  alpha = beta > 0 => Autochoose the bounds on uniform rate
%                  NOTE: No error checking to make sure that any of these 
%                  are sensible choices!
%
% OUTPUTS:
%     ans.cpt    : Array of change point locations
%     ans.last   : Indices of last change point location for \
%                  lightcurve of length equal to array index  \__\ Diagnostic
%     ans.best   : Maximum log probability values for         /  / Purposes 
%                  lightcurve of length equal to array index /     Only
%     ans.cts    : Counts in each block
%     ans.rate   : Rate in each block
%     ans.err    : (Gherels) poisson error on rate in each block
%     ans.lo_t   : Time of lower boundary (>=) of block
%     ans.hi_t   : Time of upper boundary (<) of block
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     log_post() : Associated subroutine released with this package
%     reverse()  : ISIS function to reverse the order of an array
%     cumsum()   : ISIS (>= v1.0.48) function to give the cumulative 
%                  sum of an array
% 
% MODIFICATION HISTORY:
%     Based upon Jeff Scargle's MATLAB routines to perform a 
%     Bayesian blocks decomposition of an ordered array.
%
%     Paper: "Studies in Astronomical Time Series Analysis. VI.
%             Optimum Partition of the Interval: Bayesian Blocks,
%	      Histograms, and Triggers", J.D. Scargle et al., 2003
%             (to be submitted).
%		
%     Web Site:  http://astrophysics.arc.nasa.gov/~jeffrey
%
%     Version 0.1.00: 2003/02/25 Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.1.01: 2003/02/26 MAN Corrected the output for the case where
%                                no changepoints are found.
%     Version 0.1.02: 2003/04/11 MAN Changed default to first index (but not
%                                necessarily last index) is changepoint 
%     Version 0.1.03: 2003/04/15 MAN added `variable i;' initialization
%     Version 0.1.04: 2003/04/22 MAN fixed a bug wherein sometimes the
%                                the first change point wasn't set to 0
%     Version 0.2.00: 2003/05/20 MAN Cleaned up output for no changepoints
%     Version 0.3.00: 2003/05/26 MAN Added output of blocked lightcurve
%	                         counts, rates, errors, and time bins.
%                                Initial thoughts about GTI's, and some
%                                style changes.
%     Version 0.3.01: 2003/08/22 MAN Cleaned up some error checking,
%                                added glob_opt_use();
%     Version 0.3.02: 2003/12/11 MAN minor tweak in use of glob_opt_use();
%                                and on error check of ncp_prior
%     Version 0.4.00: 2004/03/30 MAN alpha,beta now also affect event mode,
%                                as described above
%     Version 0.4.01: 2004/04/06 MAN alpha=beta > 0 defaults changed

   % Some really basic error checking for the input

   variable cell, ev_type, ncp_prior, first, alpha, beta, log_prob, lo;
   switch ( _NARGS )
   {
    case 0:
    glob_opt_use();
    return;
   }
   {
    case 3:
      ( cell, ncp_prior, ev_type ) = ();
      first = 0;
      alpha = NULL;
      beta = NULL;
   }

   {
    case 4:
      ( cell, ncp_prior, ev_type, first ) = ();
      alpha = NULL;
      beta = NULL;
   }

   {
    case 5:
      ( cell, ncp_prior, ev_type, first, alpha ) = ();
      beta = NULL;
   }

   {
    case 6:
      ( cell, ncp_prior, ev_type, first, alpha, beta ) = ();
   }

   {
      message(" ");
      message("Incorrect number of arguments.");
      glob_opt_use();
      _pop_n(_NARGS);
      return;
   }

   if ( is_struct_type(cell) )
   {
      variable fields = get_struct_field_names(cell);
      if ( length( where( fields == "pops" or fields == "size" or
                          fields == "lo_t" or fields == "hi_t" or
                          fields == "dtcor" ) ) != 5 )
      {
         message(" ");
         message("Fields on cell structure not properly set.");
         message(" ");
         return;
      }
   }
   else
   {
      message(" ");
      message("`cell' was not a structure");
      message(" ");
      return;
   }

   if (
        typeof(cell.pops) != Array_Type or 
        typeof(cell.size) != Array_Type or 
        typeof(cell.lo_t) != Array_Type or 
        typeof(cell.hi_t) != Array_Type or 
        typeof(cell.dtcor)  != Array_Type 
      ) 
   {
      message(" ");
      message("Error: All fields of structure 'cell' must be Array_Type");
      message(" ");
      return;
   }

   variable num_cells = length(cell.pops);

   if ( num_cells < 3 ) 
   {
      vmessage("Error: length of lightcurve is %d. Needs >= 3.", num_cells);
      return;
   }

   if (
        num_cells != length(cell.size) or 
        num_cells != length(cell.lo_t) or
        num_cells != length(cell.hi_t) or
        num_cells != length(cell.dtcor) 
      ) 
   {
      message(" ");
      message("Error: All arrays in cell structure must have same length");
      message(" ");
      return;
   }

   if ( typeof(ncp_prior) != Double_Type and typeof(ncp_prior) != Float_Type ) 
   {
      if ( typeof(ncp_prior) == Integer_Type )
      {
         ncp_prior = 1. * ncp_prior;
      }
      else
      {
         message(" ");
         vmessage("Error: 'ncp_prior' is type %s",string(typeof(ncp_prior)));
         message("Should be: Float_Type (suggested trial value: 6.)");
         message(" ");
         return;
      }
   }

   if ( ev_type != 1 and ev_type != 2 and ev_type != 3 ) 
   {
      message(" ");
      vmessage("Error: 'ev_type' set to %d, but must be 1, 2, or 3", ev_type);
      message(" ");
      return;
   }

   if (first == NULL) { first = 0; }

   if (typeof(first) != Double_Type or typeof(first) != Float_Type ) 
   {
      first = typecast(first,Long_Type);
   }

   if (typeof(first) != Integer_Type and typeof(first) != Long_Type ) 
   {
      message(" ");
      vmessage("Error: 'first' has type of %s",string(typeof(first)));
      message("Should be: Integer_Type (=0 or >0)");
      message(" ");
      return;
   }

   if ( first != 0 and first < 0 )
   {
      message(" ");
      vmessage("Error: 'first' set to %d, but must be=0, or >0",first);
      message(" ");
      return;
   }

   if(alpha == NULL) 
   { 
      if( ev_type != 3 )
      {
         alpha = 0.; 
      }
      else
      {
         alpha = 1.;
      }
   }

   if(beta == NULL)  
   { 
      if( ev_type != 3 )
      {
         beta = 0.; 
      }
      else
      {
         beta = 1.;
      }
   }

   if (typeof(alpha) != Integer_Type or typeof(alpha) != Long_Type ) 
   {
      alpha = typecast(alpha,Float_Type);
   }

   if (typeof(beta) != Integer_Type or typeof(beta) != Long_Type ) 
   {
      beta = typecast(beta,Float_Type);
   }

   if ( (typeof(alpha) != Double_Type and typeof(alpha) != Float_Type
         and alpha < 0. ) or
        (typeof(beta) != Double_Type and typeof(beta) != Float_Type
         and beta < 0. ) )
   {
      message(" ");
      message("Error: alpha and beta must be real numbers >= 0.");
      message(" ");
      return;
   }

   if( ev_type != 3 and alpha == 0. and beta != 0. )
   {
      lo = sum(cell.pops)/sum(cell.size);
   }
   else
   {
      lo = 1.;
   }

   if( ev_type != 3 and alpha == beta and beta > 0.)
   {
      alpha = sum(cell.pops)/sum(cell.size)/3.;
      beta = alpha*9.;
   }

   variable ans = struct{cpt,last,best,cts,rate,err,lo_t,hi_t};
   variable iw, intvl;
   variable best=[0.], last=[0];

   % s-lang doesn't do NULL arrays quite the same as MATLAB, so we have to 
   % start with the first bin.

   best[0] =  max( log_post( cell.pops[[0:0]],cell.size[[0:0]],
                             ncp_prior, ev_type, alpha, beta, lo ) );

   % Now onto the meat of the job. Slightly longer than the one line in MATLAB.

   variable besttest;
   variable i, irange;
   for( i=1; i<= num_cells-1 ; i++ )
   {
      irange = [i:0:-1];
      besttest =  [0.,best] +
                  reverse(
                            log_post( 
                                      cumsum( cell.pops[irange] ),
                                      cumsum( cell.size[irange] ),
                                      ncp_prior, ev_type, alpha, beta, lo
                                    )
                         );

      best = [best, max(besttest)];
      last = [last, min( where( besttest == best[i] ) )];

      % Is this trigger mode?

      if ( first and last[i] ) 
      {
         ans.best = best;
         ans.last = last;

         ans.cpt = [ 0, last[i] ];
	 ans.cts = Double_Type[2];
	 ans.rate = @ans.cts;
	 ans.err = @ans.cts;
	 ans.lo_t = @ans.cts;
	 ans.hi_t = @ans.cts;

         variable ia = Array_Type[2];
         ia[0] = [ 0: last[i]-1 ];
         ia[1] = [ last[i]: i ];
         i = 0;
	 loop(2)
         {
            ans.cts[i] = sum( cell.pops[ ia[0] ] );
            ans.err[i] = ( 1. + sqrt(0.75 + ans.cts[i] ) ); 

            ans.lo_t[i] = cell.lo_t[ min( ia[0] ) ];
            ans.hi_t[i] = cell.hi_t[ max( ia[0] ) ];

            intvl = sum( ( cell.hi_t[ ia[0] ] - cell.lo_t [ ia[0] ] ) 
                         * cell.dtcor[ ia[0] ] );
            ans.rate[i] = ans.cts[i] / intvl;
            ans.err[i] = ans.err[i] / intvl;

	    i++;
         }
         ans.cts = typecast(ans.cts,Integer_Type);
         return ans;
      }
   }

   % Find the optimum partition

   variable index=last[num_cells-1];
   ans.cpt = [index];
   index = last[index-1];	

   % Since s-lang doesn't handle NULL arrays the same as MATLAB,
   % need to do this first check

   while ( index )
   {
      ans.cpt = [ index, ans.cpt ];
      index = last[index-1];
   }

   % Always return the first changepoint index as 0, if not already 0
   % (as opposed to MATLAB, the latter happens because of how we 
   % initialized 'last').

   if ( ans.cpt[0] )
   {
      ans.cpt = [ 0, ans.cpt ];
   }
   ans.best=best;
   ans.last=last;

   % Make those useful things for plotting, etc.

   variable cpt_length = length(ans.cpt);
   ans.rate = Double_Type[cpt_length];
   ans.err = @ans.rate;
   ans.cts = @ans.rate;
   ans.lo_t = @ans.rate;
   ans.hi_t = @ans.rate;

   ans.lo_t[0] = cell.lo_t[0];
   ans.hi_t[cpt_length-1]=cell.hi_t[length(cell.hi_t)-1];
   i = 0;
   loop ( cpt_length-1 )
   {
      iw = [ ans.cpt[i]:ans.cpt[i+1] - 1 ];

      ans.cts[i] = sum( cell.pops[iw] );
      ans.err[i] = ( 1. + sqrt(0.75 + ans.cts[i] ) ); 

      ans.hi_t[i] = cell.hi_t[ ans.cpt[i+1] - 1 ];
      ans.lo_t[i+1] = cell.lo_t[ ans.cpt[i+1] ];
%
      intvl = sum( ( cell.hi_t[iw] - cell.lo_t [iw] ) * cell.dtcor[iw] );
      ans.rate[i] = ans.cts[i] / intvl;
      ans.err[i] = ans.err[i] / intvl;

      i++;
   }

   % Mop up that endpoint

   iw = [ ans.cpt[i]:length(cell.hi_t) - 1 ];

   ans.cts[i] = sum( cell.pops[iw] );
   ans.err[i] = ( 1. + sqrt(0.75 + ans.cts[i] ) ); 

   intvl = sum( ( cell.hi_t[iw] - cell.lo_t [iw] ) * cell.dtcor[iw] );
   ans.rate[i] = ans.cts[i] / intvl;
   ans.err[i] = ans.err[i] / intvl;

   ans.cts = typecast(ans.cts,Integer_Type);
   return ans;
}

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

private define readasm_use()
{

%       Version 2.0.00: 2003/12/08 MAN

   message(" ");
   message("Usage: ");
   message(" ");
   message(" event = sitar_readasm( file [, tstart, tstop, maxchi2, chnl, mjd] );");
   message(" ");
   message("   Variables omitted or set to 0 take on default values.");
   message("   Variables in [] are optional, but are order specific.");
   message(" ");
   message(" Inputs:");
   message("   file:  Root of asm filename, e.g. 'xa_cygx1_d1'");
   message("   tstart: Start time to be read (units of MET or MJD)");
   message("   tstop: Stop time to be read (units of MET or MJD)");
   message("   maxchi2: The maximum reduced chi2 for an 'ASM solution' which");
   message("            will be accepted to consider a data point");
   message("   chnl: =0 for total band (source.lc files = default)");
   message("        !=0 for three ASM colors + total (source.col files)"); 
   message("   mjd: !=0 changes from default of Mission Elapsed Time (days)");
   message("        to MJD = JD - 2,400,000.5");
   message(" ");
   message(" Outputs:");
   message("   event.time: Time of each event");
   message("   event.rate: Total ASM count-rate");
   message("   event.err: Uncertainty in count rate");
   message("   event.chi2: Chi2 of ASM solution");
   message(" ");
   message(" Optional Outputs:");
   message("   event.[ch1,ch2,ch3]_rate: ASM count rate in each channel");
   message("   event.[ch1,ch2,ch3]_err: ASM count rate error in each channel");
   message("   event.[ch1,ch2,ch3]_chi2: Chi2 for ASM solution in each channel");
   message("         (Overrides event.chi2, which won't be output)");
   message(" ");

   return;
}

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

public define sitar_readasm( )
{

% NAME:
%     sitar_readasm
%
% PURPOSE:
%     Read an ASM 1-dwell lightcurve
%
% CALLING SEQUENCE:
%     event  = sitar_readasm (file,tstart,tstop,maxchi2,chnl,mjd)
% 
% INPUTS:
%     file     : Root of asm filename, e.g. 'xa_cygx1_d1', leaving off
%                the '.lc' or '.col'
%
% OPTIONAL INPUTS:
%     tstart   : Start of time-interval to be read (units of MET or MJD)
%                (Default = file beginning)
%     tstop    : Stop of time-interval to be read (units of MET or MJD)
%                (Default = file ending)
%     maxchi2  : The maximum reduced chi2 for an "ASM solution"
%                which will be accepted to consider a data point
%     chnl     : 0 for total band (source.lc files = default)
%                !=0 for three ASM colors + total (source.col files)
%                Energy bands are as follows:
%                   ch1: ASM channel  410...1188 = 1.3 -  3.0 keV
%                   ch2: ASM channel 1189...1860 = 3.0 -  5.0 keV
%                   ch3: ASM channel 1861...4750 = 5.0 - 12.2 keV
%                Note: One does *not* have to pre-process the 
%                the ASM color file into separate channel files
%                using the ftools 'asmchannel'.  This will work directly
%                on the 'source.col' files.
%     mjd      : !=0 changes from default of Mission Elapsed Time (days)
%                to MJD = JD - 2,400,000.5 
%
% OUTPUT:
%     event.time : Time of each event
%     event.rate : Total ASM count-rate 
%     event.err  : Uncertainty in count rate
%     event.chi2 : Chi2 of ASM solution
%
% OPTIONAL OUTPUTS:
%     event.[ch1,ch2,ch3]_rate : ASM count rate in each channel
%     event.[ch1,ch2,ch3]_err  : ASM count rate error in each channel
%     event.[ch1,ch2,ch3]_chi2 : Chi2 for ASM solution in each channel
%                         (Overrides event.chi2, which won't be output)
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     None        
%
% MODIFICATION HISTORY:
%     Version 0.1.00, 2003/09/04, Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.2.00, 2006/01/21, MAN, Fixed bug in parsing of min/max times,
%                                 pointed out by Elisa Resconi, and minor changes
%                                 to error messages.

   variable file;
   variable tstart=-1.e16, tstop=1.e16, maxchi2 = 1.e6, chnl = 0, mjd = 0;
   variable event;
   switch (_NARGS)
   {
    case 1:
      ( file ) = ();
   }

   {
    case 2:
      ( file, tstart ) = ();
   }

   {
    case 3:
      ( file, tstart, tstop ) = ();
   }

   {
    case 4:
      ( file, tstart, tstop, maxchi2 ) = ();
   }

   {
    case 5:
      ( file, tstart, tstop, maxchi2, chnl ) = ();
   }

   {
    case 6:
      ( file, tstart, tstop, maxchi2, chnl, mjd ) = ();
   }

   {
     readasm_use();
      _pop_n(_NARGS);
      return NULL;
   }

   if ( orelse{tstart == NULL}{tstart == 0.})
   {
      tstart = -1.e16;
   }
   if ( orelse{tstop == NULL}{tstop == 0.})
   {
      tstop = 1.e16;
   }
   if ( orelse{maxchi2 == NULL}{maxchi2 == 0.})
   {
      maxchi2 = 1.e6;
   }

   if(chnl == NULL)
   {
      chnl = 0;
   }

   if ( mjd == NULL )
   {
      mjd = 0;
   }

   file = strtrans(file," ","");

   % Define filnenames

   if ( chnl != 0 ) 
   {
      file = file + ".col";
   }
   else
   {
      file = file + ".lc";
   }

   variable t, rate, error, chi2, minchan, maxchan;

   (t,rate,error,chi2,minchan,maxchan) = 
      fits_read_col(file,"TIME","RATE","ERROR","RDCHI_SQ","MINCHAN","MAXCHAN");

   variable mjdrefi = fits_read_key(file,"MJDREFI");
   variable mjdrefr = fits_read_key(file,"MJDREFF");
   mjdrefr = mjdrefr + mjdrefi;

   if ( mjd != 0 )     % Convert Mission Elapsed time to MJD if Flag Set
   {
      t = t + mjdrefr;
   }

   variable mint = min(t), maxt = max(t);
   if( maxt < tstart or mint > tstop )
   {
      message("\n    No times specified within [tstart,tstop]");
      vmessage("       [tstart,tstop] = [%e,%e]",tstart,tstop);   
      vmessage("       [min(t),max(t)] = [%e,%e]\n",mint,maxt); 
      return NULL;
   }

   variable mincolor = [ 410, 410,1189,1861];
   variable maxcolor = [4750,1188,1860,4750];
   variable ndx;

   if ( chnl != 0 )
   {
      event = struct{time, rate, ch1_rate, ch2_rate, ch3_rate, 
                     err, ch1_err, ch2_err, ch3_err, 
                     ch1_chi2, ch2_chi2, ch3_chi2};

      variable ndxa=where(minchan==mincolor[1] and maxchan==maxcolor[1]);
      variable ndxb=where(minchan==mincolor[2] and maxchan==maxcolor[2]);
      variable ndxc=where(minchan==mincolor[3] and maxchan==maxcolor[3]);

      variable ta = t[ndxa], tb = t[ndxb], tc = t[ndxc];
      variable iw = where( ta == tb and tb == tc );

      if( length(iw) != length(ta) )
      {
	message("\n    Color channels not of equal length\n");
	return NULL;
      }

      event.time   = t[ndxa];

      event.ch1_rate =  rate[ndxa];
      event.ch1_err  = error[ndxa];
      event.ch1_chi2 =  chi2[ndxa];

      event.ch2_rate =  rate[ndxb];
      event.ch2_err  = error[ndxb];
      event.ch2_chi2 =  chi2[ndxb];

      event.ch3_rate =  rate[ndxc];
      event.ch3_err  = error[ndxc];
      event.ch3_chi2 =  chi2[ndxc];

      ndx = where( event.ch1_chi2 <= maxchi2 and
                   event.ch2_chi2 <= maxchi2 and
                   event.ch3_chi2 <= maxchi2 and
		   event.time >= tstart and event.time <=tstop );
      if( length(ndx) == 0 )
      {
        vmessage("\n    No solutions with chi2 < %f found\n", maxchi2);
        return NULL;
      }

      event.time   =   event.time[ndx];
      event.ch1_rate = event.ch1_rate[ndx];
      event.ch1_err  =  event.ch1_err[ndx];
      event.ch1_chi2 = event.ch1_chi2[ndx];

      event.ch2_rate = event.ch2_rate[ndx];
      event.ch2_err  =  event.ch2_err[ndx];
      event.ch2_chi2 = event.ch2_chi2[ndx];

      event.ch3_rate = event.ch3_rate[ndx];
      event.ch3_err  =  event.ch3_err[ndx];
      event.ch3_chi2 = event.ch3_chi2[ndx];

      event.rate = event.ch1_rate + event.ch2_rate + event.ch3_rate;
      event.err = sqrt( event.ch1_err^2 + 
                        event.ch2_err^2 + 
                        event.ch3_err^2 );
   }
   else
   {
      event = struct{time, rate, err, chi2};

      ndx = where( minchan==mincolor[0] and maxchan==maxcolor[0]);

      if( length(ndx) == 0 )
      {
        message("\n    Empty lightcurve!\n");
        return NULL;
      }
      event.time =     t[ndx];
      event.rate =  rate[ndx];
      event.err  = error[ndx];
      event.chi2 =  chi2[ndx];

      ndx = where( event.chi2 <= maxchi2 and
		   event.time >= tstart and event.time <= tstop );

      if( length(ndx) == 0 )
      {
        vmessage("\n    No solutions with chi2 < %f found\n", maxchi2);
        return NULL;
      }

      event.time = event.time[ndx];
      event.rate = event.rate[ndx];
      event.err  =  event.err[ndx];
      event.chi2 = event.chi2[ndx];
   }

   return event;
}

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

private define reb_rat_use()
{

%       Version 2.1.00: 2006/10/24 MAN

   message(" ");
   message("Usage: ");
   message(" ");
   message(" lc = sitar_rebin_rate(t,rate [,err,dt,tstart,tstop,minbin,gap,cust,wght]);");
   message(" ");
   message("   Variables omitted or set to Integer_Type 0 take on defaults.)");
   message("   Variables in [] are optional, but are order specific.");
   message(" ");
   message(" Inputs:");
   message("   t:  Discrete times at which rates are measured");
   message("   rate: (Presumed GTI & Exposure Corrected) rates");
   message("   err: Error on the rates");
   message("   dt: Width of evenly spaced time bins");
   message("   tstart: Beginning of first output time bin");
   message("   tstop: End of last output time bin");
   message("   minbin: Minimum required events per bin, else the bin is");
   message("           set to 0, or ignored. Default=1.");
   message("   gap: = 0 - Empty bins are set to 0 (default)");
   message("       != 0 - Empty bins are deleted from the lightcurve");
   message("   cust.bin_lo: Lower time bounds for user defined grid");
   message("   cust.bin_hi: Upper time bounds for a user defined grid");
   message("                (Overrides dt, tstart, and tstop inputs)");
   message("   wght: Default = 0.  If != 0, and err !=0, mean is weighted by ");
   message("         the input error.  I.e., mean = (\sum rate/err^2)/(\sum 1/err^2),");
   message("         and the output variance becomes ");
   message("         (\sum (mean-rate/err^2)^2)/(\sum 1/err^2)^2");
   message(" ");
   message(" Outputs:");
   message("   lc.rate: Mean rate of resulting rebinning");
   message("   lc.bin_lo: Lower time bounds of rebinned lightcurve");
   message("   lc.bin_hi: Upper time bounds of rebinned lightcurve");
   message("   lc.num: Number of events going into a bin");
   message("   lc.var: Variance of the rate in a time bin. (Only calculated");
   message("           where lc.num >= 2, otherwise set to zero.)"); 
   message(" ");
   message(" Optional Outputs:");
   message("   lc.err: Rate errors combined in quadrature, i.e., it's ");
   message("           sqrt{\sum{ err^2 }}/N, where N is the number of points");
   message("           in that particular bin (i.e., lc.num).  Only computed if"); 
   message("           err is input [otherwise, just use sqrt(lc.var)].");
   message(" ");

   return;
}

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

public define sitar_rebin_rate( )
{

% NAME:
%     sitar_rebin_rate
%
% PURPOSE:
%     Bin a rate lightcurve to a new scale
%
% CALLING SEQUENCE:
%     lc = sitar_rebin_rate (t,rate,err,dt,tstart,tstop,minbin,gap,cust,wght)
% 
% INPUTS:
%     t           : Discrete times at which rates are measured
%     rate        : (Presumed GTI & Exposure Corrected) rates 
%
% OPTIONAL INPUTS:
%     err         : Error on the rates
%     dt          : Width of evenly spaced time bins
%                   (Defaults to 1/100 of length of lightcurve)
%     tstart      : Beginning of first time bin 
%                   (Defaults to time of first event)
%     tstop       : End of last time bin 
%                   (Defaults to time of last event)
%     minbin      : Minimum required events per bin, else the
%                   bin is set to 0, or ignored. Default=1. 
%     gap         :  = 0 - Empty bins are set to 0 (default)
%                   != 0 - Empty bins are deleted from the lightcurve
%     cust.bin_lo : Lower time bounds for a user defined binning grid 
%     cust.bin_hi : Upper time bounds for a user defined binning grid 
%                   Overrides dt, tstart inputs
%     wght        : Default = 0.  If != 0, and err !=0, mean is weighted by 
%                   the input error.  I.e., mean = (\sum rate/err^2)/(\sum 1/err^2),
%                   and the output variance becomes 
%                   (\sum (mean-rate/err^2)^2)/(\sum 1/err^2)^2
%
% OUTPUT:
%     lc.rate     : Mean rate of resulting rebinning
%     lc.bin_lo   : Lower time bounds of rebinned lightcurve
%     lc.bin_hi   : Upper time bounds of rebinned lightcurve
%     lc.num      : Number of events going into a bin
%     lc.var      : Variance of the rate in a time bin. ** Only
%                   calculated where lc.num >= 2 **, otherwise
%                   set to zero.
%
% OPTIONAL OUTPUTS:
%     lc.err      : Rate errors combined in quadrature, i.e., it's 
%                   sqrt{\sum{ err^2 }}.  No dividing by N is done! 
%                   Only computed if err is input [otherwise, just use
%                   sqrt(lc.var)].  
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     histogram   : ISIS (!! >= v1.1.5 !!) function to create a histogram 
%     linear_grid : ISIS function to create a linear grid of numbers
%     mean        : ISIS function to find the mean of an array
%     moment      : ISIS function to find statistical moments of an array
%
% MODIFICATION HISTORY:
%     Version 0.1.00, 2003/09/04, Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.2.00, 2006/10/24, Michael A. Nowak (mnowak@alum.mit.edu)
%                     Added weight parameter, and fixed bug in calculation 
%                     of lc.err
%     Version 0.3.00, 2006/11/01, Michael A. Nowak (mnowak@alum.mit.edu)
%                     Fixed a bug with default on cust variable
%     Version 0.4.00, 2007/10/15, Michael A. Nowak (mnowak@alum.mit.edu)
%                     Changed the way the grids are set up.


   variable t, rate; 
   variable err=0, dt=0., tstart=0, tstop=0, minbin=1, gap=0, custf=0,
            cust=NULL, wght=0;
   variable istck;
   variable lc;

   switch (_NARGS)
   {
    case 2:
      ( t, rate ) = ( );
   }

   {
    case 3:
      ( t, rate, err ) = ( );
   }

   {
    case 4:
      ( t, rate, err, dt ) = ( );
   }

   {
    case 5:
      ( t, rate, err, dt, tstart ) = ( );
   }

   {
    case 6:
      ( t, rate, err, dt, tstart, tstop ) = ( );
   }

   {
    case 7:
      ( t, rate, err, dt, tstart, tstop, minbin ) = ( );
   }

   {
    case 8:
      ( t, rate, err, dt, tstart, tstop, minbin, gap ) = ( );
   }

   {
    case 9:
      ( t, rate, err, dt, tstart, tstop, minbin, gap, cust ) = ( );
   }
   {
    case 10:
      ( t, rate, err, dt, tstart, tstop, minbin, gap, cust, wght ) = ( );
   }
   {
      reb_rat_use();
      _pop_n(_NARGS);
      return;
   }

   if ( is_struct_type(cust) == 1 )
   {
      variable fields = get_struct_field_names(cust);
      if ( length( where( fields == "bin_lo" or fields == "bin_hi" ) ) != 2 )
      { 
         reb_rat_use();
         message("    Incorrect fields on custom bins. Require: ");
         message("       variable cust = struct{ bin_lo, bin_hi };");
         message(" ");
         _pop_n(_NARGS);
         return;
      }
      else
      {
         custf = 1;
      }
   }
   else
   { 
      if ( cust != NULL )
      {
         reb_rat_use();
         _pop_n(_NARGS);
         return;
      }
   }

   variable lt = length(t), lr = length(rate);
   if (lt != lr )
   {
      reb_rat_use();
      message("    Incommensurate array lengths for time and rate");
      message(" ");
      return;
   }

   if( gap == NULL )
   {
      gap = 0;
   }

   if ( minbin == NULL )
   {
      minbin = 0;
   }

   if ( gap !=0 and minbin == 0 ) { minbin = 1; };

   if ( dt == NULL )
   {
      dt = 0.;
   }

   variable mint = min(t);
   variable maxt = max(t);

   if ( Integer_Type == typeof(tstart) or tstart == NULL )
   {
      if (tstart)
      {
         tstart = typecast(tstart, Double_Type);
      }
      else
      {
         tstart = mint;
      }
   }

   if ( Integer_Type == typeof(tstop) or tstop == NULL )
   {
      if (tstop)
      {
         tstop = typecast(tstop, Double_Type);
      }
      else
      {
         tstop = maxt;
      }
   }

   if ( wght == NULL )
   {
      wght = 0.;
   }

   if( maxt < tstart or mint > tstop or mint > maxt)
   {
      reb_rat_use( );
      message("    Data outside of bounds of start & stop times:");
      vmessage("       [tstart,tstop] = [%e,%e]",tstart,tstop);
      vmessage("       [min_t, max_t] = [%e,%e]",mint,maxt);
      message(" ");
      return;
   }

   variable errf = length(err);
   if( errf <= 1 )
   {
      lc = struct{ rate, bin_lo, bin_hi, num, var };
      errf = 0;
   }
   else
   {
      if (errf != lr )
      {
         reb_rat_use();
         message("    Incommensurate array length for error array");
         message(" ");
         return;
      }
      lc = struct{ rate, bin_lo, bin_hi, num, var, err };
      errf = 1;
   }

   if( dt <= 0. )
   {
      ( lc.bin_lo, lc.bin_hi ) = linear_grid(tstart,tstop,100);
   }
   else
   {
      variable nfrac = ( tstop - tstart )/dt;
      variable nbins = int(nfrac);
      lc.bin_lo = [0:nbins-1]*dt+tstart;
      lc.bin_hi = [1:nbins]*dt+tstart;
      if( nfrac-nbins > 0 )
      {
         lc.bin_lo = [lc.bin_lo,lc.bin_hi[nbins-1]];
         lc.bin_hi = [lc.bin_hi,tstop]; 
      }
   }
   
   if ( custf )
   {
      lc.bin_lo = @cust.bin_lo;
      lc.bin_hi = @cust.bin_hi;
   }

   variable rev;

   lc.num = histogram(t, lc.bin_lo, lc.bin_hi, &rev);
   lc.rate = Double_Type[length(lc.num)]; 
   lc.var = @lc.rate;
   if ( errf )
   {
      lc.err = @lc.rate;
   }

   variable i=0, mom, mome;
   loop( length(lc.num) )
   {
      if ( lc.num[i] )
      {
         if( errf !=0 and wght !=0 )
         {
            mom = moment( rate[rev[i]]/err[rev[i]] );
            mome = moment( 1/err[rev[i]] );
         }
         else
         {
            mom = moment( rate[rev[i]] );
            mome = @mom;
            mome.ave = 1.;
            mome.var = 1.;
         }

         lc.rate[i] = mom.ave/mome.ave;
	 lc.var[i] = mom.var/(mome.ave)^2;

	 if ( errf )
	 {
	    lc.err[i] = sqrt(sum( (err[rev[i]])^2 ))/lc.num[i];
	 }
      }
      i++;
   }
   
   variable idx;
   if ( minbin != 0 )
   {
      if ( gap != 0 )
      {
         idx = where( lc.num >= minbin );

	 if ( length(idx) )
         {
            lc.rate = lc.rate[idx];
	    lc.bin_lo = lc.bin_lo[idx];
	    lc.bin_hi = lc.bin_hi[idx];
	    lc.var = lc.var[idx];
	    lc.num = lc.num[idx];

            if ( errf )
            {
               lc.err = lc.err[idx];
            }
         }
	 else
	 {
            reb_rat_use( );
	    message("Minimum events requirement yields null lightcurve");
            message(" ");
	    return;
	 }
      }
      else
      {
         idx = where( lc.num < minbin );

         if ( length(idx) )
         {
            lc.rate[idx] = 0.;
            lc.var[idx] = 0.;

            if ( errf )
            {
               lc.err[idx] = 0.;
            }
         }
      }
   }

   return lc;
}

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

private define pfold_rate_use()
{

%       Version 2.0.00: 2003/12/08 MAN

   message(" ");
   message("Usage: ");
   message(" ");
   message(" profile  = sitar_pfold_rate(t,rate,p [,pdot,pddot,nphs,tzero]);");
   message(" ");
   message("   Variables omitted or set to 0 take on default values.");
   message("   Variables in [] are optional, but are order specific.");
   message(" ");
   message(" Inputs:");
   message("  ");
   message("   t: Times at which rates are measured");
   message("   rate: Lightcurve rate");
   message("   p: Period on which to fold the lightcurve");
   message("   pdot: Period derivative");
   message("   pddot: Period second derivative");
   message("   nphs: Number of phase bins in output folded lightcurve.");
   message("   tzero: Time of zero phase. Default = t[0] ");
   message(" ");
   message(" Outputs:");
   message("   profile.bin_lo: Start value of phase bin");
   message("   profile.bin_hi: Stop value of phase bin");
   message("   profile.mean: Mean rate in phase bin");
   message("   profile.var: Variance of rate in phase bin");
   message("   profile.sdm: Standard deviation of mean rate in a phase bin");
   message("   profile.num: Number of data points in phase bin");
   message(" ");
   return;
}

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

private define _sitar_pfold_profile(nphs)
{

%       Version 1.0.00: 2003/09/04 MAN
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     linear_grid : ISIS function to create a linear grid of numbers

   variable profile = struct{bin_lo, bin_hi, mean, var, sdm, num};
   (profile.bin_lo, profile.bin_hi) = linear_grid(0.,1.,nphs);
   profile.num  = Integer_Type[nphs];
   profile.mean = Double_Type[nphs];
   profile.var = Double_Type[nphs];
   profile.sdm = Double_Type[nphs];
   return profile;
}

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

private define _sitar_pfold_rate( t, rate, p, pdot, pddot, nphs, profile)
{

%       Version 1.1.00: 2006/06/27 MAN
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     histogram   : ISIS (!! >= v1.1.5 !!) function to create a histogram 
%     moment      : ISIS function to find statistical moments of an array

   % Convert time to phase

   variable phs;
   phs = p + ( pdot/2. + pddot*t/6. )*t ; 
   phs = ( t mod phs ) / phs;

   variable ndx = where( phs < 0. );
   if( length(ndx) > 0 ) 
   {
      phs[ndx] = phs[ndx]+1.;
   }

   variable rev;

   profile.num = histogram(phs, profile.bin_lo, profile.bin_hi, &rev);

   variable i, j, mom;
   _for(0, nphs-1, 1 )
   {
      i = ();
      j = rev[i];
      if ( length(j) )
      {
         mom = moment( rate[j] );
         profile.mean[i] = mom.ave;
         profile.var[i] = mom.var;
         profile.sdm[i] = mom.sdom;
      } 
   }

   return profile;
}

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

private define _sitar_pfold_rate_deriv( t, rate, p, pdot, nphs, profile)
{

%       Version 1.1.00: 2006/06/27 MAN
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     histogram   : ISIS (!! >= v1.1.5 !!) function to create a histogram 
%     moment      : ISIS function to find statistical moments of an array

   % Convert time to phase

   variable phs;
   phs = p + ( pdot/2. )*t ; 
   phs = ( t mod phs ) / phs;

   variable ndx = where( phs < 0. );
   if( length(ndx) > 0 ) 
   {
      phs[ndx] = phs[ndx]+1.;
   }

   variable rev;

   profile.num = histogram(phs, profile.bin_lo, profile.bin_hi, &rev);

   variable i, j, mom;
   _for(0, nphs-1, 1 )
   {
      i = ();
      j = rev[i];
      if ( length(j) )
      {
         mom = moment( rate[j] );
         profile.mean[i] = mom.ave;
         profile.var[i] = mom.var;
         profile.sdm[i] = mom.sdom;
      } 
   }

   return profile;
}

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

private define _sitar_pfold_rate_noderiv( t, rate, p, nphs, profile)
{

%       Version 1.0.00: 2003/09/04 MAN
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     histogram   : ISIS (!! >= v1.1.5 !!) function to create a histogram 
%     moment      : ISIS function to find statistical moments of an array

   % Convert time to phase


   t = ( t mod  p ) / p ;

   variable ndx = where( t < 0. );
   if( length(ndx) > 0 ) 
   {
      t[ndx] = t[ndx]+1.;
   }

   variable rev;

   profile.num = histogram(t, profile.bin_lo, profile.bin_hi, &rev);

   variable i, j, mom, mome;
   _for(0, nphs-1, 1 )
   {
      i = ();
      j = rev[i];
      if ( length(j) )
      {
         mom = moment( rate[j] );
         mome = @mom;
         profile.mean[i] = mom.ave;
         profile.var[i] = mom.var;
         profile.sdm[i] = mom.sdom;
      } 
   }

   return profile;
}

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

public define sitar_pfold_rate( )
{

% NAME:
%     sitar_pfold_rate
%
% PURPOSE:
%     Fold a rate lightcurve on a given period (plus first and second derivatives)
%
% CALLING SEQUENCE:
%     profile  = sitar_pfold_rate(t,rate,p,pdot,pddot,nphs,tzero);
% 
% INPUTS:
%     t        : Times at which rates are measured
%     rate     : Lightcurve rate
%     p        : Period on which to fold the lightcurve
%
% OPTIONAL INPUTS:
%     pdot     : Period derivative
%     pddot    : Period second derivative
%     nphs     : Number of phase bins in output folded lightcurve
%                Default = 20
%     tzero    : Time of zero phase
%                Default = t[0]
%
% OUTPUT:
%     profile.bin_lo : Start value of phase bin
%     profile.bin_hi : Stop value of phase bin
%     profile.mean   : Mean rate in phase bin
%     profile.var    : Variance of rate in phase bin
%     profile.sdm    : Standard deviation of the mean rate in a phase bin
%     profile.num    : Number of data points in phase bin
%
% OPTIONAL OUTPUTS:
%     NONE
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     pfold_rate_use            : Attached above; gives usage notes
%     _sitar_pfold_profile      : Attached above; creates data structure
%     _sitar_pfold_rate         : Attached above; fold with all derivatives
%     _sitar_pfold_rate_deriv   : Attached above; fold with first derivative
%     _sitar_pfold_rate_noderiv : Attached above; fold with no derivatives
%
% MODIFICATION HISTORY:
%     Version 0.1.00, 2003/09/04, Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.2.00, 2003/09/08, MAN, after suggestions by John Davis
%                                 Parameter setup occurs here, meat of
%                                 calculation occurs in subroutines

   variable profile;
   variable t, rate, p, pdot=0., pddot=0., nphs=20, tzero=-1.e18;

   switch (_NARGS)
   {
    case 3:
      ( t, rate, p ) = ( );
   }

   {
    case 4:
      ( t, rate, p, pdot ) = ( );
   }

   {
    case 5:
      ( t, rate, p, pdot, pddot ) = ( );
   }

   {
    case 6:
      ( t, rate, p, pdot, pddot, nphs ) = ( );
   }

   {
    case 7:
      ( t, rate, p, pdot, pddot, nphs, tzero ) = ( );
   }

   {
      pfold_rate_use();
      _pop_n(_NARGS);
      return;
   }

   variable ltime = length(t), lrate = length(rate);
   if( ltime < 2 or lrate < 2 or ltime != lrate )
   {
      pfold_rate_use();
      message("    Time and rate must be of equal length > 2");
      message(" ");
      return;
   }

   if ( pdot == NULL )
   {
      pdot = 0.;
   }
   if ( pddot == NULL )
   {
      pddot = 0.;
   }
   if ( orelse{ nphs == NULL }{ nphs < 1 }{nphs > ltime} )
   {
      nphs = min([20,ltime]);
   }
   else
   {
      nphs = typecast(nphs,Integer_Type);
   }
   if ( orelse{ tzero == NULL }{ tzero < -0.9e18 } )
   {
      tzero = t[0];
   }

   t = t - tzero;

   profile = _sitar_pfold_profile(nphs);

   switch ( 2*(pddot!=0)+(pdot!=0) )
   {
    case 0:
      profile = _sitar_pfold_rate_noderiv( t, rate, p, nphs, profile);
   }
   {
    case 1:
      profile = _sitar_pfold_rate_deriv( t, rate, p, pdot, nphs, profile);
   }
   {
      profile = _sitar_pfold_rate( t, rate, p, pdot, pddot, nphs, profile);
   }

   return profile;
}

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

#ifexists beta_inc
   gsl_set_error_disposition (GSL_EUNDRFLW, 0);
#endif

#ifnexists beta_inc
   private define beta_inc(x,y,z)
   {
      return 10.;
   }
#endif

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

private define fstat(v,dfn,dfd)
{

%       Version 1.0.00: 2003/09/04 MAN

   % Modelled after IDL f_pdf, which computes the probabilty (p) such that:
   % Probability(X <= v) = p, where X is a random variable from the F 
   % distribution with (dfn) and (dfd) degrees of freedom.  beta_inc, 
   % is taken from the Gnu Scientific Library

   variable x, logq=10.;
   if ( v > 0 )
   {  
      logq = beta_inc( dfd/2.0, dfn/2.0, dfd/( dfd + dfn*v ) );
      logq = log10( logq );
      return logq;
   }
   else
   {
      return 1.;
   }
}

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

private define epfold_rate_use()
{

%       Version 2.0.00: 2003/12/08 MAN

   message(" ");
   message("Usage: ");
   message(" ");
   message(" fold  = sitar_epfold_rate(t,rate,pstart,pstop [,nphs,nsrch,prds]);");
   message(" ");
   message("   Variables omitted or set to 0 take on default values.");
   message("   Variables in [] are optional, but are order specific.");
   message(" ");
   message(" Inputs:");
   message("   t: Times at which rates are measured");
   message("   rate: Lightcurve rate");
   message("   pstart: Start value of periods to search");
   message("   pstop: Stop value of periods to search");
   message("   nphs: Number of phase bins in folded lightcurve");
   message("   nsrch: abs(nsrch) = number of periods to search.");
   message("          If > 1, linear search grid");
   message("          If < 1, logarithmic search grid");
   message("   prds: Alternative array of periods to search");
   message(" ");
   message(" Outputs:");
   message("   fold.prd: Array of trial periods");
   message("   fold.lstat: Array of L-statistic for trial period");
   message("   fold.lsig: Array of log10 of single trial significance levels.");
   message("              (Note: requires use of the GNU Scientific Library"); 
   message("              package, otherwise, defaults to a value of 1.)");
   message("   fold.nphs: Array of number of phases (<= nphs) at a given");
   message("              period that had data.");
   message(" ");
   return;
}

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

public define sitar_epfold_rate( )
{

% NAME:
%     sitar_epfold_rate
%
% PURPOSE:
%     Perform an epoch fold test on a rate lightcurve
%
% CALLING SEQUENCE:
%     fold = sitar_epfold_rate(t,rate,pstart,pstop,nphs,nsrch,prds);
% 
% INPUTS:
%     t     : Times at which rates are measured
%     rate     : Lightcurve rate
%     pstart   : Starting value of periods to search
%     pstop    : Stopping value of periods to search
%
% OPTIONAL INPUTS:
%     nphs     : Number of phase bins in folded lightcurve
%                Default = 20, minimum = 2 (unless prds is set)
%     nsrch    : abs(nsrch) = number of periods to search.
%                If > 1, linear search grid
%                If < 1, logarithmic search grid
%                Default = linear, with nsrch = 
%     prds     : Alternative array of periods to search 
%
% OUTPUT:
%     fold.prd    : Array of trial periods
%     fold.lstat  : Array of L-statistic for trial period
%     fold.lsig   : Array of log10 of *single trial* significance levels
%     fold.nphs   : == nphs, number of phases searched.
%
% OPTIONAL OUTPUTS:
%     NONE
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     linear_grid               : ISIS function to create a linear grid 
%                                 of numbers
%     _sitar_pfold_profile      : Subroutine included with this package
%     _sitar_pfold_rate_noderiv : Subroutine included with this package
%
% MODIFICATION HISTORY:
%     Version 0.1.00, 2003/09/04, Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.2.00, 2003/09/10, MAN, following suggestions of John Davis
%                                 separating out calculations from 
%                                 input parameter tests

   variable fold = struct{prd,lstat,lsig,nphs};
   variable t, rate, pstart, pstop,  nphs=20, nsrch=20, prds=-1.e20;

   switch (_NARGS)
   {
    case 4:
      ( t, rate, pstart, pstop ) = ( );
   }

   {
    case 5:
      ( t, rate, pstart, pstop, nphs ) = ( );
   }

   {
    case 6:
      ( t, rate, pstart, pstop, nphs, nsrch ) = ( );
   }

   {
    case 7:
      ( t, rate, pstart, pstop, nphs, nsrch, prds ) = ( );
   }

   {
      epfold_rate_use();
      _pop_n(_NARGS);
      return;
   }

   variable ltime = length(t), lrate = length(rate);
   if ( ltime < 2 or lrate < 2 or ltime != lrate )
   {
      epfold_rate_use();
      message("    Time and rate must be of equal length > 2");
      message(" ");
      return;
   }

   if ( pstart == NULL or pstop == NULL )
   {
      epfold_rate_use();
      message("    Start and stop periods must be set.");
      message(" ");
      return;
   }

   if ( orelse{ nphs == NULL }{ nphs < 1 }{ nphs > ltime } )
   {
      nphs = min([20,ltime]);
   }
   else
   {
      nphs = typecast(nphs,Integer_Type);
   }

   if ( orelse{ nsrch == NULL }{ nsrch == 0 } )
   {
      nsrch = 20;
   }
   else
   {
      nsrch = typecast(nsrch,Integer_Type);
   }

   if ( orelse{ prds == NULL }{ min(prds) < -0.9e20 } )
   {  
      variable a,b;
      if ( sign(nsrch) < 0 )
      {
         if( sign(pstart) != sign(pstop) )
         {
            epfold_rate_use();
            message("    Cannot use logarithmic spacing if pstart/pstop have different signs");
            message(" ");
	    return;
         }
         else
         {
            (a,b) = linear_grid(0.,log10(pstop/pstart),abs(nsrch)-1);
	    prds = pstart*(10.^[a[0],b]);
         }
      }
      else
      {
         (a,b) = linear_grid(pstart,pstop,abs(nsrch)-1);
         prds = [a[0],b];
      }
   }
   else
   {
      if ( typeof(prds) != Array_Type )
      {
         epfold_rate_use();
         message("    Optional periods 'prds' must be 'Array_Type'");
         message(" ");
	 return;
      }
   }

   variable i, prof, lprds=length(prds);
   variable dltime=1.*ltime-nphs;
   variable factor = dltime / ( nphs -1. );

   % No loss of generality with subtracting the mean

   rate = rate - mean(rate);
   fold.lstat = @Double_Type[lprds];
   fold.lsig = @fold.lstat;
   fold.prd = prds;
   factor =  ( ltime - nphs ) / ( nphs -1. );

   prof = _sitar_pfold_profile(nphs);

   _for(0, lprds-1, 1 )
   {
      i = ();
      prof = _sitar_pfold_rate_noderiv(t,rate,prds[i],nphs,prof);

      fold.lstat[i] = factor * ( sum( prof.num * prof.mean^2 ) ) /
                      ( sum( prof.var * ( prof.num - 1. ) ) );

      fold.lsig[i] = fstat( fold.lstat[i], nphs-1, dltime );
   }

   fold.nphs = nphs;

   return fold;
} 

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

% If the GSL module is loaded, it's fft will overwrite ISIS's, and 
% we have to use a different normalization.

private variable gsl_fft = 0;

#ifexists _gsl_fft_complex

gsl_fft = 1;

#endif

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

public define sitar_avg_cpd()
{

% NAME: 
%     sitar_avg_cpd
%
% PURPOSE:
%     Take two evenly spaced lightcurves (presumed counts vs. time), and calculate
%     the PSD and CPD in segments of length l, averaged over the whole lightcurve.
%     Segments with data gaps are skipped.  Deadtime corrections to Poisson
%     noise are *not* made. Poisson noise is *not* subtracted from average PSD.
%
% CALLING SEQUENCE:
%     (f,psd_a,psd_b,cpd,navg,avg_cnts_a,avg_cnts_b) = 
%                                   sitar_avg_cpd(cnts_a,cnts_b,l,[tbin,times]);
%
% INPUTS:
%     cnts_a/b    : Arrays of total counts in each time bin
%     l           : Length of individual PSD segments (use a power of 2!!!)
%     tbin        : Length of evenly spaced bins (default == 1)
%     times       : Times of measurements (otherwise presumed to have no gaps)
%     norm        : Determine PSD normalization (see below)
%  
% OUTPUTS:
%     f           : PSD frequencies ( == 1/Input Time Unit)
%     psd_a/b     : Average PSD 
%     navg        : Number of data segments going into the average
%     avg_cnts_a/b: Average number of counts per segment of length l
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     fft     : ISIS implementation of an FFT (see ISIS manual for provenance)
% 
% MODIFICATION HISTORY:
%
%     Version 0.1.00: 2004/06/02 Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.2.00: 2005/08/17 MAN  Added a toggle to take care of FFT
%                                     if GSL module is installed

   variable j,s,f,l,k,alen,ntrans;
   variable transa,transb,stransa,stransb;
   variable freq,atrans=0,btrans=0,ctrans=0.+0.i,ntot=0;
   variable mcounta=0,mcountb=0;
   variable dt, wdt, nwdt, wwdt;

   variable counts_a, counts_b, len, tbin = 1., bin_time = NULL, donorm = 1.;

   switch(_NARGS)
   {
    case 3:
      (counts_a, counts_b, len) = ();
   }
   {
    case 4:
      (counts_a, counts_b, len, tbin) = ();
   }
   {
    case 5:
      (counts_a, counts_b, len, tbin, bin_time) = ();
   }
   {
    case 6:
      (counts_a, counts_b, len, tbin, bin_time, donorm) = ();
   }
   {
      variable str = 
     [" (f,psd_a,psd_b,cpd,navg,avg_a,avg_b) = sitar_avg_cpd(cnts_a,cnts_b,l,[tbin,times,norm]);"," ",
      "   Take an evenly spaced lightcurve (presumed counts vs. time), and calculate",
      "   the PSD and CPD in segments of length l, averaged over the whole lightcurve.",  
      "   Segments with data gaps are skipped.  Deadtime corrections to Poisson",
      "   noise are *not* made. Poisson noise is *not* subtracted from average PSDs."," ",
      " Inputs:",
      "   cnts_a/b: Arrays of total counts in each time bin",
      "   l       : Length of individual PSD segments (use a power of 2!!!)",
      "   tbin    : Length of evenly spaced bins (default == 1)",
      "   times   : Times of measurements (otherwise presumed to have no gaps)",
      "   norm    : Determine PSD normalization (see below)"," ",
      " Outputs:",
      "   f       : PSD frequencies ( == 1/Input Time Unit)",
      "   psd_a/b : Average PSDs (`Leahy normalization', i.e., Poisson noise == 2,",
      "             if norm=1 [default]; `rms' or `Belloni-Hasinger' or `Miyamoto'", 
      "             normalization, i.e., PSD == (rms)^2/Hz & noise== 2/Rate," ,
      "             otherwise)",
      "   cpd     : Normalized Cross Power Spectral Density (complex array)",
      "   navg    : Number of data segments going into the averages",
      "   avg_a/b : Average number of counts per segment of length l"];
      msg(str);
      return;
   }

   if( donorm == NULL ){ donorm = 1.; }

   len = typecast(len,Long_Type);

   alen = len/2;

   % Lengths of the counts and time arrays, and the number
   % of transforms we will expect to be able to do

   variable lcounts_a, lcounts_b, ltime;

   lcounts_a = length(counts_a);
   lcounts_b = length(counts_b);

   if(bin_time != NULL)
   {
      ltime = length(bin_time);
   }
   else
   {
      bin_time = [0:lcounts_a-1]*tbin;
      ltime = lcounts_a;
   }

   if(lcounts_a != ltime or lcounts_b != lcounts_b) 
   { 
      message("Arrays have unequal lengths");
      return; 
   }

   if(lcounts_a < len) 
   { 
      message("Lightcurve too short for desired segment size");
      return; 
   }

   dt = bin_time - shift(bin_time,-1);
   wdt = where( dt > 1.5*tbin );

   if(length(wdt) == 0)
   {
      ntrans = lcounts_a/len;

      _for(1,ntrans,1)
      {
         j = ();
         s = (j-1)*len;
         f = j*len -1;

         transa = counts_a[[s:f]];
         stransa = sum(transa);
         mcounta = mcounta + stransa;

         transb = counts_b[[s:f]];
         stransb = sum(transb);
         mcountb = mcountb + stransb;
         
         transa = shift( fft(transa,1) , 1 );
         atrans = atrans + 2*( abs(transa[[0:alen-1]]) )^2 / stransa;

         transb = shift( fft(transb,1) , 1 );
         btrans = btrans + 2*( abs(transb[[0:alen-1]]) )^2 / stransb;

         ctrans = ctrans + 2*transa*Conj(transb)/sqrt(stransa*stransb);
 
         ntot = ntot + 1;
      }
   }
   else
   {
      nwdt = length(wdt) + 1;
      wwdt = [0,wdt,lcounts_a-1];

      _for(1,nwdt,1)
      {
         k = ();
         ntrans = ( wwdt[k] - wwdt[k-1] )/len;
            
         _for(1,ntrans,1)
         {
            l = ();
            s = (l-1) * len + wwdt[k-1];
            f = l * len + wwdt[k-1] - 1;

            transa = counts_a[[s:f]];
            stransa = sum(transa);
            mcounta = mcounta + stransa;

            transb = counts_b[[s:f]];
            stransb = sum(transb);
            mcountb = mcountb + stransb;
         
            transa = shift( fft(transa,1) , 1 );
            atrans = atrans + 2*( abs(transa[[0:alen-1]]) )^2 / stransa;

            transb = shift( fft(transb,1) , 1 );
            btrans = btrans + 2*( abs(transb[[0:alen-1]]) )^2 / stransb;

            ctrans = ctrans + 2*transa*Conj(transb)/sqrt(stransa*stransb);

            ntot = ntot + 1;
         }
      }
   }

   if( ntot == 0 )
   {
      message("No transforms performed!");
      return;
   }

   freq = [1:alen] / (len * tbin);
   mcounta = mcounta / ntot;
   mcountb = mcountb / ntot;

   % I believe that this will give "Leahy Normalization", i.e.,
   % Poisson noise = 2, in absence of deadtime effects.

   atrans = (atrans * len) / ntot;
   btrans = (btrans * len) / ntot;
   ctrans = (ctrans * len) / ntot;

   ctrans = ctrans[[0:alen-1]];

   if( donorm != 1. )
   {
      atrans = atrans * tbin*len / mcounta;
      btrans = btrans * tbin*len / mcountb;
      ctrans = ctrans * tbin*len / sqrt(mcounta*mcountb);
   }

   if( gsl_fft == 1)
   {
      atrans = atrans/len;
      btrans = btrans/len;
      ctrans = ctrans/len;
   }

   return freq, atrans, btrans, ctrans, ntot, mcounta, mcountb;
}

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

public define sitar_avg_psd()
{

% NAME: 
%     sitar_avg_psd
%
% PURPOSE:
%     Take an evenly spaced lightcurve (presumed counts vs. time), and calculate
%     the PSD in segments of length l, averaged over the whole lightcurve.
%     Segments with data gaps are skipped.  Deadtime corrections to Poisson
%     noise are *not* made. Poisson noise is *not* subtracted from average PSD.
%
% CALLING SEQUENCE:
%     (f,psd,navg,avg_cnts) = sitar_avg_psd(cnts,l,[tbin,times]);
%
% INPUTS:
%     cnts    : Array of total counts in each time bin
%     l       : Length of individual PSD segments (use a power of 2!!!)
%     tbin    : Length of evenly spaced bins (default == 1)
%     times   : Times of measurements (otherwise presumed to have no gaps)
%     norm    : Determine PSD normalization (see below)
%  
% OUTPUTS:
%     f       : PSD frequencies ( == 1/Input Time Unit)
%     psd     : Average PSD 
%     navg    : Number of data segments going into the average
%     avg_cnts: Average number of counts per segment of length l
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     fft     : ISIS implementation of an FFT (see ISIS manual for provenance)
% 
% MODIFICATION HISTORY:
%
%     Version 0.1.00: 2004/06/02 Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.2.00: 2005/08/17 MAN  Added a toggle to take care of FFT
%                                     if GSL module is installed

   variable j,s,f,l,k,alen,ntrans,trans,strans,freq,atrans=0,ntot=0,mcount=0;
   variable dt, wdt, nwdt, wwdt;

   variable counts, len, tbin = 1., bin_time = NULL, donorm = 1.;

   switch(_NARGS)
   {
    case 2:
      (counts, len) = ();
   }
   {
    case 3:
      (counts, len, tbin) = ();
   }
   {
    case 4:
      (counts, len, tbin, bin_time) = ();
   }
   {
    case 5:
      (counts, len, tbin, bin_time, donorm) = ();
   }
   {
      variable str = 
     [" (f,psd,navg,avg_cnts) = sitar_avg_psd(cnts,l,[tbin,times,norm]);"," ",
      "   Take an evenly spaced lightcurve (presumed counts vs. time), and calculate",
      "   the PSD in segments of length l, averaged over the whole lightcurve.",  
      "   Segments with data gaps are skipped.  Deadtime corrections to Poisson",
      "   noise are *not* made. Poisson noise is *not* subtracted from average PSD."," ",
      " Inputs:",
      "   cnts    : Array of total counts in each time bin",
      "   l       : Length of individual PSD segments (use a power of 2!!!)",
      "   tbin    : Length of evenly spaced bins (default == 1)",
      "   times   : Times of measurements (otherwise presumed to have no gaps)",
      "   norm    : Determine PSD normalization (see below)"," ",
      " Outputs:",
      "   f       : PSD frequencies ( == 1/Input Time Unit)",
      "   psd     : Average PSD (`Leahy normalization', i.e., Poisson noise == 2,",
      "             if norm=1 [default]; `rms' or `Belloni-Hasinger' or `Miyamoto'", 
      "             normalization, i.e., PSD == (rms)^2/Hz & noise== 2/Rate," ,
      "             otherwise)",
      "   navg    : Number of data segments going into the average",
      "   avg_cnts: Average number of counts per segment of length l"];
      msg(str);
      return;
   }

   if( donorm == NULL ){ donorm = 1.; }

   len = typecast(len,Long_Type);

   alen = len/2;

   % Lengths of the counts and time arrays, and the number
   % of transforms we will expect to be able to do

   variable lcounts, ltime;

   lcounts = length(counts);

   if(bin_time != NULL)
   {
      ltime = length(bin_time);
   }
   else
   {
      bin_time = [0:lcounts-1]*tbin;
      ltime = lcounts;
   }

   if(lcounts != ltime ) 
   { 
      message("Arrays have unequal lengths");
      return; 
   }

   if(lcounts < len) 
   { 
      message("Lightcurve too short for desired segment size");
      return; 
   }

   dt = bin_time - shift(bin_time,-1);
   wdt = where( dt > 1.5*tbin );

   if(length(wdt) == 0)
   {
      ntrans = lcounts/len;

      _for(1,ntrans,1)
      {
         j = ();
         s = (j-1)*len;
         f = j*len -1;

         trans = counts[[s:f]];
         strans = sum(trans);
         mcount = mcount + strans;
         
         trans = shift( fft(trans,1) , 1 );
         atrans = atrans + 2*( abs(trans[[0:alen-1]]) )^2 / strans;
 
         ntot = ntot + 1;
      }
   }
   else
   {
      nwdt = length(wdt) + 1;
      wwdt = [0,wdt,lcounts-1];

      _for(1,nwdt,1)
      {
         k = ();
         ntrans = ( wwdt[k] - wwdt[k-1] )/len;
            
         _for(1,ntrans,1)
         {
            l = ();
            s = (l-1) * len + wwdt[k-1];
            f = l * len + wwdt[k-1] - 1;

            trans = counts[[s:f]];
            strans = sum(trans);
            mcount = mcount + strans;

            trans = shift( fft(trans,1) , 1 );
            atrans = atrans + 2*( abs(trans[[0:alen-1]]) )^2 / strans;

            ntot = ntot + 1;
         }
      }
   }

   if( ntot == 0 )
   {
      message("No transforms performed!");
      return;
   }

   freq = [1:alen] / (len * tbin);
   mcount = mcount / ntot;

   % I believe that this will give "Leahy Normalization", i.e.,
   % Poisson noise = 2, *** in the absence of deadtime effects ***

   atrans = (atrans * len) / ntot;

   if( donorm != 1. )
   {
      atrans = atrans * tbin*len / mcount;
   }

   if( gsl_fft == 1)
   {
      atrans = atrans/len;
   }

   return freq, atrans, ntot, mcount;
}

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

public define sitar_lags()
{

% NAME: 
%     sitar_lags
%
% PURPOSE:
%     Given two input PSD (without noise subtracted), their CPD, and their
%     associated noise levels, all as functions of Fourier frequency, calculate
%     the time lag and coherence function (and associated errors)   vs. f
%
% CALLING SEQUENCE:
%     (lag,dlag,g,dg) = sitar_lags(freq,psda,psdb,cpd,noisea,noiseb,navg);
%
% INPUTS:
%     freq    : Fourier frequency array
%     psda/b  : Power spectra array associated with two lightcurves
%               (Noise subtraction *not* applied!)
%     cpd     : Cross power spectra array associated with two lightcurves
%     noisea/b: Power spectra noise level. Array *or* a constant.
%     navg    : Number of averages (over independent data segments and frequency
%               bins) associated with each input frequency bin. Can be a constant
%               vs. f (defaulted to 1).
%  
% OUTPUTS:
%     lag     : Time lag vs. frequency. Negative values mean psdb lags psda
%     dlag    : Associated error on the lag.
%     g       : Coherence function
%     dg      : Error on the coherence function
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     ones    : ISIS function to create an array of ones
% 
% MODIFICATION HISTORY:
%
%     For background discussion, see:
%
%        Vaughan, B. A. & Nowak, M. A.  1997, ApJ, 474, L43
%
%     and references therein (especially Bendat & Piersol 1986)
%
%     Version 0.1.00: 2004/06/03 Michael A. Nowak (mnowak@alum.mit.edu)

   variable freq,psda,psdb,cpd,noisea,noiseb,navg=1;

   switch(_NARGS)
   {
    case 6:
      (freq,psda,psdb,cpd,noisea,noiseb) = ();
   }
   {
    case 7:
      (freq,psda,psdb,cpd,noisea,noiseb,navg) = ();
   }
   {
      variable str = 
     [" (lag,dlag,g,dg) = sitar_lags(freq,psda,psdb,cpd,noisea,noiseb,[navg]);",
      " ",
      "   Given two input PSD (without noise subtracted), their CPD, and their",
      "   associated noise levels, all as functions of Fourier frequency, calculate",
      "   the time lag and coherence function (and associated errors)   vs. f"," ",
      " Inputs:",
      "   freq    : Fourier frequency array",
      "   psda/b  : Power spectra array associated with two lightcurves",
      "             (Noise subtraction *not* applied!)",
      "   cpd     : Cross power spectra array associated with two lightcurves",
      "   noisea/b: Power spectra noise level. Array *or* a constant.",  
      "   navg    : Number of averages (over independent data segments and frequency",
      "             bins) associated with each input frequency bin. Can be a constant",
      "             vs. f (defaulted to 1)."," ",
      " Outputs:",
      "   lag     : Time lag vs. frequency. Negative values mean psdb lags psda",
      "   dlag    : Associated error on the lag.",
      "   g       : Coherence function",
      "   dg      : Error on the coherence function"];
      msg(str);
      return;
   }  

   variable lfreq = length(freq), lnoisea = length(noisea), lnoiseb = length(noiseb),
            lnavg = length(navg);

   if( lnoisea==1 )
   {
      noisea = ones(lfreq)*noisea;
      lnoisea = lfreq;
   }

   if( lnoiseb==1 )
   {
      noiseb = ones(lfreq)*noiseb;
      lnoiseb = lfreq;
   }

   if( lnavg == 1 )
   {
      navg = ones(lfreq)*navg;
   }

   variable lpsda = length(psda), lpsdb = length(psdb), lcpd = length(cpd);

   if( lnoisea != lnoiseb or lnoiseb != lpsda or lpsda != lpsdb or
       lpsdb != lcpd or lcpd != lnavg or lnavg != lfreq )
   {
      message("Inconsistent array sizes");
      return;
   }

   variable sa, sb, ns, gs, dgs, gerr, lag, dlag, acpd = abs(cpd);

   % Define noise subtracted PSDs

   sa = psda - noisea;
   sb = psdb - noiseb;

   ns = ( sa*noiseb + sb*noisea + noisea*noiseb ) / navg;

   % Estimate of intrinsic coherence, high signal to noise

   gs = ( acpd^2 - ns ) / sa / sb;

   % Error on intrinsic coherence, squared, sans Poisson noise

   dgs = 2. * ( 1 - gs )^2 / gs^2;

   % Overall error estimate

   gerr = sqrt( ( 2. * ns^2 * navg / ( acpd^2 - ns )^2 +
                  ((noisea+noiseb)/2)^2 * (1/sa^2 + 1/sb^2) + dgs )
                / navg );

   % Calculate the magnitude of the phase lag
   % (Note: PI is a s-lang intrinsic)

   lag = Real(cpd) / acpd;
   lag = acos( lag ) / 2 / PI / freq * Imag(cpd)/abs( Imag(cpd) );

   % Phase lag errors

   dlag = acpd^2 / psda / psdb;
   dlag = sqrt(1-dlag) / sqrt( dlag * 2 * navg ) / 2 / PI / freq;

   return lag, dlag, gs, dgs;
}

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

public define sitar_lbin_cpd()
{

% NAME: 
%     sitar_lbin_cpd
%
% PURPOSE:
%     Logarithmically rebin two PSDs (and, optionally, their associated Poisson",
%     noise levels) and the associated CPD
%
% CALLING SEQUENCE:
%     (aflo,afhi,apsda,apsdb,acpd,nf,[anoisea,anoiseb]) = 
%         sitar_lbin_cpd(f,psda,psdb,cpd,dff,[noisea,noiseb,&rev]);
%
% INPUTS:
%     f        : Array of Fourier frequencies
%     psda/b   : Arrays of PSD values
%     acpd     : Array of Cross Power Spectral Density
%     dff      : \delta f/f value for logarithmic bin spacing
%     noisea/b : Array (*or single value*) of Poisson noise levels
%     rev      : If declared (i.e., `isis> variable rev;') and input,
%                returns the reverse indices for the binning.
%  
% OUTPUTS:
%     aflo     : Lower boundary of Fourier frequency bin
%     afhi     : Upper boundary of Fourier frequency bin
%     apsda/b  : Rebinned Power Spectral Density values
%     acpd     : Rebinned Cross Power Spectral Density
%     nf       : Number of frequencies going into the bin
%     anoisea/b: Rebinned noise level (array, even for single value input)
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     histogram   : ISIS function to create a histogram
%     make_hi_grid: ISIS function to create the upper bounds of a grid
%     ones        : ISIS function to create an array of ones
% 
% MODIFICATION HISTORY:
%
%     Version 0.1.00: 2004/06/02 Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.2.00: 2005/08/16 MAN - Added afreqlo/afreqhi outputs, and 
%                                      redid the gridding to proper log10

   variable freq, psda, psdb, cpd, dff, noisea = NULL, noiseb = NULL, rrev = NULL;

   switch(_NARGS)
   {
    case 5:
      (freq, psda, psdb, cpd, dff) = ();
   }
   {
    case 6:
      (freq, psda, psdb, cpd, dff, noisea) = ();
   }
   {
    case 7:
      (freq, psda, psdb, cpd, dff, noisea, noiseb) = ();
   }
   {
    case 8:
      (freq, psda, psdb, cpd, dff, noisea, noiseb, rrev) = ();
   }
   {
      variable str = 
     [" (aflo,afhi,apsda,apsdb,acpd,nf,[anoisea,anoiseb]) = ",
      "   sitar_lbin_cpd(f,psda,psdb,cpd,dff,[noisea,noiseb,&rev]);"," ",
      "   Logarithmically rebin two PSDs (and, optionally, their associated Poisson",
      "   noise levels) and the associated Cross Power Spectral Density",
      " ",
      " Inputs:",
      "   f        : Array of Fourier frequencies",
      "   psd      : Array of PSD values",
      "   dff      : \\delta f/f value for logarithmic bin spacing",
      "   noise    : Array (*or single value*) of Poisson noise levels",
      "   rev      : If declared (i.e., `isis> variable rev;') and input,",
      "              rev returns the reverse inidices for the binning (i.e.,",
      "              (apsda[i] = mean( psda[ rev[i] ] ), etc.)"," ", 
      " Outputs:",
      "   aflo     : Lower boundary of Fourier frequency bin",
      "   afhi     : Upper boundary of Fourier frequency bin",
      "   apsda/b  : Rebinned Power Spectral Density values",
      "   acpd     : Rebinned Cross Power Spectral Density values",
      "   nf       : Number of frequencies going into the bin",
      "   anoisea/b: Rebinned noise level (array, even for single value input)"];
      msg(str);

      return;
   }

   variable lf, lpsda, lpsdb, lcpd, lnoisea, lnoiseb;
   variable fmin, fmax, nmax, flo, fhi, iw, rev;
   variable nfreq, afreq, afreqlo, afreqhi, apsda, apsdb, acpd, 
            anoisea, anoiseb;

   lf = length(freq);
   lpsda = length(psda);
   lpsdb = length(psdb);
   lcpd = length(cpd);

   if( noisea != NULL )
   {
      lnoisea = length(noisea);
      if( lnoisea == 1 )
      {
         noisea = ones(lf) * noisea;
      }
      else if( lnoisea != lpsda )
      {
         message("Noise A and PSD arrays have unequal lengths");
         return; 
      }
   }

   if( noiseb != NULL )
   {
      lnoiseb = length(noisea);
      if( lnoisea == 1 )
      {
         noiseb = ones(lf) * noiseb;
      }
      else if( lnoiseb != lpsdb )
      {
         message("Noise B and PSD arrays have unequal lengths");
         return; 
      }
   }

   if( lf != lpsda or lpsda != lpsdb or lpsdb != lcpd )
   {
      message("Inconsistent Array Sizes");
      return;
   }

   fmin = min(freq);
   fmax = max(freq);

   nmax = typecast((log10(fmax)-log10(fmin))/log10(1+dff)+0.5,Integer_Type);

   flo = fmin * 10.^([0:nmax-1]*log10(1+dff));
   fhi = make_hi_grid(flo);

   nfreq = histogram(freq, flo, fhi, &rev);

   iw = where(nfreq !=0);

   variable liw = length(iw);

   nfreq = nfreq[iw];
   afreqlo = Double_Type[liw];
   apsda = @afreqlo;
   apsdb = @afreqlo;
   acpd  = Complex_Type[liw];

   rev = rev[iw];

   variable i;

   if( noisea == NULL and noiseb == NULL)
   {
      _for (0, liw-1, 1)
      {
         i = ();
         afreqlo[i] = min(freq[rev[i]]);
         apsda[i] = sum( psda[rev[i]] );
         apsdb[i] = sum( psdb[rev[i]] );
         acpd[i] = sum( cpd[rev[i]] );
      }
   } 
   else if( noiseb == NULL )
   {
      anoisea = @afreqlo;
      anoiseb = NULL;

      _for (0, liw-1, 1)
      {
         i = ();
         afreqlo[i] = min(freq[rev[i]]);
         apsda[i] = sum( psda[rev[i]] );
         apsdb[i] = sum( psdb[rev[i]] );
         acpd[i] = sum( cpd[rev[i]] );
         anoisea[i] = sum( noisea[rev[i]] );
      }

      anoisea = anoisea / nfreq;
   }
   else if( noisea == NULL )
   {
      anoiseb = @afreqlo;
      anoisea = NULL;

      _for (0, liw-1, 1)
      {
         i = ();
         afreqlo[i] = min(freq[rev[i]]);
         apsda[i] = sum( psda[rev[i]] );
         apsdb[i] = sum( psdb[rev[i]] );
         acpd[i] = sum( cpd[rev[i]] );
         anoiseb[i] = sum( noiseb[rev[i]] );
      }

      anoiseb = anoiseb / nfreq;
   }   
   else
   {
      anoisea = @afreqlo;
      anoiseb = @afreqlo;

      _for (0, liw-1, 1)
      {
         i = ();
         afreqlo[i] = min(freq[rev[i]]);
         apsda[i] = sum( psda[rev[i]] );
         apsdb[i] = sum( psdb[rev[i]] );
         acpd[i] = sum( cpd[rev[i]] );
         anoisea[i] = sum( noisea[rev[i]] );
         anoiseb[i] = sum( noiseb[rev[i]] );
      }

      anoisea = anoisea / nfreq;
      anoiseb = anoiseb / nfreq;
   }   


   afreqhi = make_hi_grid(afreqlo);
   afreqhi[liw-1] = max(freq[rev[liw-1]]);

   apsda = apsda / nfreq;
   apsdb = apsdb / nfreq;
   acpd = acpd / nfreq;

   if(rrev != NULL)
   {
     @rrev = rev;
   }

   if( noisea == NULL and noiseb == NULL )
   {
      return afreqlo, afreqhi, apsda, apsdb, acpd, nfreq;
   } 
   else
   {
      return afreqlo, afreqhi, apsda, apsdb, acpd, nfreq, 
             anoisea, anoiseb;
   }
}

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

public define sitar_lbin_psd()
{

% NAME: 
%     sitar_lbin_psd
%
% PURPOSE:
%     Logarithmically rebin a PSD (and, optionally, its associated Poisson noise)
%
% CALLING SEQUENCE:
%     (aflo,afhi,apsd,nf,[anoise]) = sitar_lbin_psd(f,psd,dff,[noise, &rev]);
%
% INPUTS:
%     f      : Array of Fourier frequencies
%     psd    : Array of PSD values
%     dff    : \delta f/f value for logarithmic bin spacing
%     noise  : Array (*or single value*) of Poisson noise levels
%     rev    : If declared (i.e., `isis> variable rev;') and input,
%              returns the reverse indices for the binning.
%  
% OUTPUTS:
%     aflo   : Lower frequency of bin
%     afhi   : Upper frequency of bin
%     apsd   : Rebinned PSD values
%     nf     : Number of frequencies going into the bin
%     anoise : Rebinned noise level (array, even for single value input)
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     histogram   : ISIS function to create a histogram
%     make_hi_grid: ISIS function to create the upper bounds of a grid
%     ones        : ISIS function to create an array of ones
% 
% MODIFICATION HISTORY:
%
%     Version 0.1.00: 2004/06/02 Michael A. Nowak (mnowak@alum.mit.edu)
%     Version 0.2.00: 2005/08/17 MAN - Added afreqlo/afreqhi outputs, and 
%                                      redid the gridding to proper log10

   variable freq, psd, dff, noise = NULL, rrev = NULL;

   switch(_NARGS)
   {
    case 3:
      (freq, psd, dff) = ();
   }
   {
    case 4:
      (freq, psd, dff, noise) = ();
   }
   {
     case 5:
      (freq, psd, dff, noise, rrev) = ();
   }
   {
      variable str = 
     [" (aflo,afhi,apsd,nf,[anoise]) = sitar_lbin_psd(f,psd,dff,[noise,&rev]);"," ",
      "   Logarithmically rebin a PSD (and, optionally, its associated Poisson noise)",
      " ",
      " Inputs:",
      "   f      : Array of Fourier frequencies",
      "   psd    : Array of PSD values",
      "   dff    : \\delta f/f value for logarithmic bin spacing",
      "   noise  : Array (*or single value*) of Poisson noise levels",
      "   rev    : If declared (i.e., `isis> variable rev;') and input,",
      "            rev returns the reverse inidices for the binning (i.e.,",
      "            (af[i] = mean( freq[ rev[i] ] ), etc.)"," ", 
      " Outputs:",
      "   aflo   : Lower bounds of Fourier frequency bins",
      "   afhi   : Upper bounds of Fourier frequency bins",
      "   apsd   : Rebinned PSD values",
      "   nf     : Number of frequencies going into the bin",
      "   anoise : Rebinned noise level (array, even for single value input)"];
      msg(str);

      return;
   }

   variable lf, lpsd, lnoise, fmin, fmax, nmax, flo, fhi, iw, rev;
   variable nfreq, lnf, afreq, afreqlo, afreqhi, apsd, anoise;

   lf = length(freq);
   lpsd = length(psd);

   if(lf != lpsd) 
   { 
      message("Frequency and PSD arrays have unequal lengths");
      return; 
   }

   if( noise != NULL )
   {
      lnoise = length(noise);
      if( lnoise == 1 )
      {
         noise = ones(lpsd)*noise;
      }
      else if( lnoise != lpsd )
      {
         message("Noise and PSD arrays have unequal lengths");
         return; 
      }
   }

   fmin = min(freq);
   fmax = max(freq);

   nmax = typecast((log10(fmax)-log10(fmin))/log10(1+dff)+0.5,Integer_Type);

   flo = fmin * 10.^([0:nmax-1]*log10(1+dff));
   fhi = make_hi_grid(flo);

   nfreq = histogram(freq, flo, fhi, &rev);

   iw = where(nfreq !=0);

   nfreq = nfreq[iw];
   lnf = length(nfreq);

   apsd = Double_Type[lnf];
   afreqlo = @apsd;
   afreqhi = @apsd;

   rev = rev[iw];

   variable i;

   if( noise == NULL )
   {
      _for (0, lnf-1, 1)
      {
         i = ();
         afreqlo[i] = min( freq[rev[i]]);
         apsd[i] = sum( psd[rev[i]] );
      }
   } 
   else
   {
      anoise = @apsd;

      _for (0, lnf-1, 1)
      {
         i = ();
         afreqlo[i] = min( freq[rev[i]]);
         apsd[i] = sum( psd[rev[i]] );
         anoise[i] = sum( noise[rev[i]] );
      }

      anoise = anoise / nfreq;
   }

   afreqhi = make_hi_grid(afreqlo);
   afreqhi[lnf-1] = max(freq[rev[lnf-1]]);

   apsd = apsd / nfreq;

    if(rrev != NULL)
    {
      @rrev = rev;
    }

   if( noise == NULL )
   {
      return afreqlo, afreqhi, apsd, nfreq;
   } 
   else
   {
      return afreqlo, afreqhi, apsd, nfreq, anoise;
   }
}

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

#ifexists isis_exit

public define sitar_define_psd()
{

% NAME: 
%     sitar_define_psd
%
% PURPOSE:
%     Define a PSD to be a fittable dataset in ISIS, with frequency 
%     values set to keV.
%
%
% CALLING SEQUENCE:
%     id = sitar_define_psd(f_lo,f_hi,psd,err,[noise]);
%
% INPUTS:
%     f_lo   : Low value of PSD frequency bin
%     f_hi   : High value of PSD frequency bin
%     psd    : Array of PSD values
%     err    : Array of PSD Errors
%     noise  : Array (*or single value*) of Poisson noise levels
%  
% OUTPUTS:
%     id     : Dataset Index
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
%     _A            : ISIS function to convert between Angstrom and keV
%     define_counts : ISIS function to define a counts spectrum
%     _define_bgd   : ISIS function to define a background counts spectrum
%     ones          : ISIS function to create an array of ones
% 
% MODIFICATION HISTORY:
%
%     Version 0.1.00: 2005/08/16 MAN - Added afreqlo/afreqhi outputs
%     Version 0.1.01: 2005/08/18 MAN - Added header info.

   variable f_lo, f_hi, psd, err, noise=NULL;

   switch(_NARGS)
   {
    case 4:
      (f_lo, f_hi, psd, err) = ();
   }
   {
    case 5:
      (f_lo, f_hi, psd, err, noise) = ();
   }   
   {
      variable str=
     [" id = sitar_define_psd(flo,fhi,psd,err,[noise]);"," ",
      "   Define psd to be an ISIS counts data set.",
      " ",
      " Inputs:",
      "   f_lo   : Low value of PSD frequency bin (Hz -> keV)",
      "   f_hi   : High value of PSD frequency bin (Hz -> keV)",
      "   psd    : Array of PSD values (Power -> Counts/bin)",
      "   err    : Array of PSD Errors",
      "   noise  : Array (*or single value*) of Poisson noise levels.",
      "            If undefined, background is undefined.",
      "  ",
      " Outputs:",
      "   id     : Dataset Index"];
      msg(str);
      return;
   }

   variable id;

   if(noise !=NULL)
   {
      variable lnoise = length(noise); 
      if( lnoise == 1 )
      {
         noise = ones(length(psd))*noise;
         id = define_counts(_A(f_hi),_A(f_lo),reverse(psd-noise),
                            reverse(err));
      }
      else if(lnoise != length(psd))
      {
         print("Noise array length differs from psd length");
         return;
      }
      else
      {
         id = define_counts(_A(f_hi),_A(f_lo),reverse(psd-noise),
                            reverse(err));
      }
   }
   else
   {
      id = define_counts(_A(f_hi),_A(f_lo),reverse(psd),reverse(err));
   }
   return id;
}

#endif

#ifexists sherpa_version

public define sitar_define_psd()
{

% NAME: 
%     sitar_define_psd
%
% PURPOSE:
%     Define a PSD to be a fittable dataset in Sherpa, with frequency 
%     values set to keV.
%
%
% CALLING SEQUENCE:
%     id = sitar_define_psd(f_lo,f_hi,psd,err,[noise]);
%
% INPUTS:
%     id     : Dataset Index
%     f_lo   : Low value of PSD frequency bin
%     f_hi   : High value of PSD frequency bin
%     psd    : Array of PSD values
%     err    : Array of PSD Errors
%     noise  : Array (*or single value*) of Poisson noise levels
%  
% OUTPUTS:
%     None
%
% LIST OF CALLED FUNCTIONS/SUBROUTINES NOT INTRINSIC TO S-LANG:
% 
% MODIFICATION HISTORY:
%
%     Version 0.1.01s: 2005/08/18 MAN - Sherpa version

   variable id, f_lo, f_hi, psd, err, mid=0, noise=NULL;

   switch(_NARGS)
   {
    case 5:
      (id, f_lo, f_hi, psd, err) = ();
   }
   {
    case 6:
      (id, f_lo, f_hi, psd, err, mid) = ();
   }
   {
    case 7:
      (id, f_lo, f_hi, psd, err, mid, noise) = ();
   }   
   {
      variable str=
     [" sitar_define_psd(id,flo,fhi,psd,err,[mid,noise]);"," ",
      "   Define psd to be a Sherpa counts data set.",
      " ",
      " Inputs:",
      "   id     : Dataset Index",
      "   f_lo   : Low value of PSD frequency bin",
      "   f_hi   : High value of PSD frequency bin",
      "   psd    : Array of PSD values",
      "   err    : Array of PSD Errors",
      "   mid    : 0 (default) = PSD is defined as a histogram with lo/hi freq. bins,",
      "            !=0, PSD is defined at frequency midpoint f_mid = (f_lo+f_hi)/2",
      "   noise  : Array (*or single value*) of Poisson noise levels (*subtracted*).",
      "            If undefined, background is undefined, and nothing is subtracted.",
      "  ",
      " Outputs:",
      "   None"];
      msg(str);
      return;
   }

   variable da = struct{lo,hi,mid};
   if(mid ==0)
   {
      da.lo = f_lo;
      da.hi = f_hi;
      da.mid = NULL;
   }
   else
   {
      da.lo = NULL;
      da.hi = NULL;
      da.mid = (f_lo+f_hi)/2.;
   }
   () = set_axes(id,da);

   if(noise !=NULL)
   {
      variable lnoise = length(noise); 
      if( lnoise == 1 )
      {
         noise = ones(length(psd))*noise;
         psd = psd-noise;
      }
      else if(lnoise != length(psd))
      {
         print("Noise array length differs from psd length");
         return;
      }
      else
      {
         psd = psd-noise;
      }
   }

   () = set_data(id,psd);
   () = set_errors(id,err);


   return;
}

#endif