Color Landscapes (xfig example)

4 Colorful landscapes together

require("isisscripts") ;
require("rainbow_white") ;

()=evalfile("pfold_test.sl") ;
require("png") ;
require("xfig") ;
xfig_set_tmp_dir("/tmp" );

define normalize (x) {
    return (x-min(x))/max(x-min(x)) ;
    % return x/moment(x).ave ;
    % return (x-moment(x).ave)/moment(x).sdev ;
}


variable lc =  fits_read_table("/home/fuerst/data/XMM/gx301/pn_timing/lc_1sec_.3-12/light/src_sd_bary.lc") ;

variable binaryt = BinaryCor(lc.time/86400. + MJDref_satellite("XMM") ; asini = 368.3, porb = 41.498, eccentricity=0.462, t0 = 48802.79, omega = 310.4 ) * 86400. ; %binary corrected time in sec
variable t0b = min(binaryt) ;

lc = struct_combine(lc , struct {error = sqrt(lc.counts), binaryt = binaryt }) ;

variable per = 685.551 ;
variable nbins = 32 ;

%cut out some interesting parts
variable t1 = 0. ;
variable t2 =   21*per ; % 1.47e4 + 80. ;  % 
variable t3a =  46*per ; % 3.12e4 + 80. ;  % 
variable t3b =  51*per ; % 3.5e4 ;        % 
variable t3c =  62*per ; % 4.18e4 ;	% 
% variable t4 = 46735.29557657242 ;
variable t4 = 46650.0 ; % cut away last few data points as they are exactly zero.

define plotlsc (tb, tf, filename)
{
variable ndx = where(tb <= lc.time-min(lc.time)  <= tf) ;

()=printf("selected %d times\n", length(ndx)) ;
variable lcfilt = struct_filter(lc, ndx ; copy ) ;

variable nbinshist ; 
  
% variable pplc = pfold(lcfilt.binaryt, lcfilt.counts, per , lcfilt.error  ; nbins = nbins, returnlc ) ;

 variable lo, hi ;
variable mincts , maxcts ;

  if (qualifier_exists("normpp")) {
mincts =-11.8 ; maxcts =11.8 ;
}

else if (qualifier_exists("normlc")) {
mincts =0 ; maxcts =1 ;
}

else {
mincts =1 ; maxcts =200 ;
}

nbinshist =  qualifier("nbinshist", ( maxcts - mincts)/3) ;

  variable pplc = pfold(lcfilt.binaryt, lcfilt.counts, per , lcfilt.error  ; nbins = nbins, returnlc  ) ;

variable normfac =1 ;

  if (qualifier_exists("normpp"))
    {
      lcfilt.counts = (1.*lcfilt.counts- min(pplc.value))/(max(pplc.value-min(pplc.value))/2.)-1 ;
      normfac 	= (max(pplc.value-min(pplc.value))/2.) ;
      pplc.value = (1.*pplc.value- min(pplc.value))/(max(pplc.value-min(pplc.value))/2.)-1 ;
    }

  if (qualifier_exists("normlc"))
    {
      normfac 	= max(lcfilt.counts-min(lcfilt.counts)) ;
      lcfilt.counts = (1.*lcfilt.counts- min(lcfilt.counts))/normfac ;
      pplc = pfold(lcfilt.binaryt, lcfilt.counts, per , lcfilt.error  ; nbins = nbins, returnlc  ) ;
      % pplc.value = (1.*pplc.value- min(pplc.value))/normfac ;
    }

  pplc.error/=normfac ;

    
% variable nbinshist =( maxcts - mincts)/3;
% variable nbinshist = 48 ;

(lo,hi)= linear_grid(mincts, maxcts, nbinshist) ;

variable lsarr = Double_Type[nbins, nbinshist] ;

variable i ;

_for i (0, nbins-1, 1)
{
  lsarr[i,*] =  normalize(double(histogram(lcfilt.counts[pplc.revlo[i]],lo,hi)));
  % if (qualifier_exists("normpp") or qualifier_exists("normlc"))  pplc.error[i] = moment(lcfilt.counts[pplc.revlo[i]]).sdev/normfac ;
}

variable pngls = png_gray_to_rgb(lsarr, "rainbow_white") ;

png_write_flipped("pplsc_"+filename+"_bla.png",transpose( pngls)  );

variable xfls = xfig_plot_new(6,6) ;
xfls.world(0,1,mincts, maxcts) ;
  if(qualifier_exists("nozero")) {xfls.xaxis(;major=[0.5,1], minor=[0.1:1:0.1]) ;}
  xfls.xylabel(0.1, 0.9*maxcts, filename) ;
xfls.hplot([pplc.bin_lo[-1]-1, pplc.bin_lo, pplc.bin_lo[0]+1],
	   [pplc.value[-1], pplc.value, pplc.value[0]],
	   [pplc.error[-1], pplc.error, pplc.error[0]] ;; struct_combine(struct{ color="black", width = 4}, __qualifiers) ) ;
xfls.plot([0,1],[1,1] ;color="black", line=1) ;
xfls.plot([0,1],[-1,-1] ;color="black", line=1) ;
  xfls.plot_png("pplsc_"+filename+"_bla.png") ;

xfls.xlabel("Phase") ;
xfls.ylabel("Cts/sec  [normalized]") ;

xfls.render("plots/pplsc_"+filename+".eps") ;

return xfls ;
}


variable nbinsnormlc = 55 ;

variable xfall = plotlsc(t1, t4, "all"  ; normlc, eb_factor=0 , nbinshist=nbinsnormlc     ) ;
variable xfp1 = plotlsc(t1, t2,    "I"  ; normlc, eb_factor=0 , nbinshist=nbinsnormlc   ) ;
variable xfp2 = plotlsc(t2, t3a,   "II"  ; normlc, eb_factor=0, nbinshist=nbinsnormlc     ) ;
variable xfp3a = plotlsc(t3a, t3c, "III"  ; normlc, eb_factor=0, nbinshist=124/2   ) ;
% xfp3b = plotlsc(t3b, t3c, "IV" ) ;
variable xfp3c = plotlsc(t3c, t4,  "IV"  ; normlc, eb_factor=0, nbinshist=nbinsnormlc-15 , nozero) ;

xfig_multiplot( xfp1, xfp2, xfp3a, xfp3c ; cols = 2).render("plots/pplsc_comp_normlc.eps")  ;


%plot the lightcurve (again), superimpsoing the cuts we make

variable nu_lc = rebin_lc(lc, 20.; time = "binaryt", rate = "counts") ;


variable xflc = xfig_plot_new(12,6) ;
xflc.world(0.1, max(nu_lc.binaryt-t0b), 0, 224);

xflc.xlabel("Time [sec]") ;
xflc.ylabel("Counts/sec") ;

xflc.plot(nu_lc.binaryt-t0b, nu_lc.counts) ;

xflc.plot([t2,t2],[0,300] ; color = "orange") ;
xflc.plot([t3a,t3a],[0,300] ; color = "orange") ;
xflc.plot([t3c,t3c],[0,300] ; color = "orange") ;

xflc.xylabel(t2/2., 210, "I") ;
xflc.xylabel((t2+t3a)/2., 210, "II") ;
xflc.xylabel((t3a+t3c)/2., 210, "III") ;
xflc.xylabel((t3c+t4)/2., 210, "IV") ;

xflc.render("plots/xflc.eps") ;