Hello everyone.
I visulize two potentially colliding satellites at the VTS. I see them get close at TCA and I want to perform maneuver to avoid potential collision.
How can I simulate a manuver with NumericalPropagator and ConstantThrustManeuver force model in Ore-kit?
And how can I calculate required duration of maneuver?(if needed)
Hello @blacksmith16,
How can I simulate a manuver with NumericalPropagator and ConstantThrustManeuver force model in Ore-kit?
To this question, i suggest you to look at the ApogeeManeuver tutorial which will explain several mechanism involved in what you are looking to do such as:
NumericalPropagatorExtremumApproachDetector and an EventShifter for example so you can trigger the maneuver one half period before the TCA)And how can I calculate required duration of maneuver?(if needed)
A simple approach (especially with high thrust propulsion systems) is to consider putting a tangential maneuver one half-period before the estimated TCA as explained earlier and do some dichotomy on the desired miss distance / probability of collision / some other metric to find the right maneuver duration.
When considering low-thrust propulsion systems, one maneuver may not be enough and more complex algorithm must be developped which is outside of the forum’s help scope i guess.
Hope this helps !
Cheers,
Vincent
My goal is to take TLE data of satellites, calculate TLE MissDÄ°stance DeltaV then calculate a duration of maneuver to make the distance over threshold. after that i tried to implement a maneuver to one of the satellite. finally getting oem data for both satellite and one more for the one that the maneuver implemented so that i wish to see 2 satllites getting closer at tca and one imarginary satllite is maneuvering to avoid collision.
i dont get right dates and x y z positions in oem
public class Main {
public static void main(String[] args) {
try {
// Initialize Orekit
InitializeOrekit.init();
// Define TLE data for two satellites
String tleLine1Sat1 = "1 49288U 21090K 24143.89033816 .00000069 00000+0 14085-3 0 9990";
String tleLine2Sat1 = "2 49288 87.9140 95.2966 0002435 83.3835 276.7575 13.19724491127294";
TLE tleSat1 = new TLE(tleLine1Sat1, tleLine2Sat1);
String tleLine1Sat2 = "1 45318U 91056BG 24143.88247308 .00000567 00000+0 13992-2 0 9999";
String tleLine2Sat2 = "2 45318 82.5682 312.0663 0010489 174.5580 185.5659 13.19066008203960";
TLE tleSat2 = new TLE(tleLine1Sat2, tleLine2Sat2);
// Create TLE propagators
TLEPropagator propagatorSat1 = TLEPropagator.selectExtrapolator(tleSat1);
TLEPropagator propagatorSat2 = TLEPropagator.selectExtrapolator(tleSat2);
// Find TCA and miss distance
AbsoluteDate startDate = tleSat1.getDate();
AbsoluteDate endDate = startDate.shiftedBy(360 * 3600); // 15 days
double minDistance = Double.MAX_VALUE;
AbsoluteDate tcaDate = null;
AbsoluteDate currentDate = startDate;
while (currentDate.compareTo(endDate) <= 0) {
SpacecraftState stateSat1 = propagatorSat1.propagate(currentDate);
SpacecraftState stateSat2 = propagatorSat2.propagate(currentDate);
Vector3D positionSat1 = stateSat1.getPVCoordinates().getPosition();
Vector3D positionSat2 = stateSat2.getPVCoordinates().getPosition();
double distance = Vector3D.distance(positionSat2, positionSat1);
if (distance < minDistance) {
minDistance = distance;
tcaDate = currentDate;
}
currentDate = currentDate.shiftedBy(0.02); // Advance time by 60 seconds
}
if (tcaDate != null) {
System.out.println("Closest Approach Time (TCA): " + tcaDate);
System.out.println("Minimum Distance (km): " + minDistance / 1000);
SpacecraftState stateAtTcaSat1 = propagatorSat1.propagate(tcaDate);
SpacecraftState stateAtTcaSat2 = propagatorSat2.propagate(tcaDate);
Vector3D velocitySat1 = stateAtTcaSat1.getPVCoordinates().getVelocity();
Vector3D velocitySat2 = stateAtTcaSat2.getPVCoordinates().getVelocity();
Vector3D deltaV = velocitySat2.subtract(velocitySat1);
System.out.println("Delta-v at TCA (km/s): " + deltaV.getNorm() / 1000);
}
if (minDistance < 1000) {
double thrustForce = 0.1; // N
double specificImpulse = 300; // s
double satelliteMass = 500; // kg
double desiredMissDistance = 1000; // m
double initialGuess = 5; // s
double thrustDuration = calculateThrustDuration(propagatorSat1, propagatorSat2, tcaDate, thrustForce, specificImpulse, satelliteMass, desiredMissDistance, initialGuess);
System.out.println("Calculated Thrust Duration (s): " + thrustDuration);
SpacecraftState finalStateWithManeuver = applyManeuver(propagatorSat1, tcaDate, thrustDuration, thrustForce, specificImpulse, satelliteMass);
NumericalPropagator numericalPropagator = createNumericalPropagator(finalStateWithManeuver);
SpacecraftState finalStateSat1 = numericalPropagator.propagate(tcaDate);
SpacecraftState finalStateSat2 = propagatorSat2.propagate(tcaDate);
double finalDistance = Vector3D.distance(finalStateSat1.getPVCoordinates().getPosition(), finalStateSat2.getPVCoordinates().getPosition());
System.out.println("Final Miss Distance after Maneuver (m): " + finalDistance);
String outputDirectory = "C:\\Users\\emirhan\\Downloads\\";
AbsoluteDate oemStartDate = AbsoluteDate.createMJDDate(60456, 40200 * 86400.0 / 100000.0, TimeScalesFactory.getUTC());
AbsoluteDate oemEndDate = AbsoluteDate.createMJDDate(60456, 41400 * 86400.0 / 100000.0, TimeScalesFactory.getUTC());
createOEMFile(outputDirectory + "initial_sat1.txt", propagatorSat1, oemStartDate, oemEndDate);
createOEMFile(outputDirectory + "initial_sat2.txt", propagatorSat2, oemStartDate, oemEndDate);
createOEMFile(outputDirectory + "updated_sat1.txt", numericalPropagator, oemStartDate, oemEndDate);
} else {
System.out.println("No maneuver needed. Minimum distance is already greater than 1000 meters.");
}
} catch (Exception e) {
e.printStackTrace();
}
}
private static double calculateThrustDuration(TLEPropagator propagatorSat1, TLEPropagator propagatorSat2, AbsoluteDate tcaDate,
double thrustForce, double specificImpulse, double satelliteMass, double desiredMissDistance, double initialGuess) throws OrekitException {
double G = Constants.G0_STANDARD_GRAVITY;
double tolerance = 1e-6;
double maxIterations = 50;
double h = 1e-6;
double thrustDuration = initialGuess;
for (int i = 0; i < maxIterations; i++) {
double missDistance = simulateManeuverAndGetDistance(propagatorSat1, propagatorSat2, tcaDate, thrustDuration, thrustForce, specificImpulse, satelliteMass);
double missDistancePlusH = simulateManeuverAndGetDistance(propagatorSat1, propagatorSat2, tcaDate, thrustDuration + h, thrustForce, specificImpulse, satelliteMass);
double G_A = missDistance - desiredMissDistance;
double G_A_plus_h = missDistancePlusH - desiredMissDistance;
double G_prime_A = (G_A_plus_h - G_A) / h;
double nextThrustDuration = thrustDuration - G_A / G_prime_A;
if (Math.abs(nextThrustDuration - thrustDuration) < tolerance) {
return nextThrustDuration;
}
thrustDuration = nextThrustDuration;
}
return thrustDuration;
}
private static double simulateManeuverAndGetDistance(TLEPropagator propagatorSat1, TLEPropagator propagatorSat2, AbsoluteDate tcaDate,
double thrustDuration, double thrustForce, double specificImpulse, double satelliteMass) throws OrekitException {
SpacecraftState finalStateWithManeuver = applyManeuver(propagatorSat1, tcaDate, thrustDuration, thrustForce, specificImpulse, satelliteMass);
NumericalPropagator numericalPropagator = createNumericalPropagator(finalStateWithManeuver);
SpacecraftState finalStateSat1 = numericalPropagator.propagate(tcaDate);
SpacecraftState finalStateSat2 = propagatorSat2.propagate(tcaDate);
return Vector3D.distance(finalStateSat1.getPVCoordinates().getPosition(), finalStateSat2.getPVCoordinates().getPosition());
}
private static SpacecraftState applyManeuver(TLEPropagator propagatorSat1, AbsoluteDate tcaDate, double thrustDuration,
double thrustForce, double specificImpulse, double satelliteMass) throws OrekitException {
SpacecraftState initialState = propagatorSat1.propagate(tcaDate.shiftedBy(-thrustDuration));
NumericalPropagator numericalPropagator = createNumericalPropagator(initialState);
double isp = specificImpulse;
double flowRate = thrustForce / (isp * Constants.G0_STANDARD_GRAVITY);
Vector3D direction = initialState.getPVCoordinates().getVelocity().normalize();
ConstantThrustManeuver maneuver = new ConstantThrustManeuver(initialState.getDate(), thrustDuration, thrustForce, isp, direction);
numericalPropagator.addForceModel(maneuver);
return numericalPropagator.propagate(tcaDate);
}
private static NumericalPropagator createNumericalPropagator(SpacecraftState initialState) throws OrekitException {
double minStep = 0.001;
double maxStep = 1000.0;
double positionTolerance = 10.0;
OrbitType propagationType = OrbitType.CARTESIAN;
double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, initialState.getOrbit(), propagationType);
AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(propagationType);
propagator.setInitialState(initialState);
propagator.setMu(Constants.EIGEN5C_EARTH_MU);
return propagator;
}
private static void createOEMFile(String fileName, TLEPropagator propagator, AbsoluteDate startDate, AbsoluteDate endDate) throws IOException, OrekitException {
File file = new File(fileName);
file.getParentFile().mkdirs(); // Create directories if they do not exist
FileWriter writer = new FileWriter(file);
writer.write("CIC_OEM_VERS = 2.0\n");
writer.write("CREATION_DATE = " + AbsoluteDate.J2000_EPOCH + "\n");
writer.write("ORIGINATOR = Orekit\n");
writer.write("META_START\n");
writer.write("OBJECT_NAME = Satellite\n");
writer.write("OBJECT_ID = SatID\n");
writer.write("CENTER_NAME = Earth\n");
writer.write("REF_FRAME = EME2000\n");
writer.write("TIME_SYSTEM = UTC\n");
writer.write("META_STOP\n");
AbsoluteDate currentDate = startDate;
while (currentDate.compareTo(endDate) <= 0) {
SpacecraftState state = propagator.propagate(currentDate);
Vector3D position = state.getPVCoordinates().getPosition();
Vector3D velocity = state.getPVCoordinates().getVelocity();
writer.write(String.format(Locale.US, "%.5f %.10f %.3f %.3f %.3f %.6f %.6f %.6f\n",
currentDate.durationFrom(AbsoluteDate.J2000_EPOCH), currentDate.getComponents(TimeScalesFactory.getUTC()).getTime().getSecondsInUTCDay() / 86400.0,
position.getX(), position.getY(), position.getZ(),
velocity.getX(), velocity.getY(), velocity.getZ()));
currentDate = currentDate.shiftedBy(60);
}
writer.write("DATA_STOP\n");
writer.close();
}
private static void createOEMFile(String fileName, NumericalPropagator propagator, AbsoluteDate startDate, AbsoluteDate endDate) throws IOException, OrekitException {
File file = new File(fileName);
file.getParentFile().mkdirs(); // Create directories if they do not exist
FileWriter writer = new FileWriter(file);
writer.write("CIC_OEM_VERS = 2.0\n");
writer.write("CREATION_DATE = " + AbsoluteDate.J2000_EPOCH + "\n");
writer.write("ORIGINATOR = Orekit\n");
writer.write("META_START\n");
writer.write("OBJECT_NAME = Satellite\n");
writer.write("OBJECT_ID = SatID\n");
writer.write("CENTER_NAME = Earth\n");
writer.write("REF_FRAME = EME2000\n");
writer.write("TIME_SYSTEM = UTC\n");
writer.write("META_STOP\n");
AbsoluteDate currentDate = startDate;
while (currentDate.compareTo(endDate) <= 0) {
SpacecraftState state = propagator.propagate(currentDate);
Vector3D position = state.getPVCoordinates().getPosition();
Vector3D velocity = state.getPVCoordinates().getVelocity();
writer.write(String.format(Locale.US, "%.5f %.10f %.3f %.3f %.3f %.6f %.6f %.6f\n",
currentDate.durationFrom(AbsoluteDate.J2000_EPOCH), currentDate.getComponents(TimeScalesFactory.getUTC()).getTime().getSecondsInUTCDay() / 86400.0,
position.getX(), position.getY(), position.getZ(),
velocity.getX(), velocity.getY(), velocity.getZ()));
currentDate = currentDate.shiftedBy(60);
}
writer.write("DATA_STOP\n");
writer.close();
}
}
Instead of using a loop with a short time step to find the time of closest approach, I suggest to use a single propagation on one satellite with an ExtremumApproachDetector configured to use the other satellite as the secondary provider.
Also computing the velocity difference at time of closest approach is not the way you will get a collision avoidance: it is too late you are already there and the velocity difference may be really huge. You should perform a maneuver long before the time of closest approach.
There is a full ssa package that was added recently to the library, It mainly deals with computing probability, though, I don’t think it computes avoidance maneuvers.
How can i implement a maneuver?
I want to visualize maneuvered sat in VTS?