AberrationModifier.java

/* Copyright 2023 Mark Rutten
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * Mark Rutten licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.estimation.measurements.modifiers;

import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import org.hipparchus.Field;
import org.hipparchus.analysis.differentiation.Gradient;
import org.hipparchus.analysis.differentiation.GradientField;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.estimation.measurements.AngularRaDec;
import org.orekit.estimation.measurements.EstimatedMeasurement;
import org.orekit.estimation.measurements.EstimatedMeasurementBase;
import org.orekit.estimation.measurements.EstimationModifier;
import org.orekit.estimation.measurements.GroundStation;
import org.orekit.frames.FieldTransform;
import org.orekit.frames.Frame;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.Constants;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.TimeSpanMap;
import org.orekit.utils.TimeStampedFieldPVCoordinates;


/**
 * Class modifying theoretical angular measurement with (the inverse of) stellar aberration.
 * <p>
 * This class implements equation 3.252-3 from Seidelmann, "Explanatory Supplement to the Astronmical Almanac", 1992.
 *
 * @author Mark Rutten
 */
public class AberrationModifier implements EstimationModifier<AngularRaDec> {

    /** Empty constructor.
     * <p>
     * This constructor is not strictly necessary, but it prevents spurious
     * javadoc warnings with JDK 18 and later.
     * </p>
     * @since 12.0
     */
    public AberrationModifier() {
        // nothing to do
    }

    /**
     * Natural to proper correction for aberration of light.
     *
     * @param naturalRaDec the "natural" direction (in barycentric coordinates)
     * @param station      the observer ground station
     * @param date         the date of the measurement
     * @param frame        the frame of the measurement
     * @return the "proper" direction (station-relative coordinates)
     */
    @DefaultDataContext
    public static double[] naturalToProper(final double[] naturalRaDec, final GroundStation station,
                                           final AbsoluteDate date, final Frame frame) {

        // Check measurement frame is inertial
        if (!frame.isPseudoInertial()) {
            throw new OrekitException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME, frame.getName());
        }

        // Velocity of station relative to barycentre (units of c)
        final PVCoordinates baryPV = CelestialBodyFactory.getSolarSystemBarycenter()
                .getPVCoordinates(date, frame);
        final PVCoordinates stationPV = station.getBaseFrame().getPVCoordinates(date, frame);
        final Vector3D stationBaryVel = stationPV.getVelocity()
                .subtract(baryPV.getVelocity())
                .scalarMultiply(1.0 / Constants.SPEED_OF_LIGHT);

        // Delegate to private method
        return lorentzVelocitySum(naturalRaDec, stationBaryVel);
    }

    /**
     * Proper to natural correction for aberration of light.
     *
     * @param properRaDec the "proper" direction (station-relative coordinates)
     * @param station     the observer ground station
     * @param date        the date of the measurement
     * @param frame       the frame of the measurement
     * @return the "natural" direction (in barycentric coordinates)
     */
    @DefaultDataContext
    public static double[] properToNatural(final double[] properRaDec, final GroundStation station,
                                           final AbsoluteDate date, final Frame frame) {

        // Check measurement frame is inertial
        if (!frame.isPseudoInertial()) {
            throw new OrekitException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME, frame.getName());
        }

        // Velocity of barycentre relative to station (units of c)
        final PVCoordinates baryPV = CelestialBodyFactory.getSolarSystemBarycenter()
                .getPVCoordinates(date, frame);
        final PVCoordinates stationPV = station.getBaseFrame().getPVCoordinates(date, frame);
        final Vector3D baryStationVel = baryPV.getVelocity()
                .subtract(stationPV.getVelocity())
                .scalarMultiply(1.0 / Constants.SPEED_OF_LIGHT);

        // Delegate to private method
        return lorentzVelocitySum(properRaDec, baryStationVel);
    }

    /**
     * Relativistic sum of velocities.
     * This is based on equation 3.252-3 from Seidelmann, "Explanatory Supplement to the Astronmical Almanac", 1992.
     *
     * @param raDec    the direction to transform
     * @param velocity the velocity (units of c)
     * @return the transformed direction
     */
    private static double[] lorentzVelocitySum(final double[] raDec, final Vector3D velocity) {

        // Measurement as unit vector
        final Vector3D direction = new Vector3D(raDec[0], raDec[1]);

        // Coefficients for calculations
        final double inverseBeta = FastMath.sqrt(1.0 - velocity.getNormSq());
        final double velocityScale = 1.0 + direction.dotProduct(velocity) / (1.0 + inverseBeta);

        // From Seidelmann, equation 3.252-3 (unnormalised)
        final Vector3D transformDirection = (direction.scalarMultiply(inverseBeta))
                .add(velocity.scalarMultiply(velocityScale));
        return new double[] {transformDirection.getAlpha(), transformDirection.getDelta()};
    }

    /**
     * Natural to proper correction for aberration of light.
     *
     * @param naturalRaDec      the "natural" direction (in barycentric coordinates)
     * @param stationToInertial the transform from station to inertial coordinates
     * @param frame             the frame of the measurement
     * @return the "proper" direction (station-relative coordinates)
     */
    @DefaultDataContext
    public static Gradient[] fieldNaturalToProper(final Gradient[] naturalRaDec,
                                                  final FieldTransform<Gradient> stationToInertial,
                                                  final Frame frame) {

        // Check measurement frame is inertial
        if (!frame.isPseudoInertial()) {
            throw new OrekitException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME, frame.getName());
        }

        // Set up field
        final Field<Gradient> field = naturalRaDec[0].getField();
        final FieldVector3D<Gradient> zero = FieldVector3D.getZero(field);
        final FieldAbsoluteDate<Gradient> date = stationToInertial.getFieldDate();

        // Barycentre in inertial coordinates
        final PVCoordinates baryPV = CelestialBodyFactory.getSolarSystemBarycenter()
                .getPVCoordinates(date.toAbsoluteDate(), frame);

        // Station in inertial coordinates
        final TimeStampedFieldPVCoordinates<Gradient> stationPV =
                stationToInertial.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(date, zero, zero, zero));

        // Velocity of station relative to barycentre (units of c)
        final FieldVector3D<Gradient> stationBaryVel = stationPV.getVelocity()
                .subtract(baryPV.getVelocity())
                .scalarMultiply(1.0 / Constants.SPEED_OF_LIGHT);

        return fieldLorentzVelocitySum(naturalRaDec, stationBaryVel);
    }


    /**
     * Proper to natural correction for aberration of light.
     *
     * @param properRaDec       the "proper" direction (station-relative coordinates)
     * @param stationToInertial the transform from station to inertial coordinates
     * @param frame             the frame of the measurement
     * @return the "natural" direction (in barycentric coordinates)
     */
    @DefaultDataContext
    public static Gradient[] fieldProperToNatural(final Gradient[] properRaDec,
                                                  final FieldTransform<Gradient> stationToInertial,
                                                  final Frame frame) {

        // Check measurement frame is inertial
        if (!frame.isPseudoInertial()) {
            throw new OrekitException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME, frame.getName());
        }

        // Set up field
        final Field<Gradient> field = properRaDec[0].getField();
        final FieldVector3D<Gradient> zero = FieldVector3D.getZero(field);
        final FieldAbsoluteDate<Gradient> date = stationToInertial.getFieldDate();

        // Barycentre in inertial coordinates
        final PVCoordinates baryPV = CelestialBodyFactory.getSolarSystemBarycenter()
                .getPVCoordinates(date.toAbsoluteDate(), frame);

        // Station in inertial coordinates
        final TimeStampedFieldPVCoordinates<Gradient> stationPV =
                stationToInertial.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(date, zero, zero, zero));

        // Velocity of barycentre relative to station (units of c)
        final FieldVector3D<Gradient> stationBaryVel = stationPV.getVelocity()
                .negate()
                .add(baryPV.getVelocity())
                .scalarMultiply(1.0 / Constants.SPEED_OF_LIGHT);

        return fieldLorentzVelocitySum(properRaDec, stationBaryVel);
    }

    /**
     * Relativistic sum of velocities.
     * This is based on equation 3.252-3 from Seidelmann, "Explanatory Supplement to the Astronmical Almanac", 1992.
     *
     * @param raDec    the direction to transform
     * @param velocity the velocity (units of c)
     * @return the transformed direction
     */
    private static Gradient[] fieldLorentzVelocitySum(final Gradient[] raDec,
                                                      final FieldVector3D<Gradient> velocity) {

        // Measurement as unit vector
        final FieldVector3D<Gradient> direction = new FieldVector3D<>(raDec[0], raDec[1]);

        // Coefficients for calculations
        final Gradient inverseBeta = (velocity.getNormSq().negate().add(1.0)).sqrt();
        final Gradient velocityScale = (direction.dotProduct(velocity)).divide(inverseBeta.add(1.0)).add(1.0);

        // From Seidelmann, equation 3.252-3 (unnormalised)
        final FieldVector3D<Gradient> transformDirection = (direction.scalarMultiply(inverseBeta))
                .add(velocity.scalarMultiply(velocityScale));
        return new Gradient[] {transformDirection.getAlpha(), transformDirection.getDelta()};
    }


    /**
     * {@inheritDoc}
     */
    @Override
    public List<ParameterDriver> getParametersDrivers() {
        return Collections.emptyList();
    }


    @Override
    @DefaultDataContext
    public void modifyWithoutDerivatives(final EstimatedMeasurementBase<AngularRaDec> estimated) {

        // Observation date
        final AbsoluteDate date = estimated.getDate();

        // Observation station
        final GroundStation station = estimated.getObservedMeasurement().getStation();

        // Observation frame
        final Frame frame = estimated.getObservedMeasurement().getReferenceFrame();

        // Convert measurement to natural direction
        final double[] estimatedRaDec = estimated.getEstimatedValue();
        final double[] naturalRaDec = properToNatural(estimatedRaDec, station, date, frame);

        // Normalise RA
        final double[] observed           = estimated.getObservedValue();
        final double   baseRightAscension = naturalRaDec[0];
        final double   twoPiWrap          = MathUtils.normalizeAngle(baseRightAscension, observed[0]) - baseRightAscension;
        final double   rightAscension     = baseRightAscension + twoPiWrap;

        // New estimated values
        estimated.setEstimatedValue(rightAscension, naturalRaDec[1]);

    }

    @Override
    @DefaultDataContext
    public void modify(final EstimatedMeasurement<AngularRaDec> estimated) {

        // Observation date
        final AbsoluteDate date = estimated.getDate();

        // Observation frame
        final Frame frame = estimated.getObservedMeasurement().getReferenceFrame();

        // Extract RA/Dec parameters (no state derivatives)
        int nbParams = 6;
        final Map<String, Integer> indices = new HashMap<>();
        for (ParameterDriver driver : estimated.getObservedMeasurement().getParametersDrivers()) {
            if (driver.isSelected()) {
                for (TimeSpanMap.Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
                    if (!indices.containsKey(span.getData())) {
                        indices.put(span.getData(), nbParams++);
                    }
                }
            }
        }
        final Field<Gradient> field = GradientField.getField(nbParams);

        // Observation location
        final FieldTransform<Gradient> stationToInertial =
                estimated.getObservedMeasurement().getStation().getOffsetToInertial(frame, date, nbParams, indices);

        // Convert measurement to natural direction
        final double[] estimatedRaDec = estimated.getEstimatedValue();
        final Gradient[] estimatedRaDecDS = new Gradient[] {
                field.getZero().add(estimatedRaDec[0]),
                field.getZero().add(estimatedRaDec[1])
        };
        final Gradient[] naturalRaDec = fieldProperToNatural(estimatedRaDecDS, stationToInertial, frame);

        // Normalise RA
        final double[] observed = estimated.getObservedValue();
        final Gradient baseRightAscension = naturalRaDec[0];
        final double twoPiWrap = MathUtils.normalizeAngle(baseRightAscension.getReal(),
                observed[0]) - baseRightAscension.getReal();
        final Gradient rightAscension = baseRightAscension.add(twoPiWrap);

        // New estimated values
        estimated.setEstimatedValue(rightAscension.getValue(), naturalRaDec[1].getValue());

        // Derivatives (only parameter, no state)
        final double[] raDerivatives = rightAscension.getGradient();
        final double[] decDerivatives = naturalRaDec[1].getGradient();

        for (final ParameterDriver driver : estimated.getObservedMeasurement().getParametersDrivers()) {
            for (TimeSpanMap.Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
                final Integer index = indices.get(span.getData());
                if (index != null) {
                    final double[] parameterDerivative = estimated.getParameterDerivatives(driver);
                    parameterDerivative[0] += raDerivatives[index];
                    parameterDerivative[1] += decDerivatives[index];
                    estimated.setParameterDerivatives(driver, span.getStart(), parameterDerivative[0], parameterDerivative[1]);
                }
            }
        }
    }

}