/*
  Created by Jacques Bertin in 1953, this projection was the go-to choice
  of the French cartographic school when they wished to represent phenomena
  on a global scale.

  Formula designed by Philippe Rivière, 2017.
  https://visionscarto.net/bertin-projection-1953

  Port to PROJ by Philippe Rivière, 21 September 2018
*/

#include <errno.h>
#include <math.h>

#include "proj.h"
#include "proj_internal.h"

PROJ_HEAD(bertin1953, "Bertin 1953")
"\n\tMisc Sph no inv.";

namespace { // anonymous namespace
struct pj_bertin1953 {
    double cos_delta_phi, sin_delta_phi, cos_delta_gamma, sin_delta_gamma,
        deltaLambda;
};
} // anonymous namespace

static PJ_XY bertin1953_s_forward(PJ_LP lp, PJ *P) {
    PJ_XY xy = {0.0, 0.0};
    struct pj_bertin1953 *Q = static_cast<struct pj_bertin1953 *>(P->opaque);

    double fu = 1.4, k = 12., w = 1.68, d;

    /* Rotate */
    double cosphi, x, y, z, z0;
    lp.lam += PJ_TORAD(-16.5);
    cosphi = cos(lp.phi);
    x = cos(lp.lam) * cosphi;
    y = sin(lp.lam) * cosphi;
    z = sin(lp.phi);
    z0 = z * Q->cos_delta_phi + x * Q->sin_delta_phi;
    lp.lam = atan2(y * Q->cos_delta_gamma - z0 * Q->sin_delta_gamma,
                   x * Q->cos_delta_phi - z * Q->sin_delta_phi);
    z0 = z0 * Q->cos_delta_gamma + y * Q->sin_delta_gamma;
    lp.phi = asin(z0);

    lp.lam = adjlon(lp.lam);

    /* Adjust pre-projection */
    if (lp.lam + lp.phi < -fu) {
        d = (lp.lam - lp.phi + 1.6) * (lp.lam + lp.phi + fu) / 8.;
        lp.lam += d;
        lp.phi -= 0.8 * d * sin(lp.phi + M_PI / 2.);
    }

    /* Project with Hammer (1.68,2) */
    cosphi = cos(lp.phi);
    d = sqrt(2. / (1. + cosphi * cos(lp.lam / 2.)));
    xy.x = w * d * cosphi * sin(lp.lam / 2.);
    xy.y = d * sin(lp.phi);

    /* Adjust post-projection */
    d = (1. - cos(lp.lam * lp.phi)) / k;
    if (xy.y < 0.) {
        xy.x *= 1. + d;
    }
    if (xy.y > 0.) {
        xy.y *= 1. + d / 1.5 * xy.x * xy.x;
    }

    return xy;
}

PJ *PJ_PROJECTION(bertin1953) {
    struct pj_bertin1953 *Q = static_cast<struct pj_bertin1953 *>(
        calloc(1, sizeof(struct pj_bertin1953)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    P->lam0 = 0;
    P->phi0 = PJ_TORAD(-42.);

    Q->cos_delta_phi = cos(P->phi0);
    Q->sin_delta_phi = sin(P->phi0);
    Q->cos_delta_gamma = 1.;
    Q->sin_delta_gamma = 0.;

    P->es = 0.;
    P->fwd = bertin1953_s_forward;

    return P;
}