User Tools

Site Tools


isis:slxfig:monsterlightcurve

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;
isis/slxfig/monsterlightcurve.txt · Last modified: 2013-07-24 11:44 by Victoria Grinberg