Examples with DoubleDouble org.apache.sis.internal.util.DoubleDouble used on opensource projects

Example 16 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class MathFunctions method magnitude.

/**
 * Returns the magnitude of the given vector. This is defined by:
 *
 * {@preformat math
 *     sqrt(vector[0]² + vector[1]² + … + vector[length-1]²)
 * }
 *
 * <div class="section">Implementation note</div>
 * In the special case where only one element is different than zero, this method
 * returns directly the {@linkplain Math#abs(double) absolute value} of that element
 * without computing {@code sqrt(v²)}, in order to avoid rounding error. This special case
 * has been implemented because this method is often invoked for computing the length of
 * offset vectors,
 * typically aligned with the axes of a {@linkplain org.opengis.referencing.cs.CartesianCS
 * Cartesian coordinate system}.
 *
 * @param  vector  the vector for which to compute the magnitude.
 * @return the magnitude of the given vector.
 *
 * @see Math#hypot(double, double)
 */
public static double magnitude(final double... vector) {
    int i = vector.length;
    // If every elements in the array are zero, returns zero.
    double v1;
    do if (i == 0)
        return 0; while ((v1 = vector[--i]) == 0);
    // We have found a non-zero element. If it is the only one, returns it directly.
    double v2;
    do if (i == 0)
        return Math.abs(v1); while ((v2 = vector[--i]) == 0);
    // If there is exactly 2 elements, use Math.hypot which is more robust than our algorithm.
    double v3;
    do if (i == 0)
        return Math.hypot(v1, v2); while ((v3 = vector[--i]) == 0);
    // Usual magnitude computation, but using double-double arithmetic.
    final DoubleDouble sum = new DoubleDouble();
    final DoubleDouble dot = new DoubleDouble();
    sum.setToProduct(v1, v1);
    dot.setToProduct(v2, v2);
    sum.add(dot);
    dot.setToProduct(v3, v3);
    sum.add(dot);
    while (i != 0) {
        v1 = vector[--i];
        dot.setToProduct(v1, v1);
        sum.add(dot);
    }
    sum.sqrt();
    return sum.value;
}

Example 17 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class BursaWolfParameters method getPositionVectorTransformation.

/**
 * Returns the position vector transformation (geocentric domain) as an affine transform.
 * For transformations that do not depend on time, the formula is as below where {@code R}
 * is a conversion factor from arc-seconds to radians:
 *
 * <blockquote><pre> R = toRadians(1″)
 * S = 1 + {@linkplain #dS}/1000000
 * ┌    ┐    ┌                               ┐  ┌   ┐
 * │ X' │    │      S   -{@linkplain #rZ}*RS   +{@linkplain #rY}*RS   {@linkplain #tX} │  │ X │
 * │ Y' │  = │ +{@linkplain #rZ}*RS        S   -{@linkplain #rX}*RS   {@linkplain #tY} │  │ Y │
 * │ Z' │    │ -{@linkplain #rY}*RS   +{@linkplain #rX}*RS        S   {@linkplain #tZ} │  │ Z │
 * │ 1  │    │      0        0        0    1 │  │ 1 │
 * └    ┘    └                               ┘  └   ┘</pre></blockquote>
 *
 * This affine transform can be applied on <strong>geocentric</strong> coordinates.
 * This is identified as operation method 1033 in the EPSG database.
 * Those geocentric coordinates are typically converted from geographic coordinates
 * in the region or timeframe given by {@link #getDomainOfValidity()}.
 *
 * <p>If the source datum and the {@linkplain #getTargetDatum() target datum} do not use the same
 * {@linkplain DefaultGeodeticDatum#getPrimeMeridian() prime meridian}, then it is caller's responsibility
 * to apply longitude rotation before to use the matrix returned by this method.</p>
 *
 * <div class="section">Time-dependent transformation</div>
 * Some transformations use parameters that vary with time (e.g. operation method EPSG:1053).
 * Users can optionally specify a date for which the transformation is desired.
 * For transformations that do not depends on time, this date is ignored and can be null.
 * For time-dependent transformations, {@code null} values default to the transformation's
 * {@linkplain TimeDependentBWP#getTimeReference() reference time}.
 *
 * <div class="section">Inverse transformation</div>
 * The inverse transformation can be approximated by reversing the sign of the 7 parameters before to use them
 * in the above matrix. This is often considered sufficient since <cite>position vector transformations</cite>
 * are themselves approximations. However Apache SIS will rather use
 * {@link org.apache.sis.referencing.operation.matrix.MatrixSIS#inverse()} in order to increase the chances
 * that concatenation of transformations <var>A</var> → <var>B</var> followed by <var>B</var> → <var>A</var>
 * gives back the identity transform.
 *
 * @param  time  date for which the transformation is desired, or {@code null} for the transformation's reference time.
 * @return an affine transform in geocentric space created from this Bursa-Wolf parameters and the given time.
 *
 * @see DefaultGeodeticDatum#getPositionVectorTransformation(GeodeticDatum, Extent)
 */
public Matrix getPositionVectorTransformation(final Date time) {
    final DoubleDouble period = period(time);
    if (period == null && isTranslation()) {
        final Matrix4 matrix = new Matrix4();
        matrix.m03 = tX;
        matrix.m13 = tY;
        matrix.m23 = tZ;
        return matrix;
    }
    /*
         * Above was an optimization for the common case where the Bursa-Wolf parameters contain only
         * translation terms. If we have rotation or scale terms, then use double-double arithmetic.
         */
    final DoubleDouble RS = DoubleDouble.createSecondsToRadians();
    final DoubleDouble S = param(6, period);
    S.divide(PPM, 0);
    // S = 1 + dS / PPM;
    S.add(1, 0);
    // RS = toRadians(1″) * S;
    RS.multiply(S);
    final DoubleDouble X = param(3, period);
    X.multiply(RS);
    final DoubleDouble Y = param(4, period);
    Y.multiply(RS);
    final DoubleDouble Z = param(5, period);
    Z.multiply(RS);
    final DoubleDouble mX = new DoubleDouble(X);
    mX.negate();
    final DoubleDouble mY = new DoubleDouble(Y);
    mY.negate();
    final DoubleDouble mZ = new DoubleDouble(Z);
    mZ.negate();
    // Fetch Integer instance only once.
    final Integer O = 0;
    return Matrices.create(4, 4, new Number[] { S, mZ, Y, param(0, period), Z, S, mX, param(1, period), mY, X, S, param(2, period), O, O, O, 1 });
}

Example 18 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class DatumShiftGrid method getCellMean.

/**
 * Returns an average translation value for the given dimension.
 * Those average values shall provide a good "first guess" before to interpolate the actual translation value
 * at the (<var>x</var>,<var>y</var>) coordinate. This "first guess" is needed for inverse transform.
 *
 * <div class="section">Default implementation</div>
 * The default implementation computes the average of all values returned by
 * {@link #getCellValue getCellValue(dim, …)}, but subclasses may override with more specific values.
 *
 * <div class="note"><b>Example:</b>
 * In the <cite>"France geocentric interpolation"</cite> (ESPG:9655) operation method, those "average" values
 * are fixed by definition to -168, -60 and +320 metres for dimensions 0, 1 and 2 respectively
 * (geocentric <var>X</var>, <var>Y</var> and <var>Z</var>).</div>
 *
 * @param  dim  the dimension for which to get an average translation value,
 *              from 0 inclusive to {@link #getTranslationDimensions()} exclusive.
 * @return a translation value close to the average for the given dimension.
 */
public double getCellMean(final int dim) {
    final DoubleDouble sum = new DoubleDouble();
    final int nx = gridSize[0];
    final int ny = gridSize[1];
    for (int gridY = 0; gridY < ny; gridY++) {
        for (int gridX = 0; gridX < nx; gridX++) {
            sum.add(getCellValue(dim, gridX, gridY));
        }
    }
    return sum.value / (nx * ny);
}

Example 19 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class DefaultEllipsoid method getEccentricity.

/**
 * The ratio of the distance between the center and a focus of the ellipse to the length of its semi-major axis.
 * The eccentricity can alternately be computed from the equation: ℯ = √(2f - f²) where <var>f</var> is the
 * flattening factor (not inverse).
 *
 * @return ℯ, the eccentricity of this ellipsoid.
 */
public double getEccentricity() {
    final DoubleDouble e = eccentricitySquared();
    e.sqrt();
    return e.value;
}

Example 20 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class DefaultEllipsoid method semiMajorAxisDifference.

/**
 * Returns the difference between the semi-major axis length of two ellipsoids.
 * If the two ellipsoid does not use the same unit of measurement, than the axis
 * length of the other ellipsoid is converted into the units of this ellipsoid axis.
 *
 * <div class="note"><b>Example:</b>
 * {@code WGS84.semiMajorAxisDifference(ED50)} returns 251 metres. This information is a parameter of
 * {@linkplain org.apache.sis.referencing.operation.transform.MolodenskyTransform Molodensky transformations}.</div>
 *
 * @param  other  the other ellipsoid from which to get semi-major axis length difference.
 * @return (<var>other</var> ellipsoid semi-major axis) - (<var>this</var> ellipsoid semi-major axis).
 *
 * @since 0.7
 */
public double semiMajorAxisDifference(final Ellipsoid other) {
    double semiMajor = other.getSemiMajorAxis();
    // Often a no-op.
    semiMajor = other.getAxisUnit().getConverterTo(getAxisUnit()).convert(semiMajor);
    // Presumed accurate in base 10 if no unit conversion.
    final DoubleDouble a = new DoubleDouble(semiMajor);
    // Presumed accurate in base 10 (not 2) by definition.
    a.subtract(getSemiMajorAxis());
    return a.value;
}