Chandra Coverage (xfig example)
Jump to navigation
Jump to search
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);