/*
    * transform
    *
    * Parameters:
    *     from:     The geodetic position to be translated.
    *     from_a:   The semi-major axis of the "from" ellipsoid.
    *     from_f:   Flattening of the "from" ellipsoid.
    *     from_esq: Eccentricity-squared of the "from" ellipsoid.
    *     da:       Change in semi-major axis length (meters); "to" minus "from"    *     df:       Change in flattening; "to" minus "from"
    *     dx:       Change in x between "from" and "to" datum.
    *     dy:       Change in y between "from" and "to" datum.
    *     dz:       Change in z between "from" and "to" datum.
    */
    public GeodeticPosition transform (GeodeticPosition from,
                                       double from_a, double from_f,

                                       double from_esq, double da, double df,
                                       double dx, double dy, double dz)
    {
        double slat = Math.sin (from.lat);
        double clat = Math.cos (from.lat);
        double slon = Math.sin (from.lon);
        double clon = Math.cos (from.lon);
        double ssqlat = slat * slat;
        double adb = 1.0 / (1.0 - from_f);  // "a divided by b"
        double dlat, dlon, dh;

        double rn = from_a / Math.sqrt (1.0 - from_esq * ssqlat);
        double rm = from_a * (1. - from_esq) / Math.pow ((1.0 - from_esq * ssqlat), 1.5);

        dlat = (((((-dx * slat * clon - dy * slat * slon) + dz * clat)
                    + (da * ((rn * from_esq * slat * clat) / from_a)))
                + (df * (rm * adb + rn / adb) * slat * clat)))
            / (rm + from.h);

        dlon = (-dx * slon + dy * clon) / ((rn + from.h) * clat);

        dh = (dx * clat * clon) + (dy * clat * slon) + (dz * slat)
             - (da * (from_a / rn)) + ((df * rn * ssqlat) / adb);

        return new GeodeticPosition (from.lon + dlon, from.lat + dlat, from.h + dh);
    }