Lightcurve Cygx1 (xfig example)

From Remeis-Wiki
Jump to navigation Jump to search

Color-coded lightcurves

The script below contains two versions of this plot: one starting in the year 2007 (running the script as it is), one several years earlier (using some of the commented stuff instead). The reason is that the positions of some of the labels/lines in that plot needed to be adjusted when the overall plot size stayed the same but the data shown shrank.

require("isisscripts");

xfig_new_color("3815n", 0xA00000);
xfig_new_color("3815d", 0x600000);
xfig_new_color("3815s", 0xC85050);
xfig_new_color("3814", 0xA0A000);
xfig_new_color("8525", 0x00A000);
xfig_new_color("9847", 0x0000E0);
xfig_new_color("11044", 0x600080); 
xfig_new_color("gray", 0x808080);

%%%%%%%%%%%%%%%%%%%%%%%%%
%RXTE-ASM
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
%variable lc_r = RXTE_ASM_lightcurve("cygx1"; dt=1);
variable lc_r = fits_read_table("/home/miskovicova/WORK/sources/CygX1/RXTE/ASM_lc/RXTE-ASM_CygX-1.fits");
% run a get_data.sl script in that directiory to ge the most
% up-to-date light curve 

variable xxmin = min(lc_r.year);
variable xxmax = max(lc_r.year);
variable xymin = min(lc_r.rate);
variable xymax = max(lc_r.rate);

variable x2xmin = MJDofDate(1996,1,1);
%variable x2xmax = max(lc_r.mjd);
variable startyear = 2007;
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Swift
%%%%%%%%%%%%%%%%%%%%%%%%%%

variable lc_s = fits_read_table("/eu/X-ray/Swift/BAT_lightcurves/CygX-1.lc.fits");

variable yxmin = min(lc_s.time);
%variable yxmax = max(lc_s.time);
variable yymin = min(lc_s.rate);
variable yymax = max(lc_s.rate);

%%%%%%%%%%%%%%%%%%%%%%%%%
% MAXI
%%%%%%%%%%%%%%%%%%%%%%%%%

variable lc_m_time, lc_m_rate, lc_m_error, lc_m_2_4rate, lc_m_2_4err, 
         lc_m_4_10rate, lc_m_4_10err, lc_m_10_20rate, lc_m_10_20err;
(lc_m_time, lc_m_rate, lc_m_error, lc_m_2_4rate, lc_m_2_4err, lc_m_4_10rate, lc_m_4_10err, 
 lc_m_10_20rate, lc_m_10_20err) = readcol("/eu/X-ray/MAXI/lc_1day/CygX-1.txt", 1,2,3,4,5,6,7,8,9);

variable maxi = struct{ mjd = lc_m_time, rate = lc_m_rate, err = lc_m_error };

%variable zxmin = min(lc_m_time);
%variable zxmax = max(lc_m_time);
variable zymin = min(lc_m_rate);
variable zymax = max(lc_m_rate);

variable max_all = max(max(lc_r.mjd), max(lc_s.time), max(lc_m_time));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

variable q,qq, maj = Double_Type[0];
_for q (1996,2012,1)
_for qq (1,12,1)
{
   maj = [maj, MJDofDate(q,qq,1)];
}

variable mino =[MJDofDate(1996,1,1),
		MJDofDate(1997,1,1),
		MJDofDate(1998,1,1),
		MJDofDate(1999,1,1),
		MJDofDate(2000,1,1),
		MJDofDate(2001,1,1),
		MJDofDate(2002,1,1),
		MJDofDate(2003,1,1),
		MJDofDate(2004,1,1),
		MJDofDate(2005,1,1),
		MJDofDate(2006,1,1),
		MJDofDate(2007,1,1),
		MJDofDate(2008,1,1),
		MJDofDate(2009,1,1),
		MJDofDate(2010,1,1),
		MJDofDate(2011,1,1),
		MJDofDate(2012,1,1)];
variable aa = 0.25, bb = 0.15;

%variable x = xfig_plot_new(30,4.5);
variable x = xfig_plot_new(20,3.5);
%x.world(xxmin, xxmax, xymin-1, xymax+1);
%x.world1(x2xmin-50, x2xmax+200, 5, xymax+1);
%x.world2(x2xmin-50, x2xmax+200, 5, (zymin-1.5), (zymax+0.1));
x.world1(MJDofDate(startyear,1,1), max_all, xymin-1, xymax+10);
x.world2(MJDofDate(startyear,1,1), max_all, (zymin-2.5), (zymax+0.1));
%x.plot(lc_m_time, (lc_m_rate); world2, color="black", size=.5, depth=50);

%%% MAXI color array:
variable grayscale = [0:0.75:#256];

variable i, cl;
 foreach i ( [0:length(lc_m_time)-2] )
 {
    if(lc_m_time[i+1]-lc_m_time[i] >5 ) continue;
    cl = log((lc_m_rate[i]+lc_m_rate[i+1])/2/1)/log(4.5/1);
    if(cl<0)     { cl = 0;}
    if(cl>0.75)  { cl = 0.75;}
    cl = grayscale[ where(grayscale[[:-2]] <= cl <= grayscale[[1:]]) ][0];
    x.plot(lc_m_time[i+[0,1]], lc_m_rate[i+[0,1]]; world2, color=rgb2hex(cl, cl, cl;str));
 }                                                                                            %   
%message("maxi done");
%x.plot(lc_r.mjd, lc_r.rate; world1, sym="point", color="blue");

variable xtescale1 = [0:1:#100];
foreach i ( [0:length(lc_r.mjd)-2] )
{
   if(lc_r.mjd[i+1]-lc_r.mjd[i] >5 ) continue;
   $1 = log((lc_r.rate[i]+lc_r.rate[i+1])/2/13.5)/log(70/13.5);
   $2 = 0;
   if($1<0) { $2=-$1;  $1 = 0;}
   if($1>1) { $2=$1-1; $1 = 1;}
   $1 = xtescale1[ where(xtescale1[[:-2]] <= $1 <= xtescale1[[1:]] ) ][0];
   x.plot(lc_r.mjd[i+[0,1]], lc_r.rate[i+[0,1]]; world1, color=rgb2hex($1, $2, 1-$1;str));   
}
%message("rxte done");
x.plot([MJDofDate(2003,3,4), MJDofDate(2003,3,4)],   [xymin-1, xymax+10]; color="3815n");
x.plot([MJDofDate(2003,4,19), MJDofDate(2003,4,19)], [xymin-1, xymax+10]; color="3814");
x.plot([MJDofDate(2008,4,18), MJDofDate(2008,4,18)], [xymin-1, xymax+10]; color="8525");
x.plot([MJDofDate(2008,4,19)+3, MJDofDate(2008,4,19)+3], [xymin-1, xymax+10]; color="9847");
x.plot([MJDofDate(2010,1,14), MJDofDate(2010,1,14)], [xymin-1, xymax+10]; color="11044");
x.plot([MJDofDate(2011,1,6), MJDofDate(2011,1,6)], [xymin-1, xymax+10]; color="gray");
x.plot([MJDofDate(2012,2,6), MJDofDate(2012,2,6)], [xymin-1, xymax+10]; color="gray");
%variable obsid1 = xfig_new_text("ObsID\\,3815"; size="footnotesize", color="3815n", rotate=0, x0=1.45, y0=-0.2, depth=10);
%variable obsid2 = xfig_new_text("ObsID\\,3814"; size="footnotesize", color="3814", rotate=0, x0=1.75, y0=-0.5, depth=10);
%variable obsid3 = xfig_new_text("ObsID\\,8525"; size="footnotesize", color="8525", rotate=0, x0=10.5, y0=-0.2, depth=10);
%variable obsid4 = xfig_new_text("\\&\\,9847"; size="footnotesize", color="9847", rotate=0, x0=12.3, y0=-0.2, depth=10);
%variable obsid5 = xfig_new_text("ObsID\\,11044"; size="footnotesize", color="11044", rotate=0, x0=15, y0=-0.2, depth=10);
%variable obsid6 = xfig_new_text("ObsID\\,12472"; size="footnotesize", color="gray", rotate=0, x0=17.5, y0=-0.2, depth=10);
%variable obsid7 = xfig_new_text("AO13"; size="footnotesize", color="gray", rotate=0, x0=19.5, y0=-0.2, depth=10);
variable obsid3 = xfig_new_text("ObsID\\,8525"; size="footnotesize", color="8525", rotate=0, x0=3.7, y0=-0.2, depth=10);
variable obsid4 = xfig_new_text("\\&\\,9847"; size="footnotesize", color="9847", rotate=0, x0=5.5, y0=-0.2, depth=10);
variable obsid5 = xfig_new_text("11044"; size="footnotesize", color="11044", rotate=0, x0=11.7, y0=-0.2, depth=10);
variable obsid6 = xfig_new_text("12472"; size="footnotesize", color="gray", rotate=0, x0=15.45, y0=-0.2, depth=10);
variable obsid7 = xfig_new_text("AO13"; size="footnotesize", color="gray", rotate=0, x0=19.5, y0=-0.2, depth=10);
%x.add_object(obsid1);
%x.add_object(obsid2);
x.add_object(obsid3);
x.add_object(obsid4);
x.add_object(obsid5);
x.add_object(obsid6);
x.add_object(obsid7);
x.y2axis(; major=[int(zymin):int(zymax):2], major_len=0.25, minor=[int(zymin):int(zymax):0.5], minor_len=0.15);
x.y2label("\textsl{MAXI} 2--20\,keV"R);
%x.xylabel((x2xmin-380), 70, "\textsl{RXTE}-ASM"R; rotate=90);
%x.xylabel((x2xmin-280), 70, "1,5--12\,keV"R; rotate=90);
%x.xylabel(MJDofDate(startyear,1,1)-255, 70, "\textsl{RXTE}-ASM"R; rotate=90);
%x.xylabel(MJDofDate(startyear,1,1)-190, 70, "1,5--12\,keV"R; rotate=90);
x.xylabel(MJDofDate(startyear,1,1)-145, 70, "\textsl{RXTE}-ASM"R; rotate=90);
x.xylabel(MJDofDate(startyear,1,1)-105, 70, "1,5--12\,keV"R; rotate=90);
x.x1axis(; ticlabels=0, color="white");
%x.x2axis(; ticlabels=0);
x.x2axis(; major=maj, minor=mino, ticlabels=0,
	          %ticlabels=["1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012"],
	          major_len=0.15, minor_len=0.3);


%variable y = xfig_plot_new(30,3);
variable y = xfig_plot_new(20,2);
%y.world1(2003, xxmax, (yymin-0.01), (yymax+0.01));
y.world(MJDofDate(startyear,1,1), max_all, -0.08, (yymax+0.01) );%(yymin-0.01), (yymax+0.01));
%y.world2(dateOfMJD(52640), dateOfMJD(55996.5), (yymin-0.01), (yymax+0.01));
%y.plot(lc_s.time, lc_s.rate; world2, color="red");
variable j,k, swiftcl = Integer_Type[length(lc_s.time)-1];
foreach i ( [0:length(lc_s.time)-2] )
{
   if(lc_s.time[i+1]-lc_s.time[i] >5) continue;
   j = where(int(lc_r.mjd[[:-2]])+1 == lc_s.time[i]);
   if (length(j) == 0)
     {
	k = where_min ( abs( (lc_s.rate[[0:i-2]]+lc_s.rate[[0:i-2]+1] )/2. - (lc_s.rate[i]+lc_s.rate[i+1])/2.  ) );
	j = swiftcl[k];
     }
   else { j = j[0]; }
 $1 = log((lc_r.rate[j]+lc_r.rate[j+1])/2/13.5)/log(70/13.5);
 $2 = 0;
 if($1<0) {$2=-$1;  $1 = 0;}
 if($1>1) {$2=$1-1; $1 = 1;}
 swiftcl[i] = j;
   $1 = xtescale1[ where(xtescale1[[:-2]] <= $1 <= xtescale1[[1:]] ) ][0];
 y.plot(lc_s.time[i+[0,1]], lc_s.rate[i+[0,1]]; color=rgb2hex($1, $2, 1-$1; str));
}
%message("swift done");
y.plot([MJDofDate(2003,3,4), MJDofDate(2003,3,4)],   [-0.1, yymax+0.01]; color="3815n");
y.plot([MJDofDate(2003,4,19), MJDofDate(2003,4,19)], [-0.1, yymax+0.01]; color="3814");
y.plot([MJDofDate(2008,4,18), MJDofDate(2008,4,18)], [-0.1, yymax+0.01]; color="8525");
y.plot([MJDofDate(2008,4,19)+3, MJDofDate(2008,4,19)+3], [-0.1, yymax+0.01]; color="9847");
y.plot([MJDofDate(2010,1,14), MJDofDate(2010,1,14)], [-0.1, yymax+0.01]; color="11044");
y.plot([MJDofDate(2011,1,6), MJDofDate(2011,1,6)], [xymin-1, xymax+10]; color="gray");
y.plot([MJDofDate(2012,2,6), MJDofDate(2012,2,6)], [xymin-1, xymax+10]; color="gray"); 
%y.xylabel(MJDofDate(startyear,1,1)-255,0.2, "\textsl{Swift}-BAT"R; rotate=90);
%y.xylabel(MJDofDate(startyear,1,1)-190,0.2, "15--50\,keV"R; rotate=90);
y.xylabel(MJDofDate(startyear,1,1)-145,0.2, "\textsl{Swift}-BAT"R; rotate=90);
y.xylabel(MJDofDate(startyear,1,1)-105,0.2, "15--50\,keV"R; rotate=90);
y.x1axis(; world2);
variable cc = "\color{white}{a}"R;
%variable ticcs = [cc,cc,"1996",cc,cc,cc,"1997",cc,cc,cc,"1998",cc,cc,cc,"1999",cc,cc,cc,
%		  "2000",cc,cc,cc,"2001",cc,cc,cc,"2002",cc,cc,cc,"2003",cc,cc,cc,
%		  "2004",cc,cc,cc,"2005",cc,cc,cc,"2006",cc,cc,cc,"2007",cc,cc,cc,
%		  "2008",cc,cc,cc,"2009",cc,cc,cc,"2010",cc,cc,cc,"2011",cc,cc,cc,"2012",cc];
variable ticcs = String_Type[length(maj)];
ticcs[*] = cc;
_for q (0,16,1) ticcs[q*12+6] = string(1996+q);
%variable ticcss = ["1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006",
%		  "2007","2008","2009","2010","2011","2012"];
y.x1axis(; major=maj, minor=mino,
	          ticlabels=ticcs,
	          major_len=0.15, minor_len=0.3);
y.x2axis(; ticlabels=0, color="white");
y.translate(vector(0,-2.5,0));

x.add_object(y);
%message("kidiing???");
x.render("lc_cygx1.eps");