/* See zf2_flt.h */
/* Last edited on 2020-12-07 23:10:01 by jstolfi */

#define _GNU_SOURCE
#include <stdio.h>
#include <math.h>

#include <ia.h>
#include <affirm.h>
#include <flt.h>

#include <zf2_flt.h>

/*** PROTOTYPES FOR INTERNAL ROUTINES ***/

/*** IMPLEMENTATIONS ***/

void zf2_flt_grid(
    zf2_flt_eval_func_t *eval,
    zf2_flt_report_func_t *report,
    Interval xd,
    Interval yd,
    int nx,
    int ny
  )
  {
    Interval xv, yv;
    int ix, iy;
    Float xm, ym;
    Float f00, f01, f10, f11, fmm;

    ROUND_NEAR;

    for (iy = 0; iy < ny; iy++)
      {
        /* Compute the {y} interval and mudpoint of cell row {iy}: */
        ROUND_DOWN;
        yv.lo = yd.lo + ((yd.hi - yd.lo)*((float)iy))/((float)ny);
        ROUND_UP;
        yv.hi = yd.lo + ((yd.hi - yd.lo)*((float)(iy+1)))/((float)ny);
        ROUND_NEAR;
        ym = (yv.lo + yv.hi)/Two;
            
        for (ix = 0; ix < nx; ix++)
          {
            /* Compute the {x} interval and mudpoint of cell {ix} in row: */
            ROUND_DOWN;
            xv.lo = xd.lo + ((xd.hi - xd.lo)*((float)ix))/((float)nx);
            ROUND_UP;
            xv.hi = xd.lo + ((xd.hi - xd.lo)*((float)(ix+1)))/((float)nx);
            ROUND_NEAR;
            xm = (xv.lo + xv.hi)/Two;

            /* Evaluate the function at the corners and midpoint of cell: */
            f00 = eval(xv.lo, yv.lo);
            f01 = eval(xv.lo, yv.hi);
            f10 = eval(xv.hi, yv.lo);
            f11 = eval(xv.hi, yv.hi);
            fmm = eval(xm, ym);

            /* Process the four triangles: */
            Float xp[3], yp[3], zp[3];
            xp[0] = xm; yp[0] = ym;  zp[0] = fmm;
            
            xp[1] = xv.lo; yp[1] = yv.lo; zp[1] = f00;
            xp[2] = xv.hi; yp[2] = yv.lo; zp[2] = f10;
            report(xp, yp, zp);

            xp[1] = xv.hi; yp[1] = yv.lo; zp[1] = f10;
            xp[2] = xv.hi; yp[2] = yv.hi; zp[2] = f11;
            report(xp, yp, zp);

            xp[1] = xv.hi; yp[1] = yv.hi; zp[1] = f11;
            xp[2] = xv.lo; yp[2] = yv.hi; zp[2] = f01;
            report(xp, yp, zp);

            xp[1] = xv.lo; yp[1] = yv.hi; zp[1] = f01;
            xp[2] = xv.lo; yp[2] = yv.lo; zp[2] = f00;
            report(xp, yp, zp);
          }
      }
  }