%; Time-stamp: <2007-03-18 20:42:57 dph> %; Directory: ~dph/d6/Proc_data/LC_demo_44_Boo %; File: Demo.sl %; Author: D. Huenemoerder %; Original version: 2003.12.18 %;==================================================================== % version: 2 % % purpose: analyze 44 boo light curves; demonstrate light/phase curve % prototype functions. % % require("aglc"); _auto_declare = 1; % be lax about variable declarations %% utility function to find the complement of a "where()" list: % private define notwhere(idx,npts) { variable y = Integer_Type[ npts ] ; y[idx] = 1; return where( not y ) ; } % utility function to print an array of strings % (for output of demo documentation); % private define mydoc(str_array) { variable i; ()=printf("\n"); for (i=0;i>> %s\n", str_array[i]); ()=printf("\n"); } 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 ) ; } % open a plot window: putenv("PGPLOT_BACKGROUND=white"); % invert default for better view. putenv("PGPLOT_FOREGROUND=black"); pw1 = open_plot("/xwin"); charsize(1.4); resize(20,0.6); erase; %_pgsvp(0.1, 0.98, 0.15, 0.9); set_vp(0.1, 0.98, 0.15, 0.9); set_line_width(7); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mydoc(["Load pha file, +-1st orders, group, combine, and plot:"]); message("h = load_data(\"indir/pha2\");\n\n"); h = load_data("indir/pha2"); % 4-row file, for inspecting counts. group_data(h[[0,1]],8); % rebin by 8x. group_data(h[[2,3]],4); % rebin by 4x. g_meg = combine_datasets(h[[2,3]]); % dynamic summing of orders. g_heg = combine_datasets(h[[0,1]]); % define some convenient parameters: % dt = 200.0; w1 = 1.7; w2 = 8.0; w3 = 13.0; w4 = 26.0; %%% plot the spectra; counts: % xrange(w1, w4); yrange(0); title("44 Boo HETGS (binned 4,8)"); xlabel("Wavelength"); ylabel("Counts/bin"); message("hplot( get_combined(g_meg, &get_data_counts), 2);"); connect_points(1); hplot( get_combined(g_meg, &get_data_counts), 1); ohplot( get_combined(g_heg, &get_data_counts), 14); y0 = max(get_combined(g_meg, &get_data_counts).value) * 0.4; set_line_width(9); oplot([w1,w2], [y0,y0], 4); oplot([w3,w4], [y0,y0], 2); set_line_width(3); color(4); xylabel(3,y0*1.05,"\"Hard\" band"); color(2); xylabel(17,y0*1.05,"\"Soft\" band"); plot_pause; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mydoc(["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.\n"]); message("s = aglc_read_events( f_evt1a, f_evt0 );\n\n"); % read the grating coord events and the exposure (unfiltered) events: % For this demo, these two files have been filtered to include % only the columns needed here (to make them smaller). % f_evt1a = "indir/evt1a.tg"; f_evt0 = "indir/stat1.exp"; s = aglc_read_events( f_evt1a, f_evt0 ); %%% use both HEG and MEG, and orders -3 to 3 (except 0) for light curves: % g = ["H", "M"]; % grating types o = [-3,-2,-1,1,2,3]; % orders mydoc(["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.\n", "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." ] ); % compute and plot light curves in full band, and soft and hard bands: % multiplot(3,1); charsize(1.2); yrange(0,2.2); xrange(-1000); set_line_width(5); label("Time", "count_rate", "44 Boo, HETGS"); message("c= aglc(s, dt, w1, w4, g, o); % Full band\n"); c= aglc(s, dt, w1, w4, g, o); message("hplot( c.time_min, c.time_max, c.count_rate, 1);"); hplot( c.time_min, c.time_max, c.count_rate, 1); color(1); xylabel(5000, 1.5, "Full band"); message("\nch= aglc(s, dt, w1, w2, g, o); % Hard band\n"); ch=aglc(s, dt, w1, w2, g, o); message("ohplot(ch.time_min, ch.time_max, ch.count_rate, 4);"); ohplot(ch.time_min, ch.time_max, ch.count_rate, 4); color(4); xylabel(1800, 0.4, "Hard band"); message("\ncs= aglc(s, dt, w3, w4, g, o); % Soft band\n"); cs=aglc(s, dt, w3, w4, g, o); message("ohplot(cs.time_min, cs.time_max, cs.count_rate, 2);"); ohplot(cs.time_min, cs.time_max, cs.count_rate, 2); color(2); xylabel(50000, 0.35, "Soft band"); mydoc(["The light curve structure is..."]); message("print(c);\n"); print(c); plot_pause; yrange;xrange(-1000); ylabel("Hardness"); hplot( c.time_min, c.time_max, ch.count_rate/cs.count_rate, 15 ); mydoc(["Form a hardness-ratio by dividing hard counts by soft, and plot:"]); message("hplot( c.time_min, c.time_max, ch.count_rate/cs.count_rate, 15 );\n"); plot_pause; mydoc(["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();\n", "Here, we'll apply pre-computed limits."]); %%%% mark/set regions of probably-not-flaring, based on hardness: % %%% for cursor interaction, uncomment and execute: % % (x1,x2,y1,y2) = cursor_box(); % (x3,x4,y1,y2) = cursor_box(); % (x5,x6,y1,y2) = cursor_box(); %() = printf("x1=%.0f; x2=%.0f; x3=%.0f; x4=%.0f; x5=%.0f; x6=%.0f;\n", % x1, x2, x3, x4, x5, x6); % % saved values from cursor interaction: % x1=2632; x2=12844; x3=17133; x4=38851; x5=45047; x6=60910; l_lo = where( (c.time_min>x1 and c.time_minx3 and c.time_minx5 and c.time_min