Lab 2, model planetary motion

A much simpler approach is to fit the data directly to the equations (in a Cartesian coordinate system with the origin as the Earth):
	x =  
          (Subscript[A, x]) cos(ω t  +  (Subscript[α, x]));
y =  
          (Subscript[A, y]) cos(ω t  + (Subscript[α, y]));
\(TraditionalForm\`z = \ \(A\_z\) \(\(cos(ω\ t\  + \ α\_z)\)\(.\)\)\)
However, the fitting is a nonlinear problem.  So to do it properly is difficult.  I'll be agreeable if you just `read' from the figures to get the parameters.  As you can see, in general, this will not be a circle as we don't have (x^2) + (y^2) + (z^2) = 1. 

Below is a much more sophiscated way of doing things.  Hope you can learn something from it.  Note that I have used linear algebra methods several times.   Before we can reading in the data in Mathematica, I removed the column contains * and m and A, and the date/time columns to use only the Julian date.   I also make the sign + or - in a separate column.  It is important to do this to treat the negative declination correctly.  I don't explain Mathematica features, assuming you know them. Now we read in the Moon data as follows:

\(\(moonrawdata = 
        Import["D:\WangJS\teaching\CZ3102\moon.dat", \ "table"];\)\)

The raw data are in the following format (first and last 2 lines)  [Julian date, RA hr, RA min, RA sec, sign, DEC deg, DEC min, DEC sec]:

Short[moonrawdata, 10]

We convert them in more suitable format for computation,we subtract the Julian date value at 1 Jan 2002  at 00.00 (12 midnight) (2452275.16666667), so that time begin at 0 for that date and time.  We also use radian for the two angles.  First number is RA and second is DEC, both in radian, in terms of the angle 0<ϕ<2π, and  0< θ < π.  This is in standard spherical polar coordinate system.   Define sign[ ] to convert string \(TraditionalForm\`\(\(\*"""\)\(+\)\) \*"""\) and "-" to the number 1 and -1.

sign["+"] := 1;
sign["-"] := -1;

moondata  =  
        Map[ { #[[1]] - 
                2452275.16666667,  (#[[2]] + #[[3]]/
                      60.0 + #[[4]]/3600.0)*
                Pi/12.0,  (90.0  -  
                    sign[#[[5]]]*(#[[6]] + #[[7]]/
                            60.0 + #[[8]]/3600.0))* Pi/180.0} &,  
          moonrawdata];

Short[moondata, 10]

We will use the model I, geocentrical circular orbit.  It turns out the mathematics involved for other models are rather complicated.  First we need to determine the plane the orbit lays.  Note that a point on a plane (the equation for a plane) is
	A x  +  B y  +  C z  =  0
where (A, B, C) is the direction vector normal to the plane.  This direction can be determined by a least-quare method.  Each data point is a set (x, y, z).   Thus we have an over-determined system of equations to solve A,  B,  C.   The matrix is:

planeMat  =  
        Map[ {Sin[#[[3]]]*Cos[#[[2]]],  
              Sin[#[[3]]]* Sin[#[[2]]],  Cos[#[[3]]]} &,  
          moondata];

Short[planeMat, 10]

The matrix has the property
 	M (({{A

                    B}, {C}}))=0
thus we can used pseudo-inverse to find equation for A,  B,  C.  However, this is not correct.  Since  if A,  B,  C is a solution, so is λA,  λB,  λC also a solution.  So what is important is only the direction of the plane.  We can fix this  problem by requiring a unit vector. Thus we have to solve the problem of minimising 
	\(TraditionalForm\`∑\_i\ \((A\ x\_\(\(i\)\(\ \)\) + \ B\ y\_i + C\ z\_\(\(i\)\(\ \)\))\)\^2\)
subject to the constraint 
	(A^2) + (B^2) + (C^2) = 0. 
Using Lagrange multiplier method, this results to the following eigen value problem: 
 	 (M^T)M (({{A}, {B

                    C}})) =λ  (({{A

                      B}, {C}})).
  The solution is find using Eigensystem[ ].

M2  =  Transpose[planeMat] .  planeMat

res  =  Eigensystem[M2]

The three eigen vectors correspond to exactly the three directions of our new coordinates, the last one being our normal direction of the plane.  Why?  e1, e2, and e3 will be the unit vectors of three direction of the axis.  They are otho-normal.  R is precisely the needed rotation matrix to transform from the old coordinate system (x, y, z) into the new one (x', y', z').

{e1, e2, e3}  =  R =  res[[2]]

newcoord   =  Map[R . # &, planeMat];

coordXY  =  Map[ Take[#, 2] &, newcoord];

Short[coordXY, 12]

coordXY give pair of (x, y) coordinates in the new system, note that z is approximately zero.  We can make a plot of the orbit projected onto this plane.  We see the moon orbit is nearly a perfect circle. We plot the first 170 points only (about 1 revolution).

ListPlot[Take[coordXY, 170],  AspectRatio  ->  1]

We now fit the motion by  x' (t) = cos(ω t  +  α), y' (t) = sin(ω t  +  α).  This is a "uniform circular motion". To do this, we compute the angle by arctan(y/x) and fit it to a straight line.  We need to be careful when we do this as the inverse tangent function is multiple values, we need to choose the right branch.  To do this, we need to define a little funtion to remember the phase (wind is the winding number):

angle[co_] := 
      Module[{wind, a, b, dat, l  =  Length[co],  
          i}, 
wind  =  
          0; 
b =  dat  =  {ArcTan[co[[1, 1]], 
                co[[1, 2]]]}; 
For[i  = 2, 
          i ≤  
            l, ++i, 
a  =   ArcTan @@  
              co[[i]]; 
If[
            b  >  0.0  &&  
              a  <  0.0,  ++wind]; 
b  =  
            a; 
AppendTo[dat,  
            a  +  2.0*Pi*
                wind];
]; 
dat
]

angledata  =  angle[coordXY];

Short[angledata, 30]

This is the angle data vs the point number (each point is 4 hours apart for moon, or 4/24=1/6 of a Julian date).  We combine the equally spaced time data (Julian date offset from 1 Jan 2002) with the angle data to make a straight line fit.  The result is  α=2.69182+0.229975 t, see below.  This little result is the final result for (x'(t),y'(t)).   Of course, z'(t)=0.   To get (x,y,z), we need to inverse transform (R^(-1)) = (R^T).

juliavsangle  =  
        Table[ {(j - 1)*1.0/6.0,  angledata[[j]]},  {j,  1,  
            Length[angledata]}] ;

alpha[t_] = Fit[juliavsangle,  {1,  t}, t]

The last (offset) Julian date 768.0 is 0:00 clock (midnight) of Sunday, 8 Feb 2004.  Thus 9:00pm of Monday, 9 Feb 2004 is  t = 768.0 + 1.0  +  
            21.0/24.0 = 769.875.  We define the function predict[t] to give the predicted value of RA and DEC.  We have invert from Cartesian to spherical polar system:
	\(TraditionalForm\`x\  = \ \(sin(θ)\)\ \(cos(ϕ)\), 
y = \ \(sin(θ)\)\ \(sin(ϕ)\), 
z\  = \ \(\(cos(θ)\)\(.\)\)\)

predict[t_] := 
      Module[{a, xp, yp, zp, x, y, z, theta, phi, RA, RAm, RAmin, RAsec, DEC, 
          DECm, DECmin, DECsec, s}, 
a  =  
          alpha[t]; 
xp  =  
          Cos[a]; 
yp  =  
          Sin[a]; 
zp  =  
          0.0; 
{x, y, z}  =  
          Transpose[R] . {xp, yp, zp}; 
theta  =  
          ArcCos[z]; 
phi  =  
          ArcTan[x, y]; 
If[phi < 0,  
          phi  +=  2.0*Pi]; 
RA =  
          phi*12.0/Pi; 
RAhr  =  
          Floor[RA]; 
RAm  =  (RA - 
                RAhr)*60.0; 
RAmin  =  
          Floor[RAm]; 
RAsec  =  (RAm  -  
                RAmin)*60.0; 
Print["RA = ",  
          RAhr,  "hr ",  RAmin,  "min ",  
          RAsec , "sec"]; 
DEC  =  
          90.0  -  theta*180.0/Pi;  
If[
          DEC  ≥ 0,  s = "+", s = "-"; 
          DEC  =  -DEC]; 
DECdeg  =  
          Floor[DEC]; 
DECm  =  (DEC  -  
                DECdeg)*60.0; 
DECmin  =  
          Floor[DECm]; 
DECsec  =  (DECm  -  
                DECmin)*60.0; 
Print["DEC =",  
          s, " ",  DECdeg,  "deg ",  DECmin,  "min ",  
          DECsec,  "sec"];
]

predict[769.875]

According to the JPL "Ephemeris Generator" (at http://ssd.jpl.nasa.gov/cgi-bin/eph) the correct answer on 2004-Feb-09 21:00 (Local Singapore time) is
RA =  12 hr  05 min  03.94 sec
DEC = + 03 deg  46 min  48.0 sec
So my answer is off by about 15min in RA and about 1 deg off in DEC.  Not bad given that we used a very simple model.  Things will not be so lucky for Mars but I'll stop here.

Before I really quit,  I just point out that fitting the x(t),  y(t),  and z(t) to
x = x0 * cos(t+d1), y = y0 * cos(t+d2), z = z0 * cos(t+d3) 
is WRONG!  why, a plot will make this obvious:

XTdata  =  
        Map[ {#[[1]], Sin[#[[3]]]*Cos[#[[2]]]} &,  
          moondata];

Short[XTdata, 10]

Moon data x(t) vs t where x is x-component in old Cartesian coordinate system, t is Julian date offset (starting 1 Jan 2002 00:00 midnight).  Note that t increases by 1 means 1 day has past.

ListPlot[Take[XTdata, 200]]

The curve cos(t+2.0); I take x0 = 1, and d1 = 2.0.  It is obvious that the curve cannot possibly fit the data no matter what values of x0 and d1 you choose.  Why?

Plot[1.0  Cos[t + 2.0],  {t,  0, 30}]