3D orbit plots (xfig example)
Jump to navigation
Jump to search
The projection and handling of plotting in 3D is very similar to the examples Plotting in 3D (simple) and Plotting in 3D.
%-*- mode: fold; mode: slang -*-
require("isisscripts");
% we will produce two plots: one showing the orbital plane in the
% context of the tangent plane of the sky and, thus, defining the
% orbital elements. the second one defines the Be disk geometry
% with respect to the orbital plane only, i.e, *without* the
% tangent plane of the sky and the orbital elements.
% in the coordinate system used here, the tangent plane of the sky
% is defined as the xy-plane and the z-axis points away from Earth.
% that is with a camera inclination and roll angle of zero, we are
% looking directly onto the tangent plane of the sky. note that the
% x and y axis points north and east, respectively, i.e., the xy
% axis are flipped to ensure a right handed coordinate system.
xfig_add_latex_package("bm"); % for bold and italic fonts in mathmode
xfig_set_latex_preamble("\let\mathbfit\bm"R); % alias for MNRAS
xfig_set_font_style(NULL); % no bold face anywhere
%%% implement camera object %{{{
% set up the camera view port
private define _xfig_camera_setup(s, incl, zrot, roll) {
s.camera_data.incl = incl;
s.camera_data.zrot = zrot;
s.camera_data.roll = roll;
% update camera vectors
s.camera_data.cx = vector_rotate(
vector(1,0,0), vector(0,1,0), -s.camera_data.zrot*PI/180
);
s.camera_data.cy = vector_rotate(
vector_rotate(vector(0,1,0), vector(1,0,0), s.camera_data.incl*PI/180),
vector(0,1,0), -s.camera_data.zrot*PI/180
);
% optical axis of the camera
variable camaxis = crossprod(s.camera_data.cx, s.camera_data.cy);
% rotate around this vector
s.camera_data.cx = vector_rotate(s.camera_data.cx, camaxis, s.camera_data.roll*PI/180);
s.camera_data.cy = vector_rotate(s.camera_data.cy, camaxis, s.camera_data.roll*PI/180);
}
% add polygons to the plot
private define _xfig_camera_plot_polygon(s, v) {
% add starting point [0,0,0]
if (length(v.x) == 1 && not qualifier_exists("skiporigin")) {
v = struct_copy(v);
v.x = [0, v.x];
v.y = [0, v.y];
v.z = [0, v.z];
}
% add to list
list_append(s.camera_data.polygons, v);
list_append(s.camera_data.polygons_qualifiers, __qualifiers);
}
% add shaded region to the plot
private define _xfig_camera_shade_region(s, v) {
% add to the same list as for polygons, but set a qualifier
list_append(s.camera_data.polygons, v);
list_append(s.camera_data.polygons_qualifiers, struct_combine(
__qualifiers, struct { _camera_shade_region }
));
}
% add label to the plot
private define _xfig_camera_xyzlabel() {
variable s, v, l, dx = 0, dy = 0;
switch (_NARGS)
{ case 3: (s, v, l) = (); }
{ case 5: (s, v, l, dx, dy) = (); }
list_append(s.camera_data.labels, v);
list_append(s.camera_data.labels, l);
list_append(s.camera_data.labels, [dx, dy]);
list_append(s.camera_data.labels_qualifiers, __qualifiers);
}
% render polygons
private define _xfig_camera_render_polygons(s, pl) {
% loop over polygons
variable i;
_for i (0, length(s.camera_data.polygons)-1, 1) {
% shade region
if (s.camera_data.polygons_qualifiers[i] != NULL &&
struct_field_exists(s.camera_data.polygons_qualifiers[i], "_camera_shade_region")) {
pl.shade_region(
dotprod(s.camera_data.polygons[i], s.camera_data.cx),
dotprod(s.camera_data.polygons[i], s.camera_data.cy);;
s.camera_data.polygons_qualifiers[i]
);
}
% plot line
else {
pl.plot(
dotprod(s.camera_data.polygons[i], s.camera_data.cx),
dotprod(s.camera_data.polygons[i], s.camera_data.cy);;
s.camera_data.polygons_qualifiers[i]
);
}
}
% loop over labels
_for i (0, length(s.camera_data.labels)/3-1, 1) {
pl.xylabel(
dotprod(s.camera_data.labels[3*i], s.camera_data.cx),
dotprod(s.camera_data.labels[3*i], s.camera_data.cy),
s.camera_data.labels[3*i+1], __push_array(s.camera_data.labels[3*i+2]);;
s.camera_data.labels_qualifiers[i]
);
}
}
% new render function
private define _xfig_camera_render(s, file) {
% first render all polygons
s.camera_render_polygons(s, s);
% default xfig render
(@(s.camera_data.camera_xfig_render))(file);
}
% return an xfig_plot_new with a camera object included
define xfig_new_camera_plot() {
variable args = __pop_args(_NARGS);
variable s = struct_combine(
struct {
camera_setup = &_xfig_camera_setup,
camera_plot_polygon = &_xfig_camera_plot_polygon,
camera_xyzlabel = &_xfig_camera_xyzlabel,
camera_shade_region = &_xfig_camera_shade_region,
camera_render_polygons = &_xfig_camera_render_polygons,
camera_data = struct { % camera data
polygons = list_new, % list of polygons (Vector_Type[]) to plot
polygons_qualifiers = list_new, % plot qualifiers for each polygon
labels = list_new, % list of labels (Vector_Type[], String_Type)
labels_qualifiers = list_new, % qualifiers for each label
incl, zrot, roll, % camera inclination (around x), rotation (around z), and roll (around the optical axis)
cx, cy, % unit vectors
camera_xfig_render % remember xfig renderer function
}
},
xfig_plot_new(__push_args(args))
);
% default camera setup
s.camera_setup(0, 0, 0);
% remember xfig renderer function and then overwrite
s.camera_data.camera_xfig_render = s.render;
s.render = &_xfig_camera_render;
return s;
}
%}}}
%%% useful functions %{{{
% function to split a closed line given as an array of vectors using indices.
% ensure that the returned arrays of vectors define a line without breaks
define vecsplitsort(r, n) {
variable len = vector_norm(vector(diff(r.x[n]), diff(r.y[n]), diff(r.z[n])));
variable i = @n;
if (max(len) > 2*median(len)) {
i = shift(n, where_max(len)+1);
}
return vector(r.x[i], r.y[i], r.z[i]);
}
% function to shift, incline, and rotate a vector into the orbital plane
variable a, incl, omega, nOrb, ecc;
define vec2orbit(x) {
% translate according to the focal point (center of mass)
ifnot (qualifier_exists("notrans")) {
x += vector(-sqrt(a^2 - (a*(1-ecc))^2), 0, 0);
}
% incline
ifnot (qualifier_exists("noincl")) {
x = vector_rotate(x, vector(1, 0, 0), incl * PI/180);
}
% rotate
ifnot (qualifier_exists("noomega")) {
x = vector_rotate(
x, nOrb, omega * PI/180
);
}
return x;
}
% calculates an arch between the vectors v1 and v2 with radius r
define arch(v1, v2, r) {
v1 *= (1./vector_norm(v1));
v2 *= (1./vector_norm(v2));
% angle between the vectors and rotation axis, i.e, the normal
variable angle = acos(dotprod(v1, v2));
variable axis = crossprod(v1, v2);
axis *= (1./vector_norm(axis));
if (qualifier_exists("reverse")) {
axis *= -1;
angle = 2*PI - angle;
}
% calculate arch by rotating v1 around the axis
variable nphi = int(1000*angle/2/PI);
variable step = angle / nphi;
variable arch = Vector_Type[nphi];
variable i;
_for i (0, nphi-1, 1) {
arch[i] = r*vector_rotate(v1, axis, step*i);
}
return merge_struct_arrays(arch);
}
%}}}
%%% define geometry %{{{
% orbit parameters
variable ecc = .05, a = 300; % = 1 - b/a
variable omega = -50, incl = 30.;
% Be disk parameters
variable inclDisk = -50;
variable omegaDisk = 180+28;
variable radDisk = 150;
% neutron star's orbital plase
variable phiNS = .9;
% define normal vectors of the reference frame (plotted later)
% (note that x and y are flipped, see abow)
variable ex = vector(0,1,0);
variable ey = vector(1,0,0);
variable ez = vector(0,0,1);
%}}}
variable n, m;
%%% calculation of orbit, disk, and neutron star position and disk %{{{
% calculate orbit
variable nOrb = vector_rotate(vector(0, 0, 1), vector(1, 0, 0), incl * PI/180);
variable phi = [0:2*PI:#1000];
variable rOrb = vector(
Double_Type[length(phi)], Double_Type[length(phi)], Double_Type[length(phi)]
);
(rOrb.x, rOrb.y) = ellipse(a, a*(1-ecc), 0, phi);
rOrb = vec2orbit(rOrb);
% calculate Be disk
variable rDisk = vector(
Double_Type[length(phi)], Double_Type[length(phi)], Double_Type[length(phi)]
);
(rDisk.x, rDisk.y) = ellipse(radDisk, radDisk, 0, phi);
rDisk = vector_rotate(rDisk, vector(1, 0, 0), (incl+inclDisk) * PI/180);
rDisk = vector_rotate(rDisk, nOrb, omegaDisk * PI/180);
% disk normal
variable nDisk = vector_rotate(vector(0, 0, 1), vector(1, 0, 0), (incl+inclDisk) * PI/180);
nDisk = vector_rotate(nDisk, nOrb, omegaDisk * PI/180);
% neutron star position
variable n = wherefirstmin(abs(phi-phiNS*2*PI)); % orbital phase
variable rNS = vector(rOrb.x[n], rOrb.y[n], rOrb.z[n]);
%}}}
variable H = 12; % plot heights (width set individually)
%%% first xfig plot %{{{
variable pl1 = xfig_new_camera_plot(8.,H);
pl1.camera_setup(40, 40, 0);
variable xmin = -315, ymin = -350, wlen = 850;
variable WH = 1.*pl1.plot_data.plot_height/pl1.plot_data.plot_width;
pl1.world(xmin, xmin+wlen/WH, ymin, ymin+wlen);
pl1.axes(; major = 0, minor = 0, color = "white");
pl1.xylabel(.01, .99, "a)", -.5, .5; size = "small", world0);
% define colors
variable orbColor = "#ffca57";
variable tgColor = "#ffff71";
% xyz-axes
variable xyzlen = 250.; % arrow length
pl1.camera_plot_polygon(xyzlen*ex; forward_arrow, color = "red", depth = 20);
pl1.camera_plot_polygon(xyzlen*ey; forward_arrow, color = "red", depth = 20);
pl1.camera_plot_polygon(xyzlen*ez; forward_arrow, color = "red", line = 2);
pl1.camera_xyzlabel(1.05*xyzlen*ex, "x"; color = "red");
pl1.camera_xyzlabel(1.05*xyzlen*ey, "y"; color = "red");
pl1.camera_xyzlabel(1.05*xyzlen*ez, "z"; color = "red");
% plot orbit
n = where(rOrb.y < 0, &m);
pl1.camera_plot_polygon(vecsplitsort(rOrb, n); depth = 85);
pl1.camera_plot_polygon(vecsplitsort(rOrb, m); line = 2, depth = 85);
pl1.camera_shade_region(vecsplitsort(rOrb, n); color = orbColor, depth = 85);
pl1.camera_shade_region(vecsplitsort(rOrb, m); color = orbColor, depth = 95);
n = struct_filter(rOrb, int(length(rOrb.x)*.38); copy);
pl1.camera_xyzlabel(n, "orbital plane"R, .3, .2; size = "small", rotate = -25);
% plot tangent plane of the sky
variable tg = vector(
-390*[1,0,0,1,1] + 320*[0,1,1,0,0],
-250*[1,1,0,0,1] + 370*[0,0,1,1,0],
[0,0,0,0,0]
);
pl1.camera_plot_polygon(tg; depth = 90);
pl1.camera_shade_region(tg; color = tgColor, depth = 90);
n = -28;
pl1.camera_xyzlabel(struct_filter(tg, 1; copy), "tangent plane"R, .55, -.7; size = "small", rotate = n);
pl1.camera_xyzlabel(struct_filter(tg, 1; copy), "of the sky"R, .6, -.2; size = "small", rotate = n);
% intersection of the planes
n = wherefirstmin(abs(rOrb.y + rOrb.z));
pl1.camera_plot_polygon(vector(rOrb.x[n], 0, 0); line = 1);
% inclination angle
% find intersection of the orbit with the x-axis
n = where_min(abs(abs(rOrb.y) + rOrb.z)) + [-1,0,1];
m = interpol_points(0, rOrb.y[n], rOrb.x[n]);
pl1.camera_plot_polygon(vector([m,m], [0,-80], [0,0]); depth = 50);
pl1.camera_xyzlabel(vector(m, -60, 0), "$i$", .75, -.1);
% find the point on the orbit, which at a certain distance from this intersection point
n = where(rOrb.y < 0 and rOrb.x > 0);
n = n[where_min(abs(vector_norm(struct_filter(rOrb, n; copy) - vector(m, 0, 0)) - 80))];
% connect the intersection point and the orbit with an arch of the same length as the distance
pl1.camera_plot_polygon(vector(m, 0, 0) + arch(
vector(0, -1, 0), vector(rOrb.x[n] - m, rOrb.y[n], rOrb.z[n]), 80
));
% periastron and apastron position
variable pstrn = vec2orbit(vector(a, 0, 0));
variable astrn = vec2orbit(vector(-a, 0, 0));
pl1.camera_plot_polygon(merge_struct_arrays([pstrn, astrn]); line = 1, depth = 50);
pl1.camera_xyzlabel(pstrn, "P", -.5, .5);
pl1.camera_xyzlabel(astrn, "A", .3, -.8);
% prism supporting the 3D position of the periastron
pl1.camera_plot_polygon(merge_struct_arrays([pstrn, pstrn-vector(0,pstrn.y,0)]));
pl1.camera_plot_polygon(merge_struct_arrays([vector(pstrn.x,0,pstrn.z), vector(pstrn.x,0,0)]));
pl1.camera_plot_polygon(vector(pstrn.x, 0, pstrn.z));
% same for the apastron
pl1.camera_plot_polygon(merge_struct_arrays([astrn, astrn-vector(0,astrn.y,0)]); line = 2);
pl1.camera_plot_polygon(merge_struct_arrays([vector(astrn.x,0,astrn.z), vector(astrn.x,0,0)]); line = 2);
pl1.camera_plot_polygon(vector(astrn.x, 0, astrn.z); line = 2);
% omega angle
n = arch(vector(1, 0, 0), pstrn, 85);
m = length(n.x)/2; % vector to half of the arch (-> label)
pl1.camera_plot_polygon(n; line = 1);
pl1.camera_xyzlabel(.8*vector(n.x[m], n.y[m], n.z[m]), "$\omega$"R, 0, -.4);
% neutrons star position
pl1.camera_plot_polygon(rNS; forward_arrow, line = 1);
pl1.camera_xyzlabel(.5*rNS, "$\mathbfit{r}_\mathrm{ns}$"R, -.6, .5);
n = wherefirstmin(vector_norm(rOrb - rNS));
pl1.camera_plot_polygon(struct_filter(rOrb, n + int(.05*length(rOrb.x)) + [0,1]; copy); forward_arrow);
% neutron star position projected into the Be disk reference frame
pl1.camera_plot_polygon(rNS; sym = "circle", fill, size = .5, skiporigin);
variable rNSdiskZ = nDisk * dotprod(rNS, nDisk);
pl1.camera_plot_polygon(rNS-rNSdiskZ; forward_arrow, line = 3);
pl1.camera_plot_polygon(merge_struct_arrays([rNS-rNSdiskZ, rNS]); forward_arrow);
pl1.camera_xyzlabel(rNS-.5*rNSdiskZ, "$h\mathbfit{n}_\mathrm{disk}$"R, .6, -.2; rotate = -55);
pl1.camera_xyzlabel(.5*(rNS-rNSdiskZ), "$r\mathbfit{e}_\mathrm{r}$"R, .5, -.8);
% line of sight
pl1.camera_plot_polygon(rNS - vector([0,0], [0,0], [0,230]); color = "blue", forward_arrow);
pl1.camera_xyzlabel(rNS, "line of sight"R, .85, .6; color = "blue", size = "\small", rotate = 37);
pl1.camera_xyzlabel(rNS, "to Earth"R, .85, 1.2; color = "blue", size = "\small", rotate = 37);
%}}}
%%% second xfig plot %{{{
variable pl2 = xfig_new_camera_plot(10,H);
variable incCam = 20; % camera inclination with respect to the orbital plane
pl2.camera_setup(incl-90+incCam, 0, 0); % incl-90 -> edge on
% rotate the camera around the orbital normal
variable rotCam = 160;
pl2.camera_data.cx = vector_rotate(pl2.camera_data.cx, nOrb, rotCam*PI/180);
pl2.camera_data.cy = vector_rotate(pl2.camera_data.cy, nOrb, rotCam*PI/180);
variable xmin = -260, ymin = -270, wlen = 500;
variable WH = 1.*pl2.plot_data.plot_height/pl2.plot_data.plot_width;
pl2.world(xmin, xmin+wlen/WH, ymin, ymin+wlen);
pl2.axes(; major = 0, minor = 0, color = "white");
pl2.xylabel(.05, .87, "b)", -.5, .5; size = "small", world0);
xyzlen = 200;
pl2.camera_plot_polygon(xyzlen*ez; forward_arrow, color = "red", line = 2);
pl2.camera_xyzlabel(1.05*xyzlen*ez, "z"; color = "red");
% define colors
variable diskColor = "skyblue";
% z-axis projected onto the orbital plane
variable zdisk = xyzlen*(ez - (ez*nOrb)*nOrb);
%zdisk *= (1./vector_norm(zdisk));
pl2.camera_plot_polygon(zdisk; color = "red", line = 1);
pl2.camera_plot_polygon(merge_struct_arrays([xyzlen*ez, zdisk]); line = 2, color = "red");
% plot orbit
pl2.camera_plot_polygon(rOrb; depth = 75);
pl2.camera_shade_region(rOrb; color = orbColor, depth = 100);
n = where(.837 < phi/2/PI < .95); % nasty hack for the dotted line behin the Be disk
pl2.camera_plot_polygon(struct_filter(rOrb, n; copy); depth = 65, line = 2);
n = struct_filter(rOrb, int(length(rOrb.x)*.21); copy);
pl2.camera_xyzlabel(n, "orbital plane"R, 0, -.3; size = "small", rotate = -19);
pl2.camera_xyzlabel(pstrn, "P", .5, -.5);
% plot Be disk
n = where(dotprod(rDisk, nOrb) < 0, &m);
pl2.camera_plot_polygon(vecsplitsort(rDisk, n); depth = 70);
pl2.camera_plot_polygon(vecsplitsort(rDisk, m); line = 2, depth = 85);
pl2.camera_shade_region(vecsplitsort(rDisk, n); color = diskColor, depth = 70);
pl2.camera_shade_region(vecsplitsort(rDisk, m); color = diskColor, depth = 105);
n = struct_filter(rDisk, int(length(rDisk.x)*.15); copy);
pl2.camera_xyzlabel(n, "Be disk"R, 0, .5; size = "small", rotate = -23);
% highest point of the disk above the orbit
variable hline = vector_rotate(vector(0, 1, 0), vector(1, 0, 0), (incl+inclDisk)*PI/180);
hline = radDisk*vector_rotate(hline, nOrb, omegaDisk*PI/180);
pl2.camera_plot_polygon(hline; line = 3);
pl2.camera_xyzlabel(hline*1.1, "H");
pl2.camera_plot_polygon(hline - (hline*nOrb)*nOrb; line = 1);
pl2.camera_plot_polygon(merge_struct_arrays([hline, hline - (hline*nOrb)*nOrb]); line = 2);
% omega disk
n = arch(zdisk, hline - (hline*nOrb)*nOrb, 40; reverse);
pl2.camera_plot_polygon(n; line = 1);
pl2.camera_xyzlabel(.5*struct_filter(n, length(n.x)/2-1; copy), "$\theta$"R, 0, .1);
% disk inclination
n = arch(hline, hline - (hline*nOrb)*nOrb, 60);
pl2.camera_plot_polygon(n; line = 2);
pl2.camera_xyzlabel(.7*struct_filter(n, length(n.x)/2-1; copy), "$\delta$"R, .2, -.2);
% neutrons star position (new one compared to plot 1)
phiNS = .09;
n = wherefirstmin(abs(phi-phiNS*2*PI)); % orbital phase
rNS = vector(rOrb.x[n], rOrb.y[n], rOrb.z[n]);
pl2.camera_plot_polygon(rNS; sym = "circle", fill, size = .5, skiporigin);
pl2.camera_plot_polygon(rNS; forward_arrow, line = 1);
pl2.camera_xyzlabel(.6*rNS, "$\mathbfit{r}_\mathrm{ns}$"R, 1.1, .9);
n = wherefirstmin(vector_norm(rOrb - rNS));
pl2.camera_plot_polygon(struct_filter(rOrb, n + int(.05*length(rOrb.x)) + [0,1]; copy); forward_arrow);
% neutron star position projected into the Be disk reference frame
rNSdiskZ = nDisk * dotprod(rNS, nDisk);
pl2.camera_plot_polygon(rNS-rNSdiskZ; forward_arrow, line = 3);
pl2.camera_plot_polygon(merge_struct_arrays([rNS-rNSdiskZ, rNS]); forward_arrow);
pl2.camera_xyzlabel(rNS-.5*rNSdiskZ, "$h\mathbfit{n}_\mathrm{disk}$"R, 0, -.3; rotate = 39);
pl2.camera_xyzlabel(.9*(rNS-rNSdiskZ), "$r\mathbfit{e}_\mathrm{r}$"R, 0, -.7);
%}}}
% combine and render
% (note: the camera plot renderes the lines only when .render is called.
% thus we need to call these first, which is not a convenient way...)
pl1.render("../plots/reference_frames.pdf");
pl2.render("../plots/reference_frames.pdf");
xfig_new_hbox_compound(pl1,pl2).render("../plots/reference_frames.pdf");