TopocentricFrame.java

  1. /* Copyright 2002-2013 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.frames;

  18. import org.apache.commons.math3.analysis.UnivariateFunction;
  19. import org.apache.commons.math3.analysis.solvers.BracketingNthOrderBrentSolver;
  20. import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
  21. import org.apache.commons.math3.exception.TooManyEvaluationsException;
  22. import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
  23. import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
  24. import org.apache.commons.math3.util.FastMath;
  25. import org.apache.commons.math3.util.MathUtils;
  26. import org.orekit.bodies.BodyShape;
  27. import org.orekit.bodies.GeodeticPoint;
  28. import org.orekit.errors.OrekitException;
  29. import org.orekit.errors.OrekitExceptionWrapper;
  30. import org.orekit.time.AbsoluteDate;
  31. import org.orekit.utils.Constants;
  32. import org.orekit.utils.PVCoordinates;
  33. import org.orekit.utils.PVCoordinatesProvider;


  34. /** Topocentric frame.
  35.  * <p>Frame associated to a position near the surface of a body shape.</p>
  36.  * <p>
  37.  * The origin of the frame is at the defining {@link GeodeticPoint geodetic point}
  38.  * location, and the right-handed canonical trihedra is:
  39.  * </p>
  40.  * <ul>
  41.  *   <li>X axis in the local horizontal plane (normal to zenith direction) and
  42.  *   following the local parallel towards East</li>
  43.  *   <li>Y axis in the horizontal plane (normal to zenith direction) and
  44.  *   following the local meridian towards North</li>
  45.  *   <li>Z axis towards Zenith direction</li>
  46.  * </ul>
  47.  * @author V&eacute;ronique Pommier-Maurussane
  48.  */
  49. public class TopocentricFrame extends Frame implements PVCoordinatesProvider {

  50.     /** Serializable UID. */
  51.     private static final long serialVersionUID = -5997915708080966466L;

  52.     /** Body shape on which the local point is defined. */
  53.     private final BodyShape parentShape;

  54.     /** Point where the topocentric frame is defined. */
  55.     private final GeodeticPoint point;

  56.     /** Simple constructor.
  57.      * @param parentShape body shape on which the local point is defined
  58.      * @param point local surface point where topocentric frame is defined
  59.      * @param name the string representation
  60.      */
  61.     public TopocentricFrame(final BodyShape parentShape, final GeodeticPoint point,
  62.                             final String name) {

  63.         super(parentShape.getBodyFrame(),
  64.               new Transform(AbsoluteDate.J2000_EPOCH,
  65.                             new Transform(AbsoluteDate.J2000_EPOCH,
  66.                                           parentShape.transform(point).negate()),
  67.                             new Transform(AbsoluteDate.J2000_EPOCH,
  68.                                           new Rotation(point.getEast(), point.getZenith(),
  69.                                                        Vector3D.PLUS_I, Vector3D.PLUS_K),
  70.                                           Vector3D.ZERO)),
  71.               name, false);
  72.         this.parentShape = parentShape;
  73.         this.point = point;
  74.     }

  75.     /** Get the body shape on which the local point is defined.
  76.      * @return body shape on which the local point is defined
  77.      */
  78.     public BodyShape getParentShape() {
  79.         return parentShape;
  80.     }

  81.     /** Get the surface point defining the origin of the frame.
  82.      * @return surface point defining the origin of the frame
  83.      */
  84.     public GeodeticPoint getPoint() {
  85.         return point;
  86.     }

  87.     /** Get the zenith direction of topocentric frame, expressed in parent shape frame.
  88.      * <p>The zenith direction is defined as the normal to local horizontal plane.</p>
  89.      * @return unit vector in the zenith direction
  90.      * @see #getNadir()
  91.      */
  92.     public Vector3D getZenith() {
  93.         return point.getZenith();
  94.     }

  95.     /** Get the nadir direction of topocentric frame, expressed in parent shape frame.
  96.      * <p>The nadir direction is the opposite of zenith direction.</p>
  97.      * @return unit vector in the nadir direction
  98.      * @see #getZenith()
  99.      */
  100.     public Vector3D getNadir() {
  101.         return point.getNadir();
  102.     }

  103.    /** Get the north direction of topocentric frame, expressed in parent shape frame.
  104.      * <p>The north direction is defined in the horizontal plane
  105.      * (normal to zenith direction) and following the local meridian.</p>
  106.      * @return unit vector in the north direction
  107.      * @see #getSouth()
  108.      */
  109.     public Vector3D getNorth() {
  110.         return point.getNorth();
  111.     }

  112.     /** Get the south direction of topocentric frame, expressed in parent shape frame.
  113.      * <p>The south direction is the opposite of north direction.</p>
  114.      * @return unit vector in the south direction
  115.      * @see #getNorth()
  116.      */
  117.     public Vector3D getSouth() {
  118.         return point.getSouth();
  119.     }

  120.     /** Get the east direction of topocentric frame, expressed in parent shape frame.
  121.      * <p>The east direction is defined in the horizontal plane
  122.      * in order to complete direct triangle (east, north, zenith).</p>
  123.      * @return unit vector in the east direction
  124.      * @see #getWest()
  125.      */
  126.     public Vector3D getEast() {
  127.         return point.getEast();
  128.     }

  129.     /** Get the west direction of topocentric frame, expressed in parent shape frame.
  130.      * <p>The west direction is the opposite of east direction.</p>
  131.      * @return unit vector in the west direction
  132.      * @see #getEast()
  133.      */
  134.     public Vector3D getWest() {
  135.         return point.getWest();
  136.     }

  137.     /** Get the elevation of a point with regards to the local point.
  138.      * <p>The elevation is the angle between the local horizontal and
  139.      * the direction from local point to given point.</p>
  140.      * @param extPoint point for which elevation shall be computed
  141.      * @param frame frame in which the point is defined
  142.      * @param date computation date
  143.      * @return elevation of the point
  144.      * @exception OrekitException if frames transformations cannot be computed
  145.      */
  146.     public double getElevation(final Vector3D extPoint, final Frame frame,
  147.                                final AbsoluteDate date)
  148.         throws OrekitException {

  149.         // Transform given point from given frame to topocentric frame
  150.         final Transform t = frame.getTransformTo(this, date);
  151.         final Vector3D extPointTopo = t.transformPosition(extPoint);

  152.         // Elevation angle is PI/2 - angle between zenith and given point direction
  153.         return extPointTopo.getDelta();
  154.     }

  155.     /** Get the azimuth of a point with regards to the topocentric frame center point.
  156.      * <p>The azimuth is the angle between the North direction at local point and
  157.      * the projection in local horizontal plane of the direction from local point
  158.      * to given point. Azimuth angles are counted clockwise, i.e positive towards the East.</p>
  159.      * @param extPoint point for which elevation shall be computed
  160.      * @param frame frame in which the point is defined
  161.      * @param date computation date
  162.      * @return azimuth of the point
  163.      * @exception OrekitException if frames transformations cannot be computed
  164.      */
  165.     public double getAzimuth(final Vector3D extPoint, final Frame frame,
  166.                              final AbsoluteDate date)
  167.         throws OrekitException {

  168.         // Transform given point from given frame to topocentric frame
  169.         final Transform t = getTransformTo(frame, date).getInverse();
  170.         final Vector3D extPointTopo = t.transformPosition(extPoint);

  171.         // Compute azimuth
  172.         double azimuth = FastMath.atan2(extPointTopo.getX(), extPointTopo.getY());
  173.         if (azimuth < 0.) {
  174.             azimuth += MathUtils.TWO_PI;
  175.         }
  176.         return azimuth;

  177.     }

  178.     /** Get the range of a point with regards to the topocentric frame center point.
  179.      * @param extPoint point for which range shall be computed
  180.      * @param frame frame in which the point is defined
  181.      * @param date computation date
  182.      * @return range (distance) of the point
  183.      * @exception OrekitException if frames transformations cannot be computed
  184.      */
  185.     public double getRange(final Vector3D extPoint, final Frame frame,
  186.                            final AbsoluteDate date)
  187.         throws OrekitException {

  188.         // Transform given point from given frame to topocentric frame
  189.         final Transform t = frame.getTransformTo(this, date);
  190.         final Vector3D extPointTopo = t.transformPosition(extPoint);

  191.         // Compute range
  192.         return extPointTopo.getNorm();

  193.     }

  194.     /** Get the range rate of a point with regards to the topocentric frame center point.
  195.      * @param extPV point/velocity for which range rate shall be computed
  196.      * @param frame frame in which the point is defined
  197.      * @param date computation date
  198.      * @return range rate of the point (positive if point departs from frame)
  199.      * @exception OrekitException if frames transformations cannot be computed
  200.      */
  201.     public double getRangeRate(final PVCoordinates extPV, final Frame frame,
  202.                                final AbsoluteDate date)
  203.         throws OrekitException {

  204.         // Transform given point from given frame to topocentric frame
  205.         final Transform t = frame.getTransformTo(this, date);
  206.         final PVCoordinates extPVTopo = t.transformPVCoordinates(extPV);

  207.         // Compute range rate (doppler) : relative rate along the line of sight
  208.         return Vector3D.dotProduct(extPVTopo.getPosition(), extPVTopo.getVelocity()) /
  209.                extPVTopo.getPosition().getNorm();

  210.     }

  211.     /**
  212.      * Compute the limit visibility point for a satellite in a given direction.
  213.      * <p>
  214.      * This method can be used to compute visibility circles around ground stations
  215.      * for example, using a simple loop on azimuth, with either a fixed elevation
  216.      * or an elevation that depends on azimuth to take ground masks into account.
  217.      * </p>
  218.      * @param radius satellite distance to Earth center
  219.      * @param azimuth pointing azimuth from station
  220.      * @param elevation pointing elevation from station
  221.      * @return limit visibility point for the satellite
  222.      * @throws OrekitException if point cannot be found
  223.      */
  224.     public GeodeticPoint computeLimitVisibilityPoint(final double radius,
  225.                                                      final double azimuth, final double elevation)
  226.         throws OrekitException {
  227.         try {
  228.             // convergence threshold on point position: 1mm
  229.             final double deltaP = 0.001;
  230.             final UnivariateSolver solver =
  231.                     new BracketingNthOrderBrentSolver(deltaP / Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
  232.                                                       deltaP, deltaP, 5);

  233.             // find the distance such that a point in the specified direction and at the solved-for
  234.             // distance is exactly at the specified radius
  235.             final double distance = solver.solve(1000, new UnivariateFunction() {
  236.                 /** {@inheritDoc} */
  237.                 public double value(final double x) {
  238.                     try {
  239.                         final GeodeticPoint gp = pointAtDistance(azimuth, elevation, x);
  240.                         return parentShape.transform(gp).getNorm() - radius;
  241.                     } catch (OrekitException oe) {
  242.                         throw new OrekitExceptionWrapper(oe);
  243.                     }
  244.                 }
  245.             }, 0, 2 * radius);

  246.             // return the limit point
  247.             return pointAtDistance(azimuth, elevation, distance);

  248.         } catch (TooManyEvaluationsException tmee) {
  249.             throw new OrekitException(tmee);
  250.         } catch (OrekitExceptionWrapper lwe) {
  251.             throw lwe.getException();
  252.         }
  253.     }

  254.     /** Compute the point observed from the station at some specified distance.
  255.      * @param azimuth pointing azimuth from station
  256.      * @param elevation pointing elevation from station
  257.      * @param distance distance to station
  258.      * @return observed point
  259.      * @exception OrekitException if point cannot be computed
  260.      */
  261.     public GeodeticPoint pointAtDistance(final double azimuth, final double elevation,
  262.                                          final double distance)
  263.         throws OrekitException {
  264.         final double cosAz = FastMath.cos(azimuth);
  265.         final double sinAz = FastMath.sin(azimuth);
  266.         final double cosEl = FastMath.cos(elevation);
  267.         final double sinEl = FastMath.sin(elevation);
  268.         final Vector3D  observed = new Vector3D(distance * cosEl * sinAz,
  269.                                                 distance * cosEl * cosAz,
  270.                                                 distance * sinEl);
  271.         return parentShape.transform(observed, this, AbsoluteDate.J2000_EPOCH);
  272.     }

  273.     /** Get the {@link PVCoordinates} of the topocentric frame origin in the selected frame.
  274.      * @param date current date
  275.      * @param frame the frame where to define the position
  276.      * @return position/velocity of the topocentric frame origin (m and m/s)
  277.      * @exception OrekitException if position cannot be computed in given frame
  278.      */
  279.     public PVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame)
  280.         throws OrekitException {
  281.         return getTransformTo(frame, date).transformPVCoordinates(PVCoordinates.ZERO);
  282.     }

  283. }