From: martin@space.ualberta.ca (Martin Connors)
Subject: Re: Mars Observer trajectory
Reply-To: martin@space.ualberta.ca (Martin Connors)
Organization: University Of Alberta, Edmonton Canada
Date: Tue, 6 Oct 1992 15:50:23 GMT

In article <BvoKGF.GzM@well.sf.ca.us> collins@well.sf.ca.us (Steve  
Collins) writes:
> 
> I have some information on the Mars Observer trajectory.
> Now for the data. I snagged the latest trajectory and have the
> the following spacecraft state information:
> Epoch: 10/5/92 00:00:00 UTC
> Position: 1.46030475e8 3.16372935e7 1.30237608e7   km
> Velocity -7.02799688e0 3.00312953e1 1.22052289e1   km/sec
> relative to the sun in the J2000 coordinate frame.
> These can be expressed as classical elements about the sun, but I am
> not familiar with the nomenclature that the available software uses.
For those wishing to do so (perhaps for use with ephem) the following  
SHORT listing is posted (available upon request in FORTRAN too). NOTE that  
the units are AU and AU/day for the Position and Velocity so you will need  
to convert the above to those units first. This routine is based on  
Brouwer and Clemence's book (at home right now, title something like  
"Celestial Mechanics")....
#include <stdio.h>
#include <math.h>
#define pi 3.1415926536
#define twopi 6.283185308
#define radian 57.2957795
#define invrad 1.745329252e-2

main()
{
double meananomaly, eccentricity, *trueanomaly, *eccanomaly;
double lilomega,bigomega,inclination,a;
double ax,bx,px,qx,ay,by,py,qy,az,bz,pz,qz; double x,y,z,r,xd,yd,zd,n;
double obliquity=23.445785;  /* 1950.0 */
double cosl,sinl,cosb,sinb,sini,sino,cosi,coso,cose,sine;
double r2,v2,rrdot,sinu,cosu,b,u;
double k=0.01720209895; /* Gaussian constant */
/*
printf("\n enter x,y,z and derivatives");
scanf("\n %lf %lf %lf",&x,&y,&z);
scanf("\n %lf %lf %lf",&xd,&yd,&zd);
*/
x=1.074433; y=-0.483451; z=-0.214434;
xd=0.000163; yd=0.017834; zd=0.00763;
/* derivatives are in AU/day - divide by Gaussian constant */
xd /= k; yd /= k; zd /= k;
r2=x*x+y*y+z*z; r=sqrt(r2);
v2=xd*xd+yd*yd+zd*zd;
rrdot=x*xd+y*yd+z*zd;

a=1./(2./r-v2);
eccentricity=sqrt(rrdot*rrdot/a+(r*v2-1.)*(r*v2-1.));

sinu=rrdot/(eccentricity*sqrt(a));
cosu=(r*v2-1.)/eccentricity;
u=asin(sinu);
meananomaly=u-eccentricity*sinu;
b=a*sqrt(1.-eccentricity*eccentricity);
obliquity *= invrad;
coso=cos(obliquity);   sino=sin(obliquity);
px= x*cosu/r - xd*sinu*sqrt(a);
qx=(x*sinu/r + xd*sqrt(a)*(cosu-eccentricity))*a/b;
py= y*cosu/r - yd*sinu*sqrt(a);
qy=(y*sinu/r + yd*sqrt(a)*(cosu-eccentricity))*a/b;
pz= z*cosu/r - zd*sinu*sqrt(a);
qz=(z*sinu/r + zd*sqrt(a)*(cosu-eccentricity))*a/b;
n=0.01720209895/sqrt(a*a*a); /* mean motion in rad/day */
sini=sqrt(py*py*sino*sino-2.*py*pz*sino*coso+pz*pz*coso*coso
         +qy*qy*sino*sino-2.*qy*qz*sino*coso+qz*qz*coso*coso);
inclination=asin(sini);
sinl=(-py*sino+pz*coso)/sini;
cosl=(-qy*sino+qz*coso)/sini;
lilomega=asin(sinl);
bigomega=acos(px*cosl-qx*sinl);
meananomaly *= radian;
lilomega *= radian;
bigomega *= radian;
inclination *= radian;
obliquity *= radian;
printf("\n  meananomaly, eccentricity, a");
printf(" %lf %lf %lf",meananomaly,eccentricity,a);
printf("\n  lilomega, bigomega, inclination");
printf(" %lf %lf %lf",lilomega,bigomega,inclination);
}

anomaly(meananomaly,eccentricity,trueanomaly,eccanomaly)

double meananomaly, eccentricity, *trueanomaly, *eccanomaly;
{
/* This routine finds the values of the eccentric and true      */
/* anomalies in elliptical motion,                              */
/* given the mean anomaly and the eccentricity.                 */
/* The mean anomaly, eccentric anomaly, and true anomaly are    */
/* measured in radians.                                         */

double delta, m, x, y, z;

  m = meananomaly - twopi * floor(meananomaly / twopi);
  *eccanomaly = m;
  delta = *eccanomaly - (eccentricity * sin(*eccanomaly)) - m;
  while (fabs(delta) >= 1.0e-12)
    {
      delta = delta / (1.0 - (eccentricity * cos(*eccanomaly)));
      *eccanomaly = *eccanomaly - delta;
      delta = *eccanomaly - (eccentricity * sin(*eccanomaly)) - m;
    }
  x = sqrt(( 1.0 + eccentricity) / (1.0 - eccentricity));
  y = cos(*eccanomaly/2.0);
  z = sin(*eccanomaly/2.0);
  *trueanomaly = 2.0 * atan(x*tan(*eccanomaly/2.0));
}  /*end of function anomaly*/

Hope this helps, and have fun....
Martin Connors
University of Alberta
{ accuracy of listing results not vouched for. MC and UA bear no  
responsibility etc. }
Newsgroups: sci.space