Hello all,

I’ve got some questions (for Luc I guess…) about the CelestialBodyPointed class in attitudes package.

First, is there somewhere a document describing the underlying maths ? I would be interested in it.

Also, could you explain more the “compensation” computation hereafter:

Blockquote
final Vector3D v1 = Vector3D.crossProduct(rotAxisCel, phasingCel);
final Vector3D v2 = Vector3D.crossProduct(pointingP, phasingCel);
final double compensation = -Vector3D.dotProduct(v1, v2) / v2.getNormSq();
final Vector3D phasedRotAxisCel = new Vector3D(1.0, rotAxisCel, compensation, pointingP);
Blockquote

I would like to implement the same kind of attitude, but with a phasing reference that is not constant in the celestial frame (actually constant in the local orbital frame). How the rotation rate would be impacted ?
Thanks

There is no reference document for this, I’m sorry. I simply came up to these equations by myself (probably after some head banging against a wall) and used finite differences to validate everything. The first part, up to computation of `rotAxisCel` before the code snippet is trivial but I don’t think it is the one you are interested in. For the compensation, I don’t remember.

However, if I should write this again now, I would not use the same approach, and this could help you implement your attitude model. I would now use `DerivativeStructure` to compute the time derivatives, and implement by myself only the angular part. So we should start by computing `pointingP` and `pointingSat` not as instances of `Vector3D` as it is done today but as `FieldVector3D<DerivativeStructure>` instances. For `pointingP`, there is already a method for that in the class `PVCoordinates`: `toDerivativeStructureVector(order)`. For `pointingSat`, you would need to recompute it by yourself each time `getAttitude` is called instead of having it stored as a fixed vector because it will change depending on the local orbital frame. Once you have these two vectors+derivatives, you can compute the rotation as:

``````   FieldRotation<DerivativeStructure> celToSatRotationDS =
new FieldRotation<>(pointingP, phasingCel, pointingSat, phasingSat);
``````

This computes the rotation, but as the rotation elements are `DerivativeStructure` instances, you get the time derivatives for free. Then, in order to dispatch these hidden derivatives back into a regular `Transform`, you can use

``````  Transform transform = new Transform(date, new AngularCoordinates(celToSatRotationDS));
``````

This works because there is also a dedicated `AngularCoordinates` that `unpack` derivatives.

The advantages of this approach are:

• you don’t have to compute derivatives analytically
• derivatives are not limited to first order, you also get rotation acceleration

The second point (having rotation acceleration) is useful to get more accurate results from `shiftedBy`.

Note also that as an `AttitudeProvider` must have both a regular double version and a field version of getAttitude, the previous approach does work in this case too. You just need to replace `DerivativeStructure` from the regular double version by `FieldDerivativeStructure<T>` in the field version.

I’m afraid this answer my be confusing, so do not hesitate to ask for clarification if needed.

Thank you Luc for your explanation.
I will try to implement a new version of CelestialBodyPointed with a DerivativeStructure as you suggest. If it works, I will implement my own attitude based on this method.

Hello Luc,
For the conversion of CelestialBodyPointed with a DerivativeStructure, I need to convert an AbsoluteDate into FieldAbsoluteDate: how to do this?
Thanks

In the specific case of `FieldAbsoluteDate<DerivativeStructure>` where the `DerivativeStructure` is a derivative with respect to time itself, then the date must be transformed in a way it preserves its derivative with respect to time (i.e. `fDate.durationFrom(otherDate).getPartialDerivative(1)` should return 1.0, not 0.0). Assuming you have already converted some `PVCoordinates` to `FieldVector3D`, you can use the same `DSFactory` that was used during the conversion, which ensures all derivative orders are consistent: the following lines should work:

``````FieldVector3D<DerivativeStructure> pDS = pvt.toDerivativeStructureVector(someOrder);
DerivativeStructure zeroWithDerivatives = pDS.getX().getFactory().variable(0, 0.0);
FiledAbsoluteDate<T> fDate = new FieldAbsoluteDate(date, zeroWithDerivatives);
``````

The trick is to add 0.0 to the date (i.e. we don’t move it), but we most importantly introduce the dependency to the single derivation variable which is time in this very specific case.

In a more general way, you can create a `FieldAbsoluteDate` from an `AbsoluteDate` by just passing the `Field` as a first argument, field that you can retrieve from any variable by calling its `getField()` method. However, this works only if there are no fancy semantic on the field type. In your specific case, there is a dependency, hence the more complex construction.

Hi Luc,
Following your instructions, I’ve recoded CelestialBodyPointed with fields and I’ve validated the rotation and rotation rate with respect to the orginal Orekit version : I’ve got exactly the same results!
Then, I’ve implemented my own attitude class with same principles and validated with finite differences.
Thank you very much for your help

You’re welcome.

Perhaps you could contribute the revamped `CelestialBodyPointed` so it is included in the next version. It would improve the accuracy due to the second derivative and could be a good example for other attitude laws development like the one you needed.

Hi all,

I would like to pick up this thread again since I also have a question concerning the CelestialBodyPointed class. More in detail, I noticed that in the getAttitude method, the transform at the end is built as:

``````Transform transform = new Transform(date, celToSatRotation, celToSatRotation.applyTo(phasedRotAxisCel));
``````

As far as I understood this results in having for this type of Attitude providers always a zero angular acceleration. Is this correct? In case yes, could you explain why?

Thanks

I think the accelerations were not computed because they were to complex to set up and few people seemed to need them. However, since you are asking this, I guess you need the accelerations.

So as I explained earlier in this thread, the best way to compute everything simply and without arbitrary limitations would be to use `DerivativeStructure` everywhere, and to configure them with a derivation order set to 2. This would both simplify the algorithm (because with `DerivativeStructure` you write only the angular part, you don’t even bother about the angular rates) and improve the computation (because you will get the angular acceleration for free).

Do you want to give it a try and contribute the change? You can also open a feature request on the issue tracker.

Hi Luc,

Thank you for your help, I will try to use the DerivativeStructure as you suggest.
If it works, I would be more than happy to contribute and send it to you .