using System;
using UnityEngine;
using MuMech;
//*TODO, FindAN() seems to reverse the AN and DN either some, or all of the time.
namespace OrbitExtensions
{
internal static class OrbitExtensions
{
///
/// Eta to a given true anomaly.
///
///
/// Time (seconds)
///
///
/// Orbit.
///
///
/// True anomaly. (0-360)
///
public static double GetTimeToTrue(this Orbit orbit, double Tanom)
{
double t = GetTimeAtTrue(orbit, Tanom);
//t= orbit.GetTimeToMean (orbit.GetMeanAtEccentric (orbit.GetEccentricAtTrue ((360 - orbit.argumentOfPeriapsis))));
double p = orbit.period;
double now = p - orbit.timeToPe;
if (now < t)
return t - now;
else
return orbit.timeToPe + t;
}
///
/// Gets the time for a given true anomaly.
///
///
/// The time untill the target true anomaly in seconds.
///
///
/// True anomaly. (0-360)
///
public static double GetTimeAtTrue(this Orbit orbit, double TAnom)
{
return orbit.GetTimeAtMean(orbit.GetMeanAtEccentric(orbit.GetEccentricAtTrue(TAnom)));
}
///
/// Finds the coordinates of a state vector.
///
///
/// Double[] {lattitude, longitude}
///
///
/// State vector
///
public static double[] LatLonofVector(Vector3d vinput)
{
//get the geocentric latitude and longitude of a orbital element vector
double rad = Math.PI/180;
double c1 = vinput.x;
double c2 = vinput.y;
double c3 = vinput.z;
double lon = 0;
double lat = Math.Atan(c3/Math.Pow((c1*c1) + (c2*c2), .5))/rad;
if (c1 < 0)
{
lon = (Math.Atan(c2/c1)/rad) + 90;
}
else
{
lon = (Math.Atan(c2/c1)/rad) + 270;
}
var coord = new[] {lat, lon};
return coord;
}
///
/// Orbit foo, this finds the nodes of two orbits
///
///
/// The true anomaly of the ascending node(descing node is 180degrees off)
///
///
/// Orbit.
///
///
/// Target Orbit
///
public static double FindAN(this Orbit orbit, Orbit tgtorbit)
{
Vector3d normal = ARUtils.swapYZ(orbit.GetOrbitNormal());
Vector3d targetNormal = ARUtils.swapYZ(tgtorbit.GetOrbitNormal());
Vector3d vectorToAN = Vector3d.Cross(normal, targetNormal);
Vector3d vectorToPe = ARUtils.swapYZ(orbit.eccVec);
double angleFromPe = Vector3d.Angle(vectorToPe, vectorToAN);
//If the AN is actually during the infalling part of the orbit then we need to add 180 to
//angle from Pe to get the true anomaly. Test this by taking the the cross product of the
//orbit normal and vector to the periapsis. This gives a vector that points to center of the
//outgoing side of the orbit. If vectorToAN is more than 90 degrees from this vector, it occurs
//during the infalling part of the orbit.
if (Math.Abs(Vector3d.Angle(vectorToAN, Vector3d.Cross(vectorToPe, normal))) < 90)
{
return angleFromPe;
}
else
{
return 360 - angleFromPe;
}
/* 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;*/
}
///
/// Eta to this orbits ascending node relative to a second orbit.
///
///
/// The time to ascending node.
///
///
/// Orbit.
///
///
/// Target orbit
///
public static double GetTimeToRelAN(this Orbit orbit, Orbit tgtorbit)
{
return orbit.GetTimeToTrue(FindAN(orbit, tgtorbit));
}
///
/// ETA to this orbits descending node relative to a second orbit
///
///
/// The time to descending node
///
///
/// Orbit.
///
///
/// Target orbit
///
public static double GetTimeToRelDN(this Orbit orbit, Orbit tgtorbit)
{
double an = FindAN(orbit, tgtorbit);
double dn = an - 180;
if (dn < 0) dn += 360;
return orbit.GetTimeToTrue(dn);
}
///
/// ETA to ascending node (Refference Plane).
///
///
/// Remaining time in seconds.
///
public static double GetTimeToRefAN(this Orbit orbit)
{
return orbit.GetTimeToTrue(360 - orbit.argumentOfPeriapsis);
//return orbit.GetTimeToMean (orbit.GetMeanAtEccentric (orbit.GetEccentricAtTrue ((360 - orbit.argumentOfPeriapsis))));
}
///
/// ETA to descending node(Refference Plane).
///
///
/// Remaining time in seconds
///
public static double GetTimeToRefDN(this Orbit orbit)
{
if (orbit.argumentOfPeriapsis <= 180)
{
return orbit.GetTimeToTrue(180 - orbit.argumentOfPeriapsis);
}
else
return orbit.GetTimeToTrue(520 - orbit.argumentOfPeriapsis);
}
///
/// Gets the mean anomaly at a time.
///
///
/// The mean at time.
///
///
/// Orbit.
///
///
/// Time from periapsis
///
public static double GetMeanAnomalyAtTime(this Orbit orbit, double t)
{
double tau = 2*Math.PI;
double M = 0;
double p = orbit.period;
M = (t*tau)/p;
return M;
}
///
/// Takes a true anomaly(degrees from Periapsis), and computes Eccentric Anomaly at that target
///undefined if eccentricity > 1
///
///
/// Returns Eccentric Anomaly
///
/// Orbit">
///
/// True anomomaly (0-360)
///
///
public static double GetEccentricAtTrue(this Orbit orb, double TrueAnom)
{
double e = orb.eccentricity;
double E = 0f;
bool g180 = false;
if (TrueAnom > 180) g180 = true;
if (TrueAnom == 180) return Math.PI;
TrueAnom = TrueAnom*(Math.PI/180);
double cose;
double sine;
double tane;
cose = (e + Math.Cos(TrueAnom))/(1 + e*Math.Cos(TrueAnom));
sine = Math.Sqrt(1 - (cose*cose));
tane = sine/cose;
E = Math.Atan(tane);
if (E < 0.0)
{
E += 3/2*Math.PI;
}
if (g180)
{
E = (2*Math.PI - E);
}
return E;
}
///
/// You provide the eccentric anomaly, this method will bring the mean.
///
///
/// Mean anomaly
///
///
/// orbit.eccentricAnomaly (0-2pi)
///
public static double GetMeanAtEccentric(this Orbit orb, double E)
{
// M=E-e*SIN(E)
double e = orb.eccentricity;
double M = 0;
M = E - (e*Math.Sin(E));
return M;
}
///
/// Gets the time at a given mean anomaly.
///
///
/// The time from periapsis that the vessel will be at this mean anomaly
///
///
/// Orbit.
///
///
/// Mean anomaly (0-2pi)
///
///Not compatible with eccentricity >=1
public static double GetTimeAtMean(this Orbit orbit, double M)
{
double tau = 2*Math.PI;
double p = orbit.period;
return M*p/(tau);
}
///
/// Time untill vessel reaches specified mean.
///
/// Time untill vessel reaches mean anomaly in seconds
///
/// Mean anomaly
///
///Not compatible with eccentricyt >=1
public static double GetTimeToMean(this Orbit orbit, double M)
{
double t;
double p = orbit.period;
double timeFromPe;
t = GetTimeAtMean(orbit, M);
timeFromPe = p - orbit.timeToPe;
if (timeFromPe < t) return timeFromPe + t;
else return orbit.timeToPe + t;
}
///
/// Translates a true anomaly from this orbit to the argument's. Orbits should have the same reference body.
///
///
/// True anomaly on the target orbit that matches the input
///
///
/// Target Orbit
///
///
/// True Anomaly from this orbit.(0-360)
///
public static double TranslateAnomaly(this Orbit orbit, Orbit tgt_Orbit, double anomaly)
{
//a2= a1+lan1+ar1-(lan2+arg2)
double a = anomaly + orbit.argumentOfPeriapsis + orbit.LAN - (tgt_Orbit.argumentOfPeriapsis + tgt_Orbit.LAN);
if (a < 0) a += 360;
return a;
}
public static double Syncorbits(this Orbit orbit, Orbit tgt_orb, double anomaly, int orbitnum)
{
//translate true anomaly
double tgta = orbit.TranslateAnomaly(tgt_orb, anomaly);
//if t>target period subtract 1 period...
double p = tgt_orb.period;
//double tgtc = FlightGlobals.Vessels[selVessel].orbit.GetTimeAtMean(FlightGlobals.Vessels[selVessel].orbit.meanAnomaly);
//double t = this.vessel.orbit.GetTimeToTrue(anomaly);
return tgt_orb.GetTimeToTrue(tgta) + (p*orbitnum);
}
///
/// Find the Intersections of two orbits
///
///
/// 0 for no intersections found, 1 for 1, 2 for 2
/// True anomaly of intersections in out intersect1 and intersect2
///
///
/// Orbit.
///
///
/// Tgtorbit.
///
///
/// Intsect1.
///
///
/// Intsect2.
///
public static int FindIntersectSimple (this Orbit orbit, Orbit tgtorbit, out double intersect1, out double intersect2)// Doesn't work, just learning to use github
{
double anomaly = 0;
intersect1 = 361;
intersect2 = 361;
int div = 50;
double range = 360;
double lowest = 1e10;
int count = 0;
while (lowest>100) {
for (int i=0; i
30)
return 2;
}
return 0;
}
public static double f(this Orbit orbit, Orbit tgt_orbit, double x)
{
double rad = Math.PI/180;
double a1 = orbit.semiMajorAxis;
double e1 = orbit.eccentricity;
double w1 = orbit.argumentOfPeriapsis;
double l1 = orbit.LAN;
double a2 = tgt_orbit.semiMajorAxis;
double e2 = tgt_orbit.eccentricity;
double w2 = tgt_orbit.argumentOfPeriapsis;
double l2 = tgt_orbit.LAN;
return ((a2 - (e1*e1))/(e2*Math.Cos((x + w1 + l1 - (w2 + l2))*rad) + 1)) -
((a1 - (e2*e2))/(e1*Math.Cos(x*rad) + 1));
}
public static double fprime(this Orbit orbit, Orbit tgt_orbit, double x)
{
double rad = Math.PI/180;
double a1 = orbit.semiMajorAxis;
double e1 = orbit.eccentricity;
double w1 = orbit.argumentOfPeriapsis;
double l1 = orbit.LAN;
double a2 = tgt_orbit.semiMajorAxis;
double e2 = tgt_orbit.eccentricity;
double w2 = tgt_orbit.argumentOfPeriapsis;
double l2 = tgt_orbit.LAN;
return ((a2 - (e1*e1))/(e2*-Math.Sin(x))) - ((a1 - (e2*e2))/(e1*-Math.Sin(x*rad)));
}
//Orbit.getRelativePositionAtUT seems to be less precise, or something?
public static Vector3d getAbsolutePositionAtUT(this Orbit orbit, double UT)
{
Vector3d pos = orbit.getRelativePositionAtUT(UT);
pos += orbit.referenceBody.position;
return pos;
/* double timeFromNow = UT - Planetarium.GetUniversalTime();
double timeSincePe = (orbit.eccentricity > 1 ? -orbit.timeToPe : timeFromNow + (orbit.period - orbit.timeToPe));
return orbit.getPositionAtT(timeSincePe);*/
}
public static double relativeInclination(this Orbit a, Orbit b)
{
return Math.Abs(Vector3d.Angle(a.GetOrbitNormal(), b.GetOrbitNormal()));
}
}
}