Skip to content

Commit

Permalink
Feature/ans 39 (#24)
Browse files Browse the repository at this point in the history
* wcs implemented

* add rotational rate for all COs
  • Loading branch information
ATTron authored Jul 25, 2024
1 parent 96640ef commit ce406bd
Show file tree
Hide file tree
Showing 8 changed files with 272 additions and 145 deletions.
3 changes: 3 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
- [x] Orbital Details
- [x] Astronomical Coordinates
- [x] Equatorial Coordinate System
- [x] World Coordinate System
- [x] Astronomical Computation
- [x] Precession
- [x] Celestial Bodies
Expand Down Expand Up @@ -105,6 +106,8 @@ exe.root_module.addImport("astroz", astroz_mod);

- #### [Precess star to July 30, 2005](examples/precess_star.zig)

- #### [Calculate WCS values from a TLE](examples/wcs.zig)

<!-- MARKDOWN LINKS -->

[ci-shd]: https://img.shields.io/github/actions/workflow/status/ATTron/astroz/ci.yaml?branch=main&style=for-the-badge&logo=github&label=CI&labelColor=black
Expand Down
23 changes: 23 additions & 0 deletions examples/wcs.zig
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
const std = @import("std");
const astroz = @import("astroz");
const Tle = astroz.Tle;
const WorldCoordinateSystem = astroz.WorldCoordinateSystem;

pub fn main() !void {
var gpa = std.heap.GeneralPurposeAllocator(.{}){};
defer _ = gpa.deinit();
const allocator = gpa.allocator();

const test_tle =
\\1 55909U 23035B 24187.51050877 .00023579 00000+0 16099-2 0 9998
\\2 55909 43.9978 311.8012 0011446 278.6226 81.3336 15.05761711 71371
;

var tle = try Tle.parse(test_tle, allocator);
defer tle.deinit();
const wcs = WorldCoordinateSystem.fromTle(test_tle, 0.0);

std.log.debug("WCS OUTPUT: {any}", wcs);

tle.output();
}
2 changes: 1 addition & 1 deletion src/Fits.zig
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ pub fn open_and_parse(file_path: []const u8, allocator: std.mem.Allocator) !Fits
switch (hdu.content_type) {
.Image => {
image_count += 1;
std.log.debug("IMAGE FOUND NOW AT COUNT: {d}", .{image_count});
const filename = try std.fmt.bufPrint(&output_path_buffer, "image_{d}", .{image_count});
try fits_file.readImage(hdu.number, filename, .{});
},
Expand Down Expand Up @@ -289,7 +290,6 @@ fn applyStretch(self: *Fits, pixels: []f32, width: usize, height: usize, output_

var output_path_buffer: [std.fs.max_path_bytes]u8 = undefined;
if (output_path) |op| {
std.log.warn("\nOUTPUT PATH DEFINED! {s}\n", .{op});
const output = try std.fmt.bufPrint(&output_path_buffer, "{s}/{s}.png", .{ file_name, op });
try image.writeToFilePath(output, .{ .png = .{} });
} else {
Expand Down
137 changes: 5 additions & 132 deletions src/Spacecraft.zig
Original file line number Diff line number Diff line change
Expand Up @@ -148,8 +148,8 @@ pub fn propagateAttitude(self: *Spacecraft, dt: f64) void {
/// This will call the proper propagation methods based on a TLE epoch and recalculation time
/// Most of the functions this calls are private and will need code inspection to see
pub fn propagate(self: *Spacecraft, t0: f64, days: f64, h: f64, impulse_list: ?[]const Impulse) !void {
const y0_oe = self.tleToOrbitalElements();
const y0 = orbitalElementsToStateVector(y0_oe, self.orbiting_object.mu);
const y0_oe = calculations.tleToOrbitalElements(self.tle);
const y0 = calculations.orbitalElementsToStateVector(y0_oe, self.orbiting_object.mu);
var t = t0;
const tf = self.tle.first_line.epoch + days * 86400.0;
var y = y0;
Expand Down Expand Up @@ -296,7 +296,7 @@ fn cross(a: [3]f64, b: [3]f64) [3]f64 {
};
}

fn multiplyMatrices(a: [3][3]f64, b: [3][3]f64) [3][3]f64 {
pub fn multiplyMatrices(a: [3][3]f64, b: [3][3]f64) [3][3]f64 {
var result: [3][3]f64 = undefined;
for (0..3) |i| {
for (0..3) |j| {
Expand Down Expand Up @@ -429,22 +429,6 @@ fn vectorAdd(a: StateV, b: StateV) StateV {
return result;
}

fn crossProduct(a: [3]f64, b: [3]f64) [3]f64 {
return .{
a[1] * b[2] - a[2] * b[1],
a[2] * b[0] - a[0] * b[2],
a[0] * b[1] - a[1] * b[0],
};
}

fn dotProduct(a: [3]f64, b: [3]f64) f64 {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

fn magnitude(v: [3]f64) f64 {
return std.math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}

fn scalarMultiply(scalar: f64, vector: StateV) StateV {
var result: StateV = undefined;
for (0..6) |i| {
Expand Down Expand Up @@ -504,122 +488,11 @@ fn applyPlaneChange(self: *Spacecraft, y: [6]f64, plane_change: Impulse) [6]f64
const r = [3]f64{ y[0], y[1], y[2] };
const v = [3]f64{ y[3], y[4], y[5] };

var elements = stateVectorToOrbitalElements(r, v, mu);
var elements = calculations.stateVectorToOrbitalElements(r, v, mu);
elements.i += plane_change.plane_change.?.delta_inclination;
elements.raan += plane_change.plane_change.?.delta_raan;

return orbitalElementsToStateVector(elements, mu);
}

fn tleToOrbitalElements(self: Spacecraft) OrbitalElements {
const inclination = calculations.degreesToRadians(self.tle.second_line.inclination);
const raan = calculations.degreesToRadians(self.tle.second_line.right_ascension);
const eccentricity = self.tle.second_line.eccentricity;
const argument_of_perigee = calculations.degreesToRadians(self.tle.second_line.perigee);
const m_anomaly = calculations.degreesToRadians(self.tle.second_line.m_anomaly);
const m_motion = calculations.meanMotionToRadiansPerMinute(self.tle.second_line.m_motion);

const a = calculations.meanMotionToSemiMajorAxis(m_motion);

const E = solveKeplerEquation(m_anomaly, eccentricity);
const true_anomaly = 2.0 * std.math.atan(@sqrt((1 + eccentricity) / (1 - eccentricity)) * @tan(E / 2.0));

return .{
.a = a,
.e = eccentricity,
.i = inclination,
.raan = raan,
.arg_periapsis = argument_of_perigee,
.true_anomaly = true_anomaly,
};
}

/// convert orbital elements to a usable state vector
pub fn orbitalElementsToStateVector(elements: OrbitalElements, mu: f64) [6]f64 {
const p = elements.a * (1 - elements.e * elements.e);
const r = p / (1 + elements.e * @cos(elements.true_anomaly));

const r_orbital = [3]f64{
r * @cos(elements.true_anomaly),
r * @sin(elements.true_anomaly),
0,
};
const v_orbital = [3]f64{
-@sqrt(mu / p) * @sin(elements.true_anomaly),
@sqrt(mu / p) * (elements.e + @cos(elements.true_anomaly)),
0,
};

const cos_raan = @cos(elements.raan);
const sin_raan = @sin(elements.raan);
const cos_i = @cos(elements.i);
const sin_i = @sin(elements.i);
const cos_arg = @cos(elements.arg_periapsis);
const sin_arg = @sin(elements.arg_periapsis);

const rot = [3][3]f64{
.{ cos_raan * cos_arg - sin_raan * sin_arg * cos_i, -cos_raan * sin_arg - sin_raan * cos_arg * cos_i, sin_raan * sin_i },
.{ sin_raan * cos_arg + cos_raan * sin_arg * cos_i, -sin_raan * sin_arg + cos_raan * cos_arg * cos_i, -cos_raan * sin_i },
.{ sin_arg * sin_i, cos_arg * sin_i, cos_i },
};

var r_inertial: [3]f64 = undefined;
var v_inertial: [3]f64 = undefined;

r_inertial[0] = rot[0][0] * r_orbital[0] + rot[0][1] * r_orbital[1] + rot[0][2] * r_orbital[2];
r_inertial[1] = rot[1][0] * r_orbital[0] + rot[1][1] * r_orbital[1] + rot[1][2] * r_orbital[2];
r_inertial[2] = rot[2][0] * r_orbital[0] + rot[2][1] * r_orbital[1] + rot[2][2] * r_orbital[2];

v_inertial[0] = rot[0][0] * v_orbital[0] + rot[0][1] * v_orbital[1] + rot[0][2] * v_orbital[2];
v_inertial[1] = rot[1][0] * v_orbital[0] + rot[1][1] * v_orbital[1] + rot[1][2] * v_orbital[2];
v_inertial[2] = rot[2][0] * v_orbital[0] + rot[2][1] * v_orbital[1] + rot[2][2] * v_orbital[2];

return .{ r_inertial[0], r_inertial[1], r_inertial[2], v_inertial[0], v_inertial[1], v_inertial[2] };
}

pub fn stateVectorToOrbitalElements(r: [3]f64, v: [3]f64, mu: f64) OrbitalElements {
const r_mag = magnitude(r);
const v_mag = magnitude(v);
const h = crossProduct(r, v);
const h_mag = magnitude(h);
const n = crossProduct(.{ 0, 0, 1 }, h);
const n_mag = magnitude(n);

const e_vec = .{
(v_mag * v_mag - mu / r_mag) * r[0] / mu - dotProduct(r, v) * v[0] / mu,
(v_mag * v_mag - mu / r_mag) * r[1] / mu - dotProduct(r, v) * v[1] / mu,
(v_mag * v_mag - mu / r_mag) * r[2] / mu - dotProduct(r, v) * v[2] / mu,
};
const e = magnitude(e_vec);

const eps = v_mag * v_mag / 2.0 - mu / r_mag;
const a = if (@abs(e - 1.0) > 1e-10) -mu / (2.0 * eps) else std.math.inf(f64);
const i = std.math.acos(h[2] / h_mag);
const raan = std.math.atan2(n[1], n[0]);
const arg_periapsis = std.math.acos(dotProduct(n, e_vec) / (n_mag * e));
const true_anomaly = std.math.acos(dotProduct(e_vec, r) / (e * r_mag));

return .{
.a = a,
.e = e,
.i = i,
.raan = raan,
.arg_periapsis = if (r[2] >= 0) arg_periapsis else 2 * std.math.pi - arg_periapsis,
.true_anomaly = if (dotProduct(r, v) >= 0) true_anomaly else 2 * std.math.pi - true_anomaly,
};
}

fn solveKeplerEquation(m_anomaly: f64, eccentricity: f64) f64 {
var E = m_anomaly;
const tolerance: f64 = 1e-6;
var d: f64 = 1.0;

while (@abs(d) > tolerance) {
d = E - eccentricity * @sin(E) - m_anomaly;
E -= d / (1 - eccentricity * @cos(E));
}

return E;
return calculations.orbitalElementsToStateVector(elements, mu);
}

fn impulse(state: StateV, delta_v: [3]f64) StateV {
Expand Down
82 changes: 82 additions & 0 deletions src/WorldCoordinateSystem.zig
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
//! World Coordinate System is commonly used
const std = @import("std");
const constants = @import("constants.zig");
const calculations = @import("calculations.zig");
const Tle = @import("Tle.zig");
const Spacecraft = @import("Spacecraft.zig");

const Matrix3x3 = [3][3]f64;

const Vector3 = [3]f64;

const WorldCoordinateSystem = @This();

x: f64,
y: f64,
z: f64,

/// Currently this function assumes a fully parsed TLE already
pub fn fromTle(tle: Tle, t0: f64) WorldCoordinateSystem {
std.log.info("TLE PARSING, {}", .{tle});
const orbital_elements = calculations.tleToOrbitalElements(tle);

const eci = orbitalElementsToECI(orbital_elements);
const ecef = eciToECEF(eci, t0);

return .{ .x = ecef[0], .y = ecef[1], .z = ecef[2] };
}

fn orbitalElementsToECI(elements: Spacecraft.OrbitalElements) Vector3 {
const r = elements.a * (1 - elements.e * elements.e) / (1 + elements.e * @cos(elements.true_anomaly));
const x_orbit = r * @cos(elements.true_anomaly);
const y_orbit = r * @sin(elements.true_anomaly);

const R_w = Matrix3x3{
.{ @cos(elements.arg_periapsis), -@sin(elements.arg_periapsis), 0 },
.{ @sin(elements.arg_periapsis), @cos(elements.arg_periapsis), 0 },
.{ 0, 0, 1 },
};

const R_i = Matrix3x3{
.{ 1, 0, 0 },
.{ 0, @cos(elements.i), -@sin(elements.i) },
.{ 0, @sin(elements.i), @cos(elements.i) },
};

const R_o = Matrix3x3{
.{ @cos(elements.raan), -@sin(elements.raan), 0 },
.{ @sin(elements.raan), @cos(elements.raan), 0 },
.{ 0, 0, 1 },
};

const R = calculations.multiplyMatrices(R_o, calculations.multiplyMatrices(R_i, R_w));

return .{
R[0][0] * x_orbit + R[0][1] * y_orbit,
R[1][0] * x_orbit + R[1][1] * y_orbit,
R[2][0] * x_orbit + R[2][1] * y_orbit,
};
}

fn eciToECEF(eci: Vector3, time_since_epoch: f64) Vector3 {
const m = constants.earth.rotation_rate * time_since_epoch;
return .{
eci[0] * @cos(m) + eci[1] * @sin(m),
-eci[0] * @sin(m) + eci[1] * @cos(m),
eci[2],
};
}

test WorldCoordinateSystem {
const raw_tle =
\\1 55909U 23035B 24187.51050877 .00023579 00000+0 16099-2 0 9998
\\2 55909 43.9978 311.8012 0011446 278.6226 81.3336 15.05761711 71371
;
const expected_ecs: WorldCoordinateSystem = .{ .x = 4.628063569540487e3, .y = -5.164768842168279e3, .z = 7.220776206732921e0 };

var test_tle = try Tle.parse(raw_tle, std.testing.allocator);
defer test_tle.deinit();
const wcs = WorldCoordinateSystem.fromTle(test_tle, 0.0);

try std.testing.expectEqualDeep(expected_ecs, wcs);
}
Loading

0 comments on commit ce406bd

Please sign in to comment.