Orbit plot of GX 301-2 (xfig example)

Orbit plot of GX 301-2

require("isisscripts") ;
% require("/scratch2/fuerst/isisscripts/share/isisscripts") ;
% ()=evalfile("/scratch2/fuerst/isisscripts/src/plot/xfig_draw_orbit.sl") ;
()=evalfile("xfig_draw_orbit.sl") ;

% orbit parameters of GX 301-2 according to Koh et al (1997).
variable asini = 368.3 ; % lt-sec (semi major axis) 368.3
variable i = 66; % degree to rad 66
variable ecc = 0.46; %
variable omega = 310.4 ; % degrees
variable rstar = 62 ; % solar radii  in lt-sec
variable c=  2.998e8 ; %speed of light in cm/s
variable rsun=1.3914e9/2. ; %solar radius in c

%calculating the parameters of the ellipse
%(necessary to get the ASM lightcurve and its coordinate system)
% orbit parameters 
i = i /360.*2*PI ; % degree to rad
variable a  = asini/sin(i) ;
variable b = sqrt( a^2 * (1-ecc^2) ); % semi-minor axis, lt-sec
omega = omega / 360. * 2 * PI ;

%mean and eccentric anomaly for a full circle / ellipse 
variable M = [0:2*PI:#1024] ; %mean anomaly
variable eccanom = KeplerEquation (M, ecc) ;


%widht and height of th plot
variable W = 12 ;
variable H = W ;

%get the orbital phase of the XMM obs and the two Suzaku obs
variable orbphase1 = orbitalphase (55024.103 , "GX301-2_Koh1997") ;
variable orbphase2 = orbitalphase (55024.603 , "GX301-2_Koh1997") ;
variable orbphase_s1b = orbitalphase (54703.552280 , "GX301-2_Koh1997") ;
variable orbphase_s1f= orbitalphase (54704.0036689 , "GX301-2_Koh1997") ;
variable orbphase_s2b = orbitalphase (54836.439629 , "GX301-2_Koh1997") ;
variable orbphase_s2f = orbitalphase (54838.0419328  , "GX301-2_Koh1997") ;

%load the ASM lightcurve and fold into ont the orbit
variable asmlc = RXTE_ASM_lightcurve ("gx301-2") ;
variable asmpp = pfold(asmlc.time, asmlc.rate, 41.498 ;nbins = 80, t0 = 48802.79) ;

%min and max values of the ASM lightcurve
variable max4cord = 5.0 ;
variable min4cord = 0.0 ;
%normalize ASM lightcurve to these values above
variable pplc = (asmpp.value -min4cord)/(max4cord-min4cord);

%%%calculating the ellipse which is modulated by the ASM lightcurve
variable gxx2, gxy2, gxx2b, gxy2b  ;
variable gxxall = Double_Type[0];
variable gxyall = Double_Type[0];

%how far away from the orbit should the ASM lightcurve appear?
variable minscal = 0.04 ;
variable maxscal = 0.4 ;
variable j ;
_for j (0, length(pplc)-1, 1) {
   (gxx2,gxy2) = ellipse(a+asini*(minscal+maxscal*pplc[j]), b+asini*(minscal+maxscal*pplc[j]),
			 omega, eccanom[where(asmpp.bin_lo[j] < M/2/PI <= asmpp.bin_hi[j])]) ;
   gxxall = [gxxall,gxx2] ;
   gxyall = [gxyall,gxy2] ;
}
gxxall = [gxxall,gxxall[0]];
gxyall = [gxyall,gxyall[0]];


%defining eccentric and true anomaly to be used in the coordinate grid
%for the ASM lightcurve (phase indicator lines)
variable eccanom_cgrid = KeplerEquation([0.0:1:0.1]*2*PI, ecc) ;
variable theta_cgrid = 2.*atan( sqrt( (1.+ecc)/(1.-ecc) )*tan(eccanom_cgrid/2.) );  % true anomaly from eccentric anomaly

%all ellipses are shifted so that the star sits at [0,0].
variable shiftx = a*ecc*cos(omega) ;
variable shifty = a*ecc*sin(omega) ;

%define some colors:
variable darkgray = rgb2hex(0.4,0.4,0.4 ; str) ; % used for the phase labels
variable darkestgray = rgb2hex(0.15,0.15,0.15 ; str) ; % used for the ASM labels
variable asmcol = rgb2hex(0.1,0.55,0.1 ; str) ; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%opening the plot and defining the world coordinate system
variable xfgx = xfig_plot_new(W,H) ;
xfgx.world(-720,500,-500, 720) ;
variable axf, bxf ;
%plotting the orbit and highlighting the three observations
(axf, bxf)  = xfig_draw_orbit(xfgx, asini, i/2/PI*360., ecc, omega/2/PI*360., rstar ; phases={[orbphase1,orbphase2]},
		phase_color=["red", "blue", "blue"], phase_width=4, scale=1.0, trueanom,zerolines=0, width=3, depth=45,
			       forward_arrow, arrow_width=8,arrow_height=12,arrow_type=2) ;

%Earth indicating arrow
xfgx.plot([0,0],[-145,-425] ; line=1, forward_arrow, depth = 10, arrow_thickness =2) ;
xfgx.xylabel(0,-445, "To Earth" ; size ="\scriptsize"R);
%render now without ASM stuff or phase indicators.
% xfgx.render("plots/gx301_skizze2.eps") ;

%plotting the coodinate grid around the ellipse for the ASM lightcurve
variable  xlab, ylab ;
variable labscal = 0.88 ;

% _for j (0,1, 1){
%   %phase indicator lines
%   xfgx.plot([0,cos(theta_cgrid[j]+omega)]*1500,[0,sin(theta_cgrid[j]+omega)]*1500 ; color="gray", line=3) ;
%   %their labels
%   (xlab, ylab) = ellipse(a, b, omega, eccanom_cgrid[j]+0.07)  ;
%   xfgx.xylabel((xlab-shiftx)*labscal, (ylab-shifty)*labscal,sprintf("%.1f", j/10. ) ;
% 	       rotate=(theta_cgrid[j]+omega)/2/PI*360., color=darkgray, size ="\scriptsize"R ) ;
% }
% _for j (2,length(theta_cgrid)-1, 1){
%   %phase indicator lines
%   xfgx.plot([0,cos(theta_cgrid[j]+omega)]*1500,[0,sin(theta_cgrid[j]+omega)]*1500 ; color="gray", line=3, depth=81) ;
%   %their labels
%   (xlab, ylab) = ellipse(a, b, omega, eccanom_cgrid[j]-0.07)  ;
%   xfgx.xylabel((xlab-shiftx)*labscal, (ylab-shifty)*labscal,sprintf("%.1f", j/10. ) ;
% 	       rotate=(theta_cgrid[j]+omega)/2/PI*360.+180, color=darkgray, size ="\scriptsize"R,depth=20 ) ;
% }

%calculating and plotting the grid for the ASM lightcurve
variable xgrid, ygrid,g ;
%how many lines to we want to have?
variable gridpart = 1./4 ;
variable xp1, yp1, xp2, yp2 ;
for (g=0 ; g<=1 ; g+=gridpart)
{
  (xgrid, ygrid) = ellipse(a+(minscal+g*maxscal)*asini, b+(minscal+g*maxscal)*asini, omega, [0:2*PI:#256]) ;
  (xp1, yp1) = ellipse(a+(minscal+g*maxscal)*asini, b+(minscal+g*maxscal)*asini, omega, 0.3*2*PI) ;
  xfgx.plot(xgrid-shiftx, ygrid-shifty; color="gray", line=3,depth=80) ;
  %labeling the ASM coordinate system
  xfgx.xylabel(xp1-shiftx,yp1-shifty,sprintf("%.2f",g*max4cord+min4cord) ;size ="\scriptsize"R, rotate=-34, depth=7, color=darkestgray ) ;
}
%writing out  the units
(xp1, yp1) = ellipse(a+(minscal+g*maxscal)*asini, b+(minscal+g*maxscal)*asini, omega, 0.3*2*PI) ;
xfgx.xylabel(xp1-a*ecc*cos(omega),yp1-a*ecc*sin(omega), "[cts\,sec$^{-1}$]"R ;size ="\scriptsize"R, rotate=-34, depth=7, color=darkestgray ) ;

%plotting the ASM lightcurve
xfgx.plot(gxxall-shiftx,gxyall-shifty; color=asmcol, depth = 10, width=3) ;

%calculating and plotting a bezier curve pointing to the XMM observation
variable bez = bezier({[-600,-350],[-550,-120],[-200,-590],[-142,-250]}, 100) ;
xfgx.plot(bez.x,bez.y ; forward_arrow, color="red", depth = 5) ;
xfgx.xylabel(-600,-375,"XMM" ; size = "\scriptsize"R, color="red") ;
xfgx.xylabel(-600,-400,"47\\,ksec Obs." ; size = "\scriptsize"R, color="red") ;

xfgx.xylabel(0,0,"\\centering Wray\\\\ 977" ; depth=0); 

%calculating and plotting a bezier curve pointing to the Suzaku observation
% variable bezs1 = bezier({[-100,280],[-100,460],[-200,100],[-250,505]}, 100) ;
% xfgx.plot(bezs1.x,bezs1.y ; forward_arrow, color="blue", depth=5) ;
% variable bezs1 = bezier({[-80,280],[-30,460],[-50,200],[105,400]}, 100) ;
% xfgx.plot(bezs1.x,bezs1.y ; forward_arrow, color="blue", depth=5) ;
% xfgx.xylabel(-100,260,"Suzaku" ; size = "\scriptsize"R, color="blue") ;

%rendering it all
xfgx.render("plots/gx301_skizze_nosuz2.png") ;