9

I am trying to replicate a plot in Orbital Mechanics by Curtis, but I just can't quite get it. However, I have made head way by switching to np.arctan2 from np.arctan.

Maybe I am implementing arctan2 incorrectly?

import pylab
import numpy as np

e = np.arange(0.0, 1.0, 0.15).reshape(-1, 1)

nu = np.linspace(0.001, 2 * np.pi - 0.001, 50000)
M2evals = (2 * np.arctan2(1, 1 / (((1 - e) / (1 + e)) ** 0.5 * np.tan(nu / 2) -
           e * (1 - e ** 2) ** 0.5 * np.sin(nu) / (1 + e * np.cos(nu)))))

fig2 = pylab.figure()
ax2 = fig2.add_subplot(111)

for Me2, _e in zip(M2evals, e.ravel()):
    ax2.plot(nu.ravel(), Me2, label = str(_e))

pylab.legend()
pylab.xlim((0, 7.75))
pylab.ylim((0, 2 * np.pi))
pylab.show()

In the image below, there are discontinuities popping up. The function is supposed to be smooth and connect at 0 and 2 pi in the y range of (0, 2pi) not touching 0 and 2pi.

Enter image description here

Textbook plot and equation:

Enter image description here

Enter image description here

At the request of Saullo Castro, I was told that:

The problem may lie in the arctan function which gives "principle values" as output.

Thus, arctan(tan(x)) does not yield x if x is an angle in the second or third quadrant. If you plot arctan(tan(x)) from x = 0 to x = Pi, you will find that it has a discontinuous jump at x = Pi/2.

For your case, instead of writing arctan(arg), I believe you would write arctan2(1, 1/arg) where arg is the argument of your arctan function. That way, when arg becomes negative, arctan2 will yield an angle in the second quadrant rather than the fourth."

Peter Mortensen
  • 30,738
  • 21
  • 105
  • 131
dustin
  • 4,309
  • 12
  • 57
  • 79
  • 3
    To use arctan2 correctly, you need an equation for x, and an equation for y. The whole point of arctan2 is that the ratio y/x is ambiguous about the quadrant: -/- == +/+ and -/+ == +/-. If you put in y as a constant, you can still only occupy two quadrants in the result. So what you're saying doesn't make sense: this can't both be the equation and fill the quadrants you say. (Since all you're saying now is "I have this equation and it doesn't work", we don't have enough info to answer this and it's not a question that can be answered.) – tom10 May 17 '13 at 17:02
  • @tom10 it very close to working if we take out .6 - .9, the solution is correct. It is breaking down above .6 but maybe lower since the range is by .15 – dustin May 17 '13 at 17:04
  • Still, you're not giving people enough information, I'd like to know where you going, not where you started. "Close to correct" is not particularly helpful. This doesn't seem to be a question about acrtan2, but a question about what's in your textbook. – tom10 May 17 '13 at 17:07
  • @tom10 adding pics of the equation and plot then. – dustin May 17 '13 at 17:11
  • 2
    great! note that in the equation from the textbook, the arctan is only over the first term, not the whole equation. – tom10 May 17 '13 at 17:27
  • @tom10 what do you mean by that? – dustin May 17 '13 at 17:28
  • 1
    Look at the parenthesis directly right of tan^-1, where's the matching parenthesis. It's after the tan(theta/2), not at the end of the whole equation. In your equation, you're doing the arctan of the whole thing. – tom10 May 17 '13 at 17:30
  • 1
    I agree with @tom10, it seems like a case of a misplaced closing paren. Or two. –  May 17 '13 at 20:54
  • 1
    This question should just be closed, imho. It's just a parenthesis typo on the part of the OP, and sheds no light on anything else (especially not arctan vs arctan2 due to the misdirection). Therefore, closed as "too localized". – tom10 May 17 '13 at 23:40
  • @tom10 well I did learn how to moved the legend from it so it was helpful. – dustin May 17 '13 at 23:47
  • @tom10 The question should not be closed. It brings a very common problem into discussion and it suggests to use arctan2() to avoid summing 2*pi for the negatives obtained with arctan() – Saullo G. P. Castro May 18 '13 at 08:31
  • 1
    @SaultoCastro: it only brings that in tangentially (ha ha). If you want that to be what this question is about, then I suggest you answer it, with pictures of the domains, explain why arctan2(1,1/x) works, which is better, etc. Instead you're just asking others to answer that, while you've actually answered what the question is really about... making a figure that exactly matches the one in the OP's text book (ie, it's way too localized). (btw, I'm outa here... this has been a total waste of time.) – tom10 May 18 '13 at 17:25
  • This is too localized as the actual problem (and correct answer) involves neither `python`, `numpy`, nor `matplotlib`. – tacaswell May 18 '13 at 18:24

1 Answers1

12

The common practice is to sum 2pi in the negative results of arctan(), which can be done efficiently. The OP's suggestion to replace arctan(x) by arctan2(1, 1/x), also suggested by Maple 15's documentation as pointed out by Yay295, produces the same results without the need to sum 2pi. Both are shown below:

import pylab
import numpy as np
e = np.arange(0.0, 1.0, 0.15).reshape(-1, 1)
nu = np.linspace(0, 2*np.pi, 50000)
x =  ((1-e)/(1+e))**0.5 * np.tan(nu/2.)
x2 = e*(1-e**2)**0.5 * np.sin(nu)/(1 + e*np.cos(nu))
using_arctan = True
using_OP_arctan2 = False

if using_arctan:
    M2evals = 2*np.arctan(x) - x2
    M2evals[M2evals<0] += 2*np.pi
elif using_OP_arctan2:
    M2evals = 2 * np.arctan2(1,1/x) - x2

fig2 = pylab.figure()
ax2 = fig2.add_subplot(111)
for M2e, _e in zip(M2evals, e.ravel()):
    ax2.plot(nu.ravel(), M2e, label = str(_e))
pylab.legend(loc='upper left')
pylab.show()

Enter image description here

Peter Mortensen
  • 30,738
  • 21
  • 105
  • 131
Saullo G. P. Castro
  • 56,802
  • 26
  • 179
  • 234