3
$\begingroup$

All methods I've seen for computing planetary/stellar conjuctions (when two planets or a planet and a star appear close together in the sky because their angular separation, as measured from Earth, is small) are iterative: they look at the positions day by day and find the smallest angular separation.

EDIT (per comment): All of the constants (a1, b1, a2, b2, a3, c3) are known values. I'm trying to solve for t given these constants. There are, of course, an infinite number of solutions, so solutions for t in any given interval are fine too. Numerical solutions (where I plugin the constants as approximate numbers) are fine too, provided they are reasonably efficient (moreso than iterating).

I thought it would be easier with circular 2-dimensional orbits, so I tried this:

 
o[t_] = {Cos[t],Sin[t]} 
o1[t_] = a1*{Cos[a2*t+a3],Sin[a2*t+a3]} 
o2[t_] = b1*{Cos[b2*t+b3],Sin[b2*t+b3]} 

The above are the circular 2-dimensional orbit of 3 planets, where one planet's orbit has been normalized to have radius 1, period 2*Pi, and o[0] = {1,0}.

Treating "o" as Earth, the position of "o1" as viewed from Earth at time t is "o1[t]-o[t]", and "o2[t]-o[t]" for "o2".

The angles from Earth are "Apply[ArcTan,o1[t]-o[t]]" and "Apply[ArcTan,o2[t]-o[t]]". The difference between these angles is:

 
Abs[Apply[ArcTan,o1[t]-o[t]]-Apply[ArcTan,o2[t]-o[t]]] 

which is what I'm trying to minimize. In fact, because we're in 2 dimensions, the minimum value will be 0, so I just need to find the zeros of the function above (of course, the Abs[] becomes superfluous at this point).

I've tried several approaches, and nothing worked. Solve[] won't solve the equation for Reals and hangs when I don't give it a domain.

I've tried using dot products, and even tried solving the simpler problem of when two planets appear perpendicular (dot product is 0), to no avail.

Graphing with actual vales show there's no real pattern to the zeros of my function, or even a minimum distance between zeros.

I could try an iterative solution of some sort, but then I'm back to square one.

$\endgroup$
5
  • $\begingroup$ Do you want a solution in terms of the a and b variables, or is a numerical solution (for particular choices of the as and b's all right? $\endgroup$
    – march
    Aug 15, 2015 at 17:30
  • $\begingroup$ Thanks! Edited the question to reflect I'm looking for solutions in t given the constants, and that efficient numerical solutions are fine too. $\endgroup$
    – user1722
    Aug 15, 2015 at 17:44
  • $\begingroup$ I won't have time to look at this until later. In the meantime, look up NSolve and FindRoot, and probably someone else will post an answer. $\endgroup$
    – march
    Aug 15, 2015 at 18:07
  • $\begingroup$ As suggested by @march, FindRoot certainly works. For instance, with {a1 = 3.6, a2 = 1.3, a3 = 1.6, b1 = .9, b2 = 2.1, b3 = 2.6}, FindRoot[Apply[ArcTan, o1[t] - o[t]] == Apply[ArcTan, o2[t] - o[t]], {t, 0}]` yields {t -> 7.09798}. NSolve and Reduce, on the other hand, appear not to work here. $\endgroup$
    – bbgodfrey
    Aug 15, 2015 at 19:07
  • $\begingroup$ Finding multiple roots is discussed here. $\endgroup$
    – bbgodfrey
    Aug 15, 2015 at 19:13

1 Answer 1

6
$\begingroup$

Finding zeroes of the function

Apply[ArcTan, o1[t] - o[t]] - Apply[ArcTan, o2[t] - o[t]]

does not work well, because ArcTan changes discontinuously from -Pi to Pi, creating numerous spurious roots. For instance, with

{a1 = 2, a2 = 13/10, a3 = 16/10, b1 = 3, b2 = 21/10, b3 = 26/10};
Plot[Apply[ArcTan, o1[t] - o[t]] - Apply[ArcTan, o2[t] - o[t]], {t, 0, 20}]

enter image description here

A better choice is to equate slopes, even though this function does not distinguish between planets on the same side and on opposite sides of the reference planet.

Plot[First[o1[t] - o[t]] Last[o2[t] - o[t]] - First[o2[t] - o[t]] Last[o1[t] - o[t]], 
    {t, 0, 20}]

enter image description here

We now apply FindAllCrossings

FindAllCrossings[First[o1[t] - o[t]] Last[o2[t] - o[t]] - 
    First[o2[t] - o[t]] Last[o1[t] - o[t]], {t, 0, 20}]
(* {3.09583, 6.54554, 9.66901, 14.3135, 19.458} *)

and retain only those solutions for which the planets are on the same side of the reference planet.

Cases[%, t_ /; First[o1[t] - o[t]] First[o2[t] - o[t]] > 0]
(* {6.54554, 14.3135} *)

For completeness, alignment can be tested visually from

Manipulate[Graphics[{Arrow[{o[t], o1[t]}], Arrow[{o[t], o2[t]}]}, 
    PlotRange -> {{-3, 3}, {-3, 3}}], {{t, 1}, 0, 20, Appearance -> "Labeled"}]

The first alignment occurs at about t == 6.55.

enter image description here

Simpler approach

Although NSolve cannot handle the function involving ArcTan introduced in the Question, it can handle the simpler expression introduced above in this Answer.

NSolve[First[o1[t] - o[t]] Last[o2[t] - o[t]] - First[o2[t] - o[t]] 
   Last[o1[t] - o[t]] == 0 && First[o1[t] - o[t]] First[o2[t] - o[t]] > 0 && 0 < t < 20, t]
(* {{t -> 6.54554}, {t -> 14.3135}} *)

Reduce also works.

N@Reduce[First[o1[t] - o[t]] Last[o2[t] - o[t]] - First[o2[t] - o[t]] 
   Last[o1[t] - o[t]] == 0 && First[o1[t] - o[t]] First[o2[t] - o[t]] > 0 && 0 < t < 20, t]
(* t == 6.54554 || t == 14.3135 *)
$\endgroup$
1
  • 3
    $\begingroup$ Took me a second to figure out "A better choice is to equate slopes". The slope of a vector (vx,vy) is vy/vx so equating vy/vx=uy/ux for two vectors yields vy*ux = uy*vx, or vy*ux-uy*vx==0 $\endgroup$
    – user1722
    Aug 16, 2015 at 0:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.