¿Cómo calcular el tiempo de apoapsis y periapsis, dados los elementos orbitales?

Tengo los 6 elementos orbitales (excluyendo el tiempo de paso del periapsis, obviamente), 2 vectores de estado y muchos otros parámetros como el período orbital. Estoy buscando usarlos para calcular el tiempo de periapsis y el tiempo de apoapsis.

Para darle una mejor idea de los valores que ya tengo, aquí está mi código calculando todos los datos relacionados con la órbita:

 void CalculateOrbitalData(Vector3 pos1, Vector3 vel1, float m1) {
    Vector3 pos2 = planet.position;
    Vector3 vel2 = planet.velocity;
    float m2 = planet.mass;

    //Relative Position Vector
    Vector3 r = pos1 - pos2;

    //Distance between bodies
    float rmag = r.magnitude;

    //Relative Velocity Vector
    Vector3 v = vel1 - vel2;

    //Specific Angular Momentum
    Vector3 h = Vector3.Cross(r, v);

    //Standard Gravitational Parameter
    float µ = World.G * m2;

    //Eccentricity Vector
    Vector3 evec = (Vector3.Cross(v, h) / µ) - (r / Vector3.Magnitude(r));

    //Eccentricity
    float e = Vector3.Magnitude(evec);

    //Vector to Ascending Node
    Vector3 n = new Vector3(-h.x, h.z, 0);

    //True Anomaly
    float t = Mathf.Acos((Vector3.Dot(evec, r)) / (e * Vector3.Magnitude(r)));
    if (Vector3.Dot(r, v) < 0)
        t = (2 * Mathf.PI) - t;

    //Longitude of Ascending Node (2D)
    float Ω = 0;

    //Inclination (2D)
    float i = 0;

    //Argument of Periapsis
    float ω = Mathf.Atan2(evec.y, evec.x);
    float ωdegrees = ω * (180 / Mathf.PI);
    if (e == 0) {
        ω = 0;
        ωdegrees = 0;
    }

    //Eccentric Anomaly
    float E = 2 * Mathf.Atan(Mathf.Tan(t / 2) / Mathf.Sqrt((1 + e) / (1 - e)));

    //Mean Anomaly 
    float M = E - (e * Mathf.Sin(E));

    //Semi-Major Axis
    float a = 1 / ((2 / Vector3.Magnitude(r)) - (Mathf.Pow(Vector3.Magnitude(v), 2) / µ));

    //Apoapsis
    float ap = a * (1 + e);

    //Periapsis
    float pe = a * (1 - e);

    //Orbital Period
    float T = (2 * Mathf.PI) * Mathf.Sqrt(Mathf.Pow(a, 3) / µ) / 60;
    //patch orbital period to stay in line with unity timestep
    float fT = T * 1.205f;

    //Mean motion
    float m = (2 * Mathf.PI) / fT;

    //Perifocal distance
    float rp = (Vector3.Dot(h, h) / µ) / e + 1;
}