C# Orbit.getRelativePositionAtUT方法代码示例

        double tauP; //normalized parabolic transfer time

        #endregion Fields

        #region Constructors

        /// <summary>
        /// Initializes a new instance of the <see cref="ThrottleControlledAvionics.LambertSolver"/> class.
        /// </summary>
        /// <param name="orb">Starting orbit.</param>
        /// <param name="destination">Destination radius-vector.</param>
        /// <param name="UT">Starting UT.</param>
        public LambertSolver(Orbit orb, Vector3d destination, double UT)
        {
            orbit = orb;
            body = orbit.referenceBody;
            StartUT = UT;
            mu = body.gravParameter;

            r1 = orb.getRelativePositionAtUT(UT);
            var h  = Vector3d.Cross(r1, destination);
            if(h.sqrMagnitude < 0.01) h = orb.GetOrbitNormal();
            c  = destination-r1;

            cm  = c.magnitude;
            var r1m = r1.magnitude;
            var r2m = destination.magnitude;
            var rrm = r1m+r2m;
            m  = rrm+cm;
            n  = rrm-cm;
            m3 = m*m*m;

            var transfer_angle = Vector3d.Angle(r1, destination)*Mathf.Deg2Rad;
            if(h.z < 0) transfer_angle = Utils.TwoPI-transfer_angle;
            sigma = Math.Sqrt(n/m);
            if(transfer_angle > Math.PI) sigma = -sigma;
            sigma2 = sigma*sigma;
            sigma3 = sigma2*sigma;
            sigma5 = sigma2*sigma3;

            tauP = 2/3.0*(1-sigma3);
            tauME = Math.Acos(sigma)+sigma*Math.Sqrt(1-sigma2);
        }

        static void WarpShip(Vessel vessel, Orbit newOrbit)
        {
            if (newOrbit.getRelativePositionAtUT (Planetarium.GetUniversalTime ()).magnitude > newOrbit.referenceBody.sphereOfInfluence)
                throw new ArgumentException ("Destination position was above the sphere of influence");

            vessel.Landed = false;
            vessel.Splashed = false;
            vessel.landedAt = string.Empty;
            var parts = vessel.parts;
            if (parts != null) {
                var clamps = parts.Where (p => p.Modules != null && p.Modules.OfType<LaunchClamp> ().Any ()).ToList ();
                foreach (var clamp in clamps)
                    clamp.Die ();
            }

            try {
                OrbitPhysicsManager.HoldVesselUnpack (60);
            } catch (NullReferenceException) {
            }

            foreach (var v in (FlightGlobals.fetch == null ? (IEnumerable<Vessel>)new[] { vessel } : FlightGlobals.Vessels).Where(v => !v.packed))
                v.GoOnRails ();

            HardsetOrbit (vessel.orbit, newOrbit);

            vessel.orbitDriver.pos = vessel.orbit.pos.xzy;
            vessel.orbitDriver.vel = vessel.orbit.vel;
        }

        public static void SetOrbit(this Vessel vessel, Orbit newOrbit)
        {
            var destinationMagnitude = newOrbit.getRelativePositionAtUT(Planetarium.GetUniversalTime()).magnitude;
            if (destinationMagnitude > newOrbit.referenceBody.sphereOfInfluence)
            {
                View.WindowHelper.Error("Destination position was above the sphere of influence");
                return;
            }
            if (destinationMagnitude < newOrbit.referenceBody.Radius)
            {
                View.WindowHelper.Error("Destination position was below the surface");
                return;
            }

            vessel.PrepVesselTeleport();

            try
            {
                OrbitPhysicsManager.HoldVesselUnpack(60);
            }
            catch (NullReferenceException)
            {
                Extensions.Log("OrbitPhysicsManager.HoldVesselUnpack threw NullReferenceException");
            }

            var allVessels = FlightGlobals.fetch?.vessels ?? (IEnumerable<Vessel>)new[] { vessel };
            foreach (var v in allVessels)
                v.GoOnRails();

            var oldBody = vessel.orbitDriver.orbit.referenceBody;

            HardsetOrbit(vessel.orbitDriver, newOrbit);

            vessel.orbitDriver.pos = vessel.orbit.pos.xzy;
            vessel.orbitDriver.vel = vessel.orbit.vel;

            var newBody = vessel.orbitDriver.orbit.referenceBody;
            if (newBody != oldBody)
            {
                var evnt = new GameEvents.HostedFromToAction<Vessel, CelestialBody>(vessel, oldBody, newBody);
                GameEvents.onVesselSOIChanged.Fire(evnt);
            }
        }

        private Double PhaseAngleCalc(Orbit o1, Orbit o2, Double UT)
        {
            Vector3d n = o1.GetOrbitNormal();

            Vector3d p1 = o1.getRelativePositionAtUT(UT);
            Vector3d p2 = o2.getRelativePositionAtUT(UT);
            double phaseAngle = Vector3d.Angle(p1, p2);
            if (Vector3d.Angle(Quaternion.AngleAxis(90, Vector3d.forward) * p1, p2) > 90)
            {
                phaseAngle = 360 - phaseAngle;
            }

            if (o2.semiMajorAxis < o1.semiMajorAxis)
            {
                phaseAngle = phaseAngle - 360;
            }

            return LambertSolver.Deg2Rad * phaseAngle;
            //return LambertSolver.Deg2Rad * ((phaseAngle + 360) % 360);
        }

        private void initMeshFromOrbit(Vector3 bodyPosition, Mesh mesh, Orbit orbit, double startTime, double duration, Color color)
        {
            int steps = 128;

            Vector3 camPos = ScaledSpace.ScaledToLocalSpace(MapView.MapCamera.transform.position) - bodyPosition;

            double prevTA = orbit.TrueAnomalyAtUT(startTime);
            double prevTime = startTime;

            double[] stepUT = new double[steps * 4];
            int utIdx = 0;
            double maxDT = Math.Max(1.0, duration / (double)steps);
            double maxDTA = 2.0 * Math.PI / (double)steps;
            stepUT[utIdx++] = startTime;
            while(true)
            {
                double time = prevTime + maxDT;
                for (int count = 0; count < 100; ++count)
                {
                    if (count == 99)
                        Debug.Log("WARNING: infinite loop? (Trajectories.MapOverlay.initMeshFromOrbit)");
                    double ta = orbit.TrueAnomalyAtUT(time);
                    while (ta < prevTA)
                        ta += 2.0 * Math.PI;
                    if (ta - prevTA <= maxDTA)
                    {
                        prevTA = ta;
                        break;
                    }
                    time = (prevTime + time) * 0.5;
                }

                if (time > startTime + duration - (time-prevTime) * 0.5)
                    break;

                prevTime = time;

                stepUT[utIdx++] = time;
                if (utIdx >= stepUT.Length - 1)
                {
                    //Util.PostSingleScreenMessage("ut overflow", "ut overflow");
                    break; // this should never happen, but better stop than overflow if it does
                }
            }
            stepUT[utIdx++] = startTime + duration;

            var vertices = new Vector3[utIdx * 2 + 2];
            var triangles = new int[utIdx * 6];

            Vector3 prevMeshPos = Util.SwapYZ(orbit.getRelativePositionAtUT(startTime - duration / (double)steps));
            for(int i = 0; i < utIdx; ++i)
            {
                double time = stepUT[i];

                Vector3 curMeshPos = Util.SwapYZ(orbit.getRelativePositionAtUT(time));
                if (Settings.fetch.BodyFixedMode) {
                    curMeshPos = Trajectory.calculateRotatedPosition(orbit.referenceBody, curMeshPos, time);
                }

                // compute an "up" vector that is orthogonal to the trajectory orientation and to the camera vector (used to correctly orient quads to always face the camera)
                Vector3 up = Vector3.Cross(curMeshPos - prevMeshPos, camPos - curMeshPos).normalized * (lineWidth * Vector3.Distance(camPos, curMeshPos));

                // add a segment to the trajectory mesh
                vertices[i * 2 + 0] = curMeshPos - up;
                vertices[i * 2 + 1] = curMeshPos + up;

                if (i > 0)
                {
                    int idx = (i - 1) * 6;
                    triangles[idx + 0] = (i - 1) * 2 + 0;
                    triangles[idx + 1] = (i - 1) * 2 + 1;
                    triangles[idx + 2] = i * 2 + 1;

                    triangles[idx + 3] = (i - 1) * 2 + 0;
                    triangles[idx + 4] = i * 2 + 1;
                    triangles[idx + 5] = i * 2 + 0;
                }

                prevMeshPos = curMeshPos;
            }

            var colors = new Color[vertices.Length];
            for (int i = 0; i < colors.Length; ++i)
            {
                //if (color.g < 0.5)
                    colors[i] = color;
                /*else
                    colors[i] = new Color(0, (float)i / (float)colors.Length, 1.0f - (float)i / (float)colors.Length);*/
            }

            for (int i = 0; i < vertices.Length; ++i)
            {
                vertices[i] = ScaledSpace.LocalToScaledSpace(vertices[i] + bodyPosition);
            }

            mesh.Clear();
            mesh.vertices = vertices;
            mesh.colors = colors;
            mesh.triangles = triangles;
            mesh.RecalculateBounds();
//.........这里部分代码省略.........

		public static ManeuverParameters ComputeEjectionManeuver(Vector3d exit_velocity, Orbit initial_orbit, double UT_0)
		{
			double GM = initial_orbit.referenceBody.gravParameter;
			double C3 = exit_velocity.sqrMagnitude;
            double Mh = initial_orbit.referenceBody.sphereOfInfluence;
            double Rpe = initial_orbit.semiMajorAxis;

            double isma = 2 / Mh - C3 / GM; //inverted Semi-major Axis, will work for parabolic orbit
            double ecc = 1.0 - Rpe * isma;

            //double vstart = Math.Sqrt(GM * (2 / Rpe - isma)); //{ total start boost}
            //double slat = Rpe * Rpe * vstart * vstart / GM;

            //the problem here should be R for circular orbit instead of Rpe
            double slat = 2 * Rpe - isma * Rpe * Rpe; //but we don't know start point (yet) in elliptic orbit

            double theta = Math.Acos((slat / Mh - 1) / ecc);


            /*
            //old way infinity hyperbolic:
            //problems: it's not necessary hyperbolic (in case of low speed exit_velocity), 
            //and exit_velocity appears not infinite far from celestial body, but only sphereOfInfluence far
            //i.e. Mh in previous statements(theta, isma) is not infinity!

			double sma = -GM / C3;
			
			double ecc = 1 - Rpe / sma;
			double theta = Math.Acos(-1 / ecc);

            */
            Vector3d intersect_1;
			Vector3d intersect_2;

			// intersect_1 is the optimal position on the orbit assuming the orbit is circular and the sphere of influence is infinite
			IntersectConePlane(exit_velocity, theta, initial_orbit.h, out intersect_1, out intersect_2);

			Vector3d eccvec = initial_orbit.GetEccVector();

			double true_anomaly = AngleAboutAxis(eccvec, intersect_1, initial_orbit.h);
			// GetMeanAnomalyAtTrueAnomaly expects degrees and returns radians
			double mean_anomaly = initial_orbit.GetMeanAnomalyAtTrueAnomaly(true_anomaly * 180 / Math.PI);

			double delta_mean_anomaly = MuUtils.ClampRadiansPi(mean_anomaly - initial_orbit.MeanAnomalyAtUT(UT_0));

			double UT = UT_0 + delta_mean_anomaly / initial_orbit.MeanMotion();

			Vector3d V_0 = initial_orbit.getOrbitalVelocityAtUT(UT);
			Vector3d pos = initial_orbit.getRelativePositionAtUT(UT);

			double final_vel = Math.Sqrt(C3 + 2 * GM / pos.magnitude);

			Vector3d x = pos.normalized;
			Vector3d z = Vector3d.Cross(x, exit_velocity).normalized;
			Vector3d y = Vector3d.Cross(z, x);

			theta = Math.PI / 2;
			double dtheta = 0.001;
			Orbit sample = new Orbit();

			double theta_err = Double.MaxValue;

			for (int iteration = 0 ; iteration < 50 ; ++iteration)
			{
				Vector3d V_1 = final_vel * (Math.Cos(theta) * x + Math.Sin(theta) * y);
				sample.UpdateFromStateVectors(pos, V_1, initial_orbit.referenceBody, UT);
				theta_err = AngleAboutAxis(exit_velocity, sample.getOrbitalVelocityAtUT(OrbitExtensions.NextTimeOfRadius(sample, UT, sample.referenceBody.sphereOfInfluence)), z);
				if (double.IsNaN(theta_err))
					return null;
				if (Math.Abs(theta_err) <= Math.PI / 1800)
					return new ManeuverParameters((V_1 - V_0).xzy, UT);

				V_1 = final_vel * (Math.Cos(theta + dtheta) * x + Math.Sin(theta + dtheta) * y);
				sample.UpdateFromStateVectors(pos, V_1, initial_orbit.referenceBody, UT);
				double theta_err_2 = AngleAboutAxis(exit_velocity, sample.getOrbitalVelocityAtUT(OrbitExtensions.NextTimeOfRadius(sample, UT, sample.referenceBody.sphereOfInfluence)), z);

				V_1 = final_vel * (Math.Cos(theta - dtheta) * x + Math.Sin(theta - dtheta) * y);
				sample.UpdateFromStateVectors(pos, V_1, initial_orbit.referenceBody, UT);
				double theta_err_3 = AngleAboutAxis(exit_velocity, sample.getOrbitalVelocityAtUT(OrbitExtensions.NextTimeOfRadius(sample, UT, sample.referenceBody.sphereOfInfluence)), z);

				double derr = MuUtils.ClampRadiansPi(theta_err_2 - theta_err_3) / (2 * dtheta);

				theta = MuUtils.ClampRadiansTwoPi(theta - theta_err / derr);

				// if theta > pi, replace with 2pi - theta, the function ejection_angle=f(theta) is symmetrc wrt theta=pi
				if (theta > Math.PI)
					theta = 2 * Math.PI - theta;
			}

			return null;
		}

        //enum OrbitAltitudeDirection
        //{
        //    Ascending=1,
        //    Descending=2
        //}
        internal static double timeOfTargetAltitude(Orbit oOrig, double timeStart, double numDivisions, out double closestdistance, double targetAltitude)
        {
            double closestApproachTime = timeStart;
            double closestApproachAltitude = Double.MaxValue;
            double minTime = timeStart;

            //work out the end of the period to scan based on whether we are above or below the target altitude
            double UTofEndOrbitAltDirection;
            double startAltitude = oOrig.getRelativePositionAtUT(timeStart).magnitude-oOrig.referenceBody.Radius;
            if (startAltitude < targetAltitude)
                UTofEndOrbitAltDirection = oOrig.StartUT + oOrig.timeToAp;
            else
                UTofEndOrbitAltDirection = oOrig.StartUT + oOrig.timeToPe;

            double periodtoscan = UTofEndOrbitAltDirection - timeStart;
            double maxTime = UTofEndOrbitAltDirection;

            //work out iterations for precision - we only really need to within a second - so how many iterations do we actually need
            //Each iteration gets us 1/10th of the period to scan

            for (int iter = 0; iter < 8; iter++)
            {
                double dt = (maxTime - minTime) / numDivisions;
                for (int i = 0; i < numDivisions; i++)
                {
                    double t = minTime + i * dt;
                    double distance = oOrig.getRelativePositionAtUT(t).magnitude-oOrig.referenceBody.Radius;

                    //TODO: NEED TO DO SOMETHING HERE WITH PRECISION!!!! How many decimal place to compare to
                    //something about = if its not substantially closer, but is substantialy later then ignore it.
                    //ie if the later crossing is 1meter closer, but its 600 secs later then keep the earlier one
                    if (Math.Abs(distance - targetAltitude) < closestApproachAltitude)
                    {
                        closestApproachAltitude = Math.Abs(distance - targetAltitude);
                        closestApproachTime = t;
                    }
                }
                minTime = KACUtils.Clamp(closestApproachTime - dt, timeStart, timeStart + periodtoscan);
                maxTime = KACUtils.Clamp(closestApproachTime + dt, timeStart, timeStart + periodtoscan);
            }

            closestdistance = closestApproachAltitude + targetAltitude;
            return closestApproachTime;
        }

        private void initMeshFromOrbit(Vector3 bodyPosition, Mesh mesh, Orbit orbit, double startTime, double duration, Color color)
        {
            int steps = 128;

            double prevTA = orbit.TrueAnomalyAtUT(startTime);
            double prevTime = startTime;

            double[] stepUT = new double[steps * 4];
            int utIdx = 0;
            double maxDT = Math.Max(1.0, duration / (double)steps);
            double maxDTA = 2.0 * Math.PI / (double)steps;
            stepUT[utIdx++] = startTime;
            while(true)
            {
                double time = prevTime + maxDT;
                for (int count = 0; count < 100; ++count)
                {
                    if (count == 99)
                        Debug.Log("WARNING: infinite loop? (Trajectories.MapOverlay.initMeshFromOrbit)");
                    double ta = orbit.TrueAnomalyAtUT(time);
                    while (ta < prevTA)
                        ta += 2.0 * Math.PI;
                    if (ta - prevTA <= maxDTA)
                    {
                        prevTA = ta;
                        break;
                    }
                    time = (prevTime + time) * 0.5;
                }

                if (time > startTime + duration - (time-prevTime) * 0.5)
                    break;

                prevTime = time;

                stepUT[utIdx++] = time;
                if (utIdx >= stepUT.Length - 1)
                {
                    //Util.PostSingleScreenMessage("ut overflow", "ut overflow");
                    break; // this should never happen, but better stop than overflow if it does
                }
            }
            stepUT[utIdx++] = startTime + duration;

            var vertices = new Vector3[utIdx * 2 + 2];
            var uvs = new Vector2[utIdx * 2 + 2];
            var triangles = new int[utIdx * 6];

            Vector3 prevMeshPos = Util.SwapYZ(orbit.getRelativePositionAtUT(startTime - duration / (double)steps)) + bodyPosition;
            for(int i = 0; i < utIdx; ++i)
            {
                double time = stepUT[i];

                Vector3 curMeshPos = Util.SwapYZ(orbit.getRelativePositionAtUT(time));
                if (Settings.fetch.BodyFixedMode) {
                    curMeshPos = Trajectory.calculateRotatedPosition(orbit.referenceBody, curMeshPos, time);
                }
                curMeshPos += bodyPosition;

                // add a segment to the trajectory mesh
                MakeRibbonEdge(prevMeshPos, curMeshPos, lineWidth, vertices, i * 2);
                uvs[i * 2 + 0] = new Vector2(0.8f, 0);
                uvs[i * 2 + 1] = new Vector2(0.8f, 1);

                if (i > 0)
                {
                    int idx = (i - 1) * 6;
                    triangles[idx + 0] = (i - 1) * 2 + 0;
                    triangles[idx + 1] = (i - 1) * 2 + 1;
                    triangles[idx + 2] = i * 2 + 1;

                    triangles[idx + 3] = (i - 1) * 2 + 0;
                    triangles[idx + 4] = i * 2 + 1;
                    triangles[idx + 5] = i * 2 + 0;
                }

                prevMeshPos = curMeshPos;
            }

            var colors = new Color[vertices.Length];
            for (int i = 0; i < colors.Length; ++i)
            {
                //if (color.g < 0.5)
                    colors[i] = color;
                /*else
                    colors[i] = new Color(0, (float)i / (float)colors.Length, 1.0f - (float)i / (float)colors.Length);*/
            }

            mesh.Clear();
            mesh.vertices = vertices;
            mesh.uv = uvs;
            mesh.colors = colors;
            mesh.triangles = triangles;
            mesh.RecalculateBounds();
            mesh.MarkDynamic();
        }

        public void Calculate()
        {
            if (navigator.controls.showPreview) {
            // Destroy existing lines
            if (segments != null) {
            foreach(var segment in segments) {
            if (segment.line != null) {
                UnityEngine.Object.Destroy(segment.line);
            }
            }
            }
            // New segments array
            segments = new PreviewSegment[navigator.controls.ncontrols];
            // Beginning time
            double UT0 = navigator.controls.UT0;
            Orbit orbitInitial = navigator.vessel.orbit;
            // Initial mass per segment
            double m0i = navigator.vessel.GetTotalMass();
            // Calculate each segment
            for (var i = 0; i < segments.Length; i++) {
            // End time
            double UTf = UT0 + navigator.controls.controls[i].duration;
            // Calculate segment
            segments[i] = new PreviewSegment(navigator, orbitInitial, UT0, UTf, navigator.controls.controls[i], navigator.controls.controls[i].color, m0i);
            // Update initial mass for next segment
            m0i = segments[i].m1;
            // Update initial time
            UT0 = UTf;
            // Update initial orbit
            orbitInitial = segments[i].orbitf;
            }

            // Final time of trajectory
            this.UTf = UT0;

            // Draw one complete final orbit
            // Destroy existing line
            if (linef != null) {
            UnityEngine.Object.Destroy(linef);
            }
            // Create linerenderer & object
            var objf = new GameObject("Final orbit");
            linef = objf.AddComponent<LineRenderer>();
            linef.useWorldSpace = false;
            objf.layer = 10; // Map
            linef.material = MapView.fetch.orbitLinesMaterial;
            linef.SetColors(navigator.controls.colorFinal, navigator.controls.colorFinal);
            linef.SetVertexCount(361);
            // 3D points to use in linef
            linefPoints = new Vector3d[361];
            // Final orbit of trajectory
            orbitf = orbitInitial;
            // Period of final orbit
            double TPf = orbitf.period;
            // Populate points
            for(var i = 0; i <= 360; i++) {
            double UTi = this.UTf + i * TPf / 360;
            // Relative orbitf position
            Vector3d rRelOrbitf = orbitf.getRelativePositionAtUT(UTi).xzy;
            // Absolute position
            linefPoints[i] = rRelOrbitf;
            }

            // Target line
            CalculateTargetLine();
            }
        }

        public void CalculateTargetLine()
        {
            // Selected target
            var target = FlightGlobals.fetch.VesselTarget;
            // If a target is selected...
            if (target != null) {
            orbitT = target.GetOrbit(); // Target orbit
             // Spacecraft relative position at UTf
            rFinalRel = orbitf.getRelativePositionAtUT(UTf).xzy;
            // Target relative position at UTf
            rTargetFinalRel = orbitT.getRelativePositionAtUT(UTf).xzy;
            // Distance to target at UTf
            targetD = Vector3d.Distance(rFinalRel, rTargetFinalRel);
            // Relative speed to target at UTf
            targetV = Vector3d.Distance(orbitf.getOrbitalVelocityAtUT(UTf), orbitT.getOrbitalVelocityAtUT(UTf));
            // Destroy current line if present
            if (lineT != null) {
            UnityEngine.Object.Destroy(lineT);
            }
            // Make line to target at UTf
            var objT = new GameObject("Line to target");
            lineT = objT.AddComponent<LineRenderer>();
            lineT.useWorldSpace = false;
            objT.layer = 10; // Map
            lineT.material = MapView.fetch.orbitLinesMaterial;
            lineT.SetColors(Color.red, Color.red);
            lineT.SetVertexCount(2);

            // Target errors
            ApErr = orbitf.ApR - orbitT.ApR;
            PeErr = orbitf.PeR - orbitT.PeR;
            TPErr = orbitf.period - orbitT.period;
            IncErr = orbitf.inclination - orbitT.inclination;
            EccErr = orbitf.eccentricity - orbitT.eccentricity;
            LANErr = orbitf.LAN - orbitT.LAN;
            AOPErr = orbitf.argumentOfPeriapsis - orbitT.argumentOfPeriapsis;
            }
        }

		static ManeuverParameters ComputeEjectionManeuver(Vector3d exit_velocity, Orbit initial_orbit, double UT_0)
		{
			double GM = initial_orbit.referenceBody.gravParameter;
			double C3 = exit_velocity.sqrMagnitude;
			double sma = -GM / C3;
			double Rpe = initial_orbit.semiMajorAxis;
			double ecc = 1 - Rpe / sma;
			double theta = Math.Acos(-1 / ecc);

			Vector3d intersect_1;
			Vector3d intersect_2;

			// intersect_1 is the optimal position on the orbit assuming the orbit is circular and the sphere of influence is infinite
			IntersectConePlane(exit_velocity, theta, initial_orbit.h, out intersect_1, out intersect_2);

			Vector3d eccvec = initial_orbit.GetEccVector();

			double true_anomaly = AngleAboutAxis(eccvec, intersect_1, initial_orbit.h);
			// GetMeanAnomalyAtTrueAnomaly expects degrees and returns radians
			double mean_anomaly = initial_orbit.GetMeanAnomalyAtTrueAnomaly(true_anomaly * 180 / Math.PI);

			double delta_mean_anomaly = MuUtils.ClampRadiansPi(mean_anomaly - initial_orbit.MeanAnomalyAtUT(UT_0));

			double UT = UT_0 + delta_mean_anomaly / initial_orbit.MeanMotion();

			Vector3d V_0 = initial_orbit.getOrbitalVelocityAtUT(UT);
			Vector3d pos = initial_orbit.getRelativePositionAtUT(UT);

			double final_vel = Math.Sqrt(C3 + 2 * GM / pos.magnitude);

			Vector3d x = pos.normalized;
			Vector3d z = Vector3d.Cross(x, exit_velocity).normalized;
			Vector3d y = Vector3d.Cross(z, x);

			theta = Math.PI / 2;
			double dtheta = 0.001;
			Orbit sample = new Orbit();

			double theta_err = Double.MaxValue;

			for (int iteration = 0 ; iteration < 50 ; ++iteration)
			{
				Vector3d V_1 = final_vel * (Math.Cos(theta) * x + Math.Sin(theta) * y);
				sample.UpdateFromStateVectors(pos, V_1, initial_orbit.referenceBody, UT);
				theta_err = AngleAboutAxis(exit_velocity, sample.getOrbitalVelocityAtUT(OrbitExtensions.NextTimeOfRadius(sample, UT, sample.referenceBody.sphereOfInfluence)), z);
				if (double.IsNaN(theta_err))
					return null;
				if (Math.Abs(theta_err) <= Math.PI / 1800)
					return new ManeuverParameters((V_1 - V_0).xzy, UT);

				V_1 = final_vel * (Math.Cos(theta + dtheta) * x + Math.Sin(theta + dtheta) * y);
				sample.UpdateFromStateVectors(pos, V_1, initial_orbit.referenceBody, UT);
				double theta_err_2 = AngleAboutAxis(exit_velocity, sample.getOrbitalVelocityAtUT(OrbitExtensions.NextTimeOfRadius(sample, UT, sample.referenceBody.sphereOfInfluence)), z);

				V_1 = final_vel * (Math.Cos(theta - dtheta) * x + Math.Sin(theta - dtheta) * y);
				sample.UpdateFromStateVectors(pos, V_1, initial_orbit.referenceBody, UT);
				double theta_err_3 = AngleAboutAxis(exit_velocity, sample.getOrbitalVelocityAtUT(OrbitExtensions.NextTimeOfRadius(sample, UT, sample.referenceBody.sphereOfInfluence)), z);

				double derr = MuUtils.ClampRadiansPi(theta_err_2 - theta_err_3) / (2 * dtheta);

				theta = MuUtils.ClampRadiansTwoPi(theta - theta_err / derr);

				// if theta > pi, replace with 2pi - theta, the function ejection_angle=f(theta) is symmetrc wrt theta=pi
				if (theta > Math.PI)
					theta = 2 * Math.PI - theta;
			}

			return null;
		}

        private void initMeshFromOrbit(Transform meshTransform, Mesh mesh, Orbit orbit, double startTime, double duration, Color color)
        {
            int steps = 256;

            var vertices = new Vector3[steps * 2 + 2];
            var triangles = new int[steps * 6];

            Vector3 camPos = meshTransform.InverseTransformPoint(MapView.MapCamera.transform.position);

            Vector3 prevMeshPos = orbit.getRelativePositionAtUT(startTime - duration / (double)steps);
            for (int i = 0; i <= steps; ++i)
            {
                double time = startTime + duration * (double)i / (double)steps;

                Vector3 curMeshPos = orbit.getRelativePositionAtUT(time);
                if (Settings.fetch.BodyFixedMode) {
                    curMeshPos = Trajectory.calculateRotatedPosition(orbit.referenceBody, curMeshPos, time);
                }

                // compute an "up" vector that is orthogonal to the trajectory orientation and to the camera vector (used to correctly orient quads to always face the camera)
                Vector3 up = Vector3.Cross(curMeshPos - prevMeshPos, camPos - curMeshPos).normalized * (lineWidth * Vector3.Distance(camPos, curMeshPos));

                // add a segment to the trajectory mesh
                vertices[i * 2 + 0] = curMeshPos - up;
                vertices[i * 2 + 1] = curMeshPos + up;

                if (i > 0)
                {
                    int idx = (i - 1) * 6;
                    triangles[idx + 0] = (i - 1) * 2 + 0;
                    triangles[idx + 1] = (i - 1) * 2 + 1;
                    triangles[idx + 2] = i * 2 + 1;

                    triangles[idx + 3] = (i - 1) * 2 + 0;
                    triangles[idx + 4] = i * 2 + 1;
                    triangles[idx + 5] = i * 2 + 0;
                }

                prevMeshPos = curMeshPos;
            }

            var colors = new Color[vertices.Length];
            for (int i = 0; i < colors.Length; ++i)
                colors[i] = color;

            mesh.vertices = vertices;
            mesh.colors = colors;
            mesh.triangles = triangles;
            mesh.RecalculateBounds();
        }

		private double FindOrbitBodyIntersection(Orbit orbit, double startTime, double endTime, double bodyAltitude)
		{
			// binary search of entry time in atmosphere
			// I guess an analytic solution could be found, but I'm too lazy to search it

			double from = startTime;
			double to = endTime;

			int loopCount = 0;
			while (to - from > 0.1)
			{
				++loopCount;
				if (loopCount > 1000)
				{
					UnityEngine.Debug.Log("WARNING: infinite loop? (Trajectories.Trajectory.AddPatch, atmosphere limit search)");
					++errorCount_;
					break;
				}
				double middle = (from + to) * 0.5;
				if (orbit.getRelativePositionAtUT(middle).magnitude < bodyAltitude)
				{
					to = middle;
				}
				else
				{
					from = middle;
				}
			}

			return to;
		}

 /// <summary>
 /// Gets the position of the orbit at the specified time, in world coordinates.
 /// </summary>
 /// <param name="orbit"></param>
 /// <param name="time">UT</param>
 /// <returns></returns>
 private static Vector3d AbsolutePositionAt(Orbit orbit, double time)
 {
     // We use .xzy because Orbit class functions appear to use a coordinate
     // system in which Y and Z axes are swapped.  Thanks to sarbian for pointing
     // this out, here:
     // https://github.com/MuMech/MechJeb2/blob/master/MechJeb2/OrbitExtensions.cs#L18-L20
     return orbit.referenceBody.position + orbit.getRelativePositionAtUT(time).xzy;
 }

 private Vector3d GetWorldPositionAtUT(Orbit orbit, double ut)
 {
     Vector3d worldPos = Util.SwapYZ(orbit.getRelativePositionAtUT(ut));
     if (orbit.referenceBody != FlightGlobals.Bodies[0])
         worldPos += GetWorldPositionAtUT(orbit.referenceBody.orbit, ut);
     return worldPos;
 }