Monsterlightcurve (xfig example)
Jump to navigation
Jump to search
File common_pars.sl, which I use to load some parameters common to all plots in a paper:
%% -- set up variables for plotting variable width = 10; variable height = 10; variable fullwidth = 25; variable mjdmin = MJDofDate(1995,12,15); variable mjdmax = MJDofDate(2012,12,30); variable yearmin = dateOfMJD(mjdmin); variable yearmax = dateOfMJD(mjdmax); variable max_asmh = 2.5; variable min_asmh = 0.; variable max_asmr = 150.; variable min_asmr = 0.; %%ends of activity periods variable e1 = 50350; variable e2 = 51000; variable e3 = 53900; variable e4 = 55375; %%define colours variable cl_hard = "blue"; variable cl_tran = "green4"; variable cl_soft = "red"; variable cl_neutral = "gray"; variable cl_hard_bg = "blue1"; variable cl_soft_bg = "pink3"; %%define data specific stuff variable gamma_hard = 2.0; variable gamma_soft = 2.5; %define mapping parameters variable asm_th = 20.; % threshold variable asm_x0 = 0.28; %x-intersection variable asm_hs = 55.; %hard slope variable asm_ss = 350.; %soft slope %%define paths to data variable maxiname = "lc_1orbit-Cyg_X-1_55058.09375-56244.21875.csv"
Actual plotting routine:
require("isisscripts");
require("../common_pars");
variable width=25;
variable hh = 3.5;
variable asmmin = 0;
variable asmmax = 200;
variable asmamin = 0;
variable asmamax = 100;
variable asmcmin = 0;
variable asmcmax = 50;
variable batmin = -0.09;
variable batmax = 0.39;
variable asmarrow = 170;
variable ptext = 185;
variable dnd = fits_read_table("../data/all_disk_nodisk.fits");
variable obs = fits_read_table("../data/all_obs.fits");
variable tt = (obs.tstart+obs.tstop)/2.;
variable pca_hard = where(dnd.bknpower_1_phoindx1_value < gamma_hard);
variable pca_tran = where(gamma_hard <= dnd.bknpower_1_phoindx1_value < gamma_soft);
variable pca_soft = where(gamma_soft <= dnd.bknpower_1_phoindx1_value);
variable asm = RXTE_ASM_lightcurve("cygx1");
struct_filter(asm,where(asm.rate_a > 0));
struct_filter(asm,where(asm.rate_b > 0));
struct_filter(asm,where(asm.rate_c > 0));
variable asm_late = struct_filter(asm,where(asm.time > 55200);copy,dim=0);
struct_filter(asm,where(asm.time < 55200));
variable hardness = asm.rate_c/asm.rate_a;
variable asm_hard = where(asm.rate <= asm_th or asm.rate <= (hardness-asm_x0)*asm_hs);
variable asm_tran = where(asm.rate > asm_th and asm.rate > (hardness-asm_x0)*asm_hs
and asm.rate <= (hardness-asm_x0)*asm_ss);
variable asm_soft = where(asm.rate > asm_th and asm.rate > (hardness-asm_x0)*asm_ss);
%%index for the times during the longhard state
variable ind_longhard = where( 53900 < asm.time < 54900);
%BAT
variable bat = fits_read_table("../data/CygX-1.orbit.lc.fits"); %
bat.time = bat.time/3600./24.+51910+0.00074287038;
struct_filter(bat,where(bat.time < 56240));
variable bat_soft = where(bat.rate <= 0.09);
variable bat_uncat = where(bat.rate > 0.09);
%%MAXI
variable maxi = ascii_read_table("../data/"+maxiname,
[{"%F","time"},{"%F","rate"},{"%F","err"},
{"%F","lo"},{"%F","loerr"},
{"%F","mi"},{"%F","mierr"},
{"%F","hi"},{"%F","hierr"}]);
struct_filter(maxi,where(maxi.time < 56240));
variable maxi_hard = where(maxi.lo <= 1.4*maxi.mi/maxi.lo);
variable maxi_tran = where(1.4*maxi.mi/maxi.lo < maxi.lo <= 8./3.*maxi.mi/maxi.lo);
variable maxi_soft = where(maxi.lo > 8./3.*maxi.mi/maxi.lo);
%%GBM
()=evalfile("../gbm/init_gbm.sl");
variable gbm = @d;
%struct_filter(gbm, where(gbm.nr_occ>10));
struct_filter(gbm,where(gbm.mjd_start < 56240));
variable gbm_soft = where(1e-3*get_struct_field (gbm, cnfg.flux[1]) < 0.6);
variable gbm_uncat = where(1e-3*get_struct_field (gbm, cnfg.flux[1]) > 0.6);
%%axis labels
variable mino = [1996:2013];
variable maj = [1996:2013]+0.5;
variable extramin = Double_Type[0];
variable q,i;
variable ticc=String_Type[length(maj)];
_for i (0,2013-1996-1,1){
ticc[i]=string(mino[i]);
_for q (1,11,1){
extramin = [extramin,1996+i+q/12.];
}
}
variable plg = xfig_plot_new(width,hh*5./7.);
plg.world1(mjdmin,mjdmax,-1,4);
plg.world2(yearmin.year+(yearmin.month-1)/12.+yearmin.day/365.,
yearmax.year+(yearmax.month-1)/12.+yearmax.day/365.,
-1,4);
plg.x1axis(;off);
plg.x2axis(;major = maj, minor = mino, ticlabels = ticc,
major_len=0,minor_len=0.2,minor_width=2
);
plg.yaxis(;ticlabels=0,major=0,minor=0);
_for i (0,length(pca_hard)-1,1){
plg.plot([tt[pca_hard[i]],tt[pca_hard[i]]],[2,3];color=cl_hard);
};
_for i (0,length(pca_tran)-1,1){
plg.plot([tt[pca_tran[i]],tt[pca_tran[i]]],[1,2];color=cl_tran);
};
_for i (0,length(pca_soft)-1,1){
plg.plot([tt[pca_soft[i]],tt[pca_soft[i]]],[0,1];color=cl_soft);
};
plg.ylabel("\begin{center}times of\\pointed \\\textsl{RXTE} obs.\end{center}"R;rotate=-90);
plg.x2label("year");
variable pl_asm = xfig_plot_new(width,hh);
pl_asm.world1(mjdmin,mjdmax,asmmin,asmmax);
pl_asm.world2(yearmin.year+(yearmin.month-1)/12.+yearmin.day/365.,
yearmax.year+(yearmax.month-1)/12.+yearmax.day/365.,
asmmin,asmmax);
pl_asm.x2axis(;off);
pl_asm.y2axis(;ticlabels=0);
pl_asm.plot(asm.time[asm_tran],asm.rate[asm_tran];sym="point",size=0.1,color=cl_tran);
pl_asm.plot(asm.time[asm_soft],asm.rate[asm_soft];sym="point",size=0.1,color=cl_soft);
pl_asm.plot(asm.time[asm_hard],asm.rate[asm_hard];sym="point",size=0.1,color=cl_hard);
pl_asm.plot(asm_late.time,asm_late.rate;sym="point",size=0.1,color="gray");
pl_asm.plot([e1,e1],[asmmin,asmmax];line=1,width=2,depth=10);
pl_asm.plot([e2,e2],[asmmin,asmmax];line=1,width=2,depth=10);
pl_asm.plot([e3,e3],[asmmin,asmmax];line=1,width=2,depth=10);
pl_asm.plot([e4,e4],[asmmin,asmmax];line=1,width=2,depth=10);
pl_asm.plot([mjdmin,e1],[asmarrow,asmarrow];
line=0, forward_arrow, arrow_thickness =1,color="gray");
pl_asm.plot([e1,e2],[asmarrow,asmarrow];
line=0, forward_arrow, backward_arrow, arrow_thickness =1,color="gray");
pl_asm.plot([e2,e3],[asmarrow,asmarrow];
line=0, forward_arrow, backward_arrow, arrow_thickness =1,color="gray");
pl_asm.plot([e3,e4],[asmarrow,asmarrow];
line=0, forward_arrow, backward_arrow, arrow_thickness =1,color="gray");
pl_asm.plot([e4,mjdmax],[asmarrow,asmarrow];
line=0, backward_arrow, arrow_thickness =1,color="gray");
pl_asm.xylabel((mjdmin+e1)/2.,ptext,"I";depth=10);
pl_asm.xylabel((e1+e2)/2.,ptext,"II";depth=10);
pl_asm.xylabel((e2+e3)/2.,ptext,"period III";depth=10);
pl_asm.xylabel((e3+e4)/2.,ptext,"period IV";depth=10);
pl_asm.xylabel((e4+mjdmax)/2.,ptext,"period V";depth=10);
pl_asm.ylabel("\begin{center}\textsl{RXTE}-ASM \\1.5--12\,keV \\$[$cps$]$\end{center}"R);
variable pl_asmh = xfig_plot_new(width,hh);
pl_asmh.world1(mjdmin,mjdmax,0.08,15;ylog);
pl_asmh.world2(yearmin.year+(yearmin.month-1)/12.+yearmin.day/365.,
yearmax.year+(yearmax.month-1)/12.+yearmax.day/365.,
0.08,10;ylog);
pl_asmh.x2axis(;major = maj, minor = mino, ticlabels = ticc,
major_len=0,minor_len=0.2,minor_width=2
);
pl_asmh.y2axis(;ticlabels=0);
%pl_asmh.y1axis(;format="$10^{%.1f}$"R);
pl_asmh.plot(asm.time[asm_tran],asm.rate_c[asm_tran]/asm.rate_a[asm_tran];
sym="point",size=0.1,color=cl_tran);
pl_asmh.plot(asm.time[asm_soft],asm.rate_c[asm_soft]/asm.rate_a[asm_soft];
sym="point",size=0.1,color=cl_soft);
pl_asmh.plot(asm.time[asm_hard],asm.rate_c[asm_hard]/asm.rate_a[asm_hard];
sym="point",size=0.1,color=cl_hard);
pl_asmh.plot(asm_late.time,asm_late.rate_c/asm_late.rate_a;
sym="point",size=0.1,color="gray");
pl_asmh.plot([e1,e1],[asmmin,asmmax];line=1,width=2,depth=10);
pl_asmh.plot([e2,e2],[asmmin,asmmax];line=1,width=2,depth=10);
pl_asmh.plot([e3,e3],[asmmin,asmmax];line=1,width=2,depth=10);
pl_asmh.plot([e4,e4],[asmmin,asmmax];line=1,width=2,depth=10);
pl_asmh.ylabel("\begin{center}\textsl{RXTE}-ASM \\ hardness \\ (C$/$A)\end{center}"R);
variable pl_bat = xfig_plot_new(width,hh);
pl_bat.world(mjdmin,mjdmax,batmin,batmax);
pl_bat.world2(yearmin.year+(yearmin.month-1)/12.+yearmin.day/365.,
yearmax.year+(yearmax.month-1)/12.+yearmax.day/365.,batmin,batmax);
pl_bat.y2axis(;ticlabels=0,tic_depth=3);
pl_bat.x2axis(;major = maj, minor = mino, ticlabels = ticc,
major_len=0,minor_len=0.2,minor_width=2
);
pl_bat.plot(bat.time[bat_soft],bat.rate[bat_soft];sym="point",size=0.1,color=cl_soft);
pl_bat.plot(bat.time[bat_uncat],bat.rate[bat_uncat];sym="point",size=0.1);
pl_bat.plot([e1,e1],[batmin,batmax];line=1,width=2,depth=10);
pl_bat.plot([e2,e2],[batmin,batmax];line=1,width=2,depth=10);
pl_bat.plot([e3,e3],[batmin,batmax];line=1,width=2,depth=10);
pl_bat.plot([e4,e4],[batmin,batmax];line=1,width=2,depth=10);
pl_bat.plot([mjdmin,mjdmax],[0,0];color="gray",line=1);
pl_bat.y1axis(;format="%.1f",tic_depth=3);
pl_bat.ylabel("\begin{center}\textsl{Swift}-BAT \\ 15--50\,keV \\$[$cps/cm$^2]$\end{center}"R);
variable pl_maxi = xfig_plot_new(width,hh);
pl_maxi.world(mjdmin,mjdmax,0,5.5);
pl_maxi.world2(yearmin.year+(yearmin.month-1)/12.+yearmin.day/365.,
yearmax.year+(yearmax.month-1)/12.+yearmax.day/365.,0,6);
pl_maxi.x2axis(;major = maj, minor = mino, ticlabels = ticc,
major_len=0,minor_len=0.2,minor_width=2
);
pl_maxi.y2axis(;ticlabels=0);
pl_maxi.plot(maxi.time[maxi_tran],maxi.rate[maxi_tran];sym="point",size=0.1,color=cl_tran);
pl_maxi.plot(maxi.time[maxi_soft],maxi.rate[maxi_soft];sym="point",size=0.1,color=cl_soft);
pl_maxi.plot(maxi.time[maxi_hard],maxi.rate[maxi_hard];sym="point",size=0.1,color=cl_hard);
pl_maxi.plot([e1,e1],[0,6];line=1,width=2,depth=10);
pl_maxi.plot([e2,e2],[0,6];line=1,width=2,depth=10);
pl_maxi.plot([e3,e3],[0,6];line=1,width=2,depth=10);
pl_maxi.plot([e4,e4],[0,6];line=1,width=2,depth=10);
pl_maxi.ylabel("\begin{center}MAXI \\ 2--20\,keV \\$[$cps/cm$^2]$\end{center}"R);
variable pl_gbm = xfig_plot_new(width,hh);
pl_gbm.world(mjdmin,mjdmax,-0.4,1.9);
pl_gbm.world2(yearmin.year+(yearmin.month-1)/12.+yearmin.day/365.,
yearmax.year+(yearmax.month-1)/12.+yearmax.day/365.,-0.4,1.9);
pl_gbm.x2axis(;major = maj, minor = mino, ticlabels = ticc,
major_len=0,minor_len=0.2,minor_width=2
);
pl_gbm.y1axis(;format="%.1f",tic_depth=3);
pl_gbm.y2axis(;ticlabels=0,tic_depth=3);
pl_gbm.plot(0.5*(gbm.mjd_start + gbm.mjd_stop)[gbm_soft],
1e-3*get_struct_field (gbm, cnfg.flux[1])[gbm_soft],
0.5* (gbm.mjd_stop - gbm.mjd_start)[gbm_soft],
1e-3*get_struct_field (gbm, cnfg.err[1])[gbm_soft];
sym = "point",color=cl_soft);
pl_gbm.plot(0.5*(gbm.mjd_start + gbm.mjd_stop)[gbm_uncat],
1e-3*get_struct_field (gbm, cnfg.flux[1])[gbm_uncat],
0.5* (gbm.mjd_stop - gbm.mjd_start)[gbm_uncat],
1e-3*get_struct_field (gbm, cnfg.err[1])[gbm_uncat];
sym = "point");
pl_gbm.plot([mjdmin,mjdmax],[0,0];color="gray",line=1);
pl_gbm.plot([e1,e1],[-0.4,1.9];line=1,width=2,depth=10);
pl_gbm.plot([e2,e2],[-0.4,1.9];line=1,width=2,depth=10);
pl_gbm.plot([e3,e3],[-0.4,1.9];line=1,width=2,depth=10);
pl_gbm.plot([e4,e4],[-0.4,1.9];line=1,width=2,depth=10);
pl_gbm.ylabel("\begin{center}GBM \\ 25--50\,keV \\ $[$Crab$]$\end{center}"R);
variable comp = xfig_new_vbox_compound(plg,
xfig_multiplot(pl_asm,pl_asmh,
pl_maxi,
pl_bat,pl_gbm;
xlabel="MJD"));
comp.render("monster.pdf");
print("STACK:");
_print_stack;
exit;