Skip to content

Commit

Permalink
test: ✅ add testing for EquinoctialElements
Browse files Browse the repository at this point in the history
  • Loading branch information
thkruz committed Jan 14, 2024
1 parent 8120b3f commit 22ea0a8
Show file tree
Hide file tree
Showing 13 changed files with 556 additions and 200 deletions.
47 changes: 20 additions & 27 deletions src/coordinate/ClassicalElements.ts
Original file line number Diff line number Diff line change
@@ -1,24 +1,16 @@
import { Degrees, Kilometers, Radians, Seconds } from 'src/main';
import { Degrees, Kilometers, Radians, Seconds } from '../main';
import { Vector3D } from '../operations/Vector3D';
import { EpochUTC } from '../time/EpochUTC';
import { earthGravityParam, RAD2DEG, sec2min, secondsPerDay, TAU } from '../utils/constants';
import { clamp, matchHalfPlane, newtonNu } from '../utils/functions';
import { EquinoctialElements } from './EquinoctialElements';
import { OrbitRegime } from './OrbitRegime';
import { PositionVelocity, StateVector } from './StateVector';
import { ClassicalElementsParams } from '../interfaces/ClassicalElementsParams';

interface ClassicalElementsParams {
epoch: EpochUTC;
semimajorAxis: Kilometers;
eccentricity: number;
inclination: Radians;
rightAscension: Radians;
argPerigee: Radians;
trueAnomaly: Radians;
mu?: number;
}
/**
* The ClassicalElements class represents the classical orbital elements of an object.
* The ClassicalElements class represents the classical orbital elements of an
* object.
*
* @example
* ```ts
Expand Down Expand Up @@ -67,10 +59,10 @@ export class ClassicalElements {

/**
* Creates a new instance of ClassicalElements from a StateVector.
* @param state The StateVector to convert.
* @param mu The gravitational parameter of the central body. Default value is Earth's gravitational parameter.
* @returns A new instance of ClassicalElements.
* @throws Error if the StateVector is not in an inertial frame.
* @param state The StateVector to convert. @param mu The gravitational
* parameter of the central body. Default value is Earth's gravitational
* parameter. @returns A new instance of ClassicalElements. @throws Error if
* the StateVector is not in an inertial frame.
*/
static fromStateVector(state: StateVector, mu = earthGravityParam): ClassicalElements {
if (!state.inertial) {
Expand Down Expand Up @@ -155,8 +147,8 @@ export class ClassicalElements {
}

/**
* Gets the perigee of the classical elements.
* The perigee is the point in an orbit that is closest to the focus (in this case, the Earth).
* Gets the perigee of the classical elements. The perigee is the point in an
* orbit that is closest to the focus (in this case, the Earth).
* @returns The perigee distance in kilometers.
*/
get perigee(): number {
Expand All @@ -181,7 +173,7 @@ export class ClassicalElements {
* @returns The mean motion in radians.
*/
get meanMotion(): Radians {
return Math.sqrt(this.mu / (this.semimajorAxis * this.semimajorAxis * this.semimajorAxis)) as Radians;
return Math.sqrt(this.mu / this.semimajorAxis ** 3) as Radians;
}

/**
Expand Down Expand Up @@ -246,15 +238,16 @@ export class ClassicalElements {
* @returns {EquinoctialElements} The equinoctial elements.
*/
toEquinoctialElements(): EquinoctialElements {
const fr = Math.abs(this.inclination - Math.PI) < 1e-9 ? -1 : 1;
const af = this.eccentricity * Math.cos(this.argPerigee + fr * this.rightAscension);
const ag = this.eccentricity * Math.sin(this.argPerigee + fr * this.rightAscension);
const l = this.argPerigee + fr * this.rightAscension + newtonNu(this.eccentricity, this.trueAnomaly).m;
const n = this.meanMotion;
const chi = Math.tan(0.5 * this.inclination) ** fr * Math.sin(this.rightAscension);
const psi = Math.tan(0.5 * this.inclination) ** fr * Math.cos(this.rightAscension);
const I = this.inclination > Math.PI / 2 ? -1 : 1;
const h = this.eccentricity * Math.sin(this.argPerigee + I * this.rightAscension);
const k = this.eccentricity * Math.cos(this.argPerigee + I * this.rightAscension);
const meanAnomaly = newtonNu(this.eccentricity, this.trueAnomaly).m;
const lambda = (meanAnomaly + this.argPerigee + I * this.rightAscension) as Radians;
const a = this.semimajorAxis;
const p = Math.tan(0.5 * this.inclination) ** I * Math.sin(this.rightAscension);
const q = Math.tan(0.5 * this.inclination) ** I * Math.cos(this.rightAscension);

return new EquinoctialElements(this.epoch, af, ag, l, n, chi, psi, { mu: this.mu, fr });
return new EquinoctialElements({ epoch: this.epoch, k, h, lambda, a, p, q, mu: this.mu, I });
}

/**
Expand Down
189 changes: 143 additions & 46 deletions src/coordinate/EquinoctialElements.ts
Original file line number Diff line number Diff line change
@@ -1,77 +1,171 @@
import { Kilometers, Radians } from 'src/main';
import { Kilometers, Radians, Seconds } from '../main';
import { EpochUTC } from '../time/EpochUTC';
import { earthGravityParam, secondsPerDay, TAU } from '../utils/constants';
import { newtonM } from '../utils/functions';
import { ClassicalElements } from './ClassicalElements';
import { PositionVelocity } from './StateVector';
/** Equinoctial element set. */
import { EquinoctialElementsParams } from '../interfaces/EquinoctialElementsParams';

/**
* Equinoctial elements are a set of orbital elements used to describe the
* orbits of celestial bodies, such as satellites around a planet. They provide
* an alternative to the traditional Keplerian elements and are especially
* useful for avoiding singularities and numerical issues in certain types of
* orbits.
*
* Unlike Keplerian elements, equinoctial elements don't suffer from
* singularities at zero eccentricity (circular orbits) or zero inclination
* (equatorial orbits). This makes them more reliable for numerical simulations
* and analytical studies, especially in these edge cases.
*
* Reference: https://faculty.nps.edu/dad/orbital/th0.pdf
*/
export class EquinoctialElements {
epoch: EpochUTC;
/**
* The semi-major axis of the orbit in kilometers.
*/
a: Kilometers;
/**
* The h component of the eccentricity vector.
*/
h: number;
/**
* The k component of the eccentricity vector.
*/
k: number;
/**
* The p component of the ascending node vector.
*/
p: number;
/**
* The q component of the ascending node vector.
*/
q: number;
/**
* The mean longitude of the orbit in radians.
*/
lambda: Radians;
/**
* The gravitational parameter of the central body in km³/s².
*/
mu: number;
fr: number;
/** Create a new [EquinoctialElements] object from orbital elements. */
constructor(
public epoch: EpochUTC,
public af: number,
public ag: number,
public l: number,
public n: number,
public chi: number,
public psi: number,
{ mu = earthGravityParam, fr = 1 }: { mu?: number; fr?: number } = {},
) {
this.mu = mu;
this.fr = fr;
/**
* The retrograde factor. 1 for prograde orbits, -1 for retrograde orbits.
*/
I: 1 | -1;
constructor({ epoch, h, k, lambda, a, p, q, mu, I }: EquinoctialElementsParams) {
this.epoch = epoch;
this.h = h;
this.k = k;
this.lambda = lambda;
this.a = a;
this.p = p;
this.q = q;
this.mu = mu ?? earthGravityParam;
this.I = I ?? 1;
}

/**
* Returns a string representation of the EquinoctialElements object.
* @returns A string representation of the EquinoctialElements object.
*/
toString(): string {
return [
'[EquinoctialElements]',
` Epoch: ${this.epoch}`,
` af: ${this.af}`,
` ag: ${this.ag}`,
` l: ${this.l} rad`,
` n: ${this.n} rad/s`,
` chi: ${this.chi}`,
` psi: ${this.psi}`,
` a: ${this.a} km`,
` h: ${this.h}`,
` k: ${this.k}`,
` p: ${this.p}`,
` q: ${this.q}`,
` lambda: ${this.lambda} rad`,
].join('\n');
}

/** Return the orbit semimajor-axis _(km)_. */
semimajorAxis(): number {
return Math.cbrt(this.mu / (this.n * this.n));
/**
* Gets the semimajor axis.
* @returns The semimajor axis in kilometers.
*/
get semimajorAxis(): Kilometers {
return this.a;
}

/**
* Gets the mean longitude.
* @returns The mean longitude in radians.
*/
get meanLongitude(): Radians {
return this.lambda;
}

/**
* Calculates the mean motion of the celestial object.
*
* @returns The mean motion in units of radians per second.
*/
get meanMotion(): number {
return Math.sqrt(this.mu / this.a ** 3);
}

/**
* Gets the retrograde factor.
*
* @returns The retrograde factor.
*/
get retrogradeFactor(): number {
return this.I;
}

/**
* Checks if the orbit is prograde.
* @returns {boolean} True if the orbit is prograde, false otherwise.
*/
isPrograde(): boolean {
return this.I === 1;
}
/** Compute the mean motion _(rad/s)_ of this orbit. */
meanMotion(): number {
const a = this.semimajorAxis();

return Math.sqrt(this.mu / (a * a * a));
/**
* Checks if the orbit is retrograde.
* @returns {boolean} True if the orbit is retrograde, false otherwise.
*/
isRetrograde(): boolean {
return this.I === -1;
}

/** Compute the period _(seconds)_ of this orbit. */
period(): number {
return TAU * Math.sqrt(this.semimajorAxis() ** 3 / this.mu);
/**
* Gets the period of the orbit.
* @returns The period in seconds.
*/
get period(): Seconds {
return (TAU * Math.sqrt(this.semimajorAxis ** 3 / this.mu)) as Seconds;
}

// / Compute the number of revolutions completed per day for this orbit.
revsPerDay(): number {
return secondsPerDay / this.period();
/**
* Gets the number of revolutions per day.
*
* @returns The number of revolutions per day.
*/
get revsPerDay(): number {
return secondsPerDay / this.period;
}

// / Convert this to [ClassicalElements].
/**
* Converts the equinoctial elements to classical elements.
* @returns The classical elements.
*/
toClassicalElements(): ClassicalElements {
const a = this.semimajorAxis();
const e = Math.sqrt(this.af * this.af + this.ag * this.ag);
const i =
Math.PI * ((1.0 - this.fr) * 0.5) +
2.0 * this.fr * Math.atan(Math.sqrt(this.chi * this.chi + this.psi * this.psi));
const o = Math.atan2(this.chi, this.psi);
const w = Math.atan2(this.ag, this.af) - this.fr * Math.atan2(this.chi, this.psi);
const m = this.l - this.fr * o - w;
const a = this.semimajorAxis;
const e = Math.sqrt(this.k * this.k + this.h * this.h);
const i = Math.PI * ((1.0 - this.I) * 0.5) + 2.0 * this.I * Math.atan(Math.sqrt(this.p * this.p + this.q * this.q));
const o = Math.atan2(this.p, this.q);
const w = Math.atan2(this.h, this.k) - this.I * Math.atan2(this.p, this.q);
const m = this.lambda - this.I * o - w;
const v = newtonM(e, m).nu;

return new ClassicalElements({
epoch: this.epoch,
semimajorAxis: a as Kilometers,
semimajorAxis: a,
eccentricity: e,
inclination: i as Radians,
rightAscension: o as Radians,
Expand All @@ -81,7 +175,10 @@ export class EquinoctialElements {
});
}

/** Convert this to inertial position and velocity vectors. */
/**
* Converts the equinoctial elements to position and velocity.
* @returns The position and velocity in classical elements.
*/
toPositionVelocity(): PositionVelocity {
return this.toClassicalElements().toPositionVelocity();
}
Expand Down
3 changes: 2 additions & 1 deletion src/coordinate/J2000.ts
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,8 @@ export class J2000 extends StateVector {
}

/**
* Converts the coordinates from J2000 to the International Terrestrial Reference Frame (ITRF).
* Converts the coordinates from J2000 to the International Terrestrial
* Reference Frame (ITRF).
*
* This is an ECI to ECF transformation.
*/
Expand Down
6 changes: 5 additions & 1 deletion src/coordinate/StateVector.ts
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,11 @@ export type PositionVelocity = {

// / Base class for state vectors.
export abstract class StateVector {
constructor(public epoch: EpochUTC, public position: Vector3D, public velocity: Vector3D) {
constructor(
public epoch: EpochUTC,
public position: Vector3D,
public velocity: Vector3D,
) {
// Nothing to do here.
}

Expand Down
15 changes: 9 additions & 6 deletions src/coordinate/TEME.ts
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,16 @@ import { J2000 } from './J2000';
import { StateVector } from './StateVector';

/**
* True Equator Mean Equinox (TEME) is a coordinate system commonly used in satellite tracking and orbit prediction.
* It is a reference frame that defines the position and orientation of an object relative to the Earth's equator
* and equinox.
* True Equator Mean Equinox (TEME) is a coordinate system commonly used in
* satellite tracking and orbit prediction. It is a reference frame that defines
* the position and orientation of an object relative to the Earth's equator and
* equinox.
*
* By using the True Equator Mean Equinox (TEME) coordinate system, we can accurately describe the position and motion
* of satellites relative to the Earth's equator and equinox. This is particularly useful for tracking and predicting
* satellite orbits in various applications, such as satellite communication, navigation, and remote sensing.
* By using the True Equator Mean Equinox (TEME) coordinate system, we can
* accurately describe the position and motion of satellites relative to the
* Earth's equator and equinox. This is particularly useful for tracking and
* predicting satellite orbits in various applications, such as satellite
* communication, navigation, and remote sensing.
*/
export class TEME extends StateVector {
// / Create a new [TEME] object from a [ClassicalElements] object.
Expand Down
Loading

0 comments on commit 22ea0a8

Please sign in to comment.