Hello everyone. 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();
}
}