WGS84 to ETRS-TM35FIN coordinate system converter in Java

WGS84_to_ETRS89_TM35FIN.java

public class WGS84_to_ETRS89_TM35FIN {
    private static final double F = 1d / 298.257222101d; // Ellipsoidin litistyssuhde
    private static final double N = F / (2d - F);
    private static final double A1 = (6378137d / (1d + N)) * (1d + (Math.pow(N, 2) / 4d) + (Math.pow(N, 4) / 64d)); // Isoakselin puolikas
    private static final double E_POW_SQRT = Math.sqrt((2d * F) - Math.pow(F, 2));

    private static final double H1 = (1d / 2d) * N - (2d / 3d) * Math.pow(N, 2) + (5d / 16d) * Math.pow(N, 3) + (41d / 180d) * Math.pow(N, 4);
    private static final double H2 = (13d / 48d) * Math.pow(N, 2) - (3d / 5d) * Math.pow(N, 3) + (557d / 1440d) * Math.pow(N, 4);
    private static final double H3 = (61d / 240d) * Math.pow(N, 3) - (103d / 140d) * Math.pow(N, 4);
    private static final double H4 = (49561d / 161280d) * Math.pow(N, 4);

    private static final double LAMBDA_ZERO = 0.471238898d; // Keskimeridiaani (rad), 27 astetta
    private static final double K_ZERO = 0.9996d; // Mittakaavakerroin
    private static final double E_ZERO = 500000d; // Itäkoordinaatti

    public static void main(final String[] args) {
        // Tests
        System.out.println("%s => %s".formatted(
            "61.690427,27.267360",
            convert(61.690427d, 27.267360d))
        );

        System.out.println("%s => %s".formatted(
            "61.680366,27.258336",
            convert(61.680366d, 27.258336d))
        );

        System.out.println("%s => %s".formatted(
            "65.017913,25.475675",
            convert(65.017913d, 25.475675d))
        );
    }

    /**
     * Converts WGS84 to ETRS89_TM35FIN
     * <p>Based on http://www.loukko.net/koord_proj/ and https://www.suomidigi.fi/ohjeet-ja-tuki/jhs-suositukset/jhs-197-euref-fin-koordinaattijarjestelmat-niihin-liittyvat-muunnokset-ja-karttalehtijako</p>
     * @author JHaiko 2023
     */
    public static String convert(final double latitude, final double longitude) {
        final double fii = Math.toRadians(latitude);
        final double beta = Math.atan(Math.sinh(asinh(Math.tan(fii)) - E_POW_SQRT * atanh(E_POW_SQRT * Math.sin(fii))));
        final double etaDot = atanh(Math.cos(beta) * Math.sin(Math.toRadians(longitude) - LAMBDA_ZERO));
        final double zetaDot = Math.asin(Math.sin(beta) / (1d / Math.cosh(etaDot)));

        final double zeta =
            zetaDot
            + (H1 * Math.sin(2d * zetaDot) * Math.cosh(2d * etaDot))
            + (H2 * Math.sin(4d * zetaDot) * Math.cosh(4d * etaDot))
            + (H3 * Math.sin(6d * zetaDot) * Math.cosh(6d * etaDot))
            + (H4 * Math.sin(8d * zetaDot) * Math.cosh(8d * etaDot));

        final double eta =
            etaDot
            + (H1 * Math.cos(2d * zetaDot) * Math.sinh(2d * etaDot))
            + (H2 * Math.cos(4d * zetaDot) * Math.sinh(4d * etaDot))
            + (H3 * Math.cos(6d * zetaDot) * Math.sinh(6d * etaDot))
            + (H4 * Math.cos(8d * zetaDot) * Math.sinh(8d * etaDot));

        return "%s,%s".formatted(
            A1 * zeta * K_ZERO,
            A1 * eta * K_ZERO + E_ZERO
        );
    }

    // FIXME: Use Apache Commons Math
    private static double asinh(final double _double) {
        final double sign;
        final double a;

        if (Double.doubleToRawLongBits(_double) < 0) {
            a = Math.abs(_double);
            sign = -1.0d;
        } else {
            a = _double;
            sign = 1.0d;
        }

        return sign * Math.log(Math.sqrt(a * a + 1.0d) + a);
    }

    // FIXME: Use Apache Commons Math
    private static double atanh(final double _double) {
        final double mult;
        final double a;

        if (Double.doubleToRawLongBits(_double) < 0) {
            a = Math.abs(_double);
            mult = -0.5d;
        } else {
            a = _double;
            mult = 0.5d;
        }
        return mult * Math.log((1.0d + a) / (1.0d - a));
    }
}