/* decode.c : « austen decode <rep_grib> <sortie.grd> »
 * Lit tous les *.grib2 du répertoire, garde la pression au niveau de la mer (hPa) et la rafale à
 * 10 m (km/h), découpe le cadre Atlantique-Europe et range latitudes et longitudes en ordre croissant.
 * Paramètres OMM : pression réduite (0,3,1) ou pression au niveau moyen de la mer (0,3,0, ECMWF),
 * rafale (0,2,22) ; une rafale maxi (traitement statistique 2) l'emporte sur une valeur instantanée. */
#define _DEFAULT_SOURCE
#include <dirent.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "grib2.h"
#include "util.h"

static const double BOX[4] = {-40.0, 25.0, 33.0, 66.0};   /* lon mini, lon maxi, lat mini, lat maxi */

enum { V_MSL, V_GUST };

typedef struct {
    int64_t t;
    int stat;
    float *v;
} Slot;

typedef struct {
    Slot s[2][MAX_TIMES];
    int n[2];
    int set;                     /* grille fixée par le premier champ */
    int ni, nj, scan;
    double la1, lo1, di, dj;
    int *ix, *jy;                /* indices source retenus, déjà triés */
    Grid g;
} Acc;

static int want(int disc, int cat, int num, void *ctx)
{
    (void)ctx;
    return disc == 0 && ((cat == 3 && (num == 0 || num == 1)) || (cat == 2 && num == 22));
}

static int cmp_str(const void *a, const void *b) { return strcmp(*(char *const *)a, *(char *const *)b); }

typedef struct { double lon; int i; } LonIdx;
static int cmp_lon(const void *a, const void *b)
{
    double x = ((const LonIdx *)a)->lon, y = ((const LonIdx *)b)->lon;
    return x < y ? -1 : x > y;
}

static void set_grid(Acc *a, const G2Field *f)
{
    a->ni = f->ni;
    a->nj = f->nj;
    a->scan = f->scan;
    a->la1 = f->la1;
    a->lo1 = f->lo1;
    a->di = f->di;
    a->dj = f->dj;
    if (f->scan & 0x80)
        die("balayage est → ouest non géré");
    /* longitudes ramenées dans [-180, 180), triées */
    LonIdx *li = xmalloc((size_t)f->ni, sizeof *li);
    int nl = 0;
    for (int i = 0; i < f->ni; i++) {
        double lo = fmod(f->lo1 + i * f->di, 360.0);
        if (lo >= 180.0)
            lo -= 360.0;
        if (lo < -180.0)
            lo += 360.0;
        if (lo >= BOX[0] - 1e-9 && lo <= BOX[1] + 1e-9)
            li[nl++] = (LonIdx){lo, i};
    }
    qsort(li, (size_t)nl, sizeof *li, cmp_lon);
    a->ix = xmalloc((size_t)nl, sizeof(int));
    for (int k = 0; k < nl; k++)
        a->ix[k] = li[k].i;
    /* latitudes : sens selon le bit « j positif » du mode de balayage */
    double sgn = (f->scan & 0x40) ? 1.0 : -1.0;
    a->jy = xmalloc((size_t)f->nj, sizeof(int));
    int nlat = 0;
    for (int j = 0; j < f->nj; j++) {
        int jj = sgn > 0 ? j : f->nj - 1 - j;            /* parcours du sud au nord */
        double la = f->la1 + sgn * jj * f->dj;
        if (la >= BOX[2] - 1e-9 && la <= BOX[3] + 1e-9)
            a->jy[nlat++] = jj;
    }
    if (nl < 2 || nlat < 2)
        die("le cadre ne recoupe pas la grille du modèle");
    a->g.nx = nl;
    a->g.ny = nlat;
    a->g.lon0 = li[0].lon;
    a->g.dlon = f->di;
    a->g.lat0 = f->la1 + sgn * a->jy[0] * f->dj;
    a->g.dlat = f->dj;
    free(li);
    a->set = 1;
}

static void add(Acc *a, const G2Field *f)
{
    if (!a->set)
        set_grid(a, f);
    if (f->ni != a->ni || f->nj != a->nj || f->scan != a->scan || fabs(f->la1 - a->la1) > 1e-6 ||
        fabs(f->lo1 - a->lo1) > 1e-6 || fabs(f->di - a->di) > 1e-9 || fabs(f->dj - a->dj) > 1e-9)
        die("grille différente d'un message à l'autre");
    int var = f->category == 3 ? V_MSL : V_GUST;
    double scale = var == V_MSL ? 0.01 : 3.6;
    /* même heure de validité : la rafale maxi l'emporte, sinon le dernier message */
    int slot = -1;
    for (int k = 0; k < a->n[var]; k++)
        if (a->s[var][k].t == f->valid_time)
            slot = k;
    if (slot >= 0 && var == V_GUST && a->s[var][slot].stat == 2 && f->stat_proc != 2)
        return;
    if (slot < 0) {
        if (a->n[var] >= MAX_TIMES)
            die("trop d'échéances");
        slot = a->n[var]++;
    } else {
        free(a->s[var][slot].v);
    }
    size_t np = (size_t)a->g.nx * (size_t)a->g.ny;
    float *v = xmalloc(np, sizeof(float));
    for (int y = 0; y < a->g.ny; y++)
        for (int x = 0; x < a->g.nx; x++) {
            double val = f->values[(size_t)a->jy[y] * (size_t)f->ni + (size_t)a->ix[x]];
            v[(size_t)y * (size_t)a->g.nx + (size_t)x] = isnan(val) ? NAN : (float)(val * scale);
        }
    a->s[var][slot] = (Slot){f->valid_time, f->stat_proc, v};
}

static int cmp_slot(const void *a, const void *b)
{
    int64_t x = ((const Slot *)a)->t, y = ((const Slot *)b)->t;
    return x < y ? -1 : x > y;
}

int cmd_decode(int argc, char **argv)
{
    if (argc != 2)
        die("usage : austen decode <rep_grib> <sortie.grd>");
    DIR *d = opendir(argv[0]);
    if (!d)
        die("répertoire illisible : %s", argv[0]);
    char **names = NULL;
    size_t nn = 0;
    struct dirent *e;
    while ((e = readdir(d))) {
        size_t l = strlen(e->d_name);
        if (e->d_name[0] != '.' && l > 6 && !strcmp(e->d_name + l - 6, ".grib2")) {
            names = xrealloc(names, nn + 1, sizeof *names);
            names[nn++] = xstrdup(e->d_name);
        }
    }
    closedir(d);
    qsort(names, nn, sizeof *names, cmp_str);
    Acc *a = xcalloc(1, sizeof *a);
    int nerr = 0;
    for (size_t k = 0; k < nn; k++) {
        char path[4096];
        if (snprintf(path, sizeof path, "%s/%s", argv[0], names[k]) >= (int)sizeof path)
            die("chemin trop long");
        size_t len;
        uint8_t *buf = read_file(path, (size_t)2 << 30, &len);
        size_t pos = 0;
        char err[256];
        for (;;) {
            G2Field f;
            int rc = g2_next(buf, len, &pos, &f, want, NULL, err, sizeof err);
            if (rc == 0)
                break;
            if (rc < 0) {
                fprintf(stderr, "austen decode : %s : %s\n", names[k], err);
                nerr++;
                continue;
            }
            /* pression : seulement au niveau moyen de la mer (surface 101), jamais la pression au sol */
            if (f.values && !(f.category == 3 && f.surface != 101))
                add(a, &f);
            g2_free(&f);
        }
        free(buf);
        free(names[k]);
    }
    free(names);
    if (nerr)
        die("%d message(s) illisible(s) : run refusé", nerr);
    if (!a->set || a->n[V_MSL] == 0)
        die("aucune pression trouvée");
    Grd out = {0};
    out.g = a->g;
    size_t np = grid_points(&out.g);
    for (int var = 0; var < 2; var++) {
        int n = a->n[var];
        qsort(a->s[var], (size_t)n, sizeof(Slot), cmp_slot);
        int64_t *t = xmalloc((size_t)n, sizeof(int64_t));
        float *v = xmalloc((size_t)n * np, sizeof(float));
        for (int k = 0; k < n; k++) {
            t[k] = a->s[var][k].t;
            memcpy(v + (size_t)k * np, a->s[var][k].v, np * sizeof(float));
            free(a->s[var][k].v);
        }
        if (var == V_MSL) {
            out.nt_msl = n;
            out.t_msl = t;
            out.msl = v;
        } else {
            out.nt_gust = n;
            out.t_gust = t;
            out.gust = v;
        }
    }
    grd_write(argv[1], &out);
    printf("{\"nx\":%d,\"ny\":%d,\"msl\":%d,\"gust\":%d}\n", out.g.nx, out.g.ny, out.nt_msl, out.nt_gust);
    free(a->ix);
    free(a->jy);
    free(a);
    free(out.t_msl);
    free(out.t_gust);
    free(out.msl);
    free(out.gust);
    return 0;
}