Python not generating multiple polar curves

Question

I am trying to produce one plot that contains multiple polar curves; however, only one polar curve is generating. This is my first time using polar plots in matplotlib, but I assumed generating multiple curves on the same plot would work the same as generating multiple lines on a regular graph. I have seen examples of code that demonstrate exactly what I need, and even though my code looks the same, only one curve generates. The curves are also not generating one at a time in separate plots. Only one plot generates.

I have made sure that plt.show() is outside of any while or for loops, and I have even tried moving the plotting functions to an entirely new method with no success. I am unsure of what I am doing wrong and would appreciate any help I could get.

Thank you.

CODE: (Plotting function towards the bottom)

import math
import numpy as np
import matplotlib.pyplot as plt


def rFuncA(a, ec, v):
    # calculates current orbit radius from a, eccentricity, and angle
    return a * (1 - ec ** 2) / (1 + ec * math.cos(v))


def wFunc(a, ec, v, u):
    # calculates current orbit radius from a, eccentricity, and angle
    block1 = a * (1 - ec ** 2)
    block2 = 1 + ec * math.cos(v)
    return (math.sqrt(u) / block1 ** (3 / 2)) * block2 ** 2


def VoFunc(ec, v, r, w):
    return r * w * (ec * math.sin(v)) / (1 + ec * math.cos(v)), \
           r * w


def dXFunc(Fl, Fd, Ft, gamma, alpha, u, r, m, Isp, Vr, Vv, g):
    theta = alpha + gamma
    dX1 = (Vv ** 2) / r + (Fl * math.cos(gamma) - Fd * math.sin(gamma) + Ft * math.sin(theta)) / m - u / (r ** 2)
    dX2 = -(Vr * Vv / r + (Fl * math.sin(gamma) + Fd * math.cos(gamma) - Ft * math.cos(theta)) / m)
    dX3 = Vr
    dX4 = Vv / r
    dX5 = -Ft / (g * Isp)

    return dX1, dX2, dX3, dX4, dX5


def aFuncR(r, u, Vr, Vv):
    block1 = 2 * u / r
    block2 = Vr ** 2 + Vv ** 2

    return u / (block1 - block2)


def ecFuncR(r, u,  Vr, Vv):

    return (r / u) * math.sqrt((Vv ** 2 - u / r) ** 2 + (Vv * Vr) ** 2)


class Properties:
    def __init__(self):
        # Preallocate
        n = 1000000
        self.ec = np.zeros([n])
        self.a = np.zeros([n])
        self.r = np.zeros([n])
        self.gamma = np.zeros([n])
        self.X = np.zeros([n, 5])
        self.Tt = np.zeros([n])
        self.v = np.zeros([n])
        self.theta = np.zeros([n])

        # Constants
        self.Ft = 0.01  # Thrust Force, kN
        self.Fl = 0  # Lift Force, kN
        self.Fd = 0  # Drag Force, kN
        self.g = 0.009806  # gravity, km/s^2
        self.u = 3.986E5  # gravitational parameter, km^3/s^2
        self.alpha = 0  # Angle of Attack (0 for ballistic trajectory), radians
        self.Isp = 2000  # Specific impulse, sec
        self.dt = 2000  # Iterator step size, sec
        self.r_e = 6561  # Mean earth radius, km
        self.r_a = 13200  # Apogee radius at final orbit, km

        # Progressive variables initial values
        self.a[0] = 8530
        self.r[0] = rFuncA(self.a[0], self.ec[0], self.v[0])    # Changing radius, km
        self.m = np.zeros([n])
        self.m[0] = 1000 + 7.743 + 71.25    # Inital Mass
        self.w = np.zeros([n])
        self.w[0] = wFunc(self.a[0], self.ec[0], self.v[0], self.u)     # Angular velocity, rad/s

        Vo1, Vo2 = VoFunc(self.ec[0], self.v[0], self.r[0], self.w[0])
        self.Vr = np.zeros([n])     # Creates array for radial velocities
        self.Vv = np.zeros(([n]))   # Creates array for tangential velocities
        self.Vr[0] = Vo1    # Initial radial velocity, km/s
        self.Vv[0] = Vo2    # Initial tangential velocity, km/s

        # Plotting
        self.pt = np.zeros([n])
        self.pr = np.zeros([n])
        self.cr = np.zeros([n])
        self.cf = np.zeros([n])
        self.i = 0

    def runge_kutta(self):
        self.i = 0
        self.gamma[0] = math.atan(self.Vr[0] / self.Vv[0])  # Initial angle, rad
        self.X[0][:] = np.array([self.Vr[0], self.Vv[0], self.r[0], self.v[0], self.m[0]])
        while self.a[self.i] * (1 + self.ec[self.i]) < self.r_a:
            # self.Tt[i+1] = self.Tt[0] + self.dt

            dX1, dX2, dX3, dX4, dX5 = dXFunc(self.Fl, self.Fd, self.Ft, self.gamma[self.i], self.alpha, self.u, self.X[self.i][2],
                                             self.X[self.i][4], self.Isp, self.X[self.i][0], self.X[self.i][1], self.g)
            k1 = np.array([dX1, dX2, dX3, dX4, dX5])
            X2 = self.X[self.i] + k1 * self.dt / 2
            dX1, dX2, dX3, dX4, dX5 = dXFunc(self.Fl, self.Fd, self.Ft, self.gamma[self.i], self.alpha, self.u, X2[2],
                                             X2[4], self.Isp, X2[0], X2[1], self.g)
            k2 = np.array([dX1, dX2, dX3, dX4, dX5])
            X3 = X2 + k2 * self.dt / 2
            dX1, dX2, dX3, dX4, dX5 = dXFunc(self.Fl, self.Fd, self.Ft, self.gamma[self.i], self.alpha, self.u, X3[2],
                                             X3[4], self.Isp, X3[0], X3[1], self.g)
            k3 = np.array([dX1, dX2, dX3, dX4, dX5])
            X4 = X3 + k3 * self.dt
            dX1, dX2, dX3, dX4, dX5 = dXFunc(self.Fl, self.Fd, self.Ft, self.gamma[self.i], self.alpha, self.u, X4[2],
                                             X4[4], self.Isp, X4[0], X4[1], self.g)
            k4 = np.array([dX1, dX2, dX3, dX4, dX5])

            self.X[self.i+1][:] = self.X[self.i][:] + (1 / 6) * (k1 + k2 + k3 + k4 * self.dt)

            self.a[self.i+1] = aFuncR(self.X[self.i+1][2], self.u, self.X[self.i+1][0], self.X[self.i+1][1])
            self.ec[self.i+1] = ecFuncR(self.X[self.i+1][2], self.u, self.X[self.i+1][0], self.X[self.i+1][1])

            self.i += 1

        print(self.X[self.i][:])    # Final Values for Vr, Vv, r, v, and m

    def plot(self):
        i = self.i - 1
        k = 0

        # Orbit period
        while self.pt[k] < 2 * math.pi:
            self.pt[k+1] = self.pt[k] + (2 * math.pi / i)
            k += 1

        # Orbit location
        for j in range(i):
            # Transferring Orbit
            self.pr[j] = self.a[j] * (1 + self.ec[j] ** 2) / (1 + self.ec[j] * math.cos(self.pt[j]))
            # Low Earth Orbit
            self.cr[j] = self.a[0]

        plt.polar(self.pt[:i], self.pr[:i], 'r')
        plt.polar(self.pt[:i], self.cr[:i], 'g',)
        plt.show()


def main():

    P = Properties()
    P.runge_kutta()
    P.plot()


if __name__ == '__main__':
    main()

Plot: First polar curve

Second polar curve (only generates when first plot is commented out)