/* analyze.c : « austen analyze <conf.tsv> <run.grd> <sortie.json> <gmax.grd> »
 *
 * 1. Suivi du centre de la dépression :
 *    - accrochage au premier instant >= première position de référence, minimum de la pression lissée
 *      dans 350 km (puis 500 km) autour de la référence ;
 *    - suivi heure par heure, vers l'avant puis vers l'arrière : position prolongée à la vitesse
 *      récente, recherche dans 140 km + 45 km par heure d'écart ;
 *    - arrêt si le minimum touche le bord du cercle (pas de centre fermé), s'il sort de l'intérieur
 *      des données valides, ou si la pression saute de plus de 2 hPa + 1,5 hPa par heure.
 *    La pression retenue est la valeur brute la plus basse à moins de 70 km du minimum lissé.
 * 2. Rafale maxi de l'épisode (ep0, ep1] en chaque point → gmax.grd.
 * 3. Séries de rafales aux stations (interpolation bilinéaire) et maxi sur l'épisode.
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "field.h"
#include "json.h"
#include "util.h"

typedef struct {
    const Grid *g;
    int nt;
    const float *raw;          /* [nt][np] pression brute */
    float *smooth;             /* [nt][np] pression lissée */
    uint8_t *interior;         /* [np] */
    double dkm;
} Tracker;

typedef struct {
    int ok;
    double lat, lon, p;
} Hit;

static double ref_lat_lon(const Conf *c, int64_t t, double *lon)
{
    if (t <= c->ref_t[0]) {
        *lon = c->ref_lon[0];
        return c->ref_lat[0];
    }
    for (int k = 0; k + 1 < c->nref; k++)
        if (t >= c->ref_t[k] && t <= c->ref_t[k + 1]) {
            double f = (double)(t - c->ref_t[k]) / (double)(c->ref_t[k + 1] - c->ref_t[k]);
            *lon = c->ref_lon[k] + f * (c->ref_lon[k + 1] - c->ref_lon[k]);
            return c->ref_lat[k] + f * (c->ref_lat[k + 1] - c->ref_lat[k]);
        }
    *lon = c->ref_lon[c->nref - 1];
    return c->ref_lat[c->nref - 1];
}

static void tracker_init(Tracker *tr, const Grid *g, const float *raw, int nt)
{
    size_t np = grid_points(g);
    tr->g = g;
    tr->nt = nt;
    tr->raw = raw;
    tr->dkm = fabs(g->dlat) * KM_PER_DEG;
    double sigma = fmax(0.8, 45.0 / tr->dkm);
    tr->smooth = xmalloc((size_t)nt * np, sizeof(float));
    float *filled = xmalloc(np, sizeof(float));
    for (int t = 0; t < nt; t++) {
        const float *f = raw + (size_t)t * np;
        /* trous comblés au plus proche avant lissage : une valeur fixe créerait un faux minimum */
        fill_nearest(f, filled, g->nx, g->ny);
        float mx = -INFINITY;
        for (size_t k = 0; k < np; k++)
            if (!isnan(filled[k]) && filled[k] > mx)
                mx = filled[k];
        for (size_t k = 0; k < np; k++)
            if (isnan(filled[k]))
                filled[k] = mx;
        gauss_nearest(filled, tr->smooth + (size_t)t * np, g->nx, g->ny, sigma);
    }
    free(filled);
    /* intérieur : à plus de « edge » points (≈ 60 km) des bords de grille et de toute zone sans donnée */
    int edge = (int)fmax(2.0, round(60.0 / tr->dkm));
    uint8_t *fin = xmalloc(np, 1);
    for (int y = 0; y < g->ny; y++)
        for (int x = 0; x < g->nx; x++) {
            size_t k = (size_t)y * (size_t)g->nx + (size_t)x;
            fin[k] = !isnan(raw[k]) && x >= edge && y >= edge && x < g->nx - edge && y < g->ny - edge;
        }
    int *dist = xmalloc(np, sizeof(int));
    city_block(fin, dist, g->nx, g->ny);
    tr->interior = xmalloc(np, 1);
    for (size_t k = 0; k < np; k++)
        tr->interior[k] = dist[k] > edge;
    free(dist);
    free(fin);
}

static Hit find(const Tracker *tr, int t, double la0, double lo0, double R)
{
    const Grid *g = tr->g;
    size_t np = (size_t)g->nx * (size_t)g->ny;
    const float *s = tr->smooth + (size_t)t * np, *raw = tr->raw + (size_t)t * np;
    Hit h = {0, 0, 0, 0};
    /* boîte englobante du disque pour ne pas parcourir toute la grille */
    double dlat = R / KM_PER_DEG + fabs(g->dlat), dlon = R / (KM_PER_DEG * fmax(0.1, cos(la0 * DEG))) + fabs(g->dlon);
    int y0 = (int)floor((la0 - dlat - g->lat0) / g->dlat), y1 = (int)ceil((la0 + dlat - g->lat0) / g->dlat);
    int x0 = (int)floor((lo0 - dlon - g->lon0) / g->dlon), x1 = (int)ceil((lo0 + dlon - g->lon0) / g->dlon);
    y0 = y0 < 0 ? 0 : y0;
    x0 = x0 < 0 ? 0 : x0;
    y1 = y1 > g->ny - 1 ? g->ny - 1 : y1;
    x1 = x1 > g->nx - 1 ? g->nx - 1 : x1;
    double best = INFINITY;
    int bj = -1, bk = -1;
    for (int j = y0; j <= y1; j++)
        for (int k = x0; k <= x1; k++) {
            double la = g->lat0 + j * g->dlat, lo = g->lon0 + k * g->dlon;
            if (dist_km(la, lo, la0, lo0) > R)
                continue;
            double v = s[(size_t)j * (size_t)g->nx + (size_t)k];
            if (v < best) {            /* premier minimum dans l'ordre des lignes, comme numpy.argmin */
                best = v;
                bj = j;
                bk = k;
            }
        }
    if (bj < 0)
        return h;
    if (!tr->interior[(size_t)bj * (size_t)g->nx + (size_t)bk])
        return h;   /* bord du domaine : la dépression n'est pas (ou plus) dans le modèle */
    double la = g->lat0 + bj * g->dlat, lo = g->lon0 + bk * g->dlon;
    if (dist_km(la, lo, la0, lo0) > R - 2.5 * tr->dkm)
        return h;   /* minimum collé au bord du cercle : pente, pas centre fermé */
    /* pression : valeur brute la plus basse à moins de 70 km */
    double p = INFINITY;
    int r = (int)ceil(70.0 / tr->dkm) + 2;
    double coslat = fmax(0.1, cos(la * DEG));
    int rx = (int)ceil(70.0 / (fabs(g->dlon) * KM_PER_DEG * coslat)) + 2;
    for (int j = bj - r; j <= bj + r; j++)
        for (int k = bk - rx; k <= bk + rx; k++) {
            if (j < 0 || k < 0 || j >= g->ny || k >= g->nx)
                continue;
            if (dist_km(g->lat0 + j * g->dlat, g->lon0 + k * g->dlon, la, lo) > 70.0)
                continue;
            double v = raw[(size_t)j * (size_t)g->nx + (size_t)k];
            if (!isnan(v) && v < p)
                p = v;
        }
    if (!isfinite(p))
        return h;
    /* affinage parabolique de la position sur le champ lissé */
    double fla = la, flo = lo;
    if (bj > 0 && bj < g->ny - 1 && bk > 0 && bk < g->nx - 1) {
        size_t c = (size_t)bj * (size_t)g->nx + (size_t)bk;
        double a = s[c - (size_t)g->nx], b = s[c], d = s[c + (size_t)g->nx];
        double den = a - 2 * b + d;
        if (den > 0)
            fla += 0.5 * (a - d) / den * g->dlat;
        a = s[c - 1];
        d = s[c + 1];
        den = a - 2 * b + d;
        if (den > 0)
            flo += 0.5 * (a - d) / den * g->dlon;
    }
    h.ok = 1;
    h.lat = fla;
    h.lon = flo;
    h.p = p;
    return h;
}

static void walk(const Tracker *tr, const int64_t *times, Hit *found, int seed, int dir)
{
    int n = seed;
    int have_prev = 0;
    double pla = 0, plo = 0, pdt = 1;
    for (;;) {
        int m = n + dir;
        if (m < 0 || m >= tr->nt)
            return;
        double dth = fabs((double)(times[m] - times[n])) / 3600.0;
        double la = found[n].lat, lo = found[n].lon;
        if (have_prev) {
            la += (la - pla) * dth / pdt;
            lo += (lo - plo) * dth / pdt;
        }
        double R = 140.0 + 45.0 * dth;
        Hit r = find(tr, m, la, lo, R);
        if (!r.ok || fabs(r.p - found[n].p) > 2.0 + 1.5 * dth)
            return;
        pla = found[n].lat;
        plo = found[n].lon;
        pdt = dth;
        have_prev = 1;
        found[m] = r;
        n = m;
    }
}

int cmd_analyze(int argc, char **argv)
{
    if (argc != 4)
        die("usage : austen analyze <conf.tsv> <run.grd> <sortie.json> <gmax.grd>");
    Conf c;
    conf_read(argv[0], &c);
    Grd *d = grd_read(argv[1]);
    const Grid *g = &d->g;
    size_t np = grid_points(g);

    /* ---- trajectoire, sur les instants de la fenêtre [tr0, tr1] */
    int first = -1, nt = 0;
    for (int k = 0; k < d->nt_msl; k++)
        if (d->t_msl[k] >= c.tr0 && d->t_msl[k] <= c.tr1) {
            if (first < 0)
                first = k;
            nt++;
        }
    Hit *found = xcalloc((size_t)(nt ? nt : 1), sizeof(Hit));
    const int64_t *times = first >= 0 ? d->t_msl + first : NULL;
    Tracker tr = {0};
    if (nt > 0) {
        tracker_init(&tr, g, d->msl + (size_t)first * np, nt);
        int seed = 0;
        for (int k = 0; k < nt; k++)
            if (times[k] >= c.ref_t[0]) {
                seed = k;
                break;
            }
        double lo0, la0 = ref_lat_lon(&c, times[seed], &lo0);
        Hit h = find(&tr, seed, la0, lo0, 350.0);
        if (!h.ok)
            h = find(&tr, seed, la0, lo0, 500.0);
        if (!h.ok)   /* centre pas encore dans le domaine (AROME) : on essaie plus tard */
            for (int k = seed + 1; k < nt && !h.ok; k++) {
                la0 = ref_lat_lon(&c, times[k], &lo0);
                h = find(&tr, k, la0, lo0, 300.0);
                if (h.ok)
                    seed = k;
            }
        if (h.ok) {
            found[seed] = h;
            walk(&tr, times, found, seed, +1);
            walk(&tr, times, found, seed, -1);
        }
    }

    Json j;
    js_init(&j);
    js_obj(&j);
    js_key(&j, "track");
    js_arr(&j);
    char iso[20];
    for (int k = 0; k < nt; k++)
        if (found[k].ok) {
            fmt_iso(times[k], iso);
            js_obj(&j);
            js_key(&j, "t");
            js_str(&j, iso);
            js_key(&j, "lat");
            js_num(&j, found[k].lat, 3);
            js_key(&j, "lon");
            js_num(&j, found[k].lon, 3);
            js_key(&j, "p");
            js_num(&j, found[k].p, 1);
            js_end_obj(&j);
        }
    js_end_arr(&j);

    /* ---- rafales */
    if (d->nt_gust > 0) {
        float *gmax = xmalloc(np, sizeof(float));
        for (size_t k = 0; k < np; k++)
            gmax[k] = NAN;
        for (int t = 0; t < d->nt_gust; t++) {
            if (!(d->t_gust[t] > c.ep0 && d->t_gust[t] <= c.ep1))
                continue;
            const float *f = d->gust + (size_t)t * np;
            for (size_t k = 0; k < np; k++)
                if (!isnan(f[k]) && (isnan(gmax[k]) || f[k] > gmax[k]))
                    gmax[k] = f[k];
        }
        Grd out = {0};
        out.g = *g;
        out.nt_gust = 1;
        int64_t tend = c.ep1;
        out.t_gust = &tend;
        out.gust = gmax;
        out.t_msl = NULL;
        grd_write(argv[3], &out);
        double dom = -INFINITY;
        for (size_t k = 0; k < np; k++)
            if (!isnan(gmax[k]) && gmax[k] > dom)
                dom = gmax[k];

        /* séries aux stations sur (tr0, tr1] */
        js_key(&j, "gust_times");
        js_arr(&j);
        for (int t = 0; t < d->nt_gust; t++)
            if (d->t_gust[t] > c.tr0 && d->t_gust[t] <= c.tr1) {
                fmt_iso(d->t_gust[t], iso);
                js_str(&j, iso);
            }
        js_end_arr(&j);
        js_key(&j, "stations");
        js_obj(&j);
        for (int s = 0; s < c.nst; s++) {
            double mx = NAN;
            int64_t tmx = 0;
            js_key(&j, c.st_id[s]);
            js_obj(&j);
            js_key(&j, "series");
            js_arr(&j);
            for (int t = 0; t < d->nt_gust; t++) {
                if (!(d->t_gust[t] > c.tr0 && d->t_gust[t] <= c.tr1))
                    continue;
                double v = bilinear(g, d->gust + (size_t)t * np, c.st_lat[s], c.st_lon[s]);
                v = isnan(v) ? NAN : round(v * 10.0) / 10.0;
                js_num(&j, v, 1);
                if (!isnan(v) && d->t_gust[t] > c.ep0 && d->t_gust[t] <= c.ep1 && (isnan(mx) || v > mx)) {
                    mx = v;
                    tmx = d->t_gust[t];
                }
            }
            js_end_arr(&j);
            js_key(&j, "max");
            js_num(&j, mx, 1);
            js_key(&j, "tmax");
            if (isnan(mx))
                js_null(&j);
            else {
                fmt_iso(tmx, iso);
                js_str(&j, iso);
            }
            js_end_obj(&j);
        }
        js_end_obj(&j);
        js_key(&j, "gmax_domain");
        js_num(&j, dom, 1);
        free(gmax);
    }
    js_end_obj(&j);
    js_save(&j, argv[2]);
    js_free(&j);
    free(found);
    free(tr.smooth);
    free(tr.interior);
    grd_free(d);
    return 0;
}