Translation of orbital plane into ecliptic co-ords

[bdm]

I'm attempting to create a fictional asteroid belt, and I have a small mathematics problem with the orbital co-ordinates.

I have created a spreadsheet that randomly creates asteroids. The orbital parameters that are angles are Inclination, Ascending node, Mean Longitude and Periastron.

To make the orbits realistic, the original reference plane for this fictional asteroid belt needs to be in reference to the orbital plane of a large planet that is nearest to the asteroid belt. (I call this the Local Reference Plane.) This makes sense, because this planet would be the major source of perturbations for these asteroids, and if the asteroids had no orbital inclination they would orbit in roughly the same plane as the large planet.

The trouble is Celestia uses a global co-ordinate system, and it is necessary to convert the orbital parameters from one reference plane to another. I started with the idea of converting the Inclination and Ascending Node of the local reference plane to a vector, but then I got stuck.

What I need are expressions that can convert the angular orbital terms from the Local Reference Plane to Celestia's ecliptic coordinates. I expect such expressions to require the Inclination and Ascending Node of the Local Reference Plane as terms.

(It may also be possible for Celestia to specify a local reference plane in an SSC file, but I haven't found this in the documentation. Adding local reference planes to STC/SSC files may be a useful feature.)

[bdm]

This problem is very similar to the following problem on the Earth's surface:

Given a point on the Earth's surface with a specified longitude and latitude, find the longitude and latitude of another point x degrees away in direction a.

[Spaceman Spiff]

I've not got time to type in equations, but this is Spherical Geometry. I think we dealt with this before a long time ago converting ecliptic co-ordinates to galactic co-ordinates. If you can track down the topic, it's a matter of understanding each term, and adapting the equations.

Spiff.

[Apollonian]

There may be a fancier way to do it, but the simple way would be to compute a position and velocity and rotate the frame.

http://en.wikipedia.org/wiki/Orbital_state_vectors
http://en.wikipedia.org/wiki/Coordinate_rotation

The following book is a fairly good reference, though not the best
http://www.celestrak.com/software/vallado-sw.asp

You can look at his source to get the algorithm for rv2oe (radius & velocity to orbital elements). Otherwise, I can go through the calculations here if you prefer.

[bdm]

I had a quick look at the co-ordinate rotation pages that Apollonian linked in. The two-dimensional rotations I remember from twenty years ago, but it's obvious now why the three-dimensional rotations weren't taught to me. That 3 x 3 matrix would have been difficult to memorise for an exam.

First I'll restate the problem - convert an inclination and ascending node from one co-ordinate system to another. It's not important to convert mean longitudes and perihelion co-ords because for fictional solar systems one set of starting co-ordinates is as good as any other.

I've been giving the matter some thought. An inclination and ascending node is like a vector drawn on a great circle of a sphere. The inclination becomes the angular length of the vector on the great circle and the ascending node becomes the direction from the zero point of the sphere (such as a north or south pole). The tail of the vector rests on the pole of the sphere. Now if we rotate the sphere using another similar vector, where does the head of the first vector end up?

When visualised in this way, it is clear that the axis of rotation is perpendicular to the second vector, and must therefore pass through the equator of the sphere at right angles to the second vector.

It is tough trying to derive it from first principles. I'll see what I can make of Apollonian's linked pages. I think the biggest problem I have now is converting the problem to fit the supplied formulae.

What is the accepted convention for mapping Inclination and Ascending Node to an x, y, z unit vector?

[t00fri]

bdm,

your problem is conceptionally very straightforward and amounts to what Apollonian wrote. However in practice (e.g. in the context of Celestia) one has to be most careful with the partly unusual definition of the various coordinate systems involved.

In many of my catalog preparations (binaries and my recent deepsky.dsc galaxy catalog) I had to do such type of rotations explicitly in my PERL code.

Find below my simple PERL subroutine as an example that rotates (by means of quaternions, of course) from the "face-on" galaxy frame with origin in the galaxy center to the standard ecliptic frame characterized by [axis, angle] in Celestia.

The rotation angles involved are RA,DEC, the ecliptic angle epsilon= 23.4392911 deg, the position angle PA and inclination i of the galaxy relative to face-on. The tricky part is to locate the respective axes around which the individual rotations have to be performed, keep track of ALL SIGNS (!) and to remember that rotations are NOT commuting!!

Since PERL code is rather self-explanatory, here is the routine:

Code: Select all

sub orientation {
  my($ra, $dec, $Type, $i, $PA) = @_;
  my $decrot = Math::Quaternion::rotation(deg2rad(90-$dec),1,0,0);
  my $rarot  = Math::Quaternion::rotation(deg2rad(90-$ra*15),0,1,0);
  my $radecrot = $decrot*$rarot;
  my $incrot = Math::Quaternion::rotation(deg2rad($i),0,0,1);
  my $pa =  Math::Quaternion::rotation(deg2rad($PA),0,1,0);
  my $eclipticsrot = Math::Quaternion::rotation(deg2rad($epsilon),1,0,0);

  #
  if ($Type =~ /E/) {
      $rot = $pa*$radecrot*$eclipticsrot;
  } else {
    $rot = $incrot*$pa*$radecrot*$eclipticsrot;
  }

  $angle=$rot->rotation_angle*180/pi;
  @v=$rot->rotation_axis;
}


All quaternionic rotations are parametrized in an intuitive form via rotation( angle, axis-unitvector). The amazing thing is that this little routine rotates 10000+ galaxies in a fraction of a second ...

It outputs [axis,angle] given the input rotation angles $ra,$dec,$i,$pa from the published galaxy catalogs (along with some elliptic/spiral type info, since our elliptical galaxy templates are spheroidal and the spirals are spherical "pancakes")

If one has to rotate whole orbits characterized by the standard orbital parameters, the procedure is analogous but somewhat more lengthy, since almost all of the orbital parameters get transformed.

I have written a corresponding PERL routine as well.

Good luck,

Bye Fridger

[Apollonian]

Hmmm, if you want to just map inclination and ascending node, you could use my method, but you would need the following points:

- Set true anomaly to zero.
- Set argument of periapse to zero.
- Set semi-major axis and eccentricity to correspond to the following two conditions

- You need a velocity corresponding to the given inclination
- You need a radius corresponding to the ascending node location

However, for what you are doing, this isn't the best method.

I didn't really follow what you are doing with vectors, but here is what I would suggest:

Basically, what you are doing is a simple frame transformation. The inclination and ascending node can define a right-handing frame if you treat them like vectors. The ascending node is a unit vector in the direction of the ascending node. The inclination is a unit vector orthogonal to the orbit reference plane.

In the case of Earth, these vectors would correspond to
- Ascending node toward the "First point in Ares" (by convention)
- Inclination vector in the direction of the north pole

From there, all you need is the orientation of the new frame. So, if you are converting to ecliptic coordinates, this would correspond to the "Mean Obliquity of the Ecliptic". Essentially, "Mean Obliquity of the Ecliptic" is the same as "The inclination of the Ecliptic with respect to the Equator" (to first approximation).

In the case of Equator-Ecliptic, the convention for right ascension of the ascending node is the same (to first approximation). So then, it becomes a coordinate transformation as I previously mentioned - only you don't have to convert to radius and velocity.

The result is a coordinate transformation which will generalize. If you take inclination = 0 and RAAN (Right ascension of ascending node) =0 you can find the generalized transformation between the two frames. Then, you can use that transformation (matrix, quaternion, your choice) to convert any orbit from one to the other.

The following sources may be of greater help as the Euler Angle transformation correspond more obviously to orbital elements.

Looking at the mathworld article...
Phi = RAAN
Theta = inclination
Psi = Argument of periapse (zero in this case)

http://mathworld.wolfram.com/EulerAngles.html
http://en.wikipedia.org/wiki/Euler's_angle

I hope that helps.

[t00fri]

bdm,

the solution of your problem should be completely contained in my PERL routine that transforms the standard orbit parameters [from published catalogs] as measured in the sky-plane (input) into the orbit parameters in Celestia's ecliptic frame as needed for Celestia's *.stc files. Most importantly, Celestia's tricky conventions are taken care of...

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) = @_;
$del_rad = -$del_deg*$pi/180.0;
$ra_rad = $ra_deg*$pi/180.0 - $pi;
$eps = $pi/180.0*23.4392911;
$ii = $pi/180.0*(90.0 - $i);
$om = $pi/180.0*($PA_of_Node - 270.0);
$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 };
$delta=asin(cos($ii)*cos($del_rad)*sin($pi/180.0*$PA_of_Node)+sin($ii)*sin($del_rad));
$lambda=atan((sin($alpha)*cos($eps)+tan($delta)*sin($eps))/cos($alpha));
if( cos($alpha) < 0 ) { $lambda = $lambda + $pi };
$beta = asin(sin($delta)*cos($eps) - cos($delta)*sin($eps)*sin($alpha));
$alphaOm = atan(cos($om)/(-sin($del_rad))/sin($om)) + $ra_rad;
if( -sin($del_rad)*sin($om) < 0 ) { $alphaOm = $alphaOm + $pi };
$deltaOm = asin(cos($del_rad)*sin($om));
$lambdaOm = atan((sin($alphaOm)*cos($eps) + tan($deltaOm)*sin($eps))/cos($alphaOm));
if( cos($alphaOm) < 0 ) { $lambdaOm = $lambdaOm + $pi };
$betaOm = asin(sin($deltaOm)*cos($eps) - cos($deltaOm)*sin($eps)*sin($alphaOm));
$sign = $betaOm > 0? 1.0:-1.0;
$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));
}


If you want to use the routine, I can tell you any time the precise units of the input quantities in case of doubt.

That's all and straightforward to derive. It is close to trivial to translate the code into any other computer language.

The 10 input orbit parameters (and most of their their units) are enumerated in the argument list, the output parameters have the standard Celestia syntax and may be easily identified towards the end of the routine.

Of course this routine has been extensively tested.

Hence, with my routine it would be exceedingly simple to read out huge amounts of standard asteroid orbit data from the usual catalogs with PERL and print out a corresponding *.stc file with thousands of asteroids to be read in by Celestia. I guess that's precisely what you are up to. Of course you may as well follow the qualitative argumentations of Apollonian above and try to rederive the above results yourself .

Bye Fridger

[bdm]

This thread is rapidly turning into an excellent reference on rotations, thanks for your help

I'm going to need a few days to experiment with Excel. The Euler angles will probably be most useful, but the Perl code will also come in handy later.

(It helps to know that they're called Euler angles. Google is of no help if one does not know the correct terminology.)

[bdm]

An update -

I've successfully implemented a quick Excel spreadsheet that rotates co-ordinates from one plane to another.

Imagine an extrasolar planetary system with an asteroid belt in the same general plane as a large planet. The asteroids have random nodes and inclinations. To express them properly in Celestia, the coordinate frames of the asteroids must be rotated into the ecliptic plane. The spreadsheet I have created does this.

Here is a screenshot of a test image. The 16 bodies in this illustration have ascending nodes that vary in steps of 22.5 degrees, but the same inclination. The inclination and nodes are not in reference to the ecliptic (which runs horizontally in this image), but a different reference plane that is inclined 45 degrees to the ecliptic. The coordinates of the 16 planes are then rotated into the ecliptic plane so Celestia can display them.



Thanks to all who assisted with this project. Sometime in the next few weeks I will release a spreadsheet that generates random asteroid belts in arbitrary planes. The asteroids will have random orbits within a defined region. This spreadsheet will be of great assistance to the creators of fictional solar systems.

[selden]

That's great!

Hopefully you'll make the plane-changing spreadsheet available, too.

[bdm]

UPDATE

The thread on the forthcoming release of Asteroid Maker is here:

http://www.celestiaproject.net/forum/viewtopic.php?t=8089