%; Time-stamp: <2005-09-02 13:26:41 dph> 
%; Directory: ~dph/libisis
%; File:  zofind.sl
%; Author: D. Huenemoerder
%; Original version: 2003.11.10
%;====================================================================

% purpose:
%
%% prototype zo finding w/ vwhere() and isis array_fit().

% side effects:  alters current plot ranges, point and line styles,
%                labels.

% example:
%
%
% > require( "zofind" );
%
% > ( x0, y0 ) = zofind( "evt1" );   % read events, call vwhere; (follow instructures:)
%
% %  Box a trace (e.g., for MEG, two polygons, excluding zero order)......
% %  Mark another (e.g., ditto for the streak) ...
%
% > (x1,y1) = rd_zoxy( "reg1a", 1 );
%
% > (dx,dy,dr)=diff_zoxy(x0,y0,x1,y1,1);
% % dx = 0.15  dy = -0.10  dr = 0.18 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%% You must have the gtk module installed.
%%% It can be obtained from
%%% http://space.mit.edu/CXC/software/slang/modules/
%%%
%%% If your system does not have paths configured to automatically
%%% find it, you will need to  set appropriate environment variables, 
%%% or set paths in the slang interpreter via something like:
%%% 
%%%  set_import_module_path( "/nfs/cxc/a1/i686/opt/lib/slang/modules:" + get_import_module_path);
%%%  set_slang_load_path( "/nfs/cxc/a1/i686/opt/share/slsh/local-packages:" + get_slang_load_path );
%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% revision history:
%
% version 1.3.0
%          fix isis and package loading for ciao3.2

% version:1.2.0
%
%   2004.07.22    add ciao wrapper for use in sherpa.
%
%
% version: 1.1.0  2005.05.06 -
%              separate event reading to a function
%              check for HRC and read events properly.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


private variable _version = [1, 3, 0 ]   ; % major, minor, patch
variable zofind_version = sum( _version * [ 10000, 100, 1 ] ) ; 
variable zofind_version_string = sprintf("%d.%d.%d", 
                                     _version[0],
                                     _version[1],
                                     _version[2] ) ; 
%


require("vwhere");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%% CIAO/sherpa wrapper...    Conditionally load the isis and pgplot
%%% modules, if in a ciao shell.  Don't need if already in isis.
%
#ifnexists isis_exit
%
require("isis");
require("pgplot");
Isis_Verbose = -1;
%
#endif
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%% define linear function for array_fit:
%%%
private define line_fun (x, p)
{
    return p[0] + x * p[1];
}
%%%

%%%%
%
private define guess_line(x,y)
{
    % guess slope and intercept from range of x,y

    % Note that a small scatter in x will change a vertical line into
    % a very steep one.

    variable a, b, lp, lm;
    lp = (where( x==max(x) ))[0] ;
    lm = (where( x==min(x) ))[0] ;

    b = ( y[lp]-y[lm] ) / ( x[lp]-x[lm] ) ;
    a = y[lm] - b * x[lm] ;

    return a, b;
}
%


%%%%%  apply a simple sigma-clipping scheme to a linear fit:
%%  fit a line; discard points nsig from the median;
%%  refit; discard w/ smaller constraint.
%%
private define clip_fit(nsig, x, y )
{
    variable a, b, s, r, m, l, p;

    nsig = double( nsig ) ; % force to real

    (a,b) = guess_line( x, y );  % estimate slope, intercept

    (p,s) = array_fit( x, y, NULL, [a,b], NULL, NULL, &line_fun );
    a = p[0];
    b = p[1];

  % sigma clip:
    r =  y - line_fun( x, p );
    m = moment( r );
    l = where( abs( r - median(r) ) / m.sdev  < nsig );
 
   (p, s) = array_fit( x[l], y[l], NULL, [a,b], NULL, NULL, &line_fun);
    a = p[0];
    b = p[1];

  % repeat sigma clip, w/ tighter constraint...
    r =  y - line_fun( x, p );
    m = moment( r );
    l = where( abs( r - median(r) ) / m.sdev  < 0.5 );
 
   (p, s) = array_fit( x[l], y[l], NULL, [a,b], NULL, NULL, &line_fun);
    a = p[0];
    b = p[1];

    return a, b;
}
    

private define rd_acis_events( f_evt1 )
{
    variable x, y, c, g, e, s, l, xx, yy ;
    ( x, y, c, g, e, s ) = fits_read_col( f_evt1,
                     "x","y","ccd_id","grade","energy","status");
    l = where( (c==7 or c==6)
              and g<7
              and g!=1
              and g!=5
              and e<7000
              and e>500
              and s==0
              and x>0          % trap bad coords
              and y>0
              and x<1.e5
              and y<1.e5
            );

    xx = x[l];
    yy = y[l];

    variable oid, obj, dob, ttl;

    oid = fits_read_key(f_evt1, "obs_id");
    obj = fits_read_key(f_evt1, "object");
    dob = fits_read_key(f_evt1, "date-obs");

    ttl = sprintf("%s (%s) %s", obj, oid, dob );

    return xx, yy, ttl ; 
}

private define rd_hrc_events( f_evt1 )
{
    variable x, y, c, s, l, xx, yy ;
    ( x, y, c, s, xx ) = fits_read_col( f_evt1,
                     "x","y","chip_id", "status", "chipy");
    l = where( c==2
              and s==0
              and x>0          % trap bad coords
              and y>0
              and x<1.e5
              and y<1.e5
              and xx < 11000    % limit chipy range to center.
              and xx > 7000
            );

    xx = x[l];
    yy = y[l];

    variable oid, obj, dob, ttl;

    oid = fits_read_key(f_evt1, "obs_id");
    obj = fits_read_key(f_evt1, "object");
    dob = fits_read_key(f_evt1, "date-obs");

    ttl = sprintf("%s (%s) %s", obj, oid, dob );

    return xx, yy, ttl ; 
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%% put it all together; read events, call vwhere twice to get two 
%%%   groups of points to regress and intersect.
%%%   Input may be evt1; will be filtered on grade, reasonable coords, 
%%%    energies. 

define zofind( )
{
%  (x0,y0) = zofind( "AB_Dor/evt1" );

    if (_NARGS != 1)
    {
	error("USAGE:  (x0, y0) = zofind( evt_file ) ;") ;
    }

    variable f_evt1 = (); 
    variable xx, yy, ttl;

    variable instrume = fits_read_key( f_evt1, "instrume");

    if (instrume == "ACIS") ( xx, yy, ttl ) = rd_acis_events( f_evt1 );
    if (instrume == "HRC")  ( xx, yy, ttl ) = rd_hrc_events( f_evt1 );

    vmessage("%% Box a trace (e.g., for MEG, two polygons, excluding zero order)......");
    variable lm = vwhere( xx, yy );

    vmessage("%% Mark another (e.g., ditto for the streak) ...");
    variable ls = vwhere(xx,yy);

    if ( (length(lm) == 0) or (length(ls) ==0 ) )
    {
	verror ("No points selected.");
    }

    variable as, bs, am, bm ;

    % fit linear coefficients for each selection:
    %
    (as,bs) = clip_fit( 1.0,  xx[ls], yy[ls] );
    (am,bm) = clip_fit( 1.0,  xx[lm], yy[lm] );

    % generate points on solution, for overplot:
    variable x1, y1;
    x1 = [ min(xx): max(xx): 1.0 ] ;
    y1 = [ min(yy): max(yy): 1.0 ] ;

    % determine intercept:   (need /0 check)
    variable x0,y0;
    x0 = (as-am) / (bm-bs);
    y0 = as + bs*x0;

    % plot zero-order region events, overplot lines and intercept:
    variable dx=10, dy=10;
    set_line_width(1);
    connect_points(0); pointstyle(-1);
    xrange(x0-dx, x0+dx);
    yrange(y0-dy, y0+dy);

    label ("x", "y", ttl);
    
    plot( xx, yy );

    connect_points(1);
    oplot( x1, line_fun(x1, [as,bs] ), 2) ;
    oplot( x1, line_fun(x1, [am,bm] ), 2) ;
    connect_points(0);

    set_line_width(5);
    pointstyle(25); oplot(x0,y0, 3);
    pointstyle(-1);
    set_line_width(1);

    return x0, y0;
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the zero-order position from the src1a file and row
% ( there are as many rows as sources detected )
%
define rd_zoxy()  % (  f_src1a, nrow )
{
    variable f_src1a, nrow ;

    if (_NARGS != 2 )
    {
	message("%% rd_zoxy: missing/bad arguments.");
	message("%% USAGE:   (x0, y0) = rd_zoxy( f_src1a, nrow )");
	message("%%    Read zero order position from src1a or pha file, row nrow");
	return NULL ;
    }

    (f_src1a, nrow) = () ; 

    variable x0, y0 ;

    ( x0, y0 ) = fits_read_col( f_src1a, "x", "y" );

    return x0[nrow-1], y0[nrow-1] ; 
}

% compute deltas:
define diff_zoxy()  % ( x0, y0, x1, y1, verbose )
{
    variable x0, y0, x1, y1, verbose ;

    switch( _NARGS )
    {
	case 5:
	(x0, y0, x1, y1, verbose) = () ;
    }
    {
	case 4:
	(x0, y0, x1, y1) = () ;
	verbose = 0 ;
    }
    {
	message("%% diff_zoxy:  bad arguments.");
	message("%% USAGE:  (dx, dy, dr) = diff_zoxy( x0, y0, x1, y1[, verbose] );");
	message("%%   compute distances between coordinates.");
	return -1 ;
    }
	
    variable dx, dy, dr ;
    dx = x0 - x1 ;
    dy = y0 - y1 ;
    dr = hypot( dx, dy ) ;

    if ( verbose > 0 )
    {
	() = printf("%%  dx = %0.2f  dy = %0.2f  dr = %0.2f\n",
	             dx, dy, dr ) ;
    }
    return dx, dy, dr ; 
}

provide("zofind");
provide("rd_zoxy");
provide("diff_zoxy");

vmessage("%% zofind version %S loaded.", zofind_version_string ) ;