% Last Updated: Dec. 12, 2008

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

% Public Functions in This File.  Usage message (almost always) when 
% function is called without arguments.

% keV_note      : Notice positive valued bins (cued from *first* data
%                 set, based upon energy cuts.  (Positive value 
%                 restriction ties into radio data functions later on.)  
%                 Also works on combined datasets.
% keV_note_sum  : Same as above, but cues off of summed data.
% chan_note     : Notice positive valued bins based upon channels
%
%  !!! The following will also group for combined datasets !!!
%
% grppha        : Rebin data from lowest keV bound to a min counts/bin
% grppha_cut    : Rebin data from specified keV bound to a min counts/bin
% grppha_min    : ... and add in a min number of channels per bin
% grppha_sn     : ... from specified bound, to a min S/N per bin
% grppha_sn_min : ... from specified bound, to a min S/N per bin, & 
%                 minimum number of channels per bin
% group         : A simplified version that combines all of the above
% group_bins    : Group the data using specified bin boundaries
%
%  !!!  These only work on one data set at a time !!!
%
% i2x_grp       : Turn an ISIS grouping into a GRPPHA/XSPEC one 
% write_x_grp   : ... take that grouping and apply it to a file
% apply_i2x_grp : Do the above two in one fell swoop
% x2i_grp       : Turn an XSPEC/GRPPHA grouping into an ISIS one
% apply_x2i_grp : ... and apply it in one fell swoop
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

public define keV_note()
{
   variable a,b,c,data,bin_list,bin_list_a,dindx;
   bin_list = Integer_Type[0];
   switch(_NARGS)
   {
    case 3:
      (a,b,c) = ();
      a = [a];
      if ( min(b) != 0. ) b = _A(b); else return 1;
      if ( min(c) != 0. ) c = _A(c); else return 1;
      b = [b];
      c = [c];
      if ( length(c) != length(b)) return 1;
      data = get_data_counts(a[0]);

      variable i=0;
      loop(length(b))
      {
         bin_list_a = where( data.value > 0. and data.bin_lo >= c[i]
                           and data.bin_hi < b[i] );
         bin_list = [bin_list,bin_list_a];
         i++;
      }

      foreach(a)
      {
         dindx = ();
         ignore(dindx);
         notice_list(dindx,bin_list);
      }
      return;
   }
   {
      variable str = [" keV_note([a],[minkev1,minkev2,...],[maxkev1,maxkev2,...]);"," ",
      "   Exclusively notice data between energy ranges [minkev1:maxkev1],",
      "   [minkev2:maxkev2],..., for data sets [a], where the data sets are",
      "   on a common grid. Note that for binned data minkev < E.bin_lo *and*",
      "   E.bin_hi <= maxkev to be noticed. Bins with zero *in the **first** data set*",
      "   cause that bin to be excluded in all data sets."];
      msg(str);
      return;
   }
} 

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

public define keV_note_sum()
{
   variable a,b,c,data,datasum=0,bin_list,bin_list_a,dindx;
   bin_list = Integer_Type[0];
   switch(_NARGS)
   {
    case 3:
      (a,b,c) = ();
      a = [a];
      if ( min(b) != 0. ) b = _A(b); else return 1;
      if ( min(c) != 0. ) c = _A(c); else return 1;
      b = [b];
      c = [c];
      if ( length(c) != length(b)) return 1;
      data = get_data_counts(a[0]);

      foreach(a)
      {
         dindx = ();
         datasum += get_data_counts(dindx).value;
      }

      variable i=0;
      loop(length(b))
      {
         bin_list_a = where( datasum > 0. and data.bin_lo >= c[i]
                           and data.bin_hi < b[i] );
         bin_list = [bin_list,bin_list_a];
         i++;
      }

      foreach(a)
      {
         dindx = ();
         ignore(dindx);
         notice_list(dindx,bin_list);
      }
      return;
   }
   {
      variable str = [" keV_note([a],[minkev1,minkev2,...],[maxkev1,maxkev2,...]);"," ",
      "   Exclusively notice data between energy ranges [minkev1:maxkev1],",
      "   [minkev2:maxkev2],..., for data sets [a], where the data sets are",
      "   on a common grid. Note that for binned data minkev < E.bin_lo *and*",
      "   E.bin_hi <= maxkev to be noticed.  Bins with 0 counts **in the sum**",
      "   are ignored"];
      msg(str);
      return;
   }
} 

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

public define chan_note()
{
   variable a,b,c,data,bin_list,bin_list_a;
   bin_list = Integer_Type[0];
   switch(_NARGS)
   {
    case 3:
      (a,b,c) = ();
      if ( min(b) < 0 ) return 1;
      if ( min(c) < 0 ) return 1;
      b = [b];
      c = [c];
      if ( length(c) != length(b)) return 1;
      data = get_data_counts(a);
      ignore(a);
      variable i;
      if(Isis_Reverse_Channels == 1)
      {
         i=length(b)-1;
         loop(length(b))
         {
            bin_list = [bin_list,[c[i]:b[i]:-1]];
            i--;
         }
         bin_list = length(data.value) -1 - bin_list;
         print("");
         print("Selecting Based on ENERGY Channels (first channel is 0)");
         print("");
      }
      else
      {
         i=0;
         loop(length(b))
         {
            bin_list = [bin_list,[b[i]:c[i]]];
            i++;
         }
         print("");
         print("Selecting Based on WAVELENGTH Channels (first channel is 0)");
         print("");
      }
      notice_list(a,bin_list);
      return;
   }
   {
      variable str = [" chan_note(a,[minchan1,minchan2,...],[maxchan1,maxchan2,...]);"," ",
      "   Exclusively notice data between channel ranges [minchan1:maxchan1],",
      "   [minchan2:maxchan2],..., for data set a.  Channel numbers are based upon",
      "   GROUPED values (i.e., as in XSPEC)"];
      msg(str);
      return;
   }
} 

alias("keV_note", "kev_note"); 
alias("keV_note", "keVnote");
alias("keV_note","kevnote");
alias("keV_note_sum", "kev_note_sum"); 
alias("keV_note_sum", "keVnote_sum");
alias("keV_note_sum","kevnote_sum");
alias("keV_note_sum","kevnotesum");
alias("chan_note","channote");
alias("chan_note","chan_not");
alias("chan_note","channot");

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

%  Begin a series of grouping functions, to mimic ftools a bit

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

public define grppha_min()
{
   variable a,minc,mincc,minkev;
   switch(_NARGS)
   {
    case 4:
      (a,minc,mincc,minkev)=();
       a = [a];
   }
   {
      variable str = [" grppha_min(a,mincounts,minchans,minkev);", " ",
       "   Rebin (possibly combined) data sets, from a specified lower keV bound,",
       "   to a minimum counts per bin and a minimum number of channels per bin.",
       "   First bins, below the minkev cutoff, are left unbinned."," ",
       "   a = indices of data sets.  Data sets with negative indices will",
       "       be grouped, but not included in the calculation of minimum counts",
       "   mincounts = minimum counts per bin",
       "   minkev = lower energy bound for grouping start"];
      msg(str);
      return;
   }

   variable la = length(a);
   variable s = Struct_Type[la];

   variable ii;
   for (ii=0; ii <= la-1; ii++)
   {
      rebin_data(abs(a[ii]),0);
      s[ii] = get_data_counts(abs(a[ii]));
   }

   variable len = length(s[0].value);
   variable idx = Integer_Type[len] + 1;
   variable sn = 1;
   variable ssum = 0.;
   variable nchan = 0;

   variable i;
   for (i=len-1; i >=0; i--)
   {
      idx[i] = sn*idx[i];

      for (ii=la-1; ii >= 0; ii--)
      {
         if(a[ii]>0)
         {
            ssum += s[ii].value[i];
         }
      }

      nchan++;

      if( (_A(s[0].bin_hi[i]) <= minkev) or (ssum >= minc and nchan >= mincc) )
      {
         sn = sn*(-1);
         ssum = 0.;
         nchan = 0;
      }
   }

   for (ii=la-1; ii >= 0; ii--)
   {
      rebin_data(abs(a[ii]),idx);
   }
   return;
}

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

public define grppha_cut()
{
   variable a,minc,minkev;
   switch(_NARGS)
   {
    case 3:
      (a,minc,minkev)=();
      grppha_min(a,minc,1,minkev);
      return;
   }
   {
      variable str = [" grppha_cut(a,mincounts,minkev);", " ",
       "   Rebin (possibly combined) data sets, from a specified lower keV bound,",
       "   to a minimum counts per bin. First bins, below the minkev cutoff,",
       "   are left unbinned."," ",
       "   a = indices of data sets.  Data sets with negative indices will",
       "       be grouped, but not included in the calculation of minimum counts",
       "   mincounts = minimum counts per bin",
       "   minkev = lower energy bound for grouping start"];
      msg(str);
      return;
   }
}

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

public define grppha()
{
   variable a,mincounts;
   switch(_NARGS)
   {
    case 2:
      (a,mincounts)=();
      grppha_min(a,mincounts,1,0.);
      return;
   }
   {
      variable str = [" grppha(a,mincounts);"," ",
       "   Rebin data set, from its lower keV bound, to a minimum counts", 
       "   per bin.","",
       "   a = indices of data sets.  Data sets with negative indices will",
       "       be grouped, but not included in the calculation of minimum counts",
       "   mincounts = minimum counts per bin"];
      msg(str);
      return;
   }
}

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



public define grppha_sn_min()
{
   variable a,minsn,minchan,minkev;
   switch(_NARGS)
   {
    case 4:
      (a,minsn,minchan,minkev)=();
      a = [a];
   }
   {
      variable str = [" grppha_sn_min([a,b,...],minsn,minchan,minkev);", " ",
       "   Rebin (possibly combined) data sets, from a specified lower keV bound,",
       "   to a minimum signal-to-noise in the combined dataset, and a minimum",
       "   number of channels per bin (number of channels for one dataset).",
       "   First bins, below the minkev cutoff, are left unbinned.",
       "   Reads the background data set, if it exists, and presumes the total",
       "   noise goes as: ",
       "      sqrt( (Total Counts) + (Back Counts)*(Back Scale*Exp. Ratio)^2 )"," ",
       "   [a,b,...] = indices of data sets. Negative indices are grouped,",
       "               but not used in calculation of S/N.  Datasets must",
       "               be on the same starting grid.",
       "   minsn = minimum signal-to-noise ratio per bin",
       "   minchan = minimum number of channels per bin",
       "   minkev = lower energy bound for grouping start"];
      msg(str);
      return;
   }
   
   variable la = length(a);
   variable s = Struct_Type[la], b = Array_Type[la], scl = Float_Type[la];
   variable bfile, bscl, bexp, dscl, dexp, nchan=0;

   variable ii;
   for (ii=0; ii <= la-1; ii++)
   {
      rebin_data(abs(a[ii]),0);
      s[ii] = get_data_counts(abs(a[ii]));
      b[ii] = get_back(abs(a[ii]));
      scl[ii] = 1.;

      if(b[ii] == NULL)
      {
         b[ii] = Float_Type[length(s[ii].value)];
      }
      else
      { 
         dscl = get_data_backscale(abs(a[ii]))[0];
         dexp = get_data_exposure(abs(a[ii]));

         bfile = get_data_info(abs(a[ii])).bgd_file;
         if(bfile != NULL and bfile != "#_define_bgd()")
         {
            bexp = fits_read_key(bfile,"EXPOSURE");
            bscl = fits_read_key(bfile,"BACKSCAL");

            if(bexp == NULL or bscl == NULL or bexp == 0 or bscl == 0)
            {
               scl[ii] = 1.;
            }
            else 
            {
               scl[ii]= (dexp*dscl)/(bexp*bscl);
            }
         }
         else
         {
            scl[ii] = 1;
         }

      }
   }

   variable len = length(s[0].value);
   variable idx = Integer_Type[len] + 1;
   variable sn = 1;
   variable ssum = 0.;
   variable bsum = 0.;
   variable bsum_scl = 0.;

   variable i;
   for (i=len-1; i >=0; i--)
   {
      idx[i] = sn*idx[i];

      for (ii=0; ii <= la-1; ii++)
      {
         if(a[ii]>0)
         {
            ssum += s[ii].value[i];
            bsum += b[ii][i];
            bsum_scl += b[ii][i]*scl[ii];
         }
      }

      nchan++;

      if (_A(s[0].bin_hi[i]) <= minkev) 
      {
         sn = sn*(-1);
         ssum = 0.;
         bsum = 0.;
         bsum_scl = 0.;
         nchan=0;
      }

      if(ssum > 0){
         if ( (ssum-bsum) / sqrt(ssum+bsum_scl) >= minsn and nchan >= minchan)
         {
            sn = sn*(-1);
            ssum = 0.;
            bsum = 0.;
            bsum_scl = 0.;
            nchan=0;
         }
      }
   }

   foreach(a)
   {
      ii = ();
      rebin_data(abs(ii),idx);
   }
   return;
}


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


public define grppha_sn_min_array()
{
   variable a,minsn,minchan,minkev,lchans,ichans=0;
   switch(_NARGS)
   {
    case 4:
      (a,minsn,minchan,minkev)=();
      a = [a];
      minchan = [minchan];
      minkev = [minkev];
      lchans = length(minchan);
      if(lchans != length(minkev))
      {
	 print("minchan and minkev arrays have different lengths");
         return;
      }
   }
   {
      variable str = [" grppha_sn_min([a,b,...],minsn,[minchan],[minkev]);", " ",
       "   Rebin (possibly combined) data sets, from a specified *array* of lower",
       "   keV bounds, to a minimum signal-to-noise in the combined dataset, and an",
       "   *array* of minimum number of channels per bin (number of channels for one",
       "   dataset).  The [minchan] array corresponds to the [minkev] array, and the",
       "   latter is assumed to be in energ ascending order.  First bins, below the",
       "   minkev cutoff, are left unbinned.  This reads the background data set, if", 
       "   it exists, and presumes the total noise goes as: ",
       "      sqrt( (Total Counts) + (Back Counts)*(Back Scale*Exp. Ratio)^2 )"," ",
       "   [a,b,...] = indices of data sets. Negative indices are grouped,",
       "               but not used in calculation of S/N.  Datasets must",
       "               be on the same starting grid.",
       "   minsn     = minimum signal-to-noise ratio per bin",
       "   [minchan] = array of minimum number of channels per bin",
       "   [minkev]  = array of lower energy bounds for grouping start/changes"];
      msg(str);
      return;
   }
   
   variable la = length(a);
   variable s = Struct_Type[la], b = Array_Type[la], scl = Float_Type[la];
   variable bfile, bscl, bexp, dscl, dexp, nchan=0;

   variable ii;
   for (ii=0; ii <= la-1; ii++)
   {
      rebin_data(abs(a[ii]),0);
      s[ii] = get_data_counts(abs(a[ii]));
      b[ii] = get_back(abs(a[ii]));
      scl[ii] = 1.;

      if(b[ii] == NULL)
      {
         b[ii] = Float_Type[length(s[ii].value)];
      }
      else
      { 
         dscl = get_data_backscale(abs(a[ii]))[0];
         dexp = get_data_exposure(abs(a[ii]));

         bfile = get_data_info(abs(a[ii])).bgd_file;
         if(bfile != NULL and bfile != "#_define_bgd()")
         {
            bexp = fits_read_key(bfile,"EXPOSURE");
            bscl = fits_read_key(bfile,"BACKSCAL");

            if(bexp == NULL or bscl == NULL or bexp == 0 or bscl == 0)
            {
               scl[ii] = 1.;
            }
            else 
            {
               scl[ii]= (dexp*dscl)/(bexp*bscl);
            }
         }
         else
         {
            scl[ii] = 1;
         }

      }
   }

   variable len = length(s[0].value);
   variable idx = Integer_Type[len] + 1;
   variable sn = 1;
   variable ssum = 0.;
   variable bsum = 0.;
   variable bsum_scl = 0.;

   variable i;
   for (i=len-1; i >=0; i--)
   {
      idx[i] = sn*idx[i];

      for (ii=0; ii <= la-1; ii++)
      {
         if(a[ii]>0)
         {
            ssum += s[ii].value[i];
            bsum += b[ii][i];
            bsum_scl += b[ii][i]*scl[ii];
         }
      }

      nchan++;

      if (_A(s[0].bin_hi[i]) <= minkev[0]) 
      {
         sn = sn*(-1);
         ssum = 0.;
         bsum = 0.;
         bsum_scl = 0.;
         nchan=0;
      }

      if(ssum > 0){

         ichans = max( where( minkev <  _A(s[0].bin_hi[i]) ) );

         if ( (ssum-bsum) / sqrt(ssum+bsum_scl) >= minsn and nchan >= minchan[ichans])
         {
            sn = sn*(-1);
            ssum = 0.;
            bsum = 0.;
            bsum_scl = 0.;
            nchan=0;
         }
      }
   }

   foreach(a)
   {
      ii = ();
      rebin_data(abs(ii),idx);
   }
   return;
}

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

public define grppha_sn()
{
   variable a,minsn,minkev;
   switch(_NARGS)
   {
    case 3:
      (a,minsn,minkev)=();
      grppha_sn_min(a,minsn,1,minkev);
      return;
   }
   {
      variable str = [" grppha_sn([a,b,...],minsn,minkev);", " ",
       "   Rebin (possibly combined) data sets, from a specified lower keV bound,",
       "   to a minimum signal-to-noise in the combined dataset. First bins, below",
       "   the minkev cutoff, are left unbinned.  Reads the background data set,",
       "   if it exists, and presumes the total noise goes as: ",
       "      sqrt( (Total Counts) + (Back Counts)*(Back Scale*Exp. Ratio)^2 )"," ",
       "   [a,b,...] = indices of data sets. Negative indices are grouped,",
       "               but are not used in calculation of S/N.  Datasets must",
       "               be on the same starting grid.",
       "   minsn = minimum signal-to-noise ratio per bin",
       "   minkev = lower energy bound for grouping start"];
      msg(str);
      return;
   }
}

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

public define group()
{
   variable ind=Integer_Type[0],kev,chan,sn;

   switch(_NARGS)
   {
    (_NARGS>=1):
      variable i, list = __pop_list(_NARGS);
      _for i (0,length(list)-1,1)
      {
         ind = [ind,list[i]];
      }
   }
   {
      variable str = [" group([a,b,...]; qualifiers);", " ",
       "   Rebin (possibly combined) data sets, from specified lower keV bounds (optional),",
       "   to a minimum signal-to-noise (default=4.5) *and* a minimum number of channels",
       "   (default=1) in each grouped bin. Any channels below the minum keV cutoff are",
       "   left unbinned.  The background data set, if it exists, is included in the S/N",
       "   calculation. Group presumes that the total noise goes as: ","",
       "      sqrt( (Total Counts) + (Back Counts)*(Back Scale*Exp. Ratio)^2 )"," ",
       " Inputs:","",
       "   [a,b,...] = indices of data sets. Negative indices are grouped but are not used",
       "               in calculation of S/N.  Datasets must be on the same starting grid.",
       " ",
       " Qualifiers:","",
       "   sn   = minimum signal-to-noise ratio per bin (default=4.5)",
       "   kev  = single value or array of lower bounds for the start of each minimum",
       "          number of channels specification (default=0)",
       "   chan = minimum number of channels (default=1) required per bin, starting at the",
       "          corresponding specified keV bound, and proceeding to higher energies"];
      msg(str);
      return;
   }
      sn = qualifier("sn",4.5);
      kev = qualifier("kev",0);
      kev = [kev];
      chan = qualifier("chan",1);
      chan = [chan];
      if(length(chan) != length(kev))
      {
         () = printf("\n keV and minimum channel array have different lengths. \n\n");
         return;
      }
      grppha_sn_min_array(ind,sn,chan,kev);
      return;
}

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

public define group_bins()
{
   variable a,blo,bhi,lo,hi,kev;
   switch(_NARGS)
   {
    case 3:
      (a,blo,bhi)=();
      a = [a];
      blo = [blo];
      bhi = [bhi];
      if(length(blo) != length(bhi))
      {
	 print("Channel boundary arrays do not have the same length");
         return;
      }
   }
   {
      variable str = [" group_bins([a,b,...],bin_lo,bin_hi [; kev=1);", " ",
       "   Rebin (possibly combined) data sets to match as closely as possible",
       "   wavelength boundaries specified by input arrays.  Whether the input",
       "   arrays are wavelength or energy, binning is *data* lower *wavelength*",
       "   boundary >= bin input lower *wavelength* boundary and < bin input",
       "   upper *wavelength* boundary. I.e., rebinned upper wavelength boundaries",
       "   can straddle the upper wavelength boundaries of the input arrays.","",
       "   Inputs:",
       "    [a,b,...] = indices of data sets. Datasets must be on the same grid.",
       "    bin_lo    = lower wavelength boundaries of input grid",
       "    bin_hi    = upper wavelength boundaries of input grid","",
       "   Optional Input (Qualifier):",
       "    kev      != 0, to specify that input grid is in keV units"];
      msg(str);
      return;
   }
   
   kev = qualifier("kev",0);
   
   if(kev!=0)
   {
      lo = _A(bhi);
      hi = _A(blo);
   }
   else
   {
      lo = blo;
      hi = bhi;
   }

   variable la = length(a);

   variable ii;
   for (ii=0; ii <= la-1; ii++)
   {
      rebin_data(abs(a[ii]),0);
   }

   variable dlo = get_data_counts(abs(a[0])).bin_lo;
   variable idx = Integer_Type[length(dlo)] + 1;
   variable sn = 1;

   variable is=0;
   variable start_lo = lo[0];
   variable start_hi = hi[0];
   variable lbin = length(lo);

   if(max(dlo)< lo[0] or min(dlo)>= hi[-1])
   {
      () = printf(" Data and input arrays do not overlap \n");
      return;
   }

   variable i;
   for (i=0; i<=length(dlo)-1; i++)
   {
      idx[i] = sn*idx[i];

      if(dlo[i] < start_lo)
      {
         sn = sn*(-1);
      }
      else if(dlo[i] >= start_hi)
      {
         sn = sn*(-1);
         is++;

         if(is<lbin)
         {
            start_lo = lo[is];
            start_hi = hi[is];
         }
      }
   }

   foreach(a)
   {
      ii = ();
      rebin_data(abs(ii),idx);
   }
   return;
}



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

public define i2x_grp()
{
   switch(_NARGS)
   {
    case 1:
      variable i, lb, ub, id, gdi, lgdi, iw, liw;

      id = ();
      gdi = [reverse(get_data_info(id).rebin),5];
      lgdi = length(gdi)-1;

      variable grp = Integer_Type[lgdi];

      iw = [where( shift(gdi,-1) != gdi )];
      liw = length(iw);
      iw[liw-1] = lgdi;

      grp[iw[[0:liw-2]]] = 1;

      foreach i ([0:liw-2])
      {
         lb = iw[i];
         ub = iw[i+1]-1;
         if(ub>lb)
         {
            grp[[lb+1:ub]] = -1;
         }
      }
      
      return grp;
   }
   {
      variable str = [" grouping = i2x_grp(id);"," ",
      "   For any dataset, id, grouped in ISIS, create an *energy ordered*",
      "   grouping vector that follows the grppha conventions.  (ISIS groups",
      "   flagged with 0, however, become *noticed* groups of 1.)"];
      msg(str);
      return;
   }
}

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

public define write_x_grp()
{
   switch(_NARGS)
   {
    case 2:
      variable id, grouping, file, fp, colnum;

      (id, grouping) = ();
      file = get_data_info(id).file;

      fp = fits_open_file(file+"[SPECTRUM]","w");
      () = _fits_get_colnum(fp, "grouping", &colnum);
      () = _fits_write_col(fp,colnum,1,1,grouping);
      fits_close_file(fp);

      return;
   }
   {
      variable str = [" write_x_grp(id,grouping);"," ",
      "   Write an XSPEC style grouping to the *file* associated with id.",
      "   *Grouping vector must follow grppha conventions and be energy ordered.*"," ",
      "   See: grouping = i2x_grp(id);"];
      msg(str);
      return;
   }
}

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

public define apply_i2x_grp()
{
   switch(_NARGS)
   {
    case 1:
      variable id, grp;

      id = ();
      grp = i2x_grp(id);
      write_x_grp(id,grp);

      return;
   }
   {
      variable str = [" apply_i2x_grp(id);"," ",
      "   Take the ISIS defined grouping from dataset id, create an XSPEC style",
      "   grouping based upon that, and then write that grouping to the",
      "   associated file"];
      msg(str);
      return;
   }
}

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

public define x2i_grp()
{
   switch(_NARGS)
   {
    case 1:
      variable id, file, grp, lgrp, iw, liw, sn, i;

      id = ();
      file = get_data_info(id).file;
      grp = fits_read_col(file,"grouping");

      lgrp = length(grp);
      iw = where(grp == 1);
      liw = length(iw);
      iw = [iw,lgrp];

      sn = 1;
      foreach i ([0:liw-1])
      {
         grp[[iw[i]:max([iw[i],iw[i+1]-1])]] = sn;
         sn = -1*sn;
      }

      return reverse(grp);
   }
   {
      variable str = [" grouping = x2i_grp(id);"," ",
      "   Read an XSPEC style *energy ordered* grouping, associated with the",
      "   dataset indicated by id, and convert it to a *wavelength ordered* ISIS",
      "   style grouping suitable for use in rebinning, i.e., ",
      "      isis> rebin_data(indx,grouping);"];
      msg(str);
      return;
   }
}

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

public define apply_x2i_grp()
{
   switch(_NARGS)
   {
    case 1:
      variable id, grp; 

      id = ();
      grp = x2i_grp(id);
      rebin_data(id,grp);

      return;
   }
   {
      variable str = [" apply_x2i_grp(id);"," ",
      "   Take the XSPEC style grouping from the *file* associated with",
      "   dataset id, convert it to an ISIS style, wavelength ordered",
      "   grouping, and apply it to dataset id."];
      msg(str);
      return;
   }
}