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. % >>> Load pha file, +-1st orders, group, combine, and plot: h = load_data("./indir/pha2"); Reading: .... hplot( get_combined(g_meg, &get_data_counts), 2); Press RETURN to continue % >>> Read events into a structure used by light curve functions % >>> evt1a holds grating coordinates, and may be filtered. % >>> evt0 is unfiltered, and is used for exposure numbers from % >>> ALL exposed frames. s = aglc_read_events( f_evt1a, f_evt0 ); % Reading files indir/evt1a.tg, indir/stat1.exp ... % Read 379932 events. % Read 110907 exposure reference records. % Applying gti filters... 4 5 6 7 8 9 % >>> Compute light curves in full band, and soft and hard bands % >>> We could name the event files directly instead of the structure, % >>> but for reapeated use, the structure is much more efficient. % >>> NOTE: aglc returns a structure with grids, counts, rates, exposure. % >>> Variables have been defined to generalize the inputs: % >>> dt = time bin; w1-4 are wavelengths % >>> (scalars were used here; can use arrays for multiple intervals); % >>> g = ["H","M"] for HEG and MEG photons; % >>> o = [-3,-2,-1,1,2,3] for orders. c= aglc(s, dt, w1, w4, g, o); % Full band % Frames per tbin = 63 % time_step = 2.042e+02 % ontime per bin = 2.016e+02 % Computing rates for ccd_id = 4 5 6 7 8 9 % min time = 11.200000 sec % max time = 59281.600000 sec hplot( c.time_min, c.time_max, c.count_rate, 1); ch= aglc(s, dt, w1, w2, g, o); % Hard band % Frames per tbin = 63 % time_step = 2.042e+02 % ontime per bin = 2.016e+02 % Computing rates for ccd_id = 4 5 6 7 8 9 % min time = 11.200000 sec % max time = 59281.600000 sec ohplot(ch.time_min, ch.time_max, ch.count_rate, 4); cs= aglc(s, dt, w3, w4, g, o); % Soft band % Frames per tbin = 63 % time_step = 2.042e+02 % ontime per bin = 2.016e+02 % Computing rates for ccd_id = 4 5 6 8 9 % min time = 11.200000 sec % max time = 59281.600000 sec ohplot(cs.time_min, cs.time_max, cs.count_rate, 2); % >>> The light curve structure is... print(c); time = Double_Type[294] time_min = Double_Type[294] time_max = Double_Type[294] counts = Integer_Type[294] count_rate = Double_Type[294] stat_err = Double_Type[294] exposure = Double_Type[294] timezero = 7.30246e+07 mjdref = 50814 revidx1a = Array_Type[294] revidx = Array_Type[294] Press RETURN to continue % >>> Form a hardness-ratio by dividing hard counts by soft, and plot: hplot( c.time_min, c.time_max, ch.count_rate/cs.count_rate, 15 ); Press RETURN to continue % >>> There are several ways to select regions in the curve. % >>> We will select mostly non-flaring times, filter the event structure, % >>> and make phased light curves from the filtered events. % >>> A GUI-selection can be done with the slgtk module and vwhere(), like this: % >>> l_lo = vwhere(c.time_min, ch.count_rate/cs.count_rate); % >>> isis/pgplot cursor selection can be done like this: % >>> (x1,x2,y1,y2) = cursor_box(); % >>> Here, we'll apply pre-computed limits. oplot( ((c.time_min+c.time_max)[l_lo])/2, (ch.count_rate/cs.count_rate)[l_lo],2 ); oplot( ((c.time_min+c.time_max)[l_hi])/2, (ch.count_rate/cs.count_rate)[l_hi],4 ); Press RETURN to continue % >>> To make phased light curves, we have to set an ephemeris: ephem = [ 2443604.5914, 0.26781662 ]; % epoch in JD; period in days % >>> To make phased light curves for the hard and soft states, % >>> we need to filter the event structures. aglc_filter() applies the % >>> selection using information in the light curve, c, (specifically, % >>> reverse indices from the histogram function), and copies the events % >>> counted in selected bins to a new structure: s_lo = aglc_filter( c, s, l_lo ); % filter low-state % >>> Phased light curves are made with the function, aglc_phased(), which is % >>> similar to aglc(). Instead of taking a time bin, it takes % >>> phase grid parameters: pgrid = [min_phase, max_phase, phase_step]. p_lo = aglc_phased( s_lo, [0.0, 1.0, 0.009], ephem, 1.7, 26, g, o); % low state % Frames per phase bin = 63 % Phase bin adjusted to 0.0087 % Computing rates for ccd_id = 4 5 6 7 8 9 % >>> The result is a structure similar to the light curve structure. % >>> Here are the fields: print(p_lo); phase = Double_Type[114] phase_min = Double_Type[114] phase_max = Double_Type[114] counts = Integer_Type[114] count_rate = Double_Type[114] stat_err = Double_Type[114] exposure = Double_Type[114] mjdref = 50814 revidx1a = Array_Type[114] revidx = Array_Type[114] epoch = 2.4436e+06 period = 0.267817 Press RETURN to continue % >>> The exposure tag is computed from the histogram of exposure numbers % >>> phase folded. Phase gaps (such as for the flare exclusions) show % >>> as low or no exposure (bottom plot), and are properly % >>> factored into the rates. hplot(p_lo.phase_min, p_lo.phase_max, p_lo.count_rate, 4); hplot(p_lo.phase_min, p_lo.phase_max, p_lo.exposure/1000., 4); % scaled exp Press RETURN to continue % >>> Another useful quantity is the unfolded phase, or rotational cycle. % >>> Time conversion functions help here. E.g.: jdmin = aglc_tobs_to_jd(c.time_min+c.timezero, c.mjdref); % time to jd rmin = aglc_jd_to_rotnum( jdmin, jd0, pd ); % jd to rotation number rref = int(min(rmin)); % compute zero-point rmin = rmin - rref; % remove zero-point hplot( rmin, rmax, c.count_rate, 1 ); % plot the result Press RETURN to continue % >>> Another view is to plot rotation mod 1 to show the phased curve, % >>> but without accumulation of repeated phases into phase bins: plot( (rmin+rmax)[l_lo]/2 mod 1, c.count_rate[l_lo], 2 ); Press RETURN to continue % >>> Make a big summary plot of all the above, and also put into % >>> a postscript file, ./outdir/LC_Demo_44_Boo.ps Press RETURN to continue % >>> Write curves to FITS files: aglc_write_curve( "outdir/44b_c.fits", s.fevt_1a, c, "full band"); aglc_write_curve( "outdir/44b_ch.fits", s.fevt_1a, ch, "hard band"); aglc_write_curve( "outdir/44b_cs.fits", s.fevt_1a, cs, "soft band"); aglc_write_curve( "outdir/44b_pc.fits", s.fevt_1a, p_lo, "lo state"); isis> exit; >