%; Time-stamp: <2005-10-07 15:03:15 dph> 
%; Directory: ~dph/libisis
%; File: tgplot.sl  (was qgsp.sl)
%; Author: D. Huenemoerder
%; Original version: 2003.10.11
%;====================================================================
% version: 1.1  remove namespace conditional; not needed in ciao3
%
% purpose: interactive plotting/measuring of chandra grating spectra
%
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%% 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
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tgplot.sl starts here
%; Time-stamp: <2003-12-10 22:39:01 dph> 
%; Directory: ~dph/libisis
%; File: tgplot.sl
%; Author: D. Huenemoerder
%; Original version: 2003.12.06
%;====================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% version: 0.1  2003.12.06
% version: 0.2  2003.12.09 while loop added to keep plot selections repeating.
% version: 0.3  2003.12.10 add harcopy feature.  Requires caching much info.
%                          add minimal error checking on cursor input.
%                          add isis version check for data .file tag.
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% purpose: quick grating spectrum plot
%
%  Given a pha file name, arf name, and rmf name, 
%   load the data and responses,
%   plot the full spectrum (counts)
%   print the flux for the full spectrum
%   interactively select a subset range and plot it
%   overplot convolved default plasma model
%   interactively select a feature, print feature flux.
%   print lines in feature region (flux listed is for the MODEL).
%   Optionally allow multi-thermal plasma model specification.
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% SIDE EFFECTS:
%   alters line style, char size
%   ignores all loaded histograms.
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% LIMITATIONS:
%  Does not (and won't) support LETGS higher orders, since would need
%    both 1st order and high order rmfs to plot properly.
%  The ad hoc model is a thermal plasma. (Sorry all you XRB, AGN, PI folks!) 
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% NOTES:
%  For pha Type I files, use row=1.
%  Using a pager can be confusing. I use PAGER=cat and all output goes to the same window. 
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% public functions:
%
%  tgplot               % main interactive loop
%  get_pha_row()      % utility to get file row given order and grating
%  tgplot_plasma()      % set plasma temperature(s) and norm(s)
%  tgplot_hardcopy()    % ps file of last plot made to screen.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
%  EXAMPLES/NOTES:
%
% evalfile("tgplot.sl");
%
% fpha = "Capella_HA_1/P_Standard/pha2";
%
%%    A grating "PHA" fits file name.
%%    May be for HETG/ACIS-S, LETG/HRC-S, or LETG/ACIS-S
%
% farf = "Capella_HA_1/P_Standard/evt1aMEG_-1_garf.fits";
%
%%    A grating ARF.  It need not be one made for the PHA file,
%%    but it must match in grating-type and in spectral order.
%
% frmf = "./Ard/meg_m1_rmf.fits";
%
%%    A grating RMF file name.  It should match the grating and order.
%%    You can use -1 interchangeably with +1.
%
% row = get_pha_row( fpha, "MEG", -1 );
%
%%  Use this function if you don't know what row to use.
%%  The row in the PHA file which selects the order of iterest.
%%  For Type I PHA files, specify row 1.
%
% tgplot(fpha, row, farf, frmf );
%
%%   Starts the main interactive display loop
%%   Some info is typed to the terminal;
%%   The top half of a two panel plot appears;
%%   A cursor awaits marking a region for a zoomed plot in the lower panel.
%%     Click on the lower-left and upper-right corners of the desired region
%%   The lower plot is drawn, and a default model (2T) overplotted.
%%     Some info is typed to the terminal.
%%     The cursor awaits marking of a rectangle bounding a feature.
%%     Click on the lower-left and upper-right corners of the desired region
%%   The region is marked, and fluxes printed
%%          ( pflux = [photons/cm^2/s]; cps = [counts/sec] )
%%   The plot pauses for input; "q" quits. Any click starts over prompting for
%%     a new region on the top panel.
%%   
% logT = [6.5, 7.0, 7.5];   % define a logT array for plasma model
% norms = [0.5, 1.0, 0.8];  % define normalizations at each logT
% tgplot_plasma( logT, norms ) % define an aped model
%%  Define a multithermal model for the purposes of brightest line
%%    evaluation, and for a provisional folded model to overplot.
%%    The logT and norms arrays must have the same number of elements,
%%    but the length is arbitrary.
%%    No arguments will default to a 2T model, 6.8 and 7.4 with equal norms.
%%    The next call of tgplot will use the new model.
%
% qgsq_hardcopy("My_File.ps");
%
%%  Make a postscript file of the last thing plotted to the screen.
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
private variable tgplot_is_init = NULL;

private variable MEG_wmin = 1.5,
                HEG_wmin = 1.5,
                LEG_wmin = 1.0,
                MEG_wmax = 26.0,
                HEG_wmax = 19.5,
                LEG_wmax = 190.0;

private variable HEG=1, MEG=2, LEG=3;
private variable hid=-1, aid=-1, rid=-1;  % indices for counts, arf, rmf
private variable target; 
private variable wmin, wmax ;            % default full wavelength range.
private variable cwmin, cwmax, cymin, cymax; % cursor output ranges.
private variable pmdl = default_plasma_state(),   pmdl_name ;
private variable det=NULL, tg=NULL;
private variable saved_fit_fun=NULL, saved_params=NULL, saved_verbose=NULL;
private variable pwin = -1;
private variable g_all, s_lbl=line_label_default_style ;
private variable elems =  [C, N, O, Ne, Mg, Al, Si, S, Ar, Ca, Fe, Ni];
private variable sgtyp = ["", "HEG", "MEG", "LEG"];
private variable ttl ;
private variable texp ;
private variable def_charsize = 2.0;
private variable def_lbl_size = 1.67;

% variable to hold hardcopy info:
%
private variable tgplot_state =
    struct { pha,
             row,
             arf,
             rmf,
             wmin,      wmax, 
             wmin_reg,  wmax_reg,
             ymin_reg,  ymax_reg,
             wmin_feat, wmax_feat,
             ymin_feat, ymax_feat,
             flux_feat,
             flux_spec,
             title,
             fun, params,
             g_lines
            };

pmdl.temperature = NULL;           % to allow check for init.

private define plt_setup()
{
    if (pwin > 0) close_plot(pwin);
    putenv("PGPLOT_BACKGROUND=white");     % invert default for better view.
    putenv("PGPLOT_FOREGROUND=black");
    pwin = open_plot("/xwin", 1, 2 );   % set up a plot window
    resize( 24, 0.65 );
    _pgsvp(0.05, 0.98, 0.25, 0.85);
    charsize(def_charsize);
    s_lbl.char_height = def_lbl_size ;
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
define get_pha_row()
{
    if (_NARGS != 3 )
    {
	message("%% row = get_pha_row( file_pha, grating_type, order )" );
	message("%%    file_pha = grating pha file name");
	message("%%    grating_type = \"H\" | \"M\" | \"L\" ");
	message("%%    order != 0 ");
	return -1;
    }

    variable ord = () ;
    variable gtyp = () ; 
    variable fpha = () ;

    if ( ord == 0 ) error(" order=0 not supported.");

    gtyp = substr( strlow( gtyp ), 1, 1 ) ;
    if ( (gtyp != "h") and (gtyp != "m") and (gtyp != "l") )
         error("Bad grating type");

    variable agtyp = Assoc_Type[Integer_Type];
    agtyp["h"] = HEG;
    agtyp["m"] = MEG;
    agtyp["l"] = LEG;
    
    variable tmp_id = load_data( fpha );
    variable dinfo = get_data_info( tmp_id ) ;
    variable p = array_struct_field( dinfo, "part" );
    variable o = array_struct_field( dinfo, "order" );

    delete_data( tmp_id );

    variable row;
    row = where( p==agtyp[gtyp] and o==ord );

    switch ( length(row) ) 

    {
	case 0:
	vmessage("%% WARNING: no match for %s, order %d, in %s.",
	   sgtyp[ agtyp[gtyp] ], ord, fpha ) ;
	row = -1 ;
    }
    {
	case 1:
	row = row[0] + 1 ;
    }
    {
	vmessage("%% WARNING: Found more than one match for %s, order %d, in %s.",
	   sgtyp[ agtyp[gtyp] ], ord, fpha ) ;
	row = row[0] + 1 ;
    }

    return row;
}

% read the pha file header to determine the grating type and detector.
%
private define get_det_tg( pha )
{
    variable d, g;
    (d, g) = fits_read_key( pha, "detnam", "grating" );
    return d,g;
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% for use with ots arf and arbitrary exposure, force arf exposure time
% to match data exposure time.
% Do this if the two times differ by more than 10%.
%
private define texp_force( nh, na )
{
    variable atexp, dtexp;

    atexp = get_arf_exposure( na );
    dtexp = get_data_exposure( nh );

    if ( abs( 1.0 - atexp/dtexp ) > 0.1 )
    {
	set_arf_exposure( na, dtexp );
	message("%% Arf exposure differs from data exposure.");
	vmessage("%% Forced arf(%d) exposure to %f seconds", na, dtexp );
    }
    return dtexp ; 
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% set a wavelength range for the full spectrum based on
% a priori defaults.
% 
private define default_xrange( id )
{
    switch( get_data_info(id).part )
    {
	case HEG:
	wmin = HEG_wmin ;
	wmax = HEG_wmax ;
    }
    {
	case MEG:
	wmin = MEG_wmin ;
	wmax = MEG_wmax ;
    }
    {
	case LEG:
	wmin = LEG_wmin ;
	wmax = LEG_wmax ;
	% handle ACIS/LETG case 
	if ( char(det[0])=="A" ) wmax = 2.5*MEG_wmax;
    }

    % handle ACIS/LETG case 
%    if ( char(det[0])=="A" and char(tg[0])=="L" ) wmax = 2.5*MEG_wmax;

    tgplot_state.wmin = wmin;
    tgplot_state.wmax = wmax;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define reg_eflux ( h, wlo, whi )
{
  variable d, elo, ehi, l;

  d = _A( get_data_flux ( h ) ) ;   % assumes loaded, flux corrected.
  (elo, ehi) = _A( wlo, whi );
  l = where( d.bin_lo > elo   and  d.bin_hi <= ehi );

 % "integrate" energy flux [ ergs/cm^2/s ]
  return sum  ( (d.value * (d.bin_lo + d.bin_hi) / 2.)[l]) * 1.609e-9 ; 
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define tgplot_reset()
{
    if (hid > 0 ) { delete_data( hid ) ; hid = -1 ; }
    if (aid > 0 ) { delete_arf( aid )  ; aid = -1 ; }
    if (rid > 0 ) { delete_rmf( rid )  ; rid = -1 ; }

    if ( saved_fit_fun != NULL )
    {
	fit_fun( saved_fit_fun );
	set_params( saved_params );
    }
    title("");
    Isis_Verbose = saved_verbose;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define calc_fluxes( id, w1, w2, y1 )
{
    variable f = struct {wmin, wmax, ymin,  pflux, eflux, crate, cts };
    f.wmin = w1 ;
    f.wmax = w2 ;
    f.ymin = y1 ;
    f.pflux = region_flux(id, w1, w2, y1, y1 ).sum ;
    f.eflux = reg_eflux(id, w1, w2) ;
    f.crate = region_counts(id, w1, w2, y1, y1 ).sum / texp ;
    f.cts = region_counts(id, w1, w2, y1, y1).sum ;

    return f;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define pr_fluxes( f )
{
    message("");
    () = printf("%% Flux(%4.2f,%4.2f, y>%.0f) = %8.2e photons/cm^2/s\n", 
                 f.wmin, f.wmax, f.ymin, f.pflux );
    () = printf("%% 		  = %8.2e ergs/cm^2/s\n", f.eflux ) ;
    () = printf("%% 		  = %8.2e counts/s\n", f.crate );
    () = printf("%% 		  = %8.0f counts\n", f.cts ) ;
    message("");
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define mk_fit_fun( nam, idx, plas )
{ 
% e.g.,   mk_fit_fun( "xaped", 1, pmdl );
%
    variable ii;
    pmdl_name = sprintf("%s(%d)", nam, idx );  % store globally.
    create_aped_fun( nam, plas );
    fit_fun( pmdl_name );
    if (length(plas.norm) > 1 )
    {
	for (ii=0; ii<length(plas.norm); ii++)
	      set_par( sprintf("%s.norm%d", pmdl_name,ii+1),
	                plas.norm[ii], 0, 0, 10.0 );
	tie(1, [2:length(plas.norm)]);
    }
    else
    {
	set_par( sprintf("%s.norm", pmdl_name), plas.norm, 0, 0, 10.0 );
    }
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
define tgplot_plasma()
{
% tgplot_plasma([ log10_Ts[, norms] ] );

    variable i,n;
    variable logT = 6.8, norm = 0.01 ;
    switch(_NARGS)
    {
	case 0:
	pmdl.temperature = 10^logT ;
	pmdl.norm = norm ;
    }
    {
	case 1:
	logT = ();
	pmdl.temperature = 10.^logT;
	pmdl.norm = logT*0.0 + norm;
    }
    {
	case 2:
	pmdl.norm = ();
	logT = ();
	pmdl.temperature = 10^logT;
	if ( length(logT) != length(pmdl.norm) )
	{
	    message("%% WARNING: array of norms different length from logT's");
	    message("%%         Defaulting norms to match logTs.");
	    pmdl.norm = logT*0.0 + 0.01 ;
	}
    }    
    {
	message("%% USAGE:  tgplot_plasma([ log10_Ts[, norms] ] );" );
	return ; 
    }
    mk_fit_fun( "xaped", 1, pmdl );
    tgplot_state.fun = get_fit_fun;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define tgplot_init()
{
    saved_fit_fun = get_fit_fun;
    saved_params = NULL ; 
    if ( saved_fit_fun != NULL ) saved_params = get_params;
    if (pmdl.temperature == NULL) tgplot_plasma([6.8, 7.4]);
    tgplot_is_init = 0 ;
    saved_verbose = Isis_Verbose;
    Isis_Verbose = -1 ;
    plt_setup;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define eval_linelist(w1, w2)
{
% eval function on some grid to get the "brightest" set, g_all, defined:
    variable g, i, x1,x2;

    (x1,x2) = linear_grid(wmin, wmax, 8192);
    () = eval_fun(x1,x2);

    g = brightest(10, where(el_ion(Fe) and wl(w1, w2/2.) ) );
    g = [g,brightest(10, where(el_ion(Fe) and wl(w2/2., w2) ) )];

    for (i=0; i<length(elems); i++)
    {
	g = [g, brightest(6, where(el_ion(elems[i]) and wl(w1, w2) ) ) ];
    }
    return g;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
%  load pha row, corresponding arf, rmf.
%  get detnam and grating name from header.
%
private define tgplot_read()
{
    if (_NARGS != 4 )
    {
	message("%% USAGE: tgplot_read( fpha, row, farf, frmf ); ");
	message("%%  Load spectrum row and corresponding arf and rmf");
	return;
    }

    variable frmf = ();
    variable farf = ();
    variable row = ();
    variable fpha = (); 

    if (tgplot_is_init == NULL) tgplot_init;
    if (hid > 0) tgplot_reset;

    hid = load_data( fpha, row ) ;
    (det, tg) = get_det_tg( fpha ) ;
    default_xrange( hid ) ;
    g_all = eval_linelist( wmin, wmax );

    aid = load_arf( farf ) ;
    rid = load_rmf( frmf ) ;

    % check order, grating type
    if (  get_data_info(hid).part != get_arf_info(aid).part )
    {
	tgplot_reset;
	error("Data and ARF grating mismatch.");
    }
    if ( (get_data_info(hid).order != get_arf_info(aid).order )

% NOTE: not supporting LETGS higher orders.  Would require 1st order rmf to plot.
%         and (char(det[0]) != "H")  
       )
    {
	tgplot_reset;
	error("Data and ARF order mismatch.");
    }

    assign_rsp( aid, rid, hid );

    ignore(all_data);
    xnotice( hid, wmin, wmax );
    texp = texp_force( hid, aid );   % set global variable.
    flux_corr( hid );

    if (_isis_version < 10107)
    {
      ttl = sprintf("File: %s  Target: %s  Order: %d  Grating: %s",
	fpha, 
        get_data_info(hid).target,
        get_data_info(hid).order,
        sgtyp[ get_data_info(hid).part ] );
    }
    else
    {
      ttl = sprintf("File: %s  Target: %s  Order: %d  Grating: %s",
        get_data_info(hid).file,
        get_data_info(hid).target,
        get_data_info(hid).order,
        sgtyp[ get_data_info(hid).part ] );
    }

    title( ttl );
    
    variable f = calc_fluxes(hid, wmin, wmax, 0.0 ); 
    pr_fluxes( f );

    tgplot_state.pha = fpha;
    tgplot_state.row = row;
    tgplot_state.arf = farf;
    tgplot_state.rmf = frmf;
    tgplot_state.flux_spec = @f;
    tgplot_state.title = ttl;
    tgplot_state.g_lines = @g_all;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define plot_full()
{
    xrange( wmin, wmax );
    yrange(0);
    xnotice( hid, wmin, wmax );
    title( ttl );
    plot_data_counts( hid, 1 ) ;
    plot_group( g_all, 15, s_lbl ) ;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define plot_reg()
{
    message("\n%% Mark lower-left, upper-right  of rectangular region for zoom: \n");
    cwmin=0;cwmax=0;cymin=0;cymax=0;
    (cwmin, cwmax, cymin, cymax) = cursor_box();

    if (cwmin+cwmax+cymin+cymax == 0 )  %  user error - typed "q" to quit...
    {
	tgplot_reset();
	return 1;
    }

    linestyle(2);
    set_line_width(3);
    oplot([cwmin,cwmax,cwmax,cwmin,cwmin], [cymin,cymin,cymax,cymax,cymin], 2);
    linestyle(1);

    xrange( cwmin, cwmax) ;
    yrange( cymin, cymax );
    xnotice( hid, cwmin, cwmax ) ;
    () = renorm_counts;

    set_line_width(5);
    title("");
    plot_data_counts( hid, 1 );
    set_line_width(1);
    oplot_model_counts( hid, 2 );
    set_line_width(1);
    color(1);
    plot_group( brightest( 18, where( wl(cwmin, cwmax) ) ), 15, s_lbl );

    tgplot_state.wmin_reg = cwmin;
    tgplot_state.wmax_reg = cwmax;
    tgplot_state.ymin_reg = cymin;
    tgplot_state.ymax_reg = cymax;
    tgplot_state.fun = get_fit_fun;
    tgplot_state.params = get_params;

    return 0;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define tgplot_plot()
{
    variable flag;
    plot_full;
    flag = plot_reg;
    return flag;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
private define tgplot_reg_flux()
{
    message("\n%% Mark lower-left, upper-right of rectangular region for flux estimate:\n");
    (cwmin, cwmax, cymin, cymax) = cursor_box();
    if (cwmin+cwmax+cymin+cymax == 0 )  %  user error - typed "q" to quit...
    {
	tgplot_reset();
	return 1 ;
    }
	
    linestyle(2);
    set_line_width(1);
    oplot([cwmin,cwmax,cwmax,cwmin,cwmin], [cymin,cymin,cymax,cymax,cymin], 3);
    linestyle(1);
    page_group(brightest(4,where(wl(cwmin,cwmax))));
    variable f = calc_fluxes(hid, cwmin, cwmax, cymin) ; 
    pr_fluxes( f );

    color(1); charsize(1.5);
    xylabel( cwmin, cymax/2., sprintf("%8.2e pflux",f.pflux), 0, 1.0 );
    xylabel( cwmax, cymax/2, sprintf("%8.2e cps", f.crate), 0, 0.0 );
    xylabel( (cwmin+cwmax)/2, cymax+(cymax-cymin)*0.1, sprintf("%.0f cts",f.cts), 0, 0.5 );
    charsize(def_charsize);

    tgplot_state.wmin_feat = cwmin;
    tgplot_state.wmax_feat = cwmax;
    tgplot_state.ymin_feat = cymin;
    tgplot_state.ymax_feat = cymax;
    tgplot_state.flux_feat = @f;

    return 0;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
define tgplot()
{
    if (_NARGS != 4 )
    {
     error("USAGE:  tgplot( file_pha, row, file_arf, file_rmf )");
    }

    variable x, y, c ;
    variable frmf, farf, row, fpha;

    frmf=();
    farf=();
    row=();
    fpha=();

    tgplot_init();
    plt_setup;
    tgplot_read( fpha, row, farf, frmf );


    c = "x";
    while ( strlow(c) != "q" )
    {
	if ( tgplot_plot())  break ;

	!if ( tgplot_reg_flux() )
	{
	    message("\n%% ------ Click anywhere to repeat.  TYPE 'q' to quit. ------\n");
	    cursor( &x, &y, &c );
	}
	else break;
    }

    tgplot_reset();
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% attempt hardcopy, using stored info.
% may happen after reset - hence, need stored copies of all 
% relevant variables.  Use global tgplot_state structure.
%  (Should consolidate repeated code (from screen plot functions)
%    into new functions)
%
define tgplot_hardcopy()   % (outfile)
{
    variable outfile ;

    if (_NARGS != 1 )
    {
	message("%% USAGE:  tgplot_hardcopy( \"psfile\");");
	return;
    }
    
    outfile = () ; 

    variable pid = dup_plot(outfile + "/cps", pwin);
    resize( 24, 0.65 );
    _pgsvp(0.05, 0.98, 0.25, 0.85);
    charsize(def_charsize);
    s_lbl.char_height = def_lbl_size ;

    tgplot_read( tgplot_state.pha, tgplot_state.row, tgplot_state.arf, tgplot_state.rmf);
    plot_full;

    cwmin = tgplot_state.wmin_reg;
    cwmax = tgplot_state.wmax_reg;
    cymin = tgplot_state.ymin_reg;
    cymax = tgplot_state.ymax_reg;

    linestyle(2);
    set_line_width(3);
    oplot([cwmin,cwmax,cwmax,cwmin,cwmin], [cymin,cymin,cymax,cymax,cymin], 2);
    linestyle(1);

    xrange( cwmin, cwmax ) ;
    yrange( cymin, cymax );
    ignore(all_data);
    xnotice( hid, cwmin, cwmax ) ;
    fit_fun( tgplot_state.fun );
    set_params( tgplot_state.params );
    () = eval_counts;

    set_line_width(5);
    title("");
    plot_data_counts( hid, 1 );
    set_line_width(1);
    oplot_model_counts( hid, 2 );
    set_line_width(1);
    color(1);
    plot_group( brightest( 18, where( wl(cwmin, cwmax) ) ), 15, s_lbl );

    cwmin = tgplot_state.wmin_feat;
    cwmax = tgplot_state.wmax_feat;
    cymin = tgplot_state.ymin_feat;
    cymax = tgplot_state.ymax_feat;

    linestyle(2);
    set_line_width(1);
    oplot([cwmin,cwmax,cwmax,cwmin,cwmin], [cymin,cymin,cymax,cymax,cymin], 3);
    linestyle(1);
    variable f = tgplot_state.flux_feat;

    color(4); charsize(1.5);
    xylabel( cwmin, cymax/2., sprintf("%8.2e pflux",f.pflux), 0, 1.0 );
    xylabel( cwmax, cymax/2, sprintf("%8.2e cps", f.crate), 0, 0.0 );
    xylabel( (cwmin+cwmax)/2, cymax+(cymax-cymin)*0.1, sprintf("%.0f cts",f.cts), 0, 0.5 );
    charsize(def_charsize);

    _pgiden;

    close_plot(pid);
    window(pwin);

    tgplot_reset() ;
}
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

provide("get_pha_row");
provide("tgplot_plasma");
provide("tgplot");
provide("tgplot_hardcopy");

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

plasma(aped);   %%% make sure database is loaded.

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

% tgplot("Capella_HA_1/P_Standard/pha2", 3, "Capella_HA_1/P_Standard/evt1aMEG_-1_garf.fits", "Ard/meg_m1_rmf.fits");

%% try different model....
% tgplot_plasma( [6.8,7.2,7.5,8.0], [0.1, 0.5, 1.0, 0.1]);
% tgplot("Capella_HA_1/P_Standard/pha2", 3, "Capella_HA_1/P_Standard/evt1aMEG_-1_garf.fits", "Ard/meg_m1_rmf.fits");

%% try "ots" arf...:
% tgplot_plasma(6.8, 0.01);
% tgplot("Capella_HA_1/P_Standard/pha2", 4, "VW_Cep_HA/evt1aMEG_1_garf.fits", "Ard/meg_p1_rmf.fits");

%% LETGS
% tgplot("VW_Cep_LH/P_Standard/pha2", 1, "VW_Cep_LH/P_Standard/evt1aLEG_-1_garf.fits", "Ard/leg_m1.rmf");

%%%%%%%%%%%%%%%% tgplot.sl ends here