Welcome to ISIS Version 1.4.4-3
Copyright (C) 1998-2007 Massachusetts Institute of Technology

isis> .load Demo
% aglc version 1.3.2 loaded.
% >>>  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!)

s = aglc_read_events("indir/evt1a.tg", "indir/stat1.exp"); 
% Reading files indir/evt1a.tg, indir/stat1.exp ...
% Read 1190671 events.
% Read 106644 exposure reference records.
% Applying gti filters... 4 5 6 7 8 9
% >>>  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.


c = aglc( s, 71, 1.5, 10.5, ["H","M"], [-1,1]);
% >>>   curve = aglc( evt_struct, time_bin, min_waves, max_waves, gratings, orders );
% >>>   min_waves and max_waves may be arrays or scalars.

% Frames per tbin = 42
% time_step = 7.312e+01
% ontime per bin = 7.140e+01
% Computing rates for ccd_id = 5 6 7 8 9 
% min time = 699.550000 sec
% max time = 30330.550000 sec
% >>>  Here are the light curves for 1.5-10.5A (black)
% >>>  hplot(  c.time_min,  c.time_max,  c.count_rate, 1 );
c1 = aglc( s, 71, 1.5, 4.5, ["h","m"], [-1,1]);
% Frames per tbin = 42
% time_step = 7.312e+01
% ontime per bin = 7.140e+01
% Computing rates for ccd_id = 6 7 8 
% min time = 699.550000 sec
% max time = 30330.550000 sec
% >>>   ....... 1.5-4.5A (red)
c2 = aglc( s, 71, 4.5, 7.5, ["h","m"], [-1,1]);
% Frames per tbin = 42
% time_step = 7.312e+01
% ontime per bin = 7.140e+01
% Computing rates for ccd_id = 5 6 7 8 
% min time = 699.550000 sec
% max time = 30330.550000 sec
% >>>  ........ 4.5-7.5A (green)
c3 = aglc( s, 71, 7.5, 10.5, ["h","m"], [-1,1]);
% Frames per tbin = 42
% time_step = 7.312e+01
% ontime per bin = 7.140e+01
% Computing rates for ccd_id = 5 6 7 8 9 
% min time = 699.550000 sec
% max time = 30330.550000 sec
% >>>  ........ and 7.5-10.5A (blue).
Press RETURN to continue

% >>>  Ratios of the curves in bands give us a "color-color" plot:

 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);
Press RETURN to continue

% >>>  We can select via an expression, 
% >>>  or use the gtk-module GUI, "vwhere()"
% >>>  Done automatically here via...


lv_1 = where( v1 > 0.6 );
lv_2 = where( v1 <= 0.4 );
Press RETURN to continue

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


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

% >>>  Loading, interpolating,  and summing ARFs, binning spectra...
% >>>  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!
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
Press RETURN to continue
% >>>  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...

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
Press RETURN to continue

% >>>  Plot all together........
% >>>  Writing hardcopy to ./outdir/VX1_demo_1.ps....
% >>>  Write curves to FITS files:
aglc_write_curve( "./outdir/vx1_c.fits",  s.fevt_1a,    c, "full band");
aglc_write_curve( "./outdir/vx1_c1.fits", s.fevt_1a,   c1, "hard band");
aglc_write_curve( "./outdir/vx1_c2.fits", s.fevt_1a,   c2, "mid band");
aglc_write_curve( "./outdir/vx1_c3.fits", s.fevt_1a,   c3, "soft band");
isis> exit;