[FIXED] Using arctan / arctan2 to plot a from 0 to 2π

Issue

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.

Textbook plot and equation:

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

Solution

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()

Answered By - Saullo G. P. Castro