rameshvarun/starsystem.rs

## starsystem.rs
use na::Vector3;

use std::f64::consts::{PI};

const TICK_TIMESTEP: f64 = 86400.0;

const GM: f64 = 1.32712440042e20;

#[derive(Clone, Debug)]
pub struct OrbitalElements {
    pub eccentricity: f64,
    pub semimajorAxis: f64,
    pub inclination: f64,
    pub longitudeOfAscendingNode: f64,
    pub argumentOfPeriapsis: f64,
    pub initialMeanAnomaly: f64,
}

#[derive(Clone, Debug)]
pub struct Body {
    pub name: String,
    pub orbit: OrbitalElements,
}

impl OrbitalElements {
    pub fn mean_anomaly(&self, time: f64) -> f64 {
        return self.initialMeanAnomaly + time * (GM / self.semimajorAxis.powf(3.0)).sqrt();
    }

    pub fn eccentricity_anomaly(&self, time: f64) -> f64 {
        let mean_anomaly = self.mean_anomaly(time);
        let tolerance = 1e-14;

        // Newton-Raphson Iteration
        let mut solution = mean_anomaly;
        loop {
            let f = solution - self.eccentricity * f64::sin(solution) - mean_anomaly;
            let f_prime = 1.0 - self.eccentricity * f64::cos(solution);
            let next_solution = solution - f / f_prime;

            let delta = next_solution - solution;
            solution = next_solution;

            if (delta <= tolerance) {
                break;
            }
        }

        return solution;
    }

    pub fn true_anomaly(&self, time: f64) -> f64 {
        let eccentricity_anomaly = self.eccentricity_anomaly(time);
        return 2.0
            * f64::atan2(
                f64::sqrt(1.0 + self.eccentricity) * f64::sin(eccentricity_anomaly / 2.0),
                f64::sqrt(1.0 - self.eccentricity) * f64::cos(eccentricity_anomaly / 2.0),
            );
    }

    pub fn radius(&self, time: f64) -> f64 {
        let eccentricity_anomaly = self.eccentricity_anomaly(time);
        return self.semimajorAxis * (1.0 - self.eccentricity * f64::cos(eccentricity_anomaly));
    }

    pub fn orbital_frame_position(&self, time: f64) -> Vector3<f64> {
        let true_anomaly = self.true_anomaly(time);
        return self.radius(time) * Vector3::new(f64::cos(true_anomaly), f64::sin(true_anomaly), 0.0);
    }
}

#[derive(Clone, Debug)]
pub struct StarSystem {
    pub name: String,
    pub bodies: Vec<Body>,
    pub time: f64,
}

impl StarSystem {
    pub fn tick(&self) -> StarSystem {
        return StarSystem {
            name: self.name.clone(),
            bodies: self.bodies.clone(),
            time: self.time + TICK_TIMESTEP,
        };
    }
}
	use na::Vector3;

	use std::f64::consts::{PI};

	const TICK_TIMESTEP: f64 = 86400.0;

	const GM: f64 = 1.32712440042e20;

	#[derive(Clone, Debug)]
	pub struct OrbitalElements {
	pub eccentricity: f64,
	pub semimajorAxis: f64,
	pub inclination: f64,
	pub longitudeOfAscendingNode: f64,
	pub argumentOfPeriapsis: f64,
	pub initialMeanAnomaly: f64,
	}

	#[derive(Clone, Debug)]
	pub struct Body {
	pub name: String,
	pub orbit: OrbitalElements,
	}

	impl OrbitalElements {
	pub fn mean_anomaly(&self, time: f64) -> f64 {
	return self.initialMeanAnomaly + time * (GM / self.semimajorAxis.powf(3.0)).sqrt();
	}

	pub fn eccentricity_anomaly(&self, time: f64) -> f64 {
	let mean_anomaly = self.mean_anomaly(time);
	let tolerance = 1e-14;

	// Newton-Raphson Iteration
	let mut solution = mean_anomaly;
	loop {
	let f = solution - self.eccentricity * f64::sin(solution) - mean_anomaly;
	let f_prime = 1.0 - self.eccentricity * f64::cos(solution);
	let next_solution = solution - f / f_prime;

	let delta = next_solution - solution;
	solution = next_solution;

	if (delta <= tolerance) {
	break;
	}
	}

	return solution;
	}

	pub fn true_anomaly(&self, time: f64) -> f64 {
	let eccentricity_anomaly = self.eccentricity_anomaly(time);
	return 2.0
	* f64::atan2(
	f64::sqrt(1.0 + self.eccentricity) * f64::sin(eccentricity_anomaly / 2.0),
	f64::sqrt(1.0 - self.eccentricity) * f64::cos(eccentricity_anomaly / 2.0),
	);
	}

	pub fn radius(&self, time: f64) -> f64 {
	let eccentricity_anomaly = self.eccentricity_anomaly(time);
	return self.semimajorAxis * (1.0 - self.eccentricity * f64::cos(eccentricity_anomaly));
	}

	pub fn orbital_frame_position(&self, time: f64) -> Vector3<f64> {
	let true_anomaly = self.true_anomaly(time);
	return self.radius(time) * Vector3::new(f64::cos(true_anomaly), f64::sin(true_anomaly), 0.0);
	}
	}

	#[derive(Clone, Debug)]
	pub struct StarSystem {
	pub name: String,
	pub bodies: Vec<Body>,
	pub time: f64,
	}

	impl StarSystem {
	pub fn tick(&self) -> StarSystem {
	return StarSystem {
	name: self.name.clone(),
	bodies: self.bodies.clone(),
	time: self.time + TICK_TIMESTEP,
	};
	}
	}