3D orbit plots (xfig example)

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");