Insets and shaded regions (xfig example)

Insets and shaded regions

require("isisscripts");

%%% first perform an epoch folding

% load best available time resolution, time is in seconds, used to do an epoch folding
variable bat = fits_read_table("GROJ1008-57.orbit.lc.fits");
% convert time to MJD
bat.time = bat.time/86400 + MJDref_satellite("Swift");
% limit to time range of flaring activity before giant outburst
struct_filter(bat, where(56169 <= bat.time <= 56235));

% perform an epoch folding search
variable epf = epfold(bat.time, bat.rate, 1, 22; nbins = 12, nsrch = 200, exact, dt = bat.timedel/86400);


%%% load 1d-binned lightcurve for the plot (already in MJD)

bat = fits_read_table("GROJ1008-57.lc.fits");
% limit to time range from 2011 to early of 2013
struct_filter(bat, where(55889 <= bat.time <= 56295));


%%% XFIG-plot: lightcurve

% some orbital parameters of the system
variable tau = 54424.6; % time of periastron passage (MJD)
variable porb = 249.46; % orbital period (d)
% define user tic marcs
variable tics = [0:0.4:0.04], tims = [0:0.3:0.02];
variable tics2 = [5:8:.2], tims2 = [5:8:.2]+.1;

% initialize the plot
variable pl = xfig_plot_new(14,8);
pl.world(bat.time[0]+1e-6, bat.time[-1], 1e-6, 0.27); % 1e-6 to "cut" the first tic label
pl.world2((bat.time[0]-tau)/porb, (bat.time[-1]-tau)/porb, 1e-6, 0.27); % x2-axis in orbital phase
pl.xlabel("Time (MJD)");
pl.x2label("Orbital Phase");
pl.ylabel("BAT Count Rate (cts\,s$^{-1}$\,cm$^{-2}$)"R);
pl.y2axis(; ticlabels = 0);
pl.x2axis(; major = tics2, minor = tims2,
  ticlabels = array_map(String_Type, &sprintf, "%.1f", tics2 - int(tics2))); % subtract integer number from orbital phase
pl.y1axis(; major = tics, minor = tims, format = "%.2f");
pl.y2axis(; major = tics, minor = tims, ticlabels = 0);
% plot the lightcurve
pl.plot(bat.time, bat.rate, bat.error; eb_color = "gray"); % gray error-bars

% shade specific regions
variable x0, y0, dxl, dxr, n; % label positions and left and right range, respectively
% 2011 type I
x0 = 55914; y0 = 0.09; dxl = 10; dxr = 14; 
n = where(x0-dxl <= bat.time <= x0+dxr);
pl.shade_region([bat.time[n[0]], bat.time[n], bat.time[n[-1]]], [0, bat.rate[n], 0]; fillcolor = "#b3e7b3");
pl.plot(x0, y0; sym = "darr", color = "black"); % arrow
pl.xylabel(x0, y0, "Type I", 0, -1.0; color = "black"); % label
% 2012 type I
x0 = 56164; y0 = 0.09; dxl = 12; dxr = 5;
n = where(x0-dxl <= bat.time <= x0+dxr);
pl.shade_region([bat.time[n[0]], bat.time[n], bat.time[n[-1]]], [0, bat.rate[n], 0]; fillcolor = "#b3e7b3");
pl.plot(x0, y0; sym = "darr", color = "black"); % arrow
pl.xylabel(x0, y0, "Type I", 0, -1.0; color = "black"); % label
% 2012 type II
x0 = 56221; y0 = 0.16; dxl = -14; dxr = 74;
n = where(x0-dxl <= bat.time <= x0+dxr);
pl.shade_region([bat.time[n[0]], bat.time[n], bat.time[n[-1]]], [0, bat.rate[n], 0]; fillcolor = "#ffa0a0");
pl.plot(x0, y0; sym = "rarr", color = "black"); % arrow (here pointing to the right)
pl.xylabel(x0, y0, "Type II", 0.7, 0; color = "black"); % label
% flaring activity
n = where(56169 <= bat.time <= 56235);
pl.shade_region([bat.time[n[0]], bat.time[n], bat.time[n[-1]]], [0, bat.rate[n], 0]; fillcolor = "#a0a0ff", width = 0);
pl.plot(bat.time[n[0]] * [1,1], [0,0.06]; color = "black", depth = 10, line = 1); % dashed line
pl.plot(bat.time[n[-1]] * [1,1], [0,0.06]; color = "black", depth = 10, line = 1); % dashed line
pl.plot(bat.time[n[0]]+10, 0.061; sym = "larr", color = "black"); % left arrow
pl.plot(bat.time[n[-1]]-10, 0.061; sym = "rarr", color = "black"); % right arrow
pl.xylabel(.5*(bat.time[n[0]]+bat.time[n[-1]]), 0.061, "Flaring"; color = "black");

% mark observation times
variable swiobs = [55915.83, 55917.05, 55918.78, 55919.91, 55932.08, 55934.49, 55938.77, 55940.42,
                   55955.80, 55958.14, 56196.99, 56244.91, 56246.18, 56248.25, 56250.18, 56252.59];
_for n (0, length(swiobs)-1, 1) { pl.plot(swiobs[n], 0.26; sym = "darr", color = "red"); }
variable suzobs = [56251.63];
_for n (0, length(suzobs)-1, 1) { pl.plot(suzobs[n], 0.24; sym = "darr", color = "green3"); }


%%% XFIG-plot: inset epoch folding

variable ple = xfig_plot_new(5,4);
ple.world(epf.p[0], epf.p[-1]-1e-6, 1e-6, max(epf.stat));
ple.xlabel("Test Period (d)"; size = "\footnotesize"R);
ple.ylabel("$\chi^2$"R; size = "\footnotesize"R);
ple.x1axis(; ticlabel_size = "\footnotesize"R);
ple.y1axis(; ticlabel_size = "\footnotesize"R);
ple.title("Flaring Epoch Folding"; size = "\footnotesize"R);
ple.plot(epf.p, epf.stat);


%%% add inset plot and render

pl.add_object(ple, 56035, 0.13);
pl.render("1008batlcgiant2012.pdf");