Examples with GeodeticPoint org.orekit.bodies.GeodeticPoint used on opensource projects

Example 66 with GeodeticPoint

use of org.orekit.bodies.GeodeticPoint in project Orekit by CS-SI.

the class OceanLoading method displacement.

/**
 * {@inheritDoc}
 */
@Override
public Vector3D displacement(final BodiesElements elements, final Frame earthFrame, final Vector3D referencePoint) throws OrekitException {
    // allocate arrays for each species splines
    final UnivariateFunction[] realZSpline = new UnivariateFunction[mainTides.length];
    final UnivariateFunction[] imaginaryZSpline = new UnivariateFunction[mainTides.length];
    final UnivariateFunction[] realWSpline = new UnivariateFunction[mainTides.length];
    final UnivariateFunction[] imaginaryWSpline = new UnivariateFunction[mainTides.length];
    final UnivariateFunction[] realSSpline = new UnivariateFunction[mainTides.length];
    final UnivariateFunction[] imaginarySSpline = new UnivariateFunction[mainTides.length];
    // prepare splines for each species
    for (int i = 0; i < mainTides.length; ++i) {
        // compute current rates
        final double[] rates = new double[mainTides[i].length];
        for (int j = 0; j < rates.length; ++j) {
            rates[j] = mainTides[i][j].tide.getRate(elements);
        }
        // set up splines for the current rates
        realZSpline[i] = spline(rates, mainTides[i], d -> d.realZ);
        imaginaryZSpline[i] = spline(rates, mainTides[i], d -> d.imaginaryZ);
        realWSpline[i] = spline(rates, mainTides[i], d -> d.realW);
        imaginaryWSpline[i] = spline(rates, mainTides[i], d -> d.imaginaryW);
        realSSpline[i] = spline(rates, mainTides[i], d -> d.realS);
        imaginarySSpline[i] = spline(rates, mainTides[i], d -> d.imaginaryS);
    }
    // evaluate all harmonics using interpolated admittances
    double dz = 0;
    double dw = 0;
    double ds = 0;
    for (final Map.Entry<Tide, Double> entry : CARTWRIGHT_EDDEN_AMPLITUDE_MAP.entrySet()) {
        final Tide tide = entry.getKey();
        final double amplitude = entry.getValue();
        final int i = tide.getTauMultiplier();
        final double rate = tide.getRate(elements);
        // apply splines for the current rate of this tide
        final double rZ = realZSpline[i].value(rate);
        final double iZ = imaginaryZSpline[i].value(rate);
        final double rW = realWSpline[i].value(rate);
        final double iW = imaginaryWSpline[i].value(rate);
        final double rS = realSSpline[i].value(rate);
        final double iS = imaginarySSpline[i].value(rate);
        // phase for the current tide, including corrections
        final double correction;
        if (tide.getTauMultiplier() == 0) {
            correction = FastMath.PI;
        } else if (tide.getTauMultiplier() == 1) {
            correction = 0.5 * FastMath.PI;
        } else {
            correction = 0.0;
        }
        final double phase = tide.getPhase(elements) + correction;
        dz += amplitude * FastMath.hypot(rZ, iZ) * FastMath.cos(phase + FastMath.atan2(iZ, rZ));
        dw += amplitude * FastMath.hypot(rW, iW) * FastMath.cos(phase + FastMath.atan2(iW, rW));
        ds += amplitude * FastMath.hypot(rS, iS) * FastMath.cos(phase + FastMath.atan2(iS, rS));
    }
    // convert to proper frame
    final GeodeticPoint gp = earth.transform(referencePoint, earthFrame, null);
    return new Vector3D(dz, gp.getZenith(), dw, gp.getWest(), ds, gp.getSouth());
}

Example 67 with GeodeticPoint

use of org.orekit.bodies.GeodeticPoint in project Orekit by CS-SI.

the class TopocentricFrame method computeLimitVisibilityPoint.

/**
 * Compute the limit visibility point for a satellite in a given direction.
 * <p>
 * This method can be used to compute visibility circles around ground stations
 * for example, using a simple loop on azimuth, with either a fixed elevation
 * or an elevation that depends on azimuth to take ground masks into account.
 * </p>
 * @param radius satellite distance to Earth center
 * @param azimuth pointing azimuth from station
 * @param elevation pointing elevation from station
 * @return limit visibility point for the satellite
 * @throws OrekitException if point cannot be found
 */
public GeodeticPoint computeLimitVisibilityPoint(final double radius, final double azimuth, final double elevation) throws OrekitException {
    try {
        // convergence threshold on point position: 1mm
        final double deltaP = 0.001;
        final UnivariateSolver solver = new BracketingNthOrderBrentSolver(deltaP / Constants.WGS84_EARTH_EQUATORIAL_RADIUS, deltaP, deltaP, 5);
        // find the distance such that a point in the specified direction and at the solved-for
        // distance is exactly at the specified radius
        final double distance = solver.solve(1000, new UnivariateFunction() {

            /**
             * {@inheritDoc}
             */
            public double value(final double x) {
                try {
                    final GeodeticPoint gp = pointAtDistance(azimuth, elevation, x);
                    return parentShape.transform(gp).getNorm() - radius;
                } catch (OrekitException oe) {
                    throw new OrekitExceptionWrapper(oe);
                }
            }
        }, 0, 2 * radius);
        // return the limit point
        return pointAtDistance(azimuth, elevation, distance);
    } catch (MathRuntimeException mrte) {
        throw new OrekitException(mrte);
    } catch (OrekitExceptionWrapper lwe) {
        throw lwe.getException();
    }
}

Example 68 with GeodeticPoint

use of org.orekit.bodies.GeodeticPoint in project Orekit by CS-SI.

the class GeoMagneticField method calculateField.

/**
 * Calculate the magnetic field at the specified geodetic point identified
 * by latitude, longitude and altitude.
 * @param latitude the WGS84 latitude in decimal degrees
 * @param longitude the WGS84 longitude in decimal degrees
 * @param height the height above the WGS84 ellipsoid in kilometers
 * @return the {@link GeoMagneticElements} at the given geodetic point
 */
public GeoMagneticElements calculateField(final double latitude, final double longitude, final double height) {
    final GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(latitude), FastMath.toRadians(longitude), height * 1000d);
    final SphericalCoordinates sph = transformToSpherical(gp);
    final SphericalHarmonicVars vars = new SphericalHarmonicVars(sph);
    final LegendreFunction legendre = new LegendreFunction(FastMath.sin(sph.phi));
    // sum up the magnetic field vector components
    final Vector3D magFieldSph = summation(sph, vars, legendre);
    // rotate the field to geodetic coordinates
    final Vector3D magFieldGeo = rotateMagneticVector(sph, gp, magFieldSph);
    // return the magnetic elements
    return new GeoMagneticElements(magFieldGeo);
}

Example 69 with GeodeticPoint

use of org.orekit.bodies.GeodeticPoint in project Orekit by CS-SI.

the class KlobucharIonoModel method pathDelay.

/**
 * {@inheritDoc}
 */
@Override
public double pathDelay(final AbsoluteDate date, final GeodeticPoint geo, final double elevation, final double azimuth) {
    // degees to semisircles
    final double rad2semi = 1. / FastMath.PI;
    final double semi2rad = FastMath.PI;
    // Earth Centered angle
    final double psi = 0.0137 / (elevation / FastMath.PI + 0.11) - 0.022;
    // Subionospheric latitude: the latitude of the IPP (Ionospheric Pierce Point)
    // in [-0.416, 0.416], semicircle
    final double latIono = FastMath.min(FastMath.max(geo.getLatitude() * rad2semi + psi * FastMath.cos(azimuth), -0.416), 0.416);
    // Subionospheric longitude: the longitude of the IPP
    // in semicircle
    final double lonIono = geo.getLongitude() * rad2semi + (psi * FastMath.sin(azimuth) / FastMath.cos(latIono * semi2rad));
    // Geomagnetic latitude, semicircle
    final double latGeom = latIono + 0.064 * FastMath.cos((lonIono - 1.617) * semi2rad);
    // day of week and tow (sec)
    // Note: Sunday=0, Monday=1, Tuesday=2, Wednesday=3, Thursday=4, Friday=5, Saturday=6
    final DateTimeComponents dtc = date.getComponents(TimeScalesFactory.getGPS());
    final int dofweek = dtc.getDate().getDayOfWeek();
    final double secday = dtc.getTime().getSecondsInLocalDay();
    final double tow = dofweek * 86400. + secday;
    final double t = 43200. * lonIono + tow;
    // Seconds of day
    final double tsec = t - FastMath.floor(t / 86400.) * 86400;
    // Slant factor, semicircle
    final double slantFactor = 1.0 + 16.0 * FastMath.pow(0.53 - elevation / FastMath.PI, 3);
    // Period of model, seconds
    final double period = FastMath.max(72000., beta[0] + (beta[1] + (beta[2] + beta[3] * latGeom) * latGeom) * latGeom);
    // Phase of the model, radians
    // (Max at 14.00 = 50400 sec local time)
    final double x = 2.0 * FastMath.PI * (tsec - 50400.0) / period;
    // Amplitude of the model, seconds
    final double amplitude = FastMath.max(0, alpha[0] + (alpha[1] + (alpha[2] + alpha[3] * latGeom) * latGeom) * latGeom);
    // Ionospheric correction (L1)
    double ionoTimeDelayL1 = slantFactor * (5. * 1e-9);
    if (FastMath.abs(x) < 1.570) {
        ionoTimeDelayL1 += slantFactor * (amplitude * (1.0 - FastMath.pow(x, 2) / 2.0 + FastMath.pow(x, 4) / 24.0));
    }
    // Ionospheric delay for the L1 frequency, in meters, with slant correction.
    return ratio * Constants.SPEED_OF_LIGHT * ionoTimeDelayL1;
}

Example 70 with GeodeticPoint

use of org.orekit.bodies.GeodeticPoint in project Orekit by CS-SI.

the class FieldOfView method getFootprint.

/**
 * Get the footprint of the field Of View on ground.
 * <p>
 * This method assumes the Field Of View is centered on some carrier,
 * which will typically be a spacecraft or a ground station antenna.
 * The points in the footprint boundary loops are all at altitude zero
 * with respect to the ellipsoid, they correspond either to projection
 * on ground of the edges of the Field Of View, or to points on the body
 * limb if the Field Of View goes past horizon. The points on the limb
 * see the carrier origin at zero elevation. If the Field Of View is so
 * large it contains entirely the body, all points will correspond to
 * points at limb. If the Field Of View looks away from body, the
 * boundary loops will be an empty list. The points within footprint
 * the loops are sorted in trigonometric order as seen from the carrier.
 * This implies that someone traveling on ground from one point to the
 * next one will have the points visible from the carrier on his left
 * hand side, and the points not visible from the carrier on his right
 * hand side.
 * </p>
 * <p>
 * The truncation of Field Of View at limb can induce strange results
 * for complex Fields Of View. If for example a Field Of View is a
 * ring with a hole and part of the ring goes past horizon, then instead
 * of having a single loop with a C-shaped boundary, the method will
 * still return two loops truncated at the limb, one clockwise and one
 * counterclockwise, hence "closing" the C-shape twice. This behavior
 * is considered acceptable.
 * </p>
 * <p>
 * If the carrier is a spacecraft, then the {@code fovToBody} transform
 * can be computed from a {@link org.orekit.propagation.SpacecraftState}
 * as follows:
 * </p>
 * <pre>
 * Transform inertToBody = state.getFrame().getTransformTo(body.getBodyFrame(), state.getDate());
 * Transform fovToBody   = new Transform(state.getDate(),
 *                                       state.toTransform().getInverse(),
 *                                       inertToBody);
 * </pre>
 * <p>
 * If the carrier is a ground station, located using a topocentric frame
 * and managing its pointing direction using a transform between the
 * dish frame and the topocentric frame, then the {@code fovToBody} transform
 * can be computed as follows:
 * </p>
 * <pre>
 * Transform topoToBody = topocentricFrame.getTransformTo(body.getBodyFrame(), date);
 * Transform topoToDish = ...
 * Transform fovToBody  = new Transform(date,
 *                                      topoToDish.getInverse(),
 *                                      topoToBody);
 * </pre>
 * <p>
 * Only the raw zone is used, the angular margin is ignored here.
 * </p>
 * @param fovToBody transform between the frame in which the Field Of View
 * is defined and body frame.
 * @param body body surface the Field Of View will be projected on
 * @param angularStep step used for boundary loops sampling (radians)
 * @return list footprint boundary loops (there may be several independent
 * loops if the Field Of View shape is complex)
 * @throws OrekitException if some frame conversion fails or if carrier is
 * below body surface
 */
List<List<GeodeticPoint>> getFootprint(final Transform fovToBody, final OneAxisEllipsoid body, final double angularStep) throws OrekitException {
    final Frame bodyFrame = body.getBodyFrame();
    final Vector3D position = fovToBody.transformPosition(Vector3D.ZERO);
    final double r = position.getNorm();
    if (body.isInside(position)) {
        throw new OrekitException(OrekitMessages.POINT_INSIDE_ELLIPSOID);
    }
    final List<List<GeodeticPoint>> footprint = new ArrayList<List<GeodeticPoint>>();
    final List<Vertex> boundary = zone.getBoundaryLoops();
    for (final Vertex loopStart : boundary) {
        int count = 0;
        final List<GeodeticPoint> loop = new ArrayList<GeodeticPoint>();
        boolean intersectionsFound = false;
        for (Edge edge = loopStart.getOutgoing(); count == 0 || edge.getStart() != loopStart; edge = edge.getEnd().getOutgoing()) {
            ++count;
            final int n = (int) FastMath.ceil(edge.getLength() / angularStep);
            final double delta = edge.getLength() / n;
            for (int i = 0; i < n; ++i) {
                final Vector3D awaySC = new Vector3D(r, edge.getPointAt(i * delta));
                final Vector3D awayBody = fovToBody.transformPosition(awaySC);
                final Line lineOfSight = new Line(position, awayBody, 1.0e-3);
                GeodeticPoint gp = body.getIntersectionPoint(lineOfSight, position, bodyFrame, null);
                if (gp != null && Vector3D.dotProduct(awayBody.subtract(position), body.transform(gp).subtract(position)) < 0) {
                    // the intersection is in fact on the half-line pointing
                    // towards the back side, it is a spurious intersection
                    gp = null;
                }
                if (gp != null) {
                    // the line of sight does intersect the body
                    intersectionsFound = true;
                } else {
                    // the line of sight does not intersect body
                    // we use a point on the limb
                    gp = body.transform(body.pointOnLimb(position, awayBody), bodyFrame, null);
                }
                // add the point in front of the list
                // (to ensure the loop will be in trigonometric orientation)
                loop.add(0, gp);
            }
        }
        if (intersectionsFound) {
            // at least some of the points did intersect the body,
            // this loop contributes to the footprint
            footprint.add(loop);
        }
    }
    if (footprint.isEmpty()) {
        // none of the Field Of View loops cross the body
        // either the body is outside of Field Of View, or it is fully contained
        // we check the center
        final Vector3D bodyCenter = fovToBody.getInverse().transformPosition(Vector3D.ZERO);
        if (zone.checkPoint(new S2Point(bodyCenter)) != Region.Location.OUTSIDE) {
            // the body is fully contained in the Field Of View
            // we use the full limb as the footprint
            final Vector3D x = bodyCenter.orthogonal();
            final Vector3D y = Vector3D.crossProduct(bodyCenter, x).normalize();
            final double sinEta = body.getEquatorialRadius() / r;
            final double sinEta2 = sinEta * sinEta;
            final double cosAlpha = (FastMath.cos(angularStep) + sinEta2 - 1) / sinEta2;
            final int n = (int) FastMath.ceil(MathUtils.TWO_PI / FastMath.acos(cosAlpha));
            final double delta = MathUtils.TWO_PI / n;
            final List<GeodeticPoint> loop = new ArrayList<GeodeticPoint>(n);
            for (int i = 0; i < n; ++i) {
                final Vector3D outside = new Vector3D(r * FastMath.cos(i * delta), x, r * FastMath.sin(i * delta), y);
                loop.add(body.transform(body.pointOnLimb(position, outside), bodyFrame, null));
            }
            footprint.add(loop);
        }
    }
    return footprint;
}