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."

13
  • 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
    Commented May 17, 2013 at 17:02
  • 2
    great! note that in the equation from the textbook, the arctan is only over the first term, not the whole equation.
    – tom10
    Commented May 17, 2013 at 17:27
  • 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
    Commented May 17, 2013 at 17:30
  • 1
    I agree with @tom10, it seems like a case of a misplaced closing paren. Or two.
    – user629132
    Commented May 17, 2013 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
    Commented May 17, 2013 at 23:40

2 Answers 2

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

11
  • 2
    Using arctan2(1, 1/x) instead of arctan(x) here is just ridiculous.
    – tom10
    Commented May 17, 2013 at 23:16
  • @tom10 I see your point, originally I suggested to use arctan from -pi to pi, but with arctan2(1,1/x) it gave the same graphic as the text book. Please, feel free to apply any corrections using arctan and getting the same results Commented May 17, 2013 at 23:25
  • @dustin Could you explain from were you took arctan(x)=arctan2(1,1/x)? Because the common practice is to use arctan(x) and sum 2pi to the negative results Commented May 18, 2013 at 8:06
  • @SaulloCastro a physicists suggested 1, 1 / x usage in arctan2
    – dustin
    Commented May 18, 2013 at 9:18
  • @dustin It would be nice to include a reference to this in the question or in the answer Commented May 18, 2013 at 9:23
1

Given x and y coordinates, the following code will evaluate correctly the angle (assuming the angle is measured from the x-axis and counter-clockwise) between the range 0 to 2pi. It just sums 2pi if the output of arctan2 is negative in a simple and effective manner.

import numpy as np
np.arctan2(y,x) + -1*(np.sign(np.arctan2(y,x))-1)*np.pi
1

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.