#define PROG_NAME "test_dg_plot"
#define PROG_DESC "test of dyadic grids: plotting"
#define PROG_VERS "1.0"

/* Last edited on 2012-07-22 14:18:04 by stolfilocal */

#define test_dg_plot_COPYRIGHT "Copyright © 2007 by the State University of Campinas (UNICAMP)."

#define PROG_HELP \
  "  " PROG_NAME ""

#define PROG_INFO \
  "NAME\n" \
  "  " PROG_NAME " - " PROG_DESC "\n" \
  "\n" \
  "SYNOPSIS\n" \
  PROG_HELP "\n" \
  "\n" \
  "DESCRIPTION\n" \
  "  This program generates a sample plot of a function" \
  " defined on a dyadic grid.  The function is a 2D Bézier patch.\n" \
  "\n" \
  "SEE ALSO\n" \
  "  test_approx(1), test_interp(1).\n" \
  "\n" \
  "AUTHOR\n" \
  "  Started 2005-02-12 or earlier by Jorge Stolfi, IC-UNICAMP.\n" \
  "\n" \
  "MODIFICATION HISTORY\n" \
  "  2007-10-14 updated to new version of dyadic grid libs, by J. Stolfi.\n" \
  "\n" \
  "WARRANTY\n" \
  argparser_help_info_NO_WARRANTY "\n" \
  "\n" \
  "RIGHTS\n" \
  "  " test_approx_COPYRIGHT ".\n" \
  "\n" \
  argparser_help_info_STANDARD_RIGHTS

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

#include <bool.h>
#include <argparser.h>
#include <intheap.h>
#include <jsmath.h>
#include <affirm.h>
#include <rn.h>
#include <pswr.h>
#include <pswr_color.h>
#include <pswr_iso.h>

#include <udg_grid.h>
#include <udg_spline.h>
#include <udg_pulse.h>

#include <dg_grid.h>
#include <dg_tree.h>
#include <dg_tent.h>
#include <dg_tent_basis.h>
#include <dg_plot.h>
#include <dg_makegrid.h>

#include <bz_patch.h>
#include <bz_patch_plot_2D.h>

/* Dimension of grid's domain: */
#define DDIM (2)  

/* Dimension of grid's range space {(X,Y,F(X,Y))}: */
#define MAX_RDIM (5)

/* Flatness tolerance (mm): */
#define TOL_MM (0.5)

/* TRUE to debug {make_shape} instead of {dg_enum_faces}: */
#define DEBUG_SHAPE FALSE

/* INTERNAL PROTOTYPES */

int main(int argc, char **argv);
  
dg_tree_t call_make_grid(mdg_dim_t d, mdg_rank_t maxDepth, bz_patch_t *b);
  /* Creates a random finite dyadic grid of dimension {d} and
    maximum depth {maxDepth}, whose root cell has shape {b}. */

double tol(double x, double y);
  /* Maximum box radius allowed near {p}. */ 

double inside(double x, double y, double rad);
  /* TRUE if a box of radius {rad} centered at {p} intersects the 
    desired domain. */ 

double dsoft(double dx, double dy, double eps);
  /* Length of {(dx,dy)}, cooked to be at least {eps}. */

void test_plot
  ( char *name, 
    bz_degree_t g, 
    bz_patch_rdim_t nb, 
    bz_patch_rdim_t nf,
    dg_plot_2D_func_t *func,
    mdg_rank_t treeDp, 
    mdg_rank_t plotDp
  );
  /* Makes a test plot for a function {F} over a random tree of
    maximum depth {treeDp}, with maximum plotting split depth
    {plotDp}.
    
    The function {F} consists of a Bézier patch of degree {g} and
    range dimension {nb}, laminated with the function {func} which has
    range dimension {nf}. The Bézier range dimension {nb} must be at
    least 2, and the sum {nb + nf} must not exceed 5. */

/* IMPLEMENTATIONS */

int main(int argc, char **argv)
  {   
    auto void tfunc(double x[], int nx, double f[], int nf);
      /* A black-box function to plot. */
    
    void tfunc(double x[], int nx, double f[], int nf)
      { 
        demand(nx == DDIM, "bad {nx}");
        if (nf > 0) { double r = hypot(x[0],x[1]); f[0] = cos(7 * 2*M_PI * r); }
        if (nf > 1) { f[1] = cos(3 * 2*M_PI * x[0]); }
        if (nf > 2) { f[2] = cos(5 * 2*M_PI * x[1]); }
      }
    
    /* test_plot("test", 1, 3, 0, &tfunc,  1,  3); */
    test_plot("test", 1, 4, 0, &tfunc,  5,  5);
    /* test_plot("test", 1, 3, 0, &tfunc, 14, 14); */
    /* test_plot("test", 2, 3, 0, &tfunc,  5,  5); */
    /* test_plot("test", 1, 4, 0, &tfunc, 14, 12); */
    /* test_plot("test", 1, 5, 0, &tfunc, 14, 12); */
    /* test_plot("test", 1, 2, 3, &tfunc, 14, 12); */
    /* test_plot("test", 1, 2, 1, &tfunc, 14, 14); */
    fprintf(stderr, "OK so far...\n");
    return(0);
  }

void test_plot
  ( char *name, 
    bz_degree_t g, 
    bz_patch_rdim_t nb, 
    bz_patch_rdim_t nf,
    dg_plot_2D_func_t *func,
    mdg_rank_t treeDp, 
    mdg_rank_t plotDp
  )
  {
    bool_t epsformat = FALSE;
    bz_patch_t b = bz_patch_make_test(DDIM, nb, g);
    dg_tree_node_t *grid_root = call_make_grid(DDIM, treeDp, &b);

    bool_t isolines = (nb + nf - DDIM == 1);
    
    affirm(b.m == DDIM, "wrong dimension");
    fprintf(stderr, "grid cells = %d\n", grid_root->nodes);
    fprintf(stderr, "root shape =\n");
    bz_patch_print(stderr, &b, "%6.2f");
    /* bz_patch_t_dev_print(stderr, sh->dev, b.m, "%.2f"); */
    /* fprintf(stderr, "\n"); */

    /* Get the bounding box of the shape. Note that {bz->n} may exceed {DDIM}: */
    interval_t bbox[MAX_RDIM];
    bz_patch_compute_bbox(&b, bbox);

    PSStream *ps = dg_plot_init(name, epsformat, "a4", "Plot test", bbox);
    
    udg_grid_size_t psz[DDIM];
    int i;
    for (i = 0; i < DDIM; i++) { psz[i] = 1; } 
    
    if (DEBUG_SHAPE)
      { dg_plot_slow(ps, &b); }
    else
      { dg_plot_2D_func(ps, &b, nf, func, grid_root, plotDp, isolines); }
    dg_plot_finish(ps);
  }

#define x_well_1 (2.5)
#define y_well_1 (4.5)

#define x_well_2 (4.0)
#define y_well_2 (5.0)

#define y_fault (2.5)
#define r_fault (0.5)
#define w_fault (1.5)

double tol(double x, double y)
  { double d1 = dsoft(x - x_well_1, y - y_well_1, 0.01);
    // double d1 = 100000.00;
    double d2 = dsoft(x - x_well_2, y - y_well_2, 0.01);
    // double d2 = 100000.00;
    double yfx = y_fault + r_fault*sin(w_fault*x);
    double d3 = dsoft(0.0, y - yfx, 0.01);
    // double d3 = 100000.00;
    double maxr = 0.25 * 3.0/(1.0/d1 + 1.0/d2 + 1.0/d3);
    return maxr;
  }

double inside(double x, double y, double rad)
  { double dctr = hypot(x - 0.5, y - 0.5);
    return dctr - rad < 0.4; 
  }

double dsoft(double dx, double dy, double eps)
  { return sqrt(dx*dx + dy*dy + eps*eps); }

dg_tree_t call_make_grid(mdg_dim_t d, mdg_rank_t treeDp, bz_patch_t *b)
  {
    interval_t bbox[MAX_RDIM];
    double ctr[MAX_RDIM];
    int n = b->n;

    auto bool_t split(mdg_cell_index_t k, mdg_rank_t r, bz_patch_t *b, dg_tree_node_t *p);
      /* Tells {dg_make_grid} when to split a cell. Basically, when the variation in
        the coordinates and function values within the cell is larger
        than some tolerance, which depends on the cell's center. */ 

    bool_t split(mdg_cell_index_t k, mdg_rank_t r, bz_patch_t *b, dg_tree_node_t *p)
      { int i; double r2 = 0.0;
        double rad, radMax;
        if (r >= treeDp) return FALSE;
        bz_patch_compute_bbox(b, bbox);
        for (i = 0; i < n; i++) 
          { interval_t bi = bbox[i];
            double ri = (HI(bi) - LO(bi))/2; /* Half-width of interval. */
            ctr[i] = (LO(bi) + HI(bi))/2;
            r2 += ri*ri;
          }
        rad = sqrt(r2);
        radMax = tol(ctr[0], ctr[1]);
        fprintf(stderr, "  split k = %lld  rad = %7.4f radMax = %7.4f\n", k, rad, radMax);
        return rad >= radMax;
      }

    return dg_make_grid(d, mdg_ROOT_CELL, 0, b, NULL, &split);
  }