I’m trying to use a FieldNumericalPropagator with Gradient to compute the derivative of the final state with respect to a force model parameter (so that I can check I’ve implemented the derivatives correctly when using NumericalPropagator). I saw using a FieldNumericalPropagator suggested for the purpose in [1]. But I can’t figure out how to treat the force model parameter as a variable. It seems that ForceModel.addContribution(FieldSpacecraftState, FieldTimeDerivativesEquationsAdder) assumes all the force model parameters are constants (Gradient.getZero().add(...)).
How should I set up the FieldNumericalPropagator to treat a force model parameter as a Gradient.variable(...)? I’m probably missing something obvious…
I’m not sure I understand what you’re trying to do. When you talk about the NumericalPropagator, do you mean the built-in STM computation? If so, the thing is that it uses the Fielded acceleration method of ForceModel, so I’m not sure how you can validate your code with it, because it would basically be comparing to itself. Sorry if I didn’t understand.
If you’re looking to check results, there’s always the finite differences I guess.
Thanks for the response. I guess I should have provided some example code. Here is an example, with the question:
Gradient cr = Gradient.variable(1, 0, 0);
// where do I use cr when setting up the FieldNumericalPropagator?
FieldNumericalPropagator p = ...;
RadiationSensitive spacecraft = new IsotropicRadiationSingleCoefficient(1, 0);
ForceModel srp = new SolarRadiationPressure(sun, earth, spacecraft);
p.addForceModel(srp);
SpacecraftState s = p.propagate(targetDate);
// here I'd like to extract d(s)/d(cr).
Hopefully the question makes sense now. Maybe it’s just a missing feature.
I’ve not tested it end to end, but here is something dirty yet light that would do the job:
package org.orekit.forces.radiation;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.analysis.differentiation.Gradient;
import org.hipparchus.analysis.differentiation.GradientField;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.util.MathArrays;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.forces.ForceModel;
import org.orekit.forces.ForceModelModifier;
import org.orekit.frames.FramesFactory;
import org.orekit.models.earth.ReferenceEllipsoid;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.utils.ParameterDriver;
import java.util.ArrayList;
import java.util.List;
public class CustomForce implements ForceModelModifier {
private final SolarRadiationPressure pressure;
private final IsotropicRadiationSingleCoefficient singleCoefficient;
public CustomForce(final IsotropicRadiationSingleCoefficient singleCoefficient) {
this.singleCoefficient = singleCoefficient;
this.pressure = new SolarRadiationPressure(CelestialBodyFactory.getSun(), ReferenceEllipsoid.getWgs84(FramesFactory.getGTOD(true)),
singleCoefficient);
}
@Override
public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
if (s.getMass().getField() instanceof GradientField) {
final double nominalCr = singleCoefficient.getRadiationParametersDrivers().get(1).getValue();
final Gradient cr = Gradient.variable(1, 0, nominalCr);
final Gradient[] overloadParameters = (Gradient[]) MathArrays.buildArray(s.getMass().getField(), 2);
overloadParameters[0] = Gradient.constant(1, singleCoefficient.getRadiationParametersDrivers().get(0).getValue());
overloadParameters[1] = cr;
final FieldSpacecraftState<Gradient> castState = (FieldSpacecraftState<Gradient>) s;
return (FieldVector3D<T>) pressure.acceleration(castState, overloadParameters);
}
return pressure.acceleration(s, parameters);
}
@Override
public List<ParameterDriver> getParametersDrivers() {
return new ArrayList<>();
}
@Override
public ForceModel getUnderlyingModel() {
return pressure;
}
}
Assuming you’ve given an input FieldSpacecraftState with only constant terms, I believe you’ll get at the end of propagation the requested derivatives.
Edit: on second thought, might be event shorter to override the getParameters method and introduce your variable there (would still involve some casting and unchecked warnings).
Thanks for the idea! I didn’t consider using the ForceModelModifier. I’ll also think a bit more if there is a good way to include this in either the base ForceModel or Propagator.