Hi @MaximeJ and @petrus.hyvonen,
thank you so much for your replies, now that line works! And I am finally able to run the entire code.
To test the code, I firstly set the estimated values in the theoreticalEvaluation method to be equal to the observed values. As we imagine, this returns the orbit we give in input. However, this happens only if the number of measurements is <= 4. In fact, if I consider > 4 measurements I get an org.orekit.errors.OrekitException: unable to solve: singular problem.
If instead I do not put the estimated values to be equal to the observed ones (i.e. considering the true estimated values) I get either the aforementioned exception or a non existing orbit (e.g. a non valid value for the true anomaly). What could be the problem here? I have followed the Satellite Laser Ranging Orbit Determination written by @yzokras so I guess the problem is in my extension of the of the PythonAbstractMeasurement class…
Here it is the extended class: the difference between my class and Range is that the measurement consists in the distance between the satellite and the Solar System Barycentre (SSB) along a given unit vector n . Since I am considering an ICRF frame (centered on the SSB), I have defined my measurement as numpy.dot(sat_position, n) = n1*x + n2*y + n3*z i.e. the dot product between the satellite position and n. Notice that the vector of derivatives of the measurement with respect to x, y, z, vx, vy and vz is [n1, n2, n3, 0, 0, 0] i.e. the three velocity components do not impact on the “istantaneous” range measurement.
def __init__(self, date, observed, sigma, base_weight, satellite):
super().__init__(date, observed, sigma, base_weight, satellite)
def theoreticalEvaluation(self, iteration, evaluation, states):
state = states
pvaDS = super().getCoordinates(state, 0, 6)
n = [np.sin(90-22.0144)*np.cos(83.6292), np.sin(90-22.0144)*np.sin(83.6292), np.cos(90-22.0144)]
sat_position = [state.getPVCoordinates().getPosition().getX(), state.getPVCoordinates().getPosition().getY(), state.getPVCoordinates().getPosition().getZ()]
pSSB = FieldVector3D(Gradient.constant(6, 0.0), Gradient.constant(6, 0.0), Gradient.constant(6, 0.0))
vSSB = FieldVector3D(Gradient.constant(6, 0.0), Gradient.constant(6, 0.0), Gradient.constant(6, 0.0))
aSSB = FieldVector3D(Gradient.constant(6, 0.0), Gradient.constant(6, 0.0), Gradient.constant(6, 0.0))
pvaSSB = TimeStampedFieldPVCoordinates(state.getDate(), pSSB, vSSB, aSSB)
# CREATE THE ESTIMATEDMEASUREMENT OBJECT
new_range = NewRange(self.getDate(), float(self.getObservedValue()), float(self.getTheoreticalStandardDeviation()), float(self.getBaseWeight()), self.getSatellites())
estimated = EstimatedMeasurement(new_range, iteration, evaluation, [state], [pvaDS.toTimeStampedPVCoordinates(), pvaSSB.toTimeStampedPVCoordinates()])
# SET THE ESTIMATED VALUE FOR THE MEASUREMENT
# SET THE DERIVATIVE OF THE ESTIMATED MEASUREMENT WITH RESPECT TO THE CURRENT STATE (CARTESIAN)
derivative = [float(n), float(n), float(n), 0.0, 0.0, 0.0]
derivative_jarray = JArray('object')(1)
derivative_jarray = JArray('double')(derivative)
Do you see something strange in my code? Another thing I was thinking about is that maybe I do not have enough measurements (the sigma here is very large, between 1 and 3 km) or maybe I have to divide the range_weight variable (used in the laser ranging example ) by the number of observations? I haven’t understood if it is done automatically or not…
Thank you in advance.
P.S. I am giving the true orbit as a guess orbit in these first experiments…