Hi, I have a question about the Yaw Compensation Attitude Provider, specifically the algorithm used to calculate the attitude.
I have been developing my own algorithm and testing it against the orekit YawCompensation provider, and have been getting an error of about 0.006 degrees.
To get the normal vector (which should be normal to the relative velocity wrt ground), the Spacecraft velocity obtained in an ECEF frame should be sufficient right? This should take into account the movement of the earth as well.
The code is as follows:
#compute earth-fixed references at the specified date
sat_PV = pvcoordinates_provider.getPVCoordinates(date, self.earth_fixed_frame)
sat_pos = sat_PV.getPosition()
sat_vel = sat_PV.getVelocity()
# Compute nadir position in celestial frame
gp_sat = self.earth_shape.transform(
sat_PV.getPosition(), self.earth_fixed_frame, date
)
gp_nadir = GeodeticPoint(
gp_sat.getLatitude(), gp_sat.getLongitude(), float(0.0)
)
nadir_pos = self.earth_shape.transform(gp_nadir)
relativePosition = nadir_pos.subtract(sat_pos)
relativeNormal = Vector3D.crossProduct(relativePosition, sat_vel)
# compute transform from earth-fixed frame to satellite frame
earthToSatRotation = Rotation(
relativePosition, relativeNormal, Vector3D.PLUS_K, Vector3D.PLUS_J
)
However, this does not produce a similar result. I am wondering where my thought process might be wrong. Would appreciate any help. Thanks!
Simple projections are not sufficient.
The behavior is really counter-intuitive, but we have to consider the ground curvature at target point. Since the Earth is an oblate spheroid and not a sphere, this curvature has two main components, one in the East-West direction and one in the North-Sourth direction, and these two curvatures are slightly different (these are called the principal curvatures of a surface). The travel direction of the target point is generally not aligned with either East-West or North-South, so we fly over the curved surface with a slanted path. Here is where this becomes highly counter-intuitive: due to these different curvatures, the projection of the velocity is not in the plane defined by spacecraft velocity, it is skewed. The effect is small, but not negligible, I guess your 0.006° corresponds to this.
The algorithm used in Yaw compensation takes this into account, but the exact method depends on the underlying ground pointing method (NadirPointing uses a 4 points finite difference, LofOffsetPointing uses a sampling interval around current date, BodyCenterPointing uses an analytical method).
I guess we could in fact rely on an analytical method for NadirPointing, because we now have one in OneAxisEllipsoid.projectToGround that does all the magic with curvatures.
So to summarize, it is (much) more complex than a simple velocity projection.
Hey @luc. I have been trying to implement yaw compensation for off-nadir pointing. I have been developing my own algorithm and testing it against the Orekit LofOffsetPointing, and I am getting an error of about 0.008 degrees (a velocity error of 20 m/s). Based on your reply, I first calculated the target intersection on the OneAxisEllipsoid and constructed the two ellipses at that point. I then used Ellipse.projectToEllipse to project the PV on those ellipses and added them (mimicking what OneAxisEllipsoid.projectToGround does) .
The code is as follows:
#satellite PV
sat_PV = pvcoordinates_provider.getPVCoordinates(date, self.earth_fixed_frame)
sat_pos = sat_PV.getPosition()
sat_vel = self.earth_shape.projectToGround(sat_PV, frame).getVelocity()
sat_w= self.att_base.getAttitude(pvcoordinates_provider, date, frame).getSpin()
sat_R= self.att_base.getAttitude(pvcoordinates_provider, date, frame).getRotation()
# target Position on ground
targetPV= self.att_base.getTargetPV(pvcoordinates_provider, date, frame)
target_pos= targetPV.getPosition()
targetR= (target_pos.getX()**2 + target_pos.getY()**2)**0.5
sat_w_fr= sat_R.applyInverseTo(sat_w)
vv= Vector3D.crossProduct(sat_w_fr, target_pos.subtract(sat_pos) )
# Calculating the meridian ellipse
meridian = Vector3D.PLUS_I if targetR==0 else Vector3D(target_pos.getX() / targetR, target_pos.getY() / targetR, 0.)
firstEll= Ellipse(Vector3D.ZERO, meridian, Vector3D.PLUS_K, self.earth_shape.getA(), self.earth_shape.getC(), self.earth_fixed_frame)
#project to first ellipse
firstPV= firstEll.projectToEllipse(sat_PV)
firstV= firstPV.getVelocity()
firstP= firstPV.getPosition()
gr= firstP.getX()*meridian.getX() + firstP.getY()*meridian.getY()
gz= firstP.getZ()
#Calculate the second ellipse
east= Vector3D(-meridian.getY(), meridian.getX(), 0.)
zenith= Vector3D(gr * self.earth_shape.getC() / self.earth_shape.getA(), meridian, gz * self.earth_shape.getA() / self.earth_shape.getC(), Vector3D.PLUS_K).normalize()
north= Vector3D.crossProduct(zenith, east)
# project on second ellipse
secondEll= self.earth_shape.getPlaneSection(firstP, north)
secV= secondEll.projectToEllipse(sat_PV).getVelocity()
# add the two velocity
velocity= firstV.add(secV)
relativePosition = target_pos.subtract(sat_pos)
relativeNormal = Vector3D.crossProduct(relativePosition, velocity)
# compute transform from earth-fixed frame to satellite frame
celToSatRotation = Rotation(
relativePosition, relativeNormal, Vector3D.PLUS_K, Vector3D.PLUS_J
)
I am not able to move forward with this thought process. I would appreciate any help. Thanks.
From what I recall, you really need to compute the principal curvatures to computes projected velocities components along the North-South direction and then along the East-West direction before adding them up and recomputing the projected velocity vector.
One way to check your computation is correct is to validate it in a unit test with respect to finite differences or with respect to a Field version using UnivariateDerivative2 with respect to time (i.e. you just implement the point projection, and you get the derivatives as a side effect).
It is really, really counter-intuitive (in fact even mind-blowing) and requires a lot of efforts to understand and to explain. One way to represents what happens in once mind is to consider what would happen if Earth was not an ellipsoid but a torus, with one gigantic radius for the torus and one tiny radius for the tube. Then, consider a point hovering close to the surface of the torus and crossing the meridian along a slanted path (say 45 degrees). In the inertial frame, the North-South and East-West velocity components would be the same (because of the 45° inclination path), but the projection of the trajectory on the torus would be far more inclined (i.e. the projection skews the velocity direction) because one radius is far different from the other radius. So its very difficult to explain, I’m sorry I cannot do any better.