Binary orbits incompatible with radial velocities

[t00fri]

Chris,

Fridger,

Once again, your tone is insulting... If you want to have a reasonable discussion with me, you'll need to restrain yourself from filling your posts with poorly targeted shots at my work.

On the other hand, given Celestia's untransparent coordinate documentation in 2005, such a sign error can perhaps not be excluded. One should also take into consideration that my checks took place before lots of coordinate-related changes have been introduced in the code by Chris, including frames and "all that". I have never seen a clearcut crosscheck that guarantees the identity of "before and after".

Sorry, if this sounded insulting to you. It was not intended to be so -- as usual!

Note that I carefully referred above to the pre 2006 time, when the documentation of the Celestia coordinates was NOT laid down in a central place, as far as I could make out. The year 2005 was crucial in the context of our radial velocity discussion, since it was in 2005 that I tested and committed my binary star data base...

With the advent and considerable progress of the Celestia Wiki book documentation things have enormously improved. This I have never contested, however.

While I don't fancy your frames implementation very much (as we discussed earlier), I acknowledge its usefulness for specific practical applications in astrodynamics/celestial mechanics.

There may have been implicit tests of your frames-related modifications after 2006 among people at NASA, ESA ... Yet, I cannot remember any systematic "before and after" tests of these massive code changes that were made accessible to the Celestia dev team. Since Andrew and I are hunting a VERY subtle issue related to Celestia's frame-definition, I made a respective remark above (merely as a note of caution, while time was lacking for going through the post-2005-code most carefully).

...
It is however rotated 90 degrees from the conventional ecliptic coordinate system which has +z up instead of +y as in Celestia. That rotation certainly has caused a lot of headaches over the years.

Here it seems to me you effectively agree that the status of clarity about Celestia's coordinate frame was far from perfect in these earlier days of Celestia...

Fridger

[ajtribick]

Lets not get confused here: I have not implemented any reference frame stuff for .stc files. The SkyPlane reference frame (which matches Grant's spreadsheet and Fridger's binary orbits) is available for .ssc objects only. And I am seriously considering backing it out of the code if it turns out that the variant coordinate system discussed in this thread is the most common version. Trust an astronomy-related discipline to use abominations like left-handed coordinates, I suppose it is to be expected in a field where cgs units run rampant.

Also I think we should perhaps return to the subject at hand rather than trying to find insults in each others' posts, which won't get us anywhere.

[chris]

Lets not get confused here: I have not implemented any reference frame stuff for .stc files. The SkyPlane reference frame (which matches Grant's spreadsheet and Fridger's binary orbits) is available for .ssc objects only.

Pardon the confusion: I'd forgotten this was ssc only. Without actually double checking the source, I thought that you'd also added code to translate from the sky plane coordinate system to the internal coordinate system at load time.

--Chris

[t00fri]

As Andrew showed, flipping the z-axis in the skyplane frame amounts to the transformation
ω→ω+180°, Ω→Ω+180° . After generating the corresponding visualbins.stc and spectbins.stc via Perl, it is indeed impressive to confirm within Celestia that the skyplane projection is not at all affected, while the radial velocity is! That is a nice cross-check at the end of the 'pipeline'.

++++++++++++++++++++++++++++++++++++++++++++++++++++++
In this post, I want to add a different twist to this discussion, by considering another transformation that leads to precisely the same orbits as ω→ω+180°, Ω→Ω+180°.
++++++++++++++++++++++++++++++++++++++++++++++++++++++

One of the initial orbital elements is the orbital inclination i . For binaries, the inclination is originally specified relative to the plane of the sky. For Celestia we need to transform i into the output Inclination wrto the ecliptical frame.



+++++++++++++++++++++++++++++++++
It is also directly related to the definition of the z-axis. Flipping the z-axix (in the skyplane frame) from positive to negative, amounts to the transformation i -> -i
+++++++++++++++++++++++++++++++++

For completeness, let me now consider my Perl routine RotOrbits line-by-line that does the transformation from the skyplane frame to Celestia's ecliptical frame. It will be instructive to see why the effect of ω→ω+180°, Ω→Ω+180° is identical to the replacement of i -> -i with Ω, ω untouched!

Here is the routine as a reference

Code: Select all

sub RotOrbits {

my($ra_deg,$del_deg,$P,$a_arcsec,$i,$PA_of_Node,$Epoch_of_peri,$e,$Arg_of_peri
,$dist_ly) = @_;
my $del_rad = -$del_deg*$pi/180.0;
my $ra_rad = $ra_deg*$pi/180.0 - $pi;
my $eps = $pi/180.0*23.4392911;
my $ii = $pi/180.0*(90.0 - $i);
my $om = $pi/180.0*($PA_of_Node - 270.0);
my $alpha = atan(cos($ii)*cos($pi/180.0*($PA_of_Node))/(sin($ii)*cos($del_rad) -
cos($ii)*sin($del_rad)*sin($pi/180.0*($PA_of_Node)))) + $ra_rad;
if( sin($ii)*cos($del_rad)-cos($ii)*sin($del_rad)*sin($pi/180.0*$PA_of_Node) < 0 ) { $alpha = $alpha + $pi };
my $delta=asin(cos($ii)*cos($del_rad)*sin($pi/180.0*$PA_of_Node)+sin($ii)*sin(
$del_rad));
my $lambda=atan((sin($alpha)*cos($eps)+tan($delta)*sin($eps))/cos($alpha));
if( cos($alpha) < 0 ) { $lambda = $lambda + $pi };
my $beta = asin(sin($delta)*cos($eps) - cos($delta)*sin($eps)*sin($alpha));
my $alphaOm = atan(cos($om)/(-sin($del_rad))/sin($om)) + $ra_rad;
if( -sin($del_rad)*sin($om) < 0 ) { $alphaOm = $alphaOm + $pi };
my $deltaOm = asin(cos($del_rad)*sin($om));
my $lambdaOm = atan((sin($alphaOm)*cos($eps) +
tan($deltaOm)*sin($eps))/cos($alphaOm));
if( cos($alphaOm) < 0 ) { $lambdaOm = $lambdaOm + $pi };
my $betaOm = asin(sin($deltaOm)*cos($eps) -
cos($deltaOm)*sin($eps)*sin($alphaOm));
my $sign = $betaOm > 0? 1.0:-1.0;
my $dd = acos(cos($betaOm)*cos($lambdaOm - $lambda - $pi/2.0))*$sign;
$Period = $P;
$SemiMajorAxis = $dist_ly*63239.7*tan($pi/180.0*$a_arcsec/3600.0);
$Eccentricity =  $e;
$Inclination = 90 - $beta/$pi*180;
$AscendingNode = $lambda/$pi*180 + 90  - floor(($lambda/$pi*180+90)/360.0)*360;
$ArgOfPeri = $Arg_of_peri + $dd/$pi*180 - floor(($Arg_of_peri + $dd/$pi*180)/360.0)*360;
$MeanAnomaly = 360*((2000.0 - $Epoch_of_peri)/$P - floor((2000.0 - $Epoch_of_peri)/$P));
}


Since we are now talking about lots of math and since LaTeX still doesn't work with the new server, I will continue with a well commented Maple worksheet. Here it is




Note, the only orbital output parameters (in the Celestia ecliptical frame) that are affected by
ω→ω+180°, Ω→Ω+180° or equivalently i -> -i are only the following three

$Inclination
$AscendingNode
$ArgOfPeri

This I have also checked numerically, of course. There are now various possible interpretations about using i -> -i. Before getting to that, I let you first digest this result, perhaps...

Fridger

[ajtribick]

The preservation of the skyplane orbits can also be seen from the form of the Thiele-Innes constants, which only contain terms in cos(i), which is unchanged under i→-i.

For the radial velocity equation you need to know that the factor K which represents the semi-amplitude contains sin(i), which picks up a factor of -1 under i→-i. (This sin(i) term leads to the well-known mass-inclination degeneracy for the majority of the exoplanets catalogue).

[t00fri]

The preservation of the skyplane orbits can also be seen from the form of the Thiele-Innes constants, which only contain terms in cos(i), which is unchanged under i→-i.

For the radial velocity equation you need to know that the factor K which represents the semi-amplitude contains sin(i), which picks up a factor of -1 under i→-i. (This sin(i) term leads to the well-known mass-inclination degeneracy for the majority of the exoplanets catalogue).

Yes, certainly. My main motivation was also to do a further careful check of my Perl routine...
After decades of complex scientific computing, I only believe a result once I see an explicit test as well. That also applies to my own results, of course

Fridger

[ajtribick]

I guess we were addressing slightly different questions there, "does the substitution have these effects" versus "should the substitution have these effects".

Incidentally I personally prefer choosing inclination such that sin(i) is positive, as this means we can always take the semi-amplitude in the radial velocity equation to be positive, hence my describing the necessary transformation as ω→ω+180°, Ω→Ω+180°. Of course provided you are consistent in following everything through it all works out in the end.

[t00fri]

At this point of understanding, I still see 3 distinct stages where a modiified definition relative to the astronomical standard definitions might have (accidentally) occured and thus led to an inconsistent radial velocity (RV).

• in the analyses for the published binary data.
• in my Perl routine that transforms the given orbital elements
  from the skyplane to Celestia's variant of the ecliptical frame
• in the Celestia code that produces the final rendering

Each step must be examined carefully in turn.

The first possibility seems least probable, since the data used represent classical sources that have been scrutinized by many scientists (notably PhD students!! ) in the course of time.

Here is a useful PhD thesis (2004) by David Ramm (Univ. Canterbury, New Zealand),
entitled

A spectroscopic study of detached binary systems using precise radial velocities.
http://www.ir.canterbury.ac.nz/bitstrea ... lltext.pdf

Unless a better reference shows up, I suggest we take his Appendix A (Definitions and Theoretical Background) as a common basis for our further discussions. I found the presentation quite careful and well-organized. And...he uses right-handed frames throughout!

What I particularly like as a theorist is the presentation of the relevant frames and their transformations based on the familiar conservation laws of the angular momentum vector J = μ (x x v), the so-called Runge-Lenz or Perihelion vector, the total energy and the linear momentum.

Fridger

[ajtribick]

Good to see the explicit definitions there. Note that the reference frame depicted in figure A.3 differs from the one in Grant's spreadsheet: while the x axis is in the same direction, both the y and z axes point in the opposite directions. (This is of course a 180 degree rotation, so the frame is still right-handed). This incidentally results in an i=0 orbit being face-on and clockwise, as opposed to face-on and anticlockwise.

One interesting point to note is in section A.11, spectroscopic binary orbits typically quote the value of ω for the primary, while visual and astrometric orbits typically quote it for the secondary. Whether this might have any bearing on the resultant coordinate system of a combined visual+spectroscopic solution I don't know.

[ajtribick]

Ok so making the transformation between the SkyPlane as described by that thesis and the SkyPlane as Celestia is currently defining it (Grant's spreadsheet, Fridger's binary code, and my SkyPlane frame in the SVN version) can be done by making the substitutions:

ω→ω+180°
i→180°-i
Ω→180°-Ω

This can be seen from the Thiele-Innes constants (where I am also giving the constants relevant for the z-coordinate):

A = a(cos Ω cos ω - sin Ω sin ω cos i)
B = a(sin Ω cos ω + cos Ω sin ω cos i)
C = a(sin i sin ω)
F = a(-cos Ω sin ω - sin Ω cos ω cos i)
G = a(-sin Ω sin ω + cos Ω cos ω cos i)
H = a(sin i cos ω)

And given that

X = Ax + Fy
Y = Bx + Gy
Z = Cx + Hy

and in transforming between the two, we need to leave X intact while also doing Y→-Y, Z→-Z, it can be seen that the substitutions of ω, i and Ω will leave the Thiele-Innes constants A and F intact while introducing a factor of -1 to B, C, G and H.

I've attached some Perl code which does the transformations both for Celestia skyplane to ecliptic and thesis skyplane to ecliptic. It uses the Math::Quaternion module to do the various rotations.

[t00fri]

Andrew,

good idea using at last the Math::Quaternion module for doing the required binary orbit transformations. It's what I wanted to do for my visualbins/spectbins Perl scripts since a long time, after I had already used Math::Quaternion in my galaxies.pl script right from the start for the same purpose.

So, if I understood your result correctly, you claim that e.g. in my visualbins.pl / spectbins.pl Perl scripts the substitution

ω→ω+180°
i→180°-i
Ω→180°-Ω

in the (ω, i, Ω) orbital elements from (standard) astronomical catalogs takes into account ALL of Celestia's pecularities. In particular, while not affecting the sky plane projections, it leads to the physically correct radial velocities in Celestia's rendering.

One aspect that perhaps still needs some careful checking, is to make sure that the standard domains for these three parameters (ω, i, Ω) are never violated after that substitution.

Finally, one may ask, why in Celestia the +Z-axis of the skyplane frame is pointing towards Earth, unlike the standard convention in the astronomical literature? The physical reason for having conventionally the Z-axis point away from the observer (as in that thesis) is that all radial velocities measured will be consistent with positive doppler shifts being due to receding bodies. Moreover, according to the astronomical definition, 0 ≤ i < π/2 when the orbital motion is direct, whereas π/2 ≤ i ≤ π when the motion is retrograde, In the inclinations used for Celestia, this would be reversed (i→180°-i) and hence confusing to astronomers.

Fridger

[ajtribick]

So, if I understood your result correctly, you claim that e.g. in my visualbins.pl / spectbins.pl Perl scripts the substitution

ω→ω+180°
i→180°-i
Ω→180°-Ω

in the (ω, i, Ω) orbital elements from (standard) astronomical catalogs takes into account ALL of Celestia's pecularities. In particular, while not affecting the sky plane projections, it leads to the physically correct radial velocities in Celestia's rendering.

No, that is the transformation needed to go from the skyplane as defined in the thesis to the skyplane as currently defined in Celestia.

I've investigated what happens if you try to use the thesis definition of the coordinate system to interpret the Pourbaix et al. results, unfortunately it leads to different skyplane plots, and since we know that the existing spectbins.stc does give the correct skyplane plot, it means that Pourbaix et al. are not using this coordinate system. Essentially, you cannot match all the data with the elements given in Pourbaix et al. without having to use a left-handed coordinate system. You have to use ω→ω+180°, Ω→Ω+180° then use the old version of the transformation.

[t00fri]

So, if I understood your result correctly, you claim that e.g. in my visualbins.pl / spectbins.pl Perl scripts the substitution

ω→ω+180°
i→180°-i
Ω→180°-Ω

in the (ω, i, Ω) orbital elements from (standard) astronomical catalogs takes into account ALL of Celestia's pecularities. In particular, while not affecting the sky plane projections, it leads to the physically correct radial velocities in Celestia's rendering.

No, that is the transformation needed to go from the skyplane as defined in the thesis to the skyplane as currently defined in Celestia.

I've investigated what happens if you try to use the thesis definition of the coordinate system to interpret the Pourbaix et al. results, unfortunately it leads to different skyplane plots,
and since we know that the existing spectbins.stc does give the correct skyplane plot, it means that Pourbaix et al. are not using this coordinate system.

In my nomenclature, the proper sky plane (xy) plots from my spectbins.stc are correct without any substitutions, while only the "transverse" sky plane plots (radial velocities) are incorrect without ω→ω+180°, Ω→Ω+180° substitution.
Agreed?

Essentially, you cannot match all the data with the elements given in Pourbaix et al. without having to use a left-handed coordinate system.

It seems crucial here to specify precisely, which Pourbaix binaries you have examined in detail and which radial velocity data you actually took? Some of Pourbaix' radial velocities from 1999 may be simply incorrect. So a test of the whole lot (Perl) would allow a much safer conclusion!

Fridger

[ajtribick]

To aid comparison I attach a Perl script that generates SVG files of the radial velocity and skyplane plots from the Pourbaix data, under the following assumptions: ω refers to the orbit of the secondary, and the reference frame is the left-handed frame I described earlier (i.e. X points north, Y points east, Z points away from Earth).

In the orbit plot on the left-hand side, the plot is generated with North pointing down and East pointing right, for easier comparison with the orbit plots in the Sixth Orbit Catalog. For the radial velocities, the red curve is that of the primary, the green curve is that of the secondary. This can be compared with the radial velocity plots in The 9th Catalogue of Spectroscopic Binary Orbits, which uses filled squares to plot the primary star RV measurements, and open squares for the secondary.

Systems I have explicitly checked: HIP 677, HIP 2941, HIP 4463, HIP 71683/1, all of which match under the assumptions I describe above.

[t00fri]

To aid comparison I attach a Perl script that generates SVG files of the radial velocity and skyplane plots from the Pourbaix data, under the following assumptions: ω refers to the orbit of the secondary, and the reference frame is the left-handed frame I described earlier (i.e. X points north, Y points east, Z points away from Earth).

In the orbit plot on the left-hand side, the plot is generated with North pointing down and East pointing right, for easier comparison with the orbit plots in the Sixth Orbit Catalog. For the radial velocities, the red curve is that of the primary, the green curve is that of the secondary. This can be compared with the radial velocity plots in The 9th Catalogue of Spectroscopic Binary Orbits, which uses filled squares to plot the primary star RV measurements, and open squares for the secondary.

Systems I have explicitly checked: HIP 677, HIP 2941, HIP 4463, HIP 71683/1, all of which match under the assumptions I describe above.

Thanks a lot Andrew,

that's excellent for further comparisons....

Fridger