public static Vector3d CartesianToSpherical(Vector3d cartesian)
 {
     double r = cartesian.magnitude;
     double theta = Math.Atan(cartesian.y / cartesian.x) * RadToDeg;
     double phi = Math.Acos(cartesian.z / r) * RadToDeg;
     return new Vector3d(r, theta, phi);
 }

 //The KSP drag law is dv/dt = -b * v^2 where b is proportional to the air density and
 //the ship's drag coefficient. In this equation b has units of inverse length. So 1/b
 //is a characteristic length: a ship that travels this distance through air will lose a significant
 //fraction of its initial velocity
 public static double DragLength(this CelestialBody body, Vector3d pos, double dragCoeffOverMass)
 {
     double airDensity = FlightGlobals.getAtmDensity(FlightGlobals.getStaticPressure(pos, body));
     if (airDensity <= 0)
         return Double.MaxValue;
     return 1.0 / (0.5 * FlightGlobals.DragMultiplier * airDensity * dragCoeffOverMass);
 }

        private void InitializeSuffixes()
        {
            AddSuffix("X", new SetSuffix<double>(() => X, value => X = value));
            AddSuffix("Y", new SetSuffix<double>(() => Y, value => Y = value));
            AddSuffix("Z", new SetSuffix<double>(() => Z, value => Z = value));
            AddSuffix("MAG", new SetSuffix<double>(Magnitude, value =>
            {
                double oldMag = new Vector3d(X, Y, Z).magnitude;

                if (oldMag == 0) return; // Avoid division by zero

                X = X / oldMag * value;
                Y = Y / oldMag * value;
                Z = Z / oldMag * value;
            }));
            AddSuffix("VEC", new Suffix<Vector>(() => new Vector(X, Y, Z)));
            AddSuffix("NORMALIZED", new Suffix<Vector>(Normalized));
            AddSuffix("SQRMAGNITUDE", new Suffix<double>(() => new Vector3d(X, Y, Z).sqrMagnitude));
            AddSuffix("DIRECTION", new SetSuffix<Direction>(ToDirection, value =>
            {
                var newMagnitude = Vector3d.forward * new Vector3d(X, Y, Z).magnitude;

                var newVector = value.Rotation * newMagnitude;

                X = newVector.x;
                Y = newVector.y;
                Z = newVector.z;
            }));
        }

        public void Load(ConfigNode node)
        {
            if (!HighLogic.LoadedSceneIsFlight || FlightGlobals.ActiveVessel == null || parent.vessel == null) {
                return;
            }

            string[] values = node.GetValues("node");
            int max = values.Length;
            for (int k = 0; k < max; k++) {
                string[] info = values[k].Split(delimiter, StringSplitOptions.RemoveEmptyEntries);
                Vector3d deltav = new Vector3d();
                double d = 0.0;

                if (info.Length == 4) {
                    double.TryParse(info[0], out d);
                    deltav.x = d;

                    d = 0.0;
                    double.TryParse(info[1], out d);
                    deltav.y = d;

                    d = 0.0;
                    double.TryParse(info[2], out d);
                    deltav.z = d;

                    d = 0.0;
                    double.TryParse(info[3], out d);

                    // at the very least it'll /act/ like a proper maneuver node.
                    nodes.Add(new NodeState(deltav, d));
                }
            }
            Debug.Log("Node Saver loaded " + max + " nodes.");
        }

        public void RunCalculations(
            Vessel vessel,
            Quaternion gymbal)
        {
            ManeuverPresent = false;

            // Calculations thanks to Mechjeb
            Vector3d CoM = vessel.findWorldCenterOfMass();
            Vector3d velocityVesselOrbit = vessel.orbit.GetVel();
            Vector3d velocityVesselOrbitUnit = velocityVesselOrbit.normalized;
            Vector3d velocityVesselSurface = velocityVesselOrbit - vessel.mainBody.getRFrmVel(CoM);
            Vector3d velocityVesselSurfaceUnit = velocityVesselSurface.normalized;

            ProgradeOrbit = gymbal * velocityVesselOrbitUnit;
            ProgradeSurface = gymbal * velocityVesselSurfaceUnit;

            PatchedConicSolver patchedConicSolver = vessel.patchedConicSolver;
            if (patchedConicSolver == null)
                return;

            if (patchedConicSolver.maneuverNodes == null)
                return;

            if (patchedConicSolver.maneuverNodes.Count > 0)
            {
                var manuever = patchedConicSolver.maneuverNodes[0];
                ManeuverApplied = manuever.DeltaV != Vector3d.zero;

                Vector3d burnVector = manuever.GetBurnVector(vessel.orbit);
                ManeuverPlus = gymbal * burnVector.normalized;
                ManeuverPresent = true;
            }
        }

 public void AddWeightedCorrection(Vector3d cor)
 {
     var cm = cor.magnitude;
     if(cm > 1e-10) cor *= Math.Sqrt(1/cm);
     if(VSL.Physics.G > 1e-10) cor *= Utils.ClampH(Utils.G0/VSL.Physics.G, HSC.MaxCorrectionWeight);
     CourseCorrections.Add(cor);
 }

        /// <summary>
        /// Copy Constructor.
        /// </summary>
        /// <param name="orig"></param>
        public SpawnVessel(SpawnVessel orig)
        {
            deferVesselCreation = orig.deferVesselCreation;
            foreach (VesselData vessel in orig.vessels)
            {
                if (vessel.pqsCity != null)
                {
                    // Generate PQS city coordinates
                    LoggingUtil.LogVerbose(this, "Generating coordinates from PQS city for Vessel " + vessel.name);

                    // Translate by the PQS offset (inverse transform of coordinate system)
                    Vector3d position = vessel.pqsCity.transform.position;
                    LoggingUtil.LogVerbose(this, "    pqs city position = " + position);
                    Vector3d v = vessel.pqsOffset;
                    Vector3d i = vessel.pqsCity.transform.right;
                    Vector3d j = vessel.pqsCity.transform.forward;
                    Vector3d k = vessel.pqsCity.transform.up;
                    LoggingUtil.LogVerbose(this, "    i, j, k = " + i + ", " + j + "," + k);
                    Vector3d offsetPos = new Vector3d(
                        (j.y * k.z - j.z * k.y) * v.x + (i.z * k.y - i.y * k.z) * v.y + (i.y * j.z - i.z * j.y) * v.z,
                        (j.z * k.x - j.x * k.z) * v.x + (i.x * k.z - i.z * k.x) * v.y + (i.z * j.x - i.x * j.z) * v.z,
                        (j.x * k.y - j.y * k.x) * v.x + (i.y * k.x - i.x * k.y) * v.y + (i.x * j.y - i.y * j.x) * v.z
                    );
                    offsetPos *= (i.x * j.y * k.z) + (i.y * j.z * k.x) + (i.z * j.x * k.y) - (i.z * j.y * k.x) - (i.y * j.x * k.z) - (i.x * j.z * k.y);
                    vessel.latitude = vessel.body.GetLatitude(position + offsetPos);
                    vessel.longitude = vessel.body.GetLongitude(position + offsetPos);
                    LoggingUtil.LogVerbose(this, "    resulting lat, lon = (" + vessel.latitude + ", " + vessel.longitude + ")");
                }

                vessels.Add(new VesselData(vessel));
            }
        }

 /// <summary>
 /// Converts a Cartesian 3D vector to latitude and longitude.
 /// </summary>
 /// <param name="vec">Vector3d of indicated position (convention: prime meridian is +x, north pole is +z)</param>
 /// <returns>Pair of lat/lon coordinates (in degrees) of indicated position. </returns>
 public static double[] Cart2LatLon(Vector3d vec)
 {
     vec = Vector3d.Normalize (vec);
     double lat = Rad2Deg (Math.Asin (vec.z));
     double lon = Rad2Deg (Math.Atan2 (vec.x, vec.y));
     return new double[] { lat, lon };
 }

        public static ManeuverNode PlaceManeuverNode(this Vessel vessel, Orbit patch, Vector3d dV, double UT)
        {
            //placing a maneuver node with bad dV values can really mess up the game, so try to protect against that
            //and log an exception if we get a bad dV vector:
            for (int i = 0; i < 3; i++)
            {
                if (double.IsNaN(dV[i]) || double.IsInfinity(dV[i]))
                {
                    throw new Exception("VesselExtensions.PlaceManeuverNode: bad dV: " + dV);
                }
            }

            if (double.IsNaN(UT) || double.IsInfinity(UT))
            {
                throw new Exception("VesselExtensions.PlaceManeuverNode: bad UT: " + UT);
            }

            //It seems that sometimes the game can freak out if you place a maneuver node in the past, so this
            //protects against that.
            UT = Math.Max(UT, Planetarium.GetUniversalTime());

            //convert a dV in world coordinates into the coordinate system of the maneuver node,
            //which uses (x, y, z) = (radial+, normal-, prograde)
            Vector3d nodeDV = patch.DeltaVToManeuverNodeCoordinates(UT, dV);
            ManeuverNode mn = vessel.patchedConicSolver.AddManeuverNode(UT);
            mn.OnGizmoUpdated(nodeDV, UT);
            return mn;
        }

 public Direction(Quaternion q)
     : this()
 {
     rotation = q;
     euler = q.eulerAngles;
     vector = rotation * Vector3d.forward;
 }

        public void UpdateForceAndTorque(Vector3d com)
        {
            PartModule gimbal;
            VesselState.GimbalExt gimbalExt = VesselState.getGimbalExt(engine.part, out gimbal);
            Vessel v = engine.vessel;

            float currentMaxThrust, currentMinThrust;

            currentMaxThrust = engine.minFuelFlow * engine.flowMultiplier * engine.realIsp * engine.g / engine.thrustTransforms.Count;
            currentMinThrust = engine.maxFuelFlow * engine.flowMultiplier * engine.realIsp * engine.g / engine.thrustTransforms.Count;

            if (engine.throttleLocked)
            {
                currentMaxThrust *= thrustRatio;
                currentMinThrust = currentMaxThrust;
            }

            for (int i = 0; i < engine.thrustTransforms.Count; i++)
            {
                Vector3d thrust_dir = v.transform.rotation.Inverse() * gimbalExt.initialRot(gimbal, engine.thrustTransforms[i], i) * Vector3d.back;

                Vector3d pos = v.transform.rotation.Inverse() * (engine.part.transform.position - com);

                _maxVariableForce += (currentMaxThrust - currentMinThrust) * thrust_dir;
                _constantForce += currentMinThrust * thrust_dir;
                _maxVariableTorque += (currentMaxThrust - currentMinThrust) * Vector3d.Cross(pos, thrust_dir);
                _constantTorque += currentMinThrust * Vector3d.Cross(pos, thrust_dir);
            }
        }

        protected override void _Drive(FlightCtrlState s)
        {
            Vessel vessel = base.Autopilot.vessel;
            AttitudeController attitudeController = base.Autopilot.AutopilotController.AttitudeController;

            if (vessel.altitude < 10000) {
                attitudeController.Attitude = Attitude.Up;
                base.Autopilot.AutopilotController.AutoStagingController.Enable ();
                s.mainThrottle = 1.0f;
            } else if (vessel.altitude < 80000) {
                attitudeController.Pitch = 90;
                attitudeController.Yaw = 45;
                attitudeController.Attitude = Attitude.UserDefined;
                s.mainThrottle = Mathf.Clamp01 (100000.0f - (float)vessel.orbit.ApA);
            } else {
                double ut = Planetarium.GetUniversalTime () + vessel.orbit.timeToAp;

                double deltaV = Math.Sqrt (vessel.mainBody.gravParameter / (100000 + vessel.mainBody.Radius)) - vessel.orbit.getOrbitalVelocityAtUT (ut).magnitude;

                Vector3d nodeDV = new Vector3d (0, 0, deltaV);
                base.Autopilot.ManeuverNode = vessel.patchedConicSolver.AddManeuverNode (ut);
                base.Autopilot.ManeuverNode.OnGizmoUpdated (nodeDV, ut);

                base.Autopilot.AutopilotController.Disable ();
                base.Autopilot.AutopilotController.AutoStagingController.Disable ();
                attitudeController.Attitude = Attitude.None;
            }
        }

        /// <summary>
        /// Orbit foo, this finds the nodes of two orbits
        /// </summary>
        /// <returns>
        /// The true anomaly of the ascending node(descing node is 180degrees off)
        /// </returns>
        /// <param name='orbit'>
        /// Orbit.
        /// </param>
        /// <param name='tgtorbit'>
        /// Target Orbit
        /// </param>
        public static double FindAN(Orbit orbit, Orbit tgtorbit)
        {
            double rad = Math.PI/180;
            double Lan1 = orbit.LAN;
            double inc1 = orbit.inclination;
            double Lan2 = tgtorbit.LAN;
            double inc2 = tgtorbit.inclination;

            //see braeunig.us/space... cross product of two orbital planes gives the node location
            var a = new Vector3d(Math.Sin(inc1*rad)*Math.Cos(Lan1*rad), Math.Sin(inc1*rad)*Math.Sin(Lan1*rad),
                                 Math.Cos(inc1*rad));
            var b = new Vector3d(Math.Sin(inc2*rad)*Math.Cos(Lan2*rad), Math.Sin(inc2*rad)*Math.Sin(Lan2*rad),
                                 Math.Cos(inc2*rad));
            var c = new Vector3d(0, 0, 0);
            c = Vector3d.Cross(a, b);
            var coord = new double[] {0, 0};
            coord = LatLonofVector(c); //get the coordinates lat/lon of the ascending node
            double lat = coord[0];
            double lon = coord[1];

            //go look at the diagrams at braeunig.us/space
            double α = lon - Lan1; //its all greek to me
            if (α < 0) α += 360;
            double λ = Math.Atan(Math.Tan(α*rad)/Math.Cos(inc1*rad))/rad;
            double x = 180 + (λ - orbit.argumentOfPeriapsis);
            if (x > 360) return 360 - x;
            else return x;
        }

        /// <summary>
        /// This function should return exactly the same value as Vessel.atmDensity, but is more generic because you don't need an actual vessel updated by KSP to get a value at the desired location.
        /// Computations are performed for the current body position, which means it's theoritically wrong if you want to know the temperature in the future, but since body rotation is not used (position is given in sun frame), you should get accurate results up to a few weeks.
        /// </summary>
        /// <param name="position"></param>
        /// <param name="body"></param>
        /// <returns></returns>
        public static double GetTemperature(Vector3d position, CelestialBody body)
        {
            if (!body.atmosphere)
                return PhysicsGlobals.SpaceTemperature;

            double altitude = (position - body.position).magnitude - body.Radius;
            if (altitude > body.atmosphereDepth)
                return PhysicsGlobals.SpaceTemperature;

            Vector3 up = (position - body.position).normalized;
            float polarAngle = Mathf.Acos(Vector3.Dot(body.bodyTransform.up, up));
            if (polarAngle > Mathf.PI / 2.0f)
            {
                polarAngle = Mathf.PI - polarAngle;
            }
            float time = (Mathf.PI / 2.0f - polarAngle) * 57.29578f;

            Vector3 sunVector = (FlightGlobals.Bodies[0].position - position).normalized;
            float sunAxialDot = Vector3.Dot(sunVector, body.bodyTransform.up);
            float bodyPolarAngle = Mathf.Acos(Vector3.Dot(body.bodyTransform.up, up));
            float sunPolarAngle = Mathf.Acos(sunAxialDot);
            float sunBodyMaxDot = (1.0f + Mathf.Cos(sunPolarAngle - bodyPolarAngle)) * 0.5f;
            float sunBodyMinDot = (1.0f + Mathf.Cos(sunPolarAngle + bodyPolarAngle)) * 0.5f;
            float sunDotCorrected = (1.0f + Vector3.Dot(sunVector, Quaternion.AngleAxis(45f * Mathf.Sign((float)body.rotationPeriod), body.bodyTransform.up) * up)) * 0.5f;
            float sunDotNormalized = (sunDotCorrected - sunBodyMinDot) / (sunBodyMaxDot - sunBodyMinDot);
            double atmosphereTemperatureOffset = (double)body.latitudeTemperatureBiasCurve.Evaluate(time) + (double)body.latitudeTemperatureSunMultCurve.Evaluate(time) * sunDotNormalized + (double)body.axialTemperatureSunMultCurve.Evaluate(sunAxialDot);
            double temperature = body.GetTemperature(altitude) + (double)body.atmosphereTemperatureSunMultCurve.Evaluate((float)altitude) * atmosphereTemperatureOffset;

            return temperature;
        }

 public Vector3d Start(Vector3d approachPosition)
 {
     Vector3d startPos = start.Position();
     Vector3d endPos = end.Position();
     if (Vector3d.Distance(startPos, approachPosition) < Vector3d.Distance(endPos, approachPosition)) return startPos;
     else return endPos;
 }

            public double GetAngleOfAttack(Vector3d position, Vector3d velocity)
            {
                if (!horizon)
                    return angle;

                return Math.Acos(Vector3d.Dot(position, velocity) / (position.magnitude * velocity.magnitude)) - Math.PI * 0.5 + angle;
            }

            public void BuildSection(Vector3d lowerInputs, Vector3d upperInputs)
            {
                //Creates cubic from x,y,dy/dx data

                double recipXDiff, recipXDiffSq;
                recipXDiff = 1 / (lowerInputs.x - upperInputs.x);
                recipXDiffSq = recipXDiff * recipXDiff;

                a = 2 * (upperInputs.y - lowerInputs.y) * recipXDiff;
                a += upperInputs.z + lowerInputs.z;
                a *= recipXDiffSq;

                b = 3 * (upperInputs.x + lowerInputs.x) * (lowerInputs.y - upperInputs.y) * recipXDiff;
                b -= (lowerInputs.x + 2 * upperInputs.x) * lowerInputs.z;
                b -= (2 * lowerInputs.x + upperInputs.x) * upperInputs.z;
                b *= recipXDiffSq;

                c = 6 * upperInputs.x * lowerInputs.x * (upperInputs.y - lowerInputs.y) * recipXDiff;
                c += (2 * lowerInputs.x * upperInputs.x + upperInputs.x * upperInputs.x) * lowerInputs.z;
                c += (2 * lowerInputs.x * upperInputs.x + lowerInputs.x * lowerInputs.x) * upperInputs.z;
                c *= recipXDiffSq;

                d = (3 * lowerInputs.x - upperInputs.x) * upperInputs.x * upperInputs.x * lowerInputs.y;
                d += (lowerInputs.x - 3 * upperInputs.x) * lowerInputs.x * lowerInputs.x * upperInputs.y;
                d *= recipXDiff;

                d -= lowerInputs.x * upperInputs.x * upperInputs.x * lowerInputs.z;
                d -= lowerInputs.x * lowerInputs.x * upperInputs.x * upperInputs.z;
                d *= recipXDiffSq;

                upperLim = upperInputs.x;
                lowerLim = lowerInputs.x;
            }

 public bool Execute(Vector3d dV, float MinDeltaV = 0.1f, ManeuverCondition condition = null)
 {
     THR.Throttle = 0;
     dVrem.Value = dV.magnitude;
     //end if below the minimum dV
     if(dVrem < MinDeltaV) return false;
     VSL.Engines.ActivateEngines();
     if(VSL.Engines.MaxThrustM.Equals(0)) return true;
     //orient along the burning vector
     if(dVrem && VSL.Controls.RCSAvailableInDirection(-dV))
         CFG.AT.OnIfNot(Attitude.KillRotation);
     else
     {
         CFG.AT.OnIfNot(Attitude.Custom);
         ATC.SetThrustDirW(-dV);
     }
     //check the condition
     if(condition != null && !condition((float)dVrem)) return true;
     if(VSL.Controls.TranslationAvailable)
     {
         if(dVrem || VSL.Controls.AttitudeError > GLB.ATCB.AttitudeErrorThreshold)
             TRA.AddDeltaV(-VSL.LocalDir(dV));
         if(dVrem && CFG.AT[Attitude.KillRotation])
         {
             var errorF = Utils.ClampL(Vector3.Dot(VSL.Engines.Thrust.normalized, -dV.normalized), 0);
             THR.DeltaV = (float)dVrem * errorF*errorF;
         }
         else THR.DeltaV = (float)dVrem;
     }
     else THR.DeltaV = (float)dVrem;
     //			Log("\ndVrem: {}\nAttitudeError {}, DeltaV: {}, Throttle {}, RCS {}",
     //			    dVrem, VSL.Controls.AttitudeError, THR.DeltaV, THR.Throttle, VSL.Controls.RCSAvailableInDirection(-dV));//debug
     return true;
 }

		void GgsUpdate(){
			if (core.attitude.attitudeAngleFromTarget() < 1)
            {
					Vector3d ang;
					Quaternion rot;
					Quaternion ori;
					Vector3d gee;
            		ang = new Vector3d(inpX, inpY, inpZ);
					gee = FlightGlobals.getGeeForceAtPosition(this.vessel.transform.position);
					rot = Quaternion.LookRotation(gee);
					ori = Quaternion.LookRotation(gee);
					Quaternion x = Quaternion.AngleAxis(inpX, Vector3.Cross(gee,Vector3.up).normalized);
					Quaternion y = Quaternion.AngleAxis(360 - inpY, gee);
				
				if(modo == Modo.Horizonte) {
					Quaternion z = Quaternion.AngleAxis(360 - inpZ, this.vessel.transform.InverseTransformDirection(this.vessel.transform.up));	
					part.vessel.SetRotation(y * x * Quaternion.LookRotation(gee) * z);
				}
				if (modo == Modo.Prograde) {
					Quaternion z = Quaternion.AngleAxis(-(inpZ-90)+180, this.vessel.transform.InverseTransformDirection(this.vessel.transform.up));
					Vector3 pro = part.vessel.orbit.GetVel();
					//Vector3d gee = FlightGlobals.getGeeForceAtPosition(this.vessel.transform.position);
					Vector3 realPro = Vector3.Cross(gee*-1,pro);
					x = Quaternion.AngleAxis(inpY, Vector3.Cross(realPro,Vector3.up).normalized);
					y = Quaternion.AngleAxis(360 - inpX, realPro);
					part.vessel.SetRotation(y * x * Quaternion.LookRotation(realPro,pro) * z);
					
				}
            }
		}

            public Vector3d velocity; // velocity in world frame relatively to the reference body

            public VesselState(Vessel vessel)
            {
                referenceBody = vessel.orbit.referenceBody;
                time = Planetarium.GetUniversalTime();
                position = vessel.GetWorldPos3D() - referenceBody.position;
                velocity = vessel.obt_velocity;
            }