#ifndef dg_plot_2D_func_H
#define dg_plot_2D_func_H

/* Plotting functions on Bézier patches. */
/* Last edited on 2009-05-18 14:06:59 by stolfi */

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

#include <stdio.h>
#include <bool.h>
#include <pswr.h>
#include <dg_tree.h>
#include <bz_basic.h>
#include <bz_patch.h>

#define dg_plot_func_DDIM (2)  
  /* Dimension of grid's domain. */

#define dg_plot_func_MAX_VALUES (3)
  /* Maximum number of function values (besides the shape coords). */

typedef void dg_plot_2D_func_t ( double x[], int nx, double f[], int nf );
  /* Type of a generic function to be plotted. It should take the
    argument vector {x[0..nx-1]} and store the corresponding function
    value into {f[0..nf-1]}. */

void dg_plot_2D_func
  ( PSStream *ps, 
    bz_patch_t *bz,
    bz_patch_rdim_t nf,
    dg_plot_2D_func_t *func,
    dg_tree_t G,
    mdg_rank_t plotDp,
    bool_t isolines
  );
  /* Plots a function {F} over a {DDIM}-dimensional adaptive grid
    {G}, with domain distortion.
    
    The domain of the function {F} is the unit square. The two
    arguments {x[0..DDIM-1]} are interpreted as the coordinates of a point
    in the root cell of {G}, relative to that cell.
    
    The function {F} is defined by the Bézier patch {bz} and/or the
    procedure {func}. More precisely, the value {F(x)} of {F} at point
    {x} is the concatenation of the value {bz(x)} of {bz} at {x}, and the
    value {func(x)} computed by {func} at {x}.
    
    The Bézier patch {bz}, if not NULL, must have domain dimension
    {bz.m == DDIM}, range dimension {bz.n >= DDIM}, and degree
    {bz.g[0..DDIM-1]==(1,..1)}. If the Bézier patch {bz} is NULL, it
    defaults to the natural affine map from the unit square to the
    canonical root cell geometry, with {bz.m == bz.n == DDIM}.
    
    The function {func} is assumed to take a {DDIM}-argument vector
    {x[0..DDIM-1]} and return a result vector {f[0..nf-1]} with {nf}
    elements. If {func} is NULL, then {nf} must be 0.
    
    It follows that the function value {F(x)} has {nt == b.n + nf}
    components. The first {DDIM} components of {F(x)} (which always come
    from the patch {bz(x)}) are interpreted as the plotting coordinates {P(x)}
    of the point {x}. The remaining coordinates (i.e., any components
    of {bz(x)} beyond the first {DDIM}, followed by the {nf} components of
    {func(x)}) are the function values proper, which are plotted using colors
    and/or isolines.
    
    If {isolines} is FALSE, the function values are plotted as various
    shades of color. In this case the function may have from 0 to 3
    values, besides the shape (i.e., {nt} must be in the range {2..5}.
    
    If {isolines} is TRUE, the function values are mapped to discrete
    colors, so that the plot will show color bands separated by 
    isolines. In this case the function must have exacty one value
    besides the shape (i.e., {nt} must be exactly 3).  
    
    The tree {G} is implicitly truncated to depth {plotDp}. */

#endif