%; Time-stamp: <2007-03-18 20:46:55 dph> 
%; Directory: aglc/Examples/VelaX1/
%; File: Demo_1.sl
%; Author: D. Huenemoerder
%; Original version: 2003.12.15
%;====================================================================
% version: 1
%
% purpose: demonstrate/test light curve functions.
%
%%% spectral light curve, select on color-color

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

require("aglc");


_auto_declare = 1;       % be sloppy about declaring variables.



private define mydoc( str_array )
{
    () = array_map( Integer_Type, &printf, "%% >>>  %s\n", str_array ) ;
}

private define prsrc( str_array )
{
    () = array_map( Integer_Type, &printf, "%s\n", str_array ) ;
}

private variable vp = get_outer_viewport ; 
private define set_vp(x1, x2, y1, y2 )
{
    vp.xmin = x1 ;
    vp.xmax = x2 ;
    vp.ymin = y1 ;
    vp.ymax = y2 ;
    set_outer_viewport( vp ) ; 
}


putenv("PGPLOT_BACKGROUND=white");     % invert default for better view.
putenv("PGPLOT_FOREGROUND=black");

pw1 = open_plot("/xwin");
resize(17, 0.6);
set_vp(0.1,0.98,0.15,0.87);
charsize(1.75);
set_frame_line_width(5);
set_line_width(5);


f_evt1a = "indir/evt1a.tg";    % filtered grating events
f_evt0  = "indir/stat1.exp";    % exposure reference - stat1 file.

pd = 283.75;             % Vela X-1 period, in sec
dt = pd / 4.;            % oversample the period

% compute light curves in 4 bands: full, w1-w2, w2-w3, w3-w4;
%  for HEG+MEG, +/-1st orders:
%
g = ["H", "M"];          % grating parts to bin
o = [-1,1];              % grating orders to bin
w1 =  1.5;               % wavelength regions
w2 =  4.5;
w3 =  7.5;
w4 = 10.5;

msg = ["Read some event file columns into a structure. This saves time",
       "for repeated use.  Two files are required: an exposure number reference",
       "(stat1 preferred; or unfiltered events)", 
       "which is used to determine the exposure.  The other is for the",
       "grating coordinates.",
       "(Level 2 filtered data cannot be used for the exposure!)\n"];
mydoc(msg);
()= printf("s = aglc_read_events(\"%s\", \"%s\"); \n", f_evt1a, f_evt0);

s = aglc_read_events(f_evt1a, f_evt0);    % read events into a structure

% make light curves in bands, using the structure to save time re-reading
%
msg = ["Make light curves in bands, for HEG+MEG, +/-1st orders.",
       "Data to bin is in the structure; arguments specify binning,",
       "wavelength regions, gratings, and orders.",
       "(NOTE: for single-instance use, evt1a and evt0 file names can be",
       "       given instead of the structure.\n\n"];

mydoc(msg);
() = printf("c = aglc( s, %.0f, %0.1f, %0.1f, [\"%s\",\"%s\"], [%d,%d]);\n",
             dt, w1, w4, g[0], g[1], o[0], o[1] ) ;

msg= [" curve = aglc( evt_struct, time_bin, min_waves, max_waves, gratings, orders );",
      " min_waves and max_waves may be arrays or scalars.\n"];
mydoc(msg);

c =  aglc( s, dt, w1, w4, g, o );       % full band

msg = ["Here are the light curves for 1.5-10.5A (black)",
       "hplot(  c.time_min,  c.time_max,  c.count_rate, 1 );"];
mydoc(msg);

yrange(0); xrange(0);    set_line_width(5);
connect_points(1);
label("Time", "Count Rate", "Vela X1 phase=0.25");
hplot(  c.time_min,  c.time_max,  c.count_rate, 1 );

() = printf("c1 = aglc( s, %.0f, %0.1f, %0.1f, [\"%s\",\"%s\"], [%d,%d]);\n",
             dt, w1, w2, g[0], g[1], o[0], o[1] ) ;

c1 = aglc( s, dt, w1, w2, g, o );       % narrower bands...

mydoc([" ....... 1.5-4.5A (red)"]);
ohplot(c1.time_min, c1.time_max, c1.count_rate, 2 );


() = printf("c2 = aglc( s, %.0f, %0.1f, %0.1f, [\"%s\",\"%s\"], [%d,%d]);\n",
             dt, w2, w3, g[0], g[1], o[0], o[1] ) ;

c2 = aglc( s, dt, w2, w3, g, o );

mydoc(["........ 4.5-7.5A (green)"]);
ohplot(c2.time_min, c2.time_max, c2.count_rate, 3 );

() = printf("c3 = aglc( s, %.0f, %0.1f, %0.1f, [\"%s\",\"%s\"], [%d,%d]);\n",
             dt, w3, w4, g[0], g[1], o[0], o[1] ) ;

c3 = aglc( s, dt, w3, w4, g, o );

mydoc(["........ and 7.5-10.5A (blue)."]);
ohplot(c3.time_min, c3.time_max, c3.count_rate, 4 );

plot_pause;


% Form color-color ratio.
%
l = where( c1.count_rate>0 and c2.count_rate>0 );   % for non-zero denominator
v1 = c2.count_rate[l] / c1.count_rate[l] ;          % 6A/3A  low=>hard
v2 = c3.count_rate[l] / c2.count_rate[l] ;          % 9A/6A  low=>hard

msg = [ "Ratios of the curves in bands give us a \"color-color\" plot:\n" ] ;
mydoc( msg );
src = [" l = where( c1.count_rate>0 and c2.count_rate>0 );  %% for non-zero denominator", 
       "v1 = c2.count_rate[l] / c1.count_rate[l];           %% 6A/3A  lower=>harder", 
       "v2 = c3.count_rate[l] / c2.count_rate[l] ;          %% 9A/6A  lower=>harder", 
       "plot(v1,v2);" ] ; 
prsrc( src ) ; 

connect_points(0); pointstyle(-4); set_line_width(5);
xrange(0); yrange(0);
label("Hardness ratio: 6A/3A", "Hardness ratio: 9A/6A", "Vela X-1 phase 0.25");
plot(v1,v2,14);
color(4);
xylabel(0.02,0.02,"HARD");
color(2);
xylabel(0.8, 0.6, "SOFT", 1);
plot_pause;

lv_1 = where( v1 > 0.6 );
lv_2 = where( v1 <= 0.4 );

msg = ["We can select via an expression, ",
       "or use the gtk-module GUI, \"vwhere()\"",
       "Done automatically here via...\n\n"];
mydoc(msg);
()=printf("lv_1 = where( v1 > 0.6 );\n");
()=printf("lv_2 = where( v1 <= 0.4 );\n");

oplot( v1[lv_1], v2[lv_1], 2);    % overplot selections in color
oplot( v1[lv_2], v2[lv_2], 4);
plot_pause;

%%%% from each set of points, make a spectrum.
%%%% light curve structure has reverse indices:
%
% c contains light-curve and additional info, like reverse indices
% from binning.  aglc_filter selects elements in s's tags using those
% reverse indices.  The indices also contain the binning region - wavelengths,
% cross-dispersion, time, etc.  (So use the full list, s, and full curve,  c).


msg = ["The selection filters are then applied to the event structure, s, ",
       "using the \"reverse indices\" from the light curve histogram ",
       "stored in the light curve structure, c.",
       "The filter returns a new event structure.\n\n"];
mydoc(msg);
src = [ "s_v1 = aglc_filter( c, s, lv_1 );  %% filter s given c and lv_1.", 
        "s_v2 = aglc_filter( c, s, lv_2 ); %% " ] ; 
prsrc( src ) ; 

s_v1 = aglc_filter( c, s, lv_1 );      % filter s given c and lv_1.
s_v2 = aglc_filter( c, s, lv_2 );


mydoc(["Loading, interpolating,  and summing ARFs, binning spectra..."]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  load vela x-1 arfs, interpolate to MEG grid; sum (since the light curve 
%   is for H+M, -1/+1);  use MEG grid to bin selected events.

ia1 = load_arf("./indir/evt1aHEG_-1_garf.fits");
ia2 = load_arf("./indir/evt1aHEG_1_garf.fits") ;
ia3 = load_arf("./indir/evt1aMEG_-1_garf.fits");
ia4 = load_arf("./indir/evt1aMEG_1_garf.fits") ;

a3 = get_arf(ia3);
a1 = get_arf(ia1);

% total arf is sum of MEG and HEG and orders -1,1
%
atot = a3.value + get_arf(ia4).value + interpol( a3.bin_lo, a1.bin_lo, a1.value+get_arf(ia2).value );

x1 = a3.bin_lo;  % define grid, using standard MEG arf
x2 = a3.bin_hi;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% bin spectra.

s1 = double( histogram( s_v1.wave, x1, x2 ) ) ;
s2 = double( histogram( s_v2.wave, x1, x2 ) ) ;

e1 = sum( c.exposure[lv_1] );  % compute exposures
e2 = sum( c.exposure[lv_2] );

l = where( atot>0) ;               % don't div by 0
f1 = s1*0;                         % create flux arrays
f2 = s2*0;
f1[l] = s1[l] / atot[l] / e1 ;     % compute flux/bin from counts
f2[l] = s2[l] / atot[l] / e2 ;

msg = ["Given the filtered data, we can now bin spectra from the events,",
       "and load and sum arfs to compute flux from the hard and soft states.",
       "You can see that the hard (blue) state really does have",
       "a slightly steeper spectrum!"];
mydoc(msg);

src = [ "s1 = histogram( s_v1.wave, x1, x2 ) ; %% bin counts", 
	"e1 = sum( c.exposure[lv_1] ) ;        %% compute exposure", 
	"f1 = s1 / atot / e1 ;                 %% compute flux/bin from counts" ] ; 
prsrc( src ) ;


connect_points(1); set_line_width(5);
label("Wavelength", latex2pg("phot/cm^2/s/bin"), "Vela X-1 phase 0.25");
xrange(1.5, 10.5); yrange(1.e-6, 0.002); ylog; 
hplot(x1, x2, f1, 2);
ohplot(x1, x2, f2, 4);
ylin;
plot_pause;


msg=["We can also look at where the hard and soft selections lie ",
     "in the light curve itself.  We dynamically filter the ",
     "light curve array with the selections made on the color-color",
     "ratios...\n"];
mydoc(msg);

src = [ "hplot(c.time_min, c.time_max, c.count_rate, 14 ); %% gray line", 
	"oplot(c.time[lv_1], c.count_rate[lv_1], 2 );  %% red dots", 
	"oplot(c.time_min[lv_2], c.count_rate[lv_2], 4 );  %% blue dots" ] ; 
prsrc( src ) ; 

% plot the full-band light curve, and color selected bins:
%
label("Time", "Counts", "");
connect_points(1); set_line_width(3); 
xrange(0); yrange(0);
hplot(c.time_min, c.time_max, c.count_rate, 14 );
set_line_width(5); connect_points(0); pointstyle(22);
oplot(c.time[lv_1], c.count_rate[lv_1],2 );
pointstyle(23);
oplot(c.time_min[lv_2], c.count_rate[lv_2], 4 );
pointstyle(-1);
plot_pause;


mydoc(["Plot all together........"]);

% plot the color-color points, and color the selections....
%
define plt_1()
{
    connect_points(0); pointstyle(15); set_line_width(5);
    xrange(0); yrange(0);
    label("v1 (6A/3A)", "v2 (9A/6A)", "Vela X-1 phase 0.25");
    plot(v1,v2, 14);
    oplot( v1[lv_1], v2[lv_1], 2);    % overplot selections in color
    oplot( v1[lv_2], v2[lv_2], 4);
    connect_points(1) ; 
}

% plot the spectrum, as count rate, for each selection:
%
define plt_2()
{
    connect_points(1); set_line_width(5);
    label("Wavelength", latex2pg("phot/cm^2/s/bin"), "");
    xrange(1.5, 10.5); yrange(1.e-6,0.002); ylog;
    hplot( x1, x2, f1, 2);
    ohplot(x1, x2, f2, 4);
    ylin;
}

% plot the full-band light curve, and color selected bins:
%
define plt_3()
{
    label("Time", "Counts", "");
    connect_points(1); set_line_width(3); 
    xrange(0); yrange(0);
    hplot(c.time_min, c.time_max, c.count_rate, 14 );
    set_line_width(7); connect_points(0); pointstyle(22);
    oplot(c.time[lv_1], c.count_rate[lv_1], 2 );
    pointstyle(23);
    oplot(c.time_min[lv_2], c.count_rate[lv_2], 4 );
    pointstyle(-1);
}


close_plot(pw1);

pw2 = open_plot("/xwin", 1, 3);
%set_vp( 0.1, 0.98, 0.15, 0.87 ) ; 
resize(16, 1.4);
charsize(1.33);
plt_1;
plt_2;
plt_3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% hardcopy - repeat plotting.

mydoc(["Writing hardcopy to ./outdir/VX1_demo_1.ps...."]);


pid = open_plot("./outdir/VX1_demo_1.ps/vcps", 1, 3);
resize(16,1.4);
charsize(1.33);
set_frame_line_width(3);
%set_vp ( 0.1, 0.95, 0.15, 0.9) ;

plt_1;
plt_2;
plt_3;
_pgiden;

close_plot(pid);
window(pw2);

connect_points(1);  % back to default


% write curves to files:

mydoc( ["Write curves to FITS files:"]);

message( "aglc_write_curve( \"./outdir/vx1_c.fits\",  s.fevt_1a,    c, \"full band\");" ) ;
message( "aglc_write_curve( \"./outdir/vx1_c1.fits\", s.fevt_1a,   c1, \"hard band\");" ) ; 
message( "aglc_write_curve( \"./outdir/vx1_c2.fits\", s.fevt_1a,   c2, \"mid band\");" ) ; 
message( "aglc_write_curve( \"./outdir/vx1_c3.fits\", s.fevt_1a,   c3, \"soft band\");" ) ; 

aglc_write_curve( "./outdir/vx1_c.fits",  s.fevt_1a,  c, "full band");
aglc_write_curve( "./outdir/vx1_ch.fits", s.fevt_1a, c1, "hard band");
aglc_write_curve( "./outdir/vx1_cs.fits", s.fevt_1a, c2, "mid  band");
aglc_write_curve( "./outdir/vx1_pc.fits", s.fevt_1a, c3, "soft band");