%
% file: get_streak_events.sl
% 

% Modified version of findzo.sl to return the coordinates of the rotated
% streak events, "vertical in new frame".
%
% Usage:  (xs, ys) = get_streak_events("filename","m",4,gx,gy);
%


define get_streak_events( )
{

   vmsg(1, "Running FINDZO version "+ string(_version[0])+
		"."+string(_version[1])+"."+string(_version[2])+"\n");
%
% here is one hack for isis 1.2.6 or later.
%

   set_fit_method("marquardt;delta=1.0e-12");

%
% Initialize Variables
%
   variable fitsfile, gtype, ftype;
   variable gx, gy;
   variable x0, y0;
   variable boxl_strc, boxw_strc, rad_strc; 
   variable boxl_garm, boxw_garm, rad_garm, bin_garm; 
   variable detnam="INSTRUME";
   variable deg_rad=3.1416/180.;
   variable en, x, y, xp, yp, xpp, ypp;
   variable sxp, syp, sxp2, syp2, cnt, cnt2; 
   variable cltmp, l, ll, mxp, myp, strx, slp, islp, zx, zy;
   variable nx, ny;
   variable disp_bx_size, disp_n_elem;
   variable n_del, n_eng_hi, n_eng_lo; 
   variable p1, p2; 
   variable strc;

%
% note on angles: positive angle == clockwise rotation angle. 
%

   variable grang, rotang; 

%
% Read input from command line
%

   if ((_NARGS < 2) or (_NARGS == 4) or (_NARGS > 5))
   {
      vmsg(1, "  \n");
      vmsg(1, "USAGE in ISIS: (x, y) = findzo(fits_file,\"l\",[guessing_mode,gX,gY]);\n");
      vmsg(1, "    where \"l\" (LEG) can be replaced with \"h\" (HEG)\n"  );
      vmsg(1, "    or \"m\" (MEG). guessing_mode, gx and gy are optional.\n"       );
      vmsg(1, "  \n");
      vmsg(1, "USAAGE in CIAO: tg_findzo infile l guessing_mode gX gY\n");
      return -1;
   }

   if (_NARGS ==2 )
   {
      (fitsfile,gtype)=(); 
      ftype = 1; 
      x0=4096.5; y0=4096.5; 		% These are dummy values. Not actually used. 
   }

   if (_NARGS ==3 )
   { 
      (fitsfile,gtype,ftype)=();        % here, ftype can be any integer [1,2,3, or 4]. 
      if (ftype == NULL) 
      {
         vmsg(1,"ftype = NULL invalid\n");
         return -1; 
      }
      x0=4096.5; y0=4096.5;             % Initial guesses. Will not be used. 
   }

   if (_NARGS ==5 ) 
   {
      (fitsfile,gtype,ftype,x0,y0)=();  % here ftype should be 4. [x0, y0] needs 
					% to be specified. 
   }



%
% Now allowing to accept the use of cfitsio extended name syntax.
%
%
% The extention block name [EVENTS] should always be used for Chandra Event files. 
%

   variable tmp_name = strchop(fitsfile,'[',0);
   variable fitsname = tmp_name[0];
   variable fitssub = "[EVENTS]"; 

   if ( length(tmp_name) > 1 ) 
   {
      tmp_name = tmp_name[[1:length(tmp_name)-1]]; % remove the first element since 
                                                   % they should be a fits filename 

%
% Dummy proof 
%
      if ((strup(tmp_name[0]) == "EVENTS]") or (strup(tmp_name[0]) == "1]")) 
      {
	 vmsg(2,"                                             \n");
         vmsg(2,"#############################################\n");
         vmsg(2,"# Please do not specify the extention name. #\n");
         vmsg(2,"#############################################\n");
	 vmsg(2,"                                             \n");
         tmp_name = tmp_name[[1:length(tmp_name)-1]]; 
      }
%
%
      vmsg(2,"                                                                  \n");
      vmsg(2,"The extention block name is always set to be [EVENTS] by default. \n");
      vmsg(2,"                                                                  \n");

      if ( length(tmp_name) > 0 ) 
      {
         fitssub  = fitssub + "[" + strjoin(tmp_name,"[");
      }
   }


%
% Initialize the required parameters specific to this configuration. 
%
% gx, gy == guesstimated coordinate position of the source (in SKY)
% x, y   == event data arrays 
% en     == event energy
% grang  == grating arm angle 
% rotang == ajusted ROLL_NOM angle 
% detnam == DETECTOR name. 
% strc   == Streak Correction Factor
% box[w,l]_strc == Region Box Size for transfer streak (for fitting)
% rad_strc == Circle of radius exclusion zone around the zeroth order image. 
% box[w,l]_garm == Region Box Size for grating arms (for fitting)
% rad_garm == Circle of radius exclusion zone around the zeroth order image.
% bin_garm == binning width for histogramming. 
% disp_bx_size == display box size to show the fidelity of fit
% disp_n_elem == a total number of events shown in the pgplot screen.
% n_del = a level of statistical significance used in filtering (e.g., clip_array).
% n_eng_[hi,low] == Upper and lower threshold levels for energy filtering. 
%
   (gx, gy, x, y, en, grang, rotang, detnam, strc, boxw_strc, boxl_strc, rad_strc,
         boxw_garm, boxl_garm, rad_garm, bin_garm, disp_bx_size, disp_n_elem, n_del, 
         n_eng_hi, n_eng_lo) = gparam_io(fitsname, fitssub, gtype, ftype, x0 ,y0);

%
% Rotate the SKY event array (X,Y) by (-rotang) degree around (gx,gy).
% This will make the transfer streak go vertical in the new frame. 
%
   (xp,yp)=rotxy(x,y,gx,gy,-rotang*deg_rad);

%
% And draw a first region box to isolate transfer streak events  
% without zero order and with energy filtering. 
%
   l=where(xp > (gx-boxw_strc/2) and xp < (gx+boxw_strc/2) and 
	   yp > (gy-boxl_strc/2) and yp < (gy+boxl_strc/2) and
	   x > 0.0            and y    > 0.0         and  
           hypot(xp-gx,yp-gy)>rad_strc and (en > n_eng_lo and en < n_eng_hi));

   % offset w.r.t. zo for the returned event locations:
   return xp[l]-gx, yp[l]-gy;
}


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