From ce406bd34778d07d3e955de39a351caac4d7f8ba Mon Sep 17 00:00:00 2001 From: Anthony Templeton Date: Thu, 25 Jul 2024 02:46:47 -0400 Subject: [PATCH] Feature/ans 39 (#24) * wcs implemented * add rotational rate for all COs --- README.md | 3 + examples/wcs.zig | 23 ++++++ src/Fits.zig | 2 +- src/Spacecraft.zig | 137 ++------------------------------ src/WorldCoordinateSystem.zig | 82 +++++++++++++++++++ src/calculations.zig | 143 ++++++++++++++++++++++++++++++++++ src/constants.zig | 26 ++++--- src/lib.zig | 1 + 8 files changed, 272 insertions(+), 145 deletions(-) create mode 100644 examples/wcs.zig create mode 100644 src/WorldCoordinateSystem.zig diff --git a/README.md b/README.md index e195cf4..91bec72 100644 --- a/README.md +++ b/README.md @@ -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 @@ -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) + [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 diff --git a/examples/wcs.zig b/examples/wcs.zig new file mode 100644 index 0000000..34de86b --- /dev/null +++ b/examples/wcs.zig @@ -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(); +} diff --git a/src/Fits.zig b/src/Fits.zig index 149063d..2305972 100644 --- a/src/Fits.zig +++ b/src/Fits.zig @@ -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, .{}); }, @@ -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 { diff --git a/src/Spacecraft.zig b/src/Spacecraft.zig index bbe3d5d..f331d4d 100644 --- a/src/Spacecraft.zig +++ b/src/Spacecraft.zig @@ -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; @@ -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| { @@ -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| { @@ -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 { diff --git a/src/WorldCoordinateSystem.zig b/src/WorldCoordinateSystem.zig new file mode 100644 index 0000000..a3cb14d --- /dev/null +++ b/src/WorldCoordinateSystem.zig @@ -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); +} diff --git a/src/calculations.zig b/src/calculations.zig index 4c4c933..df9ad19 100644 --- a/src/calculations.zig +++ b/src/calculations.zig @@ -1,4 +1,6 @@ const std = @import("std"); +const Tle = @import("Tle.zig"); +const OrbitalElements = @import("Spacecraft.zig").OrbitalElements; const constants = @import("constants.zig"); pub fn degreesToRadians(degrees: f64) f64 { @@ -20,3 +22,144 @@ pub fn meanMotionToSemiMajorAxis(m_motion: f64) f64 { 1.0 / 3.0, ); } + +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| { + result[i][j] = 0; + for (0..3) |k| { + result[i][j] += a[i][k] * b[k][j]; + } + } + } + return result; +} + +/// convert tle to orbital elements +pub fn tleToOrbitalElements(tle: Tle) OrbitalElements { + const inclination = degreesToRadians(tle.second_line.inclination); + const raan = degreesToRadians(tle.second_line.right_ascension); + const eccentricity = tle.second_line.eccentricity; + const argument_of_perigee = degreesToRadians(tle.second_line.perigee); + const m_anomaly = degreesToRadians(tle.second_line.m_anomaly); + const m_motion = meanMotionToRadiansPerMinute(tle.second_line.m_motion); + + const a = 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 dotProduct(a: [3]f64, b: [3]f64) f64 { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +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 magnitude(v: [3]f64) f64 { + return std.math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); +} + +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; +} diff --git a/src/constants.zig b/src/constants.zig index d02591b..512d1c1 100644 --- a/src/constants.zig +++ b/src/constants.zig @@ -29,8 +29,9 @@ pub const CelestialBody = struct { j2_perturbation: f64, sea_level_density: f64, scale_height: f64, + rotation_rate: f64, // rad/s - fn init(mass: f64, mu: f64, m_fraction_solar_system: f64, m_radius: ?f64, eq_radius: ?f64, p_radius: ?f64, semi_major_axis: ?f64, perihelion: ?f64, aphelion: ?f64, period: ?f64, velocity: ?f64, eccentricity: ?f64, inclination: ?f64, oblateness: ?f64, j2_perturbation: f64, sea_level_density: f64, scale_height: f64) CelestialBody { + fn init(mass: f64, mu: f64, m_fraction_solar_system: f64, m_radius: ?f64, eq_radius: ?f64, p_radius: ?f64, semi_major_axis: ?f64, perihelion: ?f64, aphelion: ?f64, period: ?f64, velocity: ?f64, eccentricity: ?f64, inclination: ?f64, oblateness: ?f64, j2_perturbation: f64, sea_level_density: f64, scale_height: f64, rotation_rate: f64) CelestialBody { return .{ .mass = mass, .mu = mu, @@ -49,21 +50,22 @@ pub const CelestialBody = struct { .j2_perturbation = j2_perturbation, .sea_level_density = sea_level_density, .scale_height = scale_height, + .rotation_rate = rotation_rate, }; } }; -pub const sun = CelestialBody.init(1.8842e30, 1.32712e11, 9.98657e-1, null, 695700, null, null, null, null, null, null, null, null, null, 0.0000002, 1e-12, 50000.0); -pub const mercury = CelestialBody.init(3.30101e23, 2.20319e4, 1.65789e-7, 2439.4, 2440.53, 2438.26, 5.79091e7, 4.60009e7, 6.98173e7, 87.97, 47.87, 0.20564, 7.01, 0.000, 0.00006, 1e-12, 200.0); -pub const venus = CelestialBody.init(4.86732e24, 3.24859e5, 2.44455e-6, 6051.8, 6051.8, 6051.8, 1.08209e8, 1.07477e8, 1.08940e8, 224.70, 35.02, 0.00676, 3.39, 0.000, 0.000027, 65.0, 15.9); -pub const earth = CelestialBody.init(5.97219e24, 3.98600e5, 2.99946e-6, 6371.0084, 6378.1366, 6356.7519, 1.49598e8, 1.47100e8, 1.52096e8, 365.26, 29.78, 0.01670, 0.00, 0.003353, 0.00108263, 1.225, 7.249); -pub const moon = CelestialBody.init(7.34581e22, 4.90280e3, 3.68934e-8, 1737.4, null, null, 3.83398e5, 3.62106e5, 4.04689e5, 27.18, 1.03, 0.05555, 23.71, 0.0012, 0.0002027, 5e-13, 100.0); -pub const mars = CelestialBody.init(6.41693e23, 4.28284e4, 3.22282e-7, 3389.50, 3396.19, 3376.20, 2.27939e8, 2.06645e8, 2.49233e8, 686.97, 24.13, 0.09342, 1.85, 0.00648, 0.001964, 0.020, 11.1); -pub const jupiter = CelestialBody.init(1.89852e27, 1.26713e8, 9.53510e-4, 69911, 71492, 66854, 7.78321e8, 7.40603e8, 8.16038e8, 4332.52, 13.06, 0.04846, 1.30, 0.06487, 0.014736, 0.16, 27.0); -pub const saturn = CelestialBody.init(5.68460e26, 3.79406e7, 2.85502e-4, 58232, 60268, 54364, 1.42910e9, 1.35096e9, 1.50724e9, 10783.05, 9.64, 0.05468, 2.49, 0.09796, 0.016298, 0.19, 59.5); -pub const uranus = CelestialBody.init(8.68192e25, 5.79456e6, 4.36039e-5, 25362, 25559, 24973, 2.87479e9, 2.73854e9, 3.01104e9, 30768.84, 6.79, 0.04739, 0.77, 0.02293, 0.003343, 0.42, 27.7); -pub const neptune = CelestialBody.init(1.02431e26, 6.83653e6, 5.14447e-5, 24622, 24764, 24341, 4.50489e9, 4.46384e9, 4.54594e9, 60357.05, 5.43, 0.00911, 1.77, 0.01708, 0.003411, 0.45, 19.7); -pub const pluto = CelestialBody.init(1.46158e22, 9.75500e2, 7.34061e-9, 1188.3, null, null, 5.91540e9, 4.44212e9, 7.38868e9, 90821.51, 4.74, 0.24906, 17.14, 0.002, 0.00039, 1e-6, 50.0); +pub const sun = CelestialBody.init(1.8842e30, 1.32712e11, 9.98657e-1, null, 695700, null, null, null, null, null, null, null, null, null, 0.0000002, 1e-12, 50000.0, 2.865e-6); +pub const mercury = CelestialBody.init(3.30101e23, 2.20319e4, 1.65789e-7, 2439.4, 2440.53, 2438.26, 5.79091e7, 4.60009e7, 6.98173e7, 87.97, 47.87, 0.20564, 7.01, 0.000, 0.00006, 1e-12, 200.0, 1.24e-6); +pub const venus = CelestialBody.init(4.86732e24, 3.24859e5, 2.44455e-6, 6051.8, 6051.8, 6051.8, 1.08209e8, 1.07477e8, 1.08940e8, 224.70, 35.02, 0.00676, 3.39, 0.000, 0.000027, 65.0, 15.9, -2.99e-7); +pub const earth = CelestialBody.init(5.97219e24, 3.98600e5, 2.99946e-6, 6371.0084, 6378.1366, 6356.7519, 1.49598e8, 1.47100e8, 1.52096e8, 365.26, 29.78, 0.01670, 0.00, 0.003353, 0.00108263, 1.225, 7.249, 7.2921150e-5); +pub const moon = CelestialBody.init(7.34581e22, 4.90280e3, 3.68934e-8, 1737.4, null, null, 3.83398e5, 3.62106e5, 4.04689e5, 27.18, 1.03, 0.05555, 23.71, 0.0012, 0.0002027, 5e-13, 100.0, 2.6617e-6); +pub const mars = CelestialBody.init(6.41693e23, 4.28284e4, 3.22282e-7, 3389.50, 3396.19, 3376.20, 2.27939e8, 2.06645e8, 2.49233e8, 686.97, 24.13, 0.09342, 1.85, 0.00648, 0.001964, 0.020, 11.1, 7.088e-5); +pub const jupiter = CelestialBody.init(1.89852e27, 1.26713e8, 9.53510e-4, 69911, 71492, 66854, 7.78321e8, 7.40603e8, 8.16038e8, 4332.52, 13.06, 0.04846, 1.30, 0.06487, 0.014736, 0.16, 27.0, 1.7585e-4); +pub const saturn = CelestialBody.init(5.68460e26, 3.79406e7, 2.85502e-4, 58232, 60268, 54364, 1.42910e9, 1.35096e9, 1.50724e9, 10783.05, 9.64, 0.05468, 2.49, 0.09796, 0.016298, 0.19, 59.5, 1.6378e-4); +pub const uranus = CelestialBody.init(8.68192e25, 5.79456e6, 4.36039e-5, 25362, 25559, 24973, 2.87479e9, 2.73854e9, 3.01104e9, 30768.84, 6.79, 0.04739, 0.77, 0.02293, 0.003343, 0.42, 27.7, -1.012e-4); +pub const neptune = CelestialBody.init(1.02431e26, 6.83653e6, 5.14447e-5, 24622, 24764, 24341, 4.50489e9, 4.46384e9, 4.54594e9, 60357.05, 5.43, 0.00911, 1.77, 0.01708, 0.003411, 0.45, 19.7, 1.083e-4); +pub const pluto = CelestialBody.init(1.46158e22, 9.75500e2, 7.34061e-9, 1188.3, null, null, 5.91540e9, 4.44212e9, 7.38868e9, 90821.51, 4.74, 0.24906, 17.14, 0.002, 0.00039, 1e-6, 50.0, -1.139e-5); test CelestialBody { try std.testing.expectEqual(5.97219e24, earth.mass); diff --git a/src/lib.zig b/src/lib.zig index 7c0cf4e..20dc56e 100644 --- a/src/lib.zig +++ b/src/lib.zig @@ -9,6 +9,7 @@ pub const Parser = @import("parsers.zig").Parser; pub const Spacecraft = @import("Spacecraft.zig"); pub const calculations = @import("calculations.zig"); pub const EquatorialCoordinateSystem = @import("EquatorialCoordinateSystem.zig"); +pub const WorldCoordinateSystem = @import("WorldCoordinateSystem.zig"); pub const Fits = @import("Fits.zig"); test {