hi @Vincent and @Serrof,
Thanks for your time and availability on this topic.
To answer chronologically to both of your questions:
-
@Serrof the covariance matrix is retrieved as-is from the CDM file. It’s triangularized by the screening provider.
-
@Vincent Yes indeed the probability of collision provided in the CDM from the above scenario is of 3.076069e-4.
-
The CDM is set on the CCSDS 508.0-B-1 standard where the probability is given as a value defined in the interval [0.0, 1.0].
I’ve noticed that the reference frame set within the extracted CDM files is the fixed ITRF one.
It causes the method compute() inputting a cdm file to crash when trying to build an orekit orbit instance requiring a non-inertial frame. See below the related code and thrown exception.
It might be interesting to set a checker with a frame transformation loop within the mentionned method.
// Patera2005 method
ShortTermEncounter2DPOCMethod method = new Patera2005();
String CDM_PATH = "collision-resources/ION_SCV8_vs_STARLINK_1233.cdm";
ClassLoader classloader = Thread.currentThread().getContextClassLoader();
URL inputUrl = classloader.getResource(CDM_PATH);
DataSource data = new DataSource(inputUrl.getPath());
Cdm cdm = new ParserBuilder().buildCdmParser().parseMessage(data);
double primaryRadius = 5.0;
double secondaryRadius = 5.0;
double PoC = method.compute(cdm, primaryRadius, secondaryRadius).getValue();
System.out.println("Computed probability with Patera2005: " + PoC);
Thrown exception:
Exception in thread “main” org.orekit.errors.OrekitIllegalArgumentException: non pseudo-inertial frame “ITRF-2020/CIO/2010-based ITRF simple EOP”
at org.orekit.orbits.Orbit.ensurePseudoInertialFrame(Orbit.java:205)
at org.orekit.orbits.Orbit.(Orbit.java:145)
at org.orekit.orbits.CartesianOrbit.(CartesianOrbit.java:100)
at org.orekit.orbits.CartesianOrbit.(CartesianOrbit.java:123)
at org.orekit.ssa.collision.shorttermencounter.probability.twod.AbstractShortTermEncounter2DPOCMethod.getPrimaryFromCdm(AbstractShortTermEncounter2DPOCMethod.java:242)
at org.orekit.ssa.collision.shorttermencounter.probability.twod.AbstractShortTermEncounter2DPOCMethod.compute(AbstractShortTermEncounter2DPOCMethod.java:95)
at org.orekit.ssa.collision.shorttermencounter.probability.twod.ShortTermEncounter2DPOCMethod.compute(ShortTermEncounter2DPOCMethod.java:64)
at org.orekit.CollisionRiskProbabilityComputationWithCDM.main(CollisionRiskProbabilityComputationWithCDM.java:28)
I confirm that defining the covariance matrices with the QSW_INERTIAL local orbital frame does increase the probability of collision to some order of magnitudes close to the reference PoC. However it doesn’t seem to make sense to virtually set it referenced to a local orbital frame when they are actually defined within the ITRF one.
I’ve decided to project the input PV coordinates and covariance matrices from the CDM (defined in the ITRF frame) respectively to the GCRF frame (allowing the definitition of satellite orbits) and the QSW_INERTIAL local orbital frame. The results are very encouraging:
- Computed probability with Laas2015: 0.001483649317441629
- Computed probability with Alfano2005: 0.0014836595274094796
- Computed probability with Alfriend1999: 0.0015020406713417813
- Computed probability with Alfriend1999Max: 0.19427013238554633
- Computed probability with Chan1997: 0.0015009131774305908
- Computed probability with Patera2005: 0.0014836493174408775
When the reference PoC provided within the CDM is 0.004450713 (Foster1992).
The order of magnitude are matching apart from Alfriend1999Max.
The residual deviation might be inherant to each model algorithms…
Here’s the reference CDM file retrieved from recent conjunctions notified by CSpOC.
ION_SCV8_vs_STARLINK_1233.txt (12.9 KB)
Let me know what do you think about it.
You can find below the related script:
// initial setup
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
Frame inertialFrame = FramesFactory.getGCRF();
LOFType LOF = LOFType.QSW_INERTIAL;
LofOffset lofOffset = new LofOffset(inertialFrame, LOF);
TimeScale timeScale = TimeScalesFactory.getUTC();
AbsoluteDate TCA = new AbsoluteDate(2023, 7, 5, 20, 31, 15.893, timeScale);
double mu = Constants.IERS2010_EARTH_MU;
double[][] jacobian = new double[6][6];
// --- define primary object --- //
// build initial state
double xPrimary = -5719153.201;
double yPrimary = -2486155.271;
double zPrimary= -3021252.701;
double vXPrimary = 2333.174842 ;
double vYPrimary = 2825.732323 ;
double vZPrimary = -6727.808538;
Vector3D positionPrimary = new Vector3D(xPrimary, yPrimary, zPrimary);
Vector3D velocityPrimary = new Vector3D(vXPrimary, vYPrimary, vZPrimary);
PVCoordinates pvaCoordinatesPrimary = new PVCoordinates(positionPrimary, velocityPrimary);
AbsolutePVCoordinates absPVACoordinatesPrimary = new AbsolutePVCoordinates(itrf, TCA, pvaCoordinatesPrimary);
Attitude primaryAttitude = lofOffset.getAttitude(absPVACoordinatesPrimary, TCA, itrf);
SpacecraftState primaryState = new SpacecraftState(absPVACoordinatesPrimary, primaryAttitude);
Transform primaryTransform = primaryState.toTransform();
primaryTransform.getJacobian(CartesianDerivativesFilter.USE_PV, jacobian);
RealMatrix primarydC2dC1 = MatrixUtils.createRealMatrix(jacobian);
// setup the covariance
double CR_R = 127.5401258877026 ;
double CT_R = -216.5469301895778;
double CT_T = 9738.762752029239;
double CN_R = 23.41184551274494;
double CN_T = -5.346171445289681;
double CN_N = 24.60870138594973;
double CRDOT_R = 0.2482739312599897;
double CRDOT_T = -10.55496325674788 ;
double CRDOT_N = 0.005950144116086071;
double CRDOT_RDOT = 0.01147849769311488 ;
double CTDOT_R = -0.1375882369099475 ;
double CTDOT_T = 0.1222189661319367 ;
double CTDOT_N = -0.02576925085939057;
double CTDOT_RDOT = -0.0001458466399925159;
double CTDOT_TDOT = 0.0001499094267380701;
double CNDOT_R = 0.008141795305978281;
double CNDOT_T = -0.05578838789977865;
double CNDOT_N = -0.0003344474012710994;
double CNDOT_RDOT = 0.00005801069817622688;
double CNDOT_TDOT = -0.000008757429553163563;
double CNDOT_DNOT = 0.00009287778622823746;
double[][] dPrimary = {{CR_R, 0, 0, 0, 0, 0},
{CT_R, CT_T, 0, 0, 0, 0},
{CN_R, CN_T, CN_N, 0, 0, 0},
{CRDOT_R, CRDOT_T, CRDOT_N, CRDOT_RDOT, 0, 0},
{CTDOT_R, CTDOT_T, CTDOT_N, CTDOT_RDOT, CTDOT_TDOT, 0},
{CNDOT_R, CNDOT_T, CNDOT_N, CNDOT_RDOT, CNDOT_TDOT, CNDOT_DNOT}};
Array2DRowRealMatrix orbitalCovariancePrimary = new Array2DRowRealMatrix(dPrimary);
OrbitType orbitType = OrbitType.CARTESIAN;
PositionAngle angleType = PositionAngle.TRUE;
// --- define secondary object --- //
// build initial state
double xSecondary = -5719163.147;
double ySecondary = -2486109.2;
double zSecondary = -3021222.87;
double vXSecondary = -415.96327;
double vYSecondary= -5206.88041;
double vZSecondary = 5081.948896;
Vector3D positionSecondary = new Vector3D(xSecondary, ySecondary, zSecondary);
Vector3D velocitySecondary = new Vector3D(vXSecondary, vYSecondary, vZSecondary);
PVCoordinates pvaCoordinatesSecondary = new PVCoordinates(positionSecondary, velocitySecondary);
AbsolutePVCoordinates absPVACoordinatesSecondary = new AbsolutePVCoordinates(itrf, TCA, pvaCoordinatesSecondary);
Attitude secondaryAttitude = lofOffset.getAttitude(absPVACoordinatesSecondary, TCA, itrf);
SpacecraftState secondaryState = new SpacecraftState(absPVACoordinatesSecondary, secondaryAttitude);
Transform secondaryTransform = secondaryState.toTransform();
secondaryTransform.getJacobian(CartesianDerivativesFilter.USE_PV, jacobian);
RealMatrix secondarydC2dC1 = MatrixUtils.createRealMatrix(jacobian);
// setup the covariance
CR_R = 964.6447977021088;
CT_R = -29096.52093664356;
CT_T = 1555885.738355947;
CN_R = -610.4622634358643;
CN_T = 7105.637820746857;
CN_N = 1325.505208766663;
CRDOT_R = 31.55505053547703;
CRDOT_T = -1687.047343882137;
CRDOT_N = -7.503566951636299;
CRDOT_RDOT = 1.829454352163062;
CTDOT_R = -0.7822616662183071;
CTDOT_T = 16.88006830847172;
CTDOT_N = 0.627700053500938;
CTDOT_RDOT = -0.01829894383644205;
CTDOT_TDOT = 0.000706074549697974;
CNDOT_R = 0.004835998593925756;
CNDOT_T = -0.1102870161755583;
CNDOT_N = -0.2277988353274012;
CNDOT_RDOT = 0.00005700360443949522;
CNDOT_TDOT = -0.00001234324761245689;
CNDOT_DNOT = 0.0003799085974184537;
double[][] dSecondary = {{CR_R, 0, 0, 0, 0, 0},
{CT_R, CT_T, 0, 0, 0, 0},
{CN_R, CN_T, CN_N, 0, 0, 0},
{CRDOT_R, CRDOT_T, CRDOT_N, CRDOT_RDOT, 0, 0},
{CTDOT_R, CTDOT_T, CTDOT_N, CTDOT_RDOT, CTDOT_TDOT, 0},
{CNDOT_R, CNDOT_T, CNDOT_N, CNDOT_RDOT, CNDOT_TDOT, CNDOT_DNOT}};
Array2DRowRealMatrix orbitalCovarianceSecondary = new Array2DRowRealMatrix(dSecondary);
// transform Covariances from ITRF to defined LOF
RealMatrix transformedPrimaryCovarianceRealMatrix = primarydC2dC1.multiply(orbitalCovariancePrimary).multiplyTransposed(primarydC2dC1);
RealMatrix transformedSecondaryCovarianceRealMatrix = secondarydC2dC1.multiply(orbitalCovarianceSecondary).multiplyTransposed(secondarydC2dC1);
StateCovariance primaryTransformedCovariance = new StateCovariance(transformedPrimaryCovarianceRealMatrix, TCA, LOF);
StateCovariance secondaryTransformedCovariance = new StateCovariance(transformedSecondaryCovarianceRealMatrix, TCA, LOF);
// transform PVCoordinates from ITRF to GCRF and build associated orbits
Transform frameTransform = itrf.getTransformTo(inertialFrame, TCA);
PVCoordinates transformedPvCoordinatesPrimary = frameTransform.transformPVCoordinates(pvaCoordinatesPrimary);
PVCoordinates transformedPvCoordinatesSecondary = frameTransform.transformPVCoordinates(pvaCoordinatesSecondary);
Orbit orbitPrimary = new CartesianOrbit(transformedPvCoordinatesPrimary, inertialFrame, TCA, mu);
Orbit orbitSecondary = new CartesianOrbit(transformedPvCoordinatesSecondary, inertialFrame, TCA, mu);
// compute collision probability
double primaryRadius = 5.000000;
double secondaryRadius = 5.000000;
double combinedRadius = primaryRadius + secondaryRadius;
// Laas2015 method
ShortTermEncounter2DPOCMethod method = new Laas2015();
ProbabilityOfCollision probabilityOfCollision = method.compute(orbitPrimary, primaryTransformedCovariance, orbitSecondary, secondaryTransformedCovariance, combinedRadius, 1e-10);
if (probabilityOfCollision.getValue() == 0.004450713) {
System.out.println("It matches with Laas2015!");
} else {
System.out.println("Computed probability with Laas2015: " + probabilityOfCollision.getValue());
}
// Alfano2005 method
method = new Alfano2005();
probabilityOfCollision = method.compute(orbitPrimary, primaryTransformedCovariance, orbitSecondary, secondaryTransformedCovariance, combinedRadius, 1e-10);
if (probabilityOfCollision.getValue() == 0.004450713) {
System.out.println("It matches with Alfano2005!");
} else {
System.out.println("Computed probability with Alfano2005: " + probabilityOfCollision.getValue());
}
// Alfriend1999 method
method = new Alfriend1999();
probabilityOfCollision = method.compute(orbitPrimary, primaryTransformedCovariance, orbitSecondary, secondaryTransformedCovariance, combinedRadius, 1e-10);
if (probabilityOfCollision.getValue() == 0.004450713) {
System.out.println("It matches with Alfriend1999!");
} else {
System.out.println("Computed probability with Alfriend1999: " + probabilityOfCollision.getValue());
}
// Alfriend1999Max method
method = new Alfriend1999Max();
probabilityOfCollision = method.compute(orbitPrimary, primaryTransformedCovariance, orbitSecondary, secondaryTransformedCovariance, combinedRadius, 1e-10);
if (probabilityOfCollision.getValue() == 0.004450713) {
System.out.println("It matches with Alfriend1999Max!");
} else {
System.out.println("Computed probability with Alfriend1999Max: " + probabilityOfCollision.getValue());
}
// Chan1997 method
method = new Chan1997();
probabilityOfCollision = method.compute(orbitPrimary, primaryTransformedCovariance, orbitSecondary, secondaryTransformedCovariance, combinedRadius, 1e-10);
if (probabilityOfCollision.getValue() == 0.004450713) {
System.out.println("It matches with Chan1997!");
} else {
System.out.println("Computed probability with Chan1997: " + probabilityOfCollision.getValue());
}
// Patera2005 method
method = new Patera2005();
probabilityOfCollision = method.compute(orbitPrimary, primaryTransformedCovariance, orbitSecondary, secondaryTransformedCovariance, combinedRadius, 1e-10);
if (probabilityOfCollision.getValue() == 0.004450713) {
System.out.println("It matches with Patera2005!");
} else {
System.out.println("Computed probability with Patera2005: " + probabilityOfCollision.getValue());
}
}
.
!
.
) ?
.