<commit_before>/*
 Copyright 2015 <NAME>, <NAME>, <NAME>

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

 http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
package com.jhlabs.map.proj;

import java.awt.geom.Point2D;

/**
 * The Hufnagel projection family was introduced by <NAME> in
 * "Hufnagel, H. 1989. Ein System unecht-zylindrischer Kartennetze fr
 * Erdkarten. Kartographische Nachrichten, 39(3), 8996." All projections are
 * equal-area. Some parameter combinations break the projection into pieces.
 * This is due to the difficulty of finding good initial guesses for the
 * iterative computation of coordinates. Implementation by <NAME> and
 * <NAME>, Oregon State University, with help by <NAME>rebe,
 * Mapthematics. November 2014 to July 2015.
 *
 * @author <NAME>, <NAME>
 */
public class HufnagelProjection extends PseudoCylindricalProjection {

    private static final double EPS = 1e-6;
    private static final int MAX_ITER = 100;
    private static final int LUT_SIZE = 101;

    private double A, B, psiMax, aspectRatio;
    private double[] latLUT, yLUT, psiLUT;
    private double K2, K, C;

    public HufnagelProjection() {
        this(2. / 3., 1. / 3., Math.PI / 4., 2.0);
    }

    public HufnagelProjection(double A, double B, double psiMax, double aspectRatio) {
        this.A = A;
        this.B = B;
        this.aspectRatio = aspectRatio;
        this.psiMax = psiMax;
        initializeHufnagel();
    }

    private void initializeHufnagel() {
        if (psiMax == 0) {
            K = Math.sqrt(aspectRatio / Math.PI);
        } else {
            Point2D.Double outPoint = new Point2D.Double();
            double width;

            K2 = (4. * Math.PI) / (2. * psiMax + (1. + A - B / 2.) * Math.sin(2. * psiMax)
                    + ((A + B) / 2.) * Math.sin(4. * psiMax) + (B / 2.) * Math.sin(6. * psiMax));
            K = Math.sqrt(K2);
            C = 1.;

            initializeLookUpTables();

            this.project(Math.PI, 0., outPoint);
            width = outPoint.x;
            this.project(0., Math.PI / 2., outPoint);
            C = Math.sqrt(aspectRatio / (width / outPoint.y));
        }
    }

    public void initializeLookUpTables() {
        latLUT = new double[LUT_SIZE];
        yLUT = new double[LUT_SIZE];
        psiLUT = new double[LUT_SIZE];

        for (int i = 0; i < LUT_SIZE; i++) {
            double psi = psiMax * i / (LUT_SIZE - 1.);
            double phi;
            if (i == 0) {
                phi = 0.0;
            } else if (i == LUT_SIZE - 1) {
                phi = Math.PI / 2.;
            } else {
                double sin2Psi = Math.sin(2. * psi);
                double sin4Psi = Math.sin(4. * psi);
                double sin6Psi = Math.sin(6. * psi);
                double sinPhi = 0.25 / Math.PI * K2 * (2 * psi + (1 + A - 0.5 * B) * sin2Psi
                        + 0.5 * (A + B) * sin4Psi + 0.5 * B * sin6Psi);
                if (Math.abs(sinPhi) > 1) {
                    phi = sinPhi > 0 ? Math.PI / 2 : -Math.PI / 2;
                } else {
                    phi = Math.asin(sinPhi);
                }
            }

            double r = Math.sqrt(1 + A * Math.cos(2 * psi) + B * Math.cos(4 * psi));
            double y = r * Math.sin(psi);
            if (i > 0) {
                if (y < yLUT[i - 1] || phi < latLUT[i - 1]) {
                    y = yLUT[i - 1];
                    phi = latLUT[i - 1];
                }
            }

            latLUT[i] = phi;
            yLUT[i] = y;
            psiLUT[i] = psi;
        }
    }

    private double approximatePsiFromTable(double v, double[] table) {
        double w, psi, vAbs = Math.abs(v);
        int imid, imin = 0, imax = table.length;

        while (true) {
            imid = (int) Math.floor((imin + imax) / 2);
            if (imid == imin) {
                // This also handles abs(phi) == latitudeTable[0] because mid must == min.
                break;
            } else if (vAbs > table[imid]) {
                imin = imid;
            } else {
                // abs(phi) < latitudeTable[mid], or abs(phi) == latitudeTable[mid] and mid  0
                imax = imid;
            }
        }

        w = (vAbs - table[imin]) / (table[imin + 1] - table[imin]);
        psi = w * (psiLUT[imin + 1] - psiLUT[imin]) + psiLUT[imin];

        return v < 0 ? -psi : psi;
    }

    private double psiFromPhi(double phi) {
        double psi = approximatePsiFromTable(phi, latLUT);
        int i = 0;
        if (Math.abs(psi) < Math.PI / 2.) {
            while (true) {
                double psi_x_2 = psi * 2.;
                double deltaPsiNominator = (.25 * K2) * (psi_x_2 + (1. + A - .5 * B)
                        * Math.sin(psi_x_2) + (0.5 * (A + B)) * Math.sin(2. * psi_x_2)
                        + (0.5 * B) * Math.sin(3. * psi_x_2)) - Math.PI * Math.sin(phi);
                if (Math.abs(deltaPsiNominator) < EPS) {
                    break;
                }
                double deltaPsiDenominator = (.5 * K2) * (1. + (1. + A - .5 * B)
                        * Math.cos(psi_x_2) + (A + B) * Math.cos(2. * psi_x_2)
                        + (3. * .5 * B) * Math.cos(3. * psi_x_2));
                double deltaPsi = deltaPsiNominator / deltaPsiDenominator;
                if (Double.isNaN(deltaPsi) || Double.isInfinite(deltaPsi) || i++ > MAX_ITER) {
                    return Double.NaN;
                }
                psi = psi - deltaPsi;
            }
        }
        return psi;
    }

    private double psiFromY(double y) {
        double psi = approximatePsiFromTable(y, yLUT);
        int i = 0;
        if (Math.abs(psi) < Math.PI / 2.) {
            while (true) {
                double cosPsi = Math.cos(psi);
                double cos2Psi = Math.cos(2 * psi);
                double cos4Psi = Math.cos(4 * psi);
                double r = Math.sqrt(1 + A * cos2Psi + B * cos4Psi);
                double deltaPsiNominator = r * Math.sin(psi) - y;
                if (Math.abs(deltaPsiNominator) < EPS) {
                    break;
                }
                double deltaPsiDenominator = (cosPsi * (1 - A + 2 * B
                        + (2 * A - 4 * B) * cos2Psi + 3 * B * cos4Psi)) / r;
                double deltaPsi = deltaPsiNominator / deltaPsiDenominator;
                if (Double.isNaN(deltaPsi) || Double.isInfinite(deltaPsi) || i++ > MAX_ITER) {
                    return Double.NaN;
                }
                psi = psi - deltaPsi;
            }
        }
        return psi;
    }

    @Override
    public Point2D.Double project(double lplam, double lpphi, Point2D.Double out) {
        if (psiMax == 0) {<commit_msg>Simplified Hufnagel<commit_after>/*
 Copyright 2015 <NAME>, <NAME>, <NAME>

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

 http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
package com.jhlabs.map.proj;

import java.awt.geom.Point2D;

/**
 * The Hufnagel projection family was introduced by <NAME> in
 * "Hufnagel, H. 1989. Ein System unecht-zylindrischer Kartennetze fr
 * Erdkarten. Kartographische Nachrichten, 39(3), 8996." All projections are
 * equal-area. Some parameter combinations break the projection into pieces.
 * This is due to the difficulty of finding good initial guesses for the
 * iterative computation of coordinates. Implementation by <NAME> and
 * <NAME>, Oregon State University, with help by <NAME>,
 * Mapthematics. November 2014 to July 2015.
 *
 * @author <NAME>, <NAME>
 */
public class HufnagelProjection extends PseudoCylindricalProjection {

    private static final double EPS = 1e-6;
    private static final int MAX_ITER = 100;
    private static final int LUT_SIZE = 101;

    private double A, B, psiMax, aspectRatio;
    private double[] latLUT, yLUT, psiLUT;
    private double K2, K, C;

    public HufnagelProjection() {
        this(2. / 3., 1. / 3., Math.PI / 4., 2.0);
    }

    public HufnagelProjection(double A, double B, double psiMax, double aspectRatio) {
        this.A = A;
        this.B = B;
        this.aspectRatio = aspectRatio;
        this.psiMax = psiMax;
        initializeHufnagel();
    }

    private void initializeHufnagel() {
        if (psiMax == 0) {
            K = Math.sqrt(aspectRatio / Math.PI);
        } else {
            Point2D.Double outPoint = new Point2D.Double();
            double width;

            K2 = (4. * Math.PI) / (2. * psiMax + (1. + A - B / 2.) * Math.sin(2. * psiMax)
                    + ((A + B) / 2.) * Math.sin(4. * psiMax) + (B / 2.) * Math.sin(6. * psiMax));
            K = Math.sqrt(K2);
            C = 1.;

            initializeLookUpTables();

            this.project(Math.PI, 0., outPoint);
            width = outPoint.x;
            this.project(0., Math.PI / 2., outPoint);
            C = Math.sqrt(aspectRatio / (width / outPoint.y));
        }
    }

    public void initializeLookUpTables() {
        latLUT = new double[LUT_SIZE];
        yLUT = new double[LUT_SIZE];
        psiLUT = new double[LUT_SIZE];

        for (int i = 0; i < LUT_SIZE; i++) {
            double psi = psiMax * i / (LUT_SIZE - 1.);
            double phi;
            if (i == 0) {
                phi = 0.0;
            } else if (i == LUT_SIZE - 1) {
                phi = Math.PI / 2.;
            } else {
                double sin2Psi = Math.sin(2. * psi);
                double sin4Psi = Math.sin(4. * psi);
                double sin6Psi = Math.sin(6. * psi);
                double sinPhi = 0.25 / Math.PI * K2 * (2 * psi + (1 + A - 0.5 * B) * sin2Psi
                        + 0.5 * (A + B) * sin4Psi + 0.5 * B * sin6Psi);
                if (Math.abs(sinPhi) > 1) {
                    phi = sinPhi > 0 ? Math.PI / 2 : -Math.PI / 2;
                } else {
                    phi = Math.asin(sinPhi);
                }
            }

            double r = Math.sqrt(1 + A * Math.cos(2 * psi) + B * Math.cos(4 * psi));
            double y = r * Math.sin(psi);
            if (i > 0) {
                if (y < yLUT[i - 1] || phi < latLUT[i - 1]) {
                    y = yLUT[i - 1];
                    phi = latLUT[i - 1];
                }
            }

            latLUT[i] = phi;
            yLUT[i] = y;
            psiLUT[i] = psi;
        }
    }

    private double approximatePsiFromTable(double v, double[] table) {
        double w, psi, vAbs = Math.abs(v);
        int imid, imin = 0, imax = table.length;

        while (true) {
            imid = (int) Math.floor((imin + imax) / 2);
            if (imid == imin) {
                // This also handles abs(phi) == latitudeTable[0] because mid must == min.
                break;
            } else if (vAbs > table[imid]) {
                imin = imid;
            } else {
                // abs(phi) < latitudeTable[mid], or abs(phi) == latitudeTable[mid] and mid  0
                imax = imid;
            }
        }

        w = (vAbs - table[imin]) / (table[imin + 1] - table[imin]);
        psi = w * (psiLUT[imin + 1] - psiLUT[imin]) + psiLUT[imin];

        return v < 0 ? -psi : psi;
    }

    private double psiFromPhi(double phi) {
        double psi = approximatePsiFromTable(phi, latLUT);
        int i = 0;
        while (true) {
            double psi_x_2 = psi * 2.;
            double deltaPsiNominator = (.25 * K2) * (psi_x_2 + (1. + A - .5 * B)
                    * Math.sin(psi_x_2) + (0.5 * (A + B)) * Math.sin(2. * psi_x_2)
                    + (0.5 * B) * Math.sin(3. * psi_x_2)) - Math.PI * Math.sin(phi);
            if (Math.abs(deltaPsiNominator) < EPS) {
                break;
            }
            double deltaPsiDenominator = (.5 * K2) * (1. + (1. + A - .5 * B)
                    * Math.cos(psi_x_2) + (A + B) * Math.cos(2. * psi_x_2)
                    + (3. * .5 * B) * Math.cos(3. * psi_x_2));
            double deltaPsi = deltaPsiNominator / deltaPsiDenominator;
            if (Double.isNaN(deltaPsi) || Double.isInfinite(deltaPsi) || i++ > MAX_ITER) {
                return Double.NaN;
            }
            psi = psi - deltaPsi;
        }
        return psi;
    }

    private double psiFromY(double y) {
        double psi = approximatePsiFromTable(y, yLUT);
        int i = 0;
        while (true) {
            double cosPsi = Math.cos(psi);
            double cos2Psi = Math.cos(2 * psi);
            double cos4Psi = Math.cos(4 * psi);
            double r = Math.sqrt(1 + A * cos2Psi + B * cos4Psi);
            double deltaPsiNominator = r * Math.sin(psi) - y;
            if (Math.abs(deltaPsiNominator) < EPS) {
                break;
            }
            double deltaPsiDenominator = (cosPsi * (1 - A + 2 * B
                    + (2 * A - 4 * B) * cos2Psi + 3 * B * cos4Psi)) / r;
            double deltaPsi = deltaPsiNominator / deltaPsiDenominator;
            if (Double.isNaN(deltaPsi) || Double.isInfinite(deltaPsi) || i++ > MAX_ITER) {
                return Double.NaN;
            }
            psi = psi - deltaPsi;
        }
        return psi;
    }

    @Override
    public Point2D.Double project(double lplam, double lpphi, Point2D.Double out) {
        if (psiMax == 0) {
{'fraction_non_alphanumeric': 0.12421383647798742, 'fraction_numerical': 0.022584333905088624, 'mean_word_length': 1.5516046681254558, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3513378', 'n_tokens_mistral': 4810, 'n_tokens_neox': 4488, 'n_words': 1273}