Chandra Coverage (xfig example)

Chandra coverage of Cyg X-1

Here you can see two different things: first how to create a nice circular diagram, and second how to calculate the Roche lobe for the companion star in isis.

We start with the script for the figure:

require("isisscripts");
require("xfig");

% define a couple of colors in hex-RGB:
xfig_new_color("3815n", 0xA00000);
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);

% encode the observing modes for Chandra:
variable TE = 0;
variable CC = 1;

% collect the information about all the Chandra observations:
% { ObsID, orbital phases (min, max), plot-radius, observing mode, ASM-rate, (plot-color) }
variable observations =
[ {  "107", 0.908, 0.940, 0.5, TE, 32.35 },  % 0
  {  "107", 0.723, 0.751, 0.5, TE, 32.35 },  % 1
  { "1511", 0.819, 0.845, 0.5, CC, 21.4  },  % 2
  { "2415", 0.734, 0.796, 0.7, CC, 29.3  },  % 3
  { "3407", 0.847, 0.909, 0.7, CC, 49.07 },  % 4
  { "2741", 0.196, 0.206, 0.6, TE, 69.75 },  % 5
  { "2742", 0.522, 0.532, 0.6, TE, 66.59 },  % 6
  { "2743", 0.704, 0.714, 0.6, TE, 43.09 },  % 7
  { "3724", 0.964, 0.019, 0.5, CC, 61.29 },  % 8
  { "3815", 0.703, 0.818, 0.8, CC, 25.79, "3815n" },  % 9
  { "3814", 0.925, 0.025, 0.8, TE, 17.55, "3814" },  % 10
  { "8525", 0.017, 0.079, 0.85, TE, 13.52, "8525" },  % 11
  { "9847", 0.17 , 0.21 , 0.85, TE, 17.52, "9847" },  % 12
  {"11044", 0.4775,0.54 , 0.85, TE, 18.38, "11044" },  % 13
];

% The radius of the plot
variable R = 5;

% We want to draw a lot of arcs, therefore we ease our work
% by creating a function that does the job for us.
% The minimum required information are center and radius of the
% underlying circle. This arc defaults to a full circle if no start 
% stop values for the angle are given.
define xfig_new_polyline_arc()
{
  variable x, y, r, phi1=0, phi2=2*PI;
  switch(_NARGS)
  { case 3: (x, y, r) = (); }
  { case 5: (x, y, r, phi1, phi2) = (); }

  variable N = qualifier("N", 720);
  variable phi = [phi1 : phi2 : #N];
  return xfig_new_polyline(x+r*cos(phi), y+r*sin(phi), 0*phi;; __qualifiers());
}

% load Roche lobe of star:
% How to create the fits file with the Roche lobe information
% is shown in the separate script below
$1 = "Roche_contours_xy_q0.567.fits+2";
vmessage("V/V[Roche] = %f", fits_read_key($1, "V_VR"));
variable Roche = fits_read_table($1);
struct_filter(Roche, [[::50]]);  % take only every 50th element of the struct's arrays

% show only the observations corresponding to these indices:
variable inds=[9:13];
{
  variable xfig = xfig_new_compound();

  % drawing the coordinate system:
  xfig.insert( xfig_new_polyline_arc(0, 0, R; N=100) );
  variable phi;
  for(phi=0; phi<0.99; phi+=0.1)
  { % tic marks:
    $1 = -[0.98, 1.02]*R;
    variable c = cos(2*PI*phi);
    variable s = sin(2*PI*phi);
    xfig.insert(xfig_new_polyline($1*c, $1*s));
    % ticlabels:
    $1 = -1.05*R;
    xfig.insert( xfig_new_text(sprintf("%.1f", phi); rotate=360*phi+90, x0=$1*c, y0=$1*s) );
  }

  % star
  $1 = .2*R;
  xfig_new_color("lightblue", 0x80C0FF);
  xfig.insert( xfig_new_polyline($1*(Roche.x-1), $1*Roche.y; fillcolor="lightblue", depth=60, closed) );
  xfig.insert( xfig_new_text("\texttt{\small x}"R) );

  % observations
  variable obs;
  variable i = 0;
  foreach obs (observations[inds])
  {
    variable phi1 = obs[1]*2*PI;
    variable phi2 = obs[2]*2*PI;
    if(phi2<phi1)  phi2 += 2*PI;
    variable r = -obs[3]*R;  % for many observationv this is just a more or less random number 
                             % to avoid overlap between distinct observations
    i++;
    variable col = obs[6]; 
% The commented code produces the same colors as in the previous example with the light curve.
% They are color coding the hardness of the observation.
%    $1 = int( log(obs[5]/13.5)/log(70/13.5)*256);
%    xfig_new_color(col, $1 * 0x010000 + 0xFF-$1);
    xfig.insert(xfig_new_polyline_arc(0, 0, r, phi1, phi2; N=50, width=7, line=obs[4], color=col));

    variable phi = [phi1:phi2:#50];
    variable x = r*cos(phi);
    variable y = r*sin(phi);

    % draw dashed lines from the center to the edges of the observation
    xfig.insert( xfig_new_polyline([x[ 0], 0], [y[ 0], 0]; line=1, color=col) );
    xfig.insert( xfig_new_polyline([x[-1], 0], [y[-1], 0]; line=1, color=col) );
    
    % and finally label the observation with its obsID
    phi = (phi1+phi2)/2;
    $1 = r*.9;
    variable col = xfig_lookup_color_rgb(obs[6]);
    xfig.insert( xfig_new_text(sprintf("\definecolor{c}{rgb}{%.3f,%.3f,%.3f}\color{c}%s"R,
                                       (col >> 16)*1./255, ((col & 0x00FF00) >> 8)*1./255, (col & 0x0000FF)*1./255, obs[0]);
                               rotate=180*phi/PI+90, x0=$1*cos(phi), y0=$1*sin(phi)) );
  }

  % we are done and can render the figure:
  xfig.render("Chandra_coverage.pdf");
}

And now: the Roche lobe contours in the x-y-plane. The x-z-plane works analoguely. The whole secret here is to know that our isisscripts provide functions to calculate the Roche potential; namely, Roche_lobe_dimension, Roche_potential, and Roche_equipotential_volume. (get_contour_lines is used for historical reasons. gcontour_compute from the gcontour module shipped with slxfig is much more efficient!) Feel free to read the help for each of them!

require("isisscripts");

variable q = 0.567; % mass ration of the two objects
variable RL = Roche_lobe_dimension(q); 
variable x = [1.01*RL.xmin : 1.01*RL.xmax : #2000];
variable y = [1.01*RL.ymin : 1.01*RL.ymax : #2000];
variable X, Y; (X,Y) = get_grid(x, y; reshape);

variable phi = Roche_potential(X, Y; q=q);
plot_image(phi, 1, x, y);

variable fp = fits_open_file(sprintf("Roche_contours_xy_q%.3f.fits", q), "c");
variable VRoche = NULL;
% fits_write_image_hdu(fp, sprintf("Roche potential(q=%.3f)", q), phi);

variable potential;
foreach potential ({
 [1., RL.potential],
 [.9, -3.177503], % star filling 90% of Roche lobe -> will be 2nd table in fits file
#iffalse
 [.95, -3.139556],
 [.8, -3.264192],
 [.7, -3.368641],
 [.6, -3.497279],
 [.5, -3.660313],
 [.4, -3.876381],
 [.3, -4.183464],
 [.2, -4.677502],
#endif
})
{
  vmessage("V/V[Roche] = %.2f,  phi = %f", potential[0], potential[1]);
  variable cont = get_contour_lines(phi, potential[1]);
  oplot(x[cont[0].x], y[cont[0].y]);
  variable V = Roche_equipotential_volume(q; phi=potential[1]);
  if(VRoche==NULL)  VRoche = V;
  fits_write_binary_table(fp, sprintf("V/V[Roche(q=%.3f)]=%.3f", q, V/VRoche),
        struct { x = x[cont[0].x], y = y[cont[0].y] },
        struct { q = q, phiRoche = potential[1], volume = V, V_VR = V/VRoche, r_rR = (V/VRoche)^(1./3) }
  );
}
fits_close_file(fp);