#!/usr/bin/env isis-script

private variable Script_Version_String = "1.0.0";

%    pile_estimate.sl      2010 April 23
%
%    This file is the pile_estimate routine and associated software.
%    Copyright (C) 2008 Massachusetts Institute of Technology
%                       Author: Michael A. Nowak, mnowak@alum.mit.edu
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.

% Program requires: ISIS, & the S-Lang cmdopt, XPA, & GTK modules.
% ("vwhere" is part of the GTK module.)
% cmdopt module is normally bundled with S-lang and/or ISIS.
% See the ISIS web page: http://space.mit.edu/cxc/isis/index.html
% And the S-lang modules page: http://space.mit.edu/cxc/software/slang/modules

% vwhere for lightcurve filtering
define vwhere();
variable nofilter = 0;
variable show_select =0;
try
{
  require("vwhere");
}
catch AnyError:
{
  nofilter = 1;
  ()=fprintf(stderr, "** Warning: Unable to import the gtk module (vwhere), assuming --nofilter\n");
}

% DS9 for displaying - requires the XPA module
autoload("ds9_new","ds9object");

% Taken from John Davis's aeattcor.sl
private define open_ds9 ()
{
  variable e, ds9 = NULL;
  try (e)
  {
    ds9 = ds9_new(sprintf("ds9.%s.%d", path_basename(__argv[0]), getpid()); verbose=1, timeout=90);
  }
  catch AnyError:
  {
    ()=fprintf (stderr, "** Warning: Unable to open ds9: %S\n", e.message);
  }
  return ds9;
}

% cmdopt module for parsing the command line
require("cmdopt");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% A simple boxcar 3x3 convolution written by John Davis

private define conv3x3 (img)
{
  % An FFT may be faster but this is simple
  variable dims = array_shape (img);
  variable n0 = dims[0];
  variable n1 = dims[1];
  variable i = [0:n0-1];
  variable j = [0:n1-1];
  variable i0 = (i-1);
  variable j0 = (j-1);
  variable i1 = (i+1) mod n0;
  variable j1 = (j+1) mod n1;

  variable img1 = @img;
  img1[i,j] += img[i0,j0]; img1[i,j] += img[i0,j]; img1[i,j] += img[i0,j1];
  img1[i,j] += img[i,j0];                          img1[i,j] += img[i,j1];
  img1[i,j] += img[i1,j0]; img1[i,j] += img[i1,j]; img1[i,j] += img[i1,j1];

  return __tmp(img1) / 9.0f;
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% A simple program written by John Davis to watch for a key input

private define getkey()
{
  ()=system("stty raw");
  variable ch;
  if(1 != read(fileno(stdin), &ch, 1))  ch=0;
  ()=system("stty sane");
  return ch;
}

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

private define exit_usage()
{
  variable fp = stderr;
  ()=fprintf (fp, "\nUsage: %s [options] evtfile\n",__argv[0]);
  ()=fprintf (fp, "\n%s\n", %{{{
+"   Create a pileup map out of the Suzaku event data in `evtfile'. A rate\n"
+"   lightcurve will be created in 3 energy bins (ra, rb, rc). Unless the\n"
+"   --nofilter option has been chosen, these data will be passed to `vwhere',\n"
+"   which can be used to select times of uniform rate and colors. Time\n"
+"   selections are printed to screen and written to a file, e.g., for xselect.\n"
+"\n"
+"   A box-car smoothed (3x3 pixel bins) image is displayed, with the image\n"
+"   rescaled to be the *estimated* pileup fraction.  Pileup fraction\n"
+"   here is defined to be the ratio of events lost via grade *or* energy\n"
+"   migration to the events expected in the absence of  pileup.  The\n"
+"   image displays discrete steps that represent the *minimum* pileup\n"
+"   fraction in the displayed region, with the exception of the highest\n"
+"   displayed value. This value corresponds to the image maximum, and is\n"
+"   shown over regions at this maximum down to half way towards the next\n"
+"   lowest displayed value. (A continuous pileup image can be displayed by\n"
+"   using the --cont-pfrac option.).\n"
+"\n"
+"   One can quit without proceeding further, or create a region file. Any\n"
+"   combination circular, elliptical, box, and annular exclude/include regions\n"
+"   are supported, so long as there is at least one include region. Include\n"
+"   regions are applied first, and then the exclude regions are applied.\n"
+"   Regions must be in physical coordinates and use the ds9/funtools format.\n"
+"   Choosing 'y' will automatically save the region to a file, and the\n"
+"   --show_select option redraws the image with the region selection applied.\n"
+"\n"
+"   Using the time and region filtered events, an estimate is made of the\n"
+"   values to which one should set pileup parameters (frame time, number\n"
+"   of regions, PSF fraction) if fitting the filtered events with the ISIS\n"
+"   pileup model. Values are printed to a file. (This is experimental, and\n"
+"   may not be 100\% quantitatively accurate.)\n"
+"\n"
+" Options (must precede input evtfile; arrays/ranges do not use () or []):\n"
+"\n"
+"   --nofilter\n"
+"         The selection of events with `vwhere' is skipped.\n"
+"         If the gtk module is not installed, --nofilter is the default.\n"
+"   --tbin=val\n"
+"         The lightcurve is created with time bins of the total frame time\n"
+"         (exposure, plus delay time, e.g., for burst mode observations),\n"
+"         multiplied by an integer binning factor, tbin.  Rates, however,\n"
+"         are per exposure time.  [Default=16]\n"
+"   --save-lightcurve=val\n"
+"          Save lightcurve into FITS file (ignored, if --nofilter is chosen).\n"
+"   --pfrac=vals\n"
+"         The minimum fractional pileup levels for the discretized pileup\n"
+"         image.  [Default=0.01,0.02,0.05,0.1]\n"
+"   --cont-pfrac\n"
+"         Show a continuous image of the fractional pileup.\n"
+"         (The --pfrac option is ignored if --cont-pfrac is used.)\n"
+"   --only-save-pileup-img=val\n"
+"         The estimated pileup fraction is only saved as a FITS image,\n"
+"         but not displayed through ds9.\n"
+"   --show-select\n"
+"         Display the pileup image with selections applied at the end.\n"
+"   --outfile=val\n"
+"         The root name for the output files, i.e., val.flt (time filter),\n"
+"         val.reg (region file), val.parms (suggested pileup parameters).\n"
+"         [Default=input file name]\n"
+"   --filter-file=val\n"
+"         (Optional input.) The name of a file containing pairs of times to\n"
+"         be considered good times.  Creates a new column, filter, that is \n"
+"         passed to vwhere.  This column will contain 1's for times within \n"
+"         these GTI values, and 0's for times outside of them.  This allows\n"
+"         one to apply previously defined filters to the data. [default=NULL]\n"
+"   --a-band=vals  [default: 0.5-1.5]\n"
+"   --b-band=vals  [default: 1.5-3]\n"
+"   --c-band=vals  [default: 3-9]\n"
+"         The energy bands (keV) for the three colors in the lightcurve,\n"
+"         specified as 'lo-hi', or, equivalently, 'lo,hi'\n"
+"   --alpha=val\n"
+"         The usual pileup alpha parameter (see the Chandra ABC Guide to\n"
+"         Pileup, or Davis 2001, ApJ, 562, 575).  Affects the estimated\n"
+"         pileup fractions and suggested model parameters [Default=0.5]\n"
+"   -v, --version\n"
+"         Print version.\n"
+"   -h, --help\n"
+"         This message.\n"); %}}}
  exit(1);
}

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

private define exit_version()
{
  ()=fprintf(stdout, "Version: %S\n", Script_Version_String);
  exit (0);
}

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

% Define the types of regions we will parse from DS9.  Find the points
% on the interior of the regions.

define rotate(xpix,ypix,xcent,ycent,angl)
{
   angl *= PI/180.;
   return cos(angl)*(xpix-xcent) + sin(angl)*(ypix-ycent),
         -sin(angl)*(xpix-xcent) + cos(angl)*(ypix-ycent);
}

define circle(xpix,ypix,xcent,ycent,rad)
{
  return where(hypot(xpix-xcent, ypix-ycent) <= rad);
}

define ellipse(xpix,ypix,xcent,ycent,radx,rady,angl)
{
  variable xtran, ytran;
  (xtran, ytran) = rotate(xpix,ypix,xcent,ycent,angl);
  variable rpix = (xtran/radx)^2+(ytran/rady)^2;
  return where(rpix <= 1);
}

define box(xpix,ypix,xcent,ycent,xlength,ylength,angl)
{
  variable xtran, ytran;
  (xtran, ytran) = rotate(xpix,ypix,xcent,ycent,angl);
  return where( abs(xtran) <= xlength/2. and abs(ytran) <= ylength/2. );
}

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

define slsh_main ()
{
  % The binning of the lightcurve, in units of the multiple of the
  % basic time step between frames.
  variable tbin = 16;

  variable save_lightcurve = NULL;

  % The pileup levels at which we'll create the image
  variable pfrac = [0.01,0.02,0.05,0.1];

  variable continuous_pfrac = 0;
  variable save_pileup_img = NULL;

  variable filename, outfile=NULL, filter_file=NULL;

  % energy bands (in keV)
  variable a_band,  ea_lo=0.5, ea_hi=1.5,
           b_band,  eb_lo=1.5, eb_hi=3.,
           c_band,  ec_lo=3.,  ec_hi=9.;

  % Usual pileup alpha parameter.
  variable alpha=qualifier("alpha",0.5);

  % Read in the user choices
  variable c = cmdopt_new ();
  c.add("nofilter", &nofilter; inc);
  c.add("tbin", &tbin; type="float");
  c.add("save-lightcurve", &save_lightcurve; type="str");
  c.add("pfrac", &pfrac; type="str");
  c.add("cont-pfrac", &continuous_pfrac; inc);
  c.add("only-save-pileup-img", &save_pileup_img; type="str");
  c.add("show-select", &show_select; inc);
  c.add("outfile", &outfile; type="str");
  c.add("filter-file", &filter_file; type="str");
  c.add("a-band", &a_band; type="str");
  c.add("b-band", &b_band; type="str");
  c.add("c-band", &c_band; type="str");
  c.add("alpha", &alpha; type="float");
  c.add("h|help", &exit_usage);
  c.add("v|version", &exit_version);
  variable i = c.process (__argv,1);
  if(i+1 != __argc)  exit_usage();
  filename = __argv[i];
  if(outfile==NULL)  outfile = filename;

  if(typeof(pfrac)==String_Type)
    pfrac = array_map(Double_Type, &atof, strtok(pfrac,"-([,;])"));

  if(typeof(a_band)==String_Type)
    a_band = array_map(Double_Type, &atof, strtok(a_band,"-([,;])"));

  if(typeof(b_band)==String_Type)
    b_band = array_map(Double_Type, &atof, strtok(b_band,"-([,;])"));

  if(typeof(c_band)==String_Type)
    c_band = array_map(Double_Type, &atof, strtok(c_band,"-([,;])"));

  % Some precautionary measures to make sure inputs are sane

  tbin = int(tbin);
  if( tbin < 1 ) tbin= 1;

  if(length(pfrac)<1)  pfrac=[0.01];
  pfrac = pfrac[array_sort(pfrac)];

  ifnot(length(a_band)==2 and 0<__is_numeric(a_band)<3)  a_band = [ea_lo, ea_hi];
  ifnot(length(b_band)==2 and 0<__is_numeric(b_band)<3)  b_band = [eb_lo, eb_hi];
  ifnot(length(c_band)==2 and 0<__is_numeric(c_band)<3)  c_band = [ec_lo, ec_hi];

  ea_lo = min(a_band);  ea_hi = max(a_band);
  eb_lo = min(b_band);  eb_hi = max(b_band);
  ec_lo = min(c_band);  ec_hi = max(c_band);

  % The base pileup level below which we set to the mean of non-zero pixels
  variable base = 3.e-3;
  if(base > pfrac[0])  base=pfrac[0]/2.;

  % event information (time, pha/pi value, x & y coordinates)
  variable evts = fits_read_table(filename, ["time","pi","x","y"]);

  % arrays for GTI times
  variable gti_lo, gti_hi;
  (gti_lo,gti_hi) = fits_read_col(filename+"[GTI]","start","stop");

  % Total exposure time, frame time, integration time/frame, dead
  % time/frame, fractional exposure, start time, stop time
  variable expos, snapti, delay, tstart, tstop;
  (expos, snapti, delay, tstart, tstop) =
    fits_read_key(filename,"EXPOSURE","SNAPTI1","DELAY1","TSTART","TSTOP");
  variable frame = snapti + delay;
  variable frac = snapti/frame;

  if(gti_lo[0] < tstart)  gti_lo[0] = tstart;
  if(gti_hi[-1] > tstop)  gti_hi[-1] = tstop;

  % We will first pre-filter the data & only keep the good times.
  variable gti_hist_bounds = Double_Type[2*length(gti_lo)];
  gti_hist_bounds[[0::2]] = gti_lo;
  gti_hist_bounds[[1::2]] = gti_hi;
  variable rtot;
  variable gti_hist = histogram(evts.time, gti_hist_bounds, make_hi_grid(gti_hist_bounds), &rtot);
  variable gti_indx_array, gti_indx = Integer_Type[0];
  foreach gti_indx_array (rtot[[0::2]])
    gti_indx = [gti_indx,gti_indx_array];
  struct_filter(evts, gti_indx);
  evts = struct_combine(evts, struct { kev = 3.65e-3*evts.pi });

  % First create a lightcurve of rate in each frame.  (We ignore any
  % slight misalignment of the histogram bins with the start/stop times
  % of each frame, and just work with the event times as being in the
  % middle of our bin.  Final GTI filtering will take care of that, and
  % lose at most one frame on any GTI boundary.)
  variable lo = [tstart:tstop:frame];
  variable hi = lo + frame;
  variable llo = length(lo);

  variable htot = histogram(evts.time, lo, hi, &rtot);
  % Exposure time in each frame - we will set to zero where there are
  % no counts.
  variable texp=Double_Type[llo];
  texp[where(htot>0)] = frame*frac;

  % Need to add logic in case there are no good times

  variable rbin=Array_Type[llo/tbin+1];
  $0 = Double_Type[llo/tbin+1];
  variable texp_tot=@$0, bin_lo=@$0, bin_hi=@$0;
  ifnot(nofilter)
  {
    variable lightcurve = struct { time=@$0, rate=@$0, error=@$0, ra=@$0, rb=@$0, rc=@$0 };
    variable ha = histogram(evts.time[where( ea_lo <= evts.kev < ea_hi )], lo, hi);
    variable hb = histogram(evts.time[where( eb_lo <= evts.kev < eb_hi )], lo, hi);
    variable hc = histogram(evts.time[where( ec_lo <= evts.kev < ec_hi )], lo, hi);
  }

  % Bin up the lightcurve by rates. Only include good time frames.
  _for i (0,llo-tbin,tbin)
  {
    variable j=i/tbin;
    variable indx = [i:min([i+tbin,llo])-1];

    variable indx_indx, rbin_indx=Integer_Type[0];
    foreach indx_indx (indx)
      rbin_indx = [rbin_indx,rtot[indx_indx]];
    rbin[j] = rbin_indx;

    texp_tot[j] = sum(texp[indx]);
    bin_lo[j]=lo[indx[0]];
    bin_hi[j]=hi[indx[-1]];

    if(texp_tot[j]>0 and not nofilter)
    {
      lightcurve.rate[j]  =      sum(htot[indx]) /texp_tot[j];
      lightcurve.error[j] = sqrt(sum(htot[indx]))/texp_tot[j];
      lightcurve.ra[j]    =      sum(ha  [indx]) /texp_tot[j];
      lightcurve.rb[j]    =      sum(hb  [indx]) /texp_tot[j];
      lightcurve.rc[j]    =      sum(hc  [indx]) /texp_tot[j];
    }
  }
  variable igd = where(texp_tot>0);
  variable ivw = Integer_Type[0], iw;
  ifnot(nofilter)
  {
    lightcurve.time = (bin_lo + bin_hi)/2.;
    struct_filter(lightcurve, igd);

    if(filter_file != NULL)
    {
       lightcurve = struct_combine(lightcurve, 
                           struct{filter=Char_Type[length(lightcurve.time)]});

       (gti_lo,gti_hi) = readcol(filter_file,1,2);
       _for i (0,length(gti_lo)-1,1)
       {
          iw = where( lightcurve.time >= gti_lo[i] and 
                      lightcurve.time <= gti_hi[i]     );
          if(length(iw)>0) lightcurve.filter[iw] = 1;
       }
    }
      
    lightcurve.time -= tstart;
    ivw = vwhere(lightcurve);
    lightcurve.time += tstart;

    if(filter_file == NULL)
    {
       lightcurve = struct_combine(lightcurve, 
                           struct{filter=Char_Type[length(lightcurve.time)]});
    }
    if(ivw != NULL)
    {
       lightcurve.filter[*] = 0;
       lightcurve.filter[ivw] = 1;
    }

    if(save_lightcurve != NULL)
      fits_write_binary_table(save_lightcurve, "light curve", lightcurve,
                              struct {
                                TUNIT1 = "s",
                                TUNIT2 = "c/s",
                                TUNIT3 = "c/s",
                                TUNIT4 = "c/s",
                                TUNIT5 = "c/s",
                                TUNIT6 = "c/s",
                                TUNIT7 = "boolean",
                                A_BAND = sprintf("%.3f-%.3f keV", ea_lo, ea_hi),
                                B_BAND = sprintf("%.3f-%.3f keV", eb_lo, eb_hi),
                                C_BAND = sprintf("%.3f-%.3f keV", ec_lo, ec_hi),
                                TELESCOP = fits_read_key(filename, "TELESCOP"),
                                INSTRUME = fits_read_key(filename, "INSTRUME"),
                                OBJECT   = fits_read_key(filename, "OBJECT"),
                                MJDREFI  = fits_read_key(filename, "MJDREFI"),
                                MJDREFF  = fits_read_key(filename, "MJDREFF"),
                              });
  }

  variable itm;
  if(length(ivw) > 0)
  {
    if(length(ivw)==length(igd))
      ()=fprintf(stdout,"\nAll data selected; resulting filter is consistent with internal GTI.\n");
    itm = igd[ivw];
  }
  else
  {
    ()=fprintf(stdout,"\nNull vwhere filter; resulting filter is consistent with internal GTI.\n");
    itm = igd;
  }

  % Take the selected points from vwhere, and figure out the
  % contiguous time intervals that they represent
  variable litm = length(itm);

  % Look for where there are gaps in the selected lightcurve.
  variable igap = where( (itm - shift(itm,-1)) > 1);
  if( length(igap) == 0 )
    igap = [0,litm];
  else
    igap = [0,igap, litm];

  % Print out results in format suitable for xselect filtering
  ()=fprintf(stdout,"\nRequired Time Filter:\n\n");
  variable fp = fopen(outfile+".flt","w");
  _for i (1,length(igap)-1,1)
  {
    ()=fprintf(stdout,"%15.10e        %15.10e\n", bin_lo[ itm[igap[i-1]] ], bin_hi[ itm[igap[i]-1] ]);
    ()=fprintf(fp    ,"%15.10e        %15.10e\n", bin_lo[ itm[igap[i-1]] ], bin_hi[ itm[igap[i]-1] ]);
  }
  ()=fclose(fp);
  ()=fprintf(stdout,"\n");

  % Generate an array of indices of the events we have kept Need some
  % logic in case the choice results in a null selection

  variable lrev = array_map(Integer_Type,&length,rbin);
  itm = itm[where(lrev[itm] > 0)];
  lrev = lrev[itm];

  variable irev = int(cumsum([0,lrev]));
  variable ikeep = Integer_Type[int(sum(lrev))];

  _for i (0,length(itm)-1,1)
    ikeep[[irev[i]:irev[i]+lrev[i]-1]] = rbin[itm[i]];

  % The coordinates of what we just kept
  variable x = evts.x[ikeep];
  variable y = evts.y[ikeep];

  % Regrid to 1x1 pixel.
  variable dx = min(x)+0.5, dy=min(y)+0.5;
  variable xhist = [1:int(max(x)+dx):1];
  variable yhist = [1:int(max(y)+dy):1];
  variable piximg = histogram2d(y,x,yhist,xhist);
  variable lx = length(xhist), ly = length(yhist);
  variable xpix = Float_Type[ly,lx];
  variable ypix = Float_Type[ly,lx];

  _for i (0,ly-1,1)  xpix[i,*] = xhist;
  _for i (0,lx-1,1)  ypix[*,i] = yhist;

  % Smooth on 3x3, multiply by 9, divide by binning factors, to get
  % counts/3x3 pixels/frame
  piximg = (9.*conv3x3(piximg)/sum(texp_tot[itm]))*(frame*frac);

  % Convert to counts we *should* have seen without pileup.  Third
  % order, alpha-dependent approximation
  piximg=(1+((2-alpha)/2.+(2*alpha^2-9*alpha+9)/6.*piximg)*piximg)*piximg;

  % Convert to pileup fraction - this is the "harshest" definition of
  % pileup fraction, f_t, from the Chandra ABC Guide.  It is the
  % fraction of incident events affected by *either* grade *or*
  % energy migration.  Other definitions of pileup fraction (of which
  % there are many) will give smaller numbers.
  piximg = 1 - exp(-piximg);

  variable pileimg=@piximg;

  ifnot(continuous_pfrac)
  {
    % Transform to discontinuous pileup fractions

    variable pmx = max(pileimg), npmx;

    i = where(pileimg < base and pileimg > 0);
    if(length(i)>0) { npmx = mean(pileimg[i]); pileimg[i] = npmx; }

    i = where(pileimg >= base and pileimg < pfrac[0]);
    if(length(i)>0) { pileimg[i] = base; npmx = base; }

    variable lpf = length(pfrac);
    _for j (1,lpf-1,1)
    {
      i = where(pileimg >= pfrac[j-1] and pileimg < pfrac[j]);
      if(length(i)>0) { pileimg[i] = pfrac[j-1]; npmx = pfrac[j-1]; }
    }

    i = where(pileimg >= pfrac[lpf-1]);
    if(length(i)>0) { pileimg[i] = pfrac[lpf-1]; npmx = pfrac[lpf-1]; }

    i = where(piximg >= npmx + (pmx-npmx)/2.);
    if(length(i)>0) pileimg[i] = pmx;
  }

  if(save_pileup_img != NULL)
  {
    fits_write_image_hdu(save_pileup_img, "estimated pileup fraction", pileimg,
                         struct {
                           CRPIX1 = fits_read_key(filename, "TCRPX10"),
                           CDELT1 = fits_read_key(filename, "TCDLT10"),
                           CTYPE1 = fits_read_key(filename, "TCTYP10"),
                           CRVAL1 = fits_read_key(filename, "TCRVL10"),
                           CUNIT1 = fits_read_key(filename, "TUNIT10"),
                           CRPIX2 = fits_read_key(filename, "TCRPX11"),
                           CDELT2 = fits_read_key(filename, "TCDLT11"),
                           CTYPE2 = fits_read_key(filename, "TCTYP11"),
                           CRVAL2 = fits_read_key(filename, "TCRVL11"),
                           CUNIT2 = fits_read_key(filename, "TUNIT11"),
                         });
    exit(0);
  }

  % Open DS9
  variable ds9 = open_ds9();
  if(ds9==NULL)  exit(1);  % probably because of "Could not launch or contact ds9, giving up"...

  ds9.put_array(pileimg);

  % Let's apply a World Coordinate System
  ds9.set_wcs(struct {
                crpix = [1,1],
                crval = [1,1],
                cdelt = [1,1],
                ctype = ["linear","linear"],
                cunit = ["pixel","pixel"]
              }; alt='p', append);

                        % Need to check how generic these keyword names are
  ds9.set_wcs(struct {
                crpix = [fits_read_key(filename,"TCRPX10","TCRPX11")],
                crval = [fits_read_key(filename,"TCRVL10","TCRVL11")],
                cdelt = [fits_read_key(filename,"TCDLT10","TCDLT11")],
                ctype = [fits_read_key(filename,"TCTYP10","TCTYP11")],
                cunit = ["deg","deg"]
              });
  ds9.set_cmap("hsv");

  % Make sure region formats and coordinate systems are the right ones
  ds9.send("regions system physical");
  ds9.send("regions format ds9");

  % Get the ds9 selected region and properly parse it
  variable dsreg=String_Type[1],ldsr=1;

  ()=fprintf(stdout,"\n  Create include/exclude regions (circle, annulus, ellipse, box) using physical coordinates.\n");

  variable reg_all = ["annulus","circle","ellipse","box"];
  variable nreg_type = length(reg_all);
  variable phys_coord, reg_type, reg_add, reg_par;
  variable key = "n";
  while(key != "q" and key != "y")
  {
    ()=fprintf(stdout,"\n  Enter \'y\' when finished, or \'q\' to finish without a region file:\n");
    key = getkey();
    if(key=="y")
    {
      dsreg = ds9.get_regions(;system="physical",format="ds9");
      ldsr = length(dsreg);

      phys_coord = Integer_Type[ldsr];
      reg_type = Integer_Type[2*ldsr];
      reg_add = Integer_Type[2*ldsr];
      reg_par = Array_Type[2*ldsr];

      % There better be at least two lines in the region file
      if(ldsr<2)
      {
        key="n";
        ()=fprintf(stdout,"\n   Hey!   Could not parse region\n");
      }
      else
      {
        _for i (0,ldsr-1,1)
        {
          % One line better say that these are physical coordinates
          phys_coord[i] = is_substr(dsreg[i],"physical");
          % Check to see what kinds of regions are present, and
          % if they are include or exclude regions
          _for j (0,nreg_type-1,1)
          {
            if(is_substr(dsreg[i],reg_all[j]))
            {
              reg_type[i] = max([1,j]);
              reg_add[i] = 1;
              reg_par[i] = strtok(dsreg[i],"[(,)]");

              % Turn the region string into just its parameters
              reg_par[i] = atof( reg_par[i][[1:length(reg_par[i])-1]] );

              if(is_substr(dsreg[i],"-"+reg_all[j]))  reg_add[i] = -1;

              % If annulus, break it up into its include and exclude regions
              if(j==0)
              {
                reg_type[i+ldsr] = reg_type[i];
                reg_add[i+ldsr] = -1;
                reg_par[i+ldsr] = reg_par[i][[0:2]];
                reg_par[i] = [ reg_par[i][[0:1]], reg_par[i][3] ];
              }

              % If ellipse or box, check to see if it's an annulus
              % and if so, break it up into its include and exclude regions
              if(j>1 and length(reg_par[i])==7)
              {
                reg_type[i+ldsr] = reg_type[i];
                reg_add[i+ldsr] = -1;
                reg_par[i+ldsr] = [ reg_par[i][[0:3]], reg_par[i][6] ];
                reg_par[i] = [ reg_par[i][[0:1]], reg_par[i][[4:6]] ];
              }
            }
          }
        }
        if(max(reg_add)<1)
        {
          key="n";
          ()=fprintf(stdout,"\n      Must be at least one include region\n");
        }
        else if(max(reg_type)<1)
        {
          key="n";
          ()=fprintf(stdout,"\n      Could not parse region(s)\n");
        }
        else if(max(phys_coord)<1)
        {
          key="n";
          ()=fprintf(stdout,"\n      File coordinates must be in physical coordinates.\n");
        }
        else if(key!="q")
        {
          fp = fopen(outfile+".reg","w");
          % We have to make sure the header comes first, plus
          % regions come next, then annulus regions, and minus
          % regions come last
          foreach j (where(reg_add[[0:ldsr-1]]==0))
            ()=fprintf(fp,"%s\n",dsreg[j]);
          foreach j ( where(reg_add[[0:ldsr-1]]==1 and
                            reg_add[[ldsr:2*ldsr-1]]!=-1) )
            ()=fprintf(fp,"%s\n",dsreg[j]);
          foreach j ( where(reg_add[[0:ldsr-1]]==1 and
                            reg_add[[ldsr:2*ldsr-1]]==-1) )
             ()=fprintf(fp,"%s\n",dsreg[j]);
          foreach j (where(reg_add[[0:ldsr-1]]==-1))
             ()=fprintf(fp,"%s\n",dsreg[j]);
          ()=fclose(fp);
        }
      }
    }
  }

  if(key=="q")
  {
    () = printf("\n");
    return;
  }

  ikeep = Integer_Type[length(xpix)];
  variable ireg;

  % Loop through the include regions
  foreach i (where(reg_add==1))
  {
    switch(reg_type[i])
    { case 1:
        ireg = circle(xpix,ypix,reg_par[i][0],reg_par[i][1],reg_par[i][2]);
    }
    { case 2:
        ireg = ellipse(xpix,ypix,reg_par[i][0],reg_par[i][1],reg_par[i][2],
                                 reg_par[i][3],reg_par[i][4]);
    }
    { case 3:
        ireg = box(xpix,ypix,reg_par[i][0],reg_par[i][1],reg_par[i][2],
                             reg_par[i][3],reg_par[i][4]);
    }
    if(length(ireg) > 0)  ikeep[ireg]=1;
  }

  % Loop through and then exclude the minus regions
  foreach i (where(reg_add==-1))
  {
    switch(reg_type[i])
    { case 1:
        ireg = circle(xpix,ypix,reg_par[i][0],reg_par[i][1],reg_par[i][2]);
    }
    { case 2:
        ireg = ellipse(xpix,ypix,reg_par[i][0],reg_par[i][1],reg_par[i][2],
                                 reg_par[i][3],reg_par[i][4]);
    }
    { case 3:
        ireg = box(xpix,ypix,reg_par[i][0],reg_par[i][1],reg_par[i][2],
                             reg_par[i][3],reg_par[i][4]);
    }
    if(length(ireg) > 0)  ikeep[ireg]=0;
  }

  % Estimate the pileup parameters for what's left, then save
  if(max(ikeep) > 0)
  {
    variable ilose = where(ikeep==0);
    ikeep = where(ikeep > 0);
    variable lm = sum(piximg[ikeep]);
    variable ls = sum(piximg[ikeep]^2);
    variable lc = sum(piximg[ikeep]^3);
    variable lavg = lm*lc^2/ls^3;
    variable fpsf = ls^2/lm/lc;
    variable neqv = ls^3/lc^2/9;

    if(length(ilose) > 0 and show_select)
    {
      pileimg[ilose] = 0.;
      ds9.put_array(pileimg);
    }

    variable text = "\n"
           +sprintf("Suggested pileup model parameters for %s,\n", filename)
           +sprintf("using region file %s.reg and filter file %s.flt:\n", outfile, outfile)
           +        "\n"
           +sprintf("   Frame Time (set_frame_time)   : %6.4f\n", snapti)
           +sprintf("   Effective piled PSF fraction  : %6.4f\n", fpsf)
           +sprintf("   Equivalent number of regions  : %6i\n", int(neqv))
           +sprintf("   Max equivalent pileup fraction: %6.4f * %6.4f = %6.4f\n", lavg, fpsf, lavg*fpsf)
           +        "\n";
    ()=fprintf(stdout, text);
    fp = fopen(outfile+".parms","w");
    ()=fprintf(fp, text);
    ()=fclose(fp);
  }
  else
    ()=printf("Empty region.\n\n");
}