/* See figs2d.h. */
/* Last edited on 2009-08-21 18:34:17 by stolfilocal */

#include "plot2d.h"

void plot_2d_tree_leaves
  ( PSStream *ps,
    PlotOptions *o, 
    dg_tree_node_t *sub,
    interval_t rootCell[]
  )
  {
    int d = 2;
    
    auto void do_plot(mdg_rank_t r, dg_tree_node_t *n, interval_t cell[]);
    
    void do_plot(mdg_rank_t r, dg_tree_node_t *n, interval_t cell[])
      { if (n == NULL)
          { /* Do nothing */ }
        else if ((LOCH(*n) == NULL) || (HICH(*n) == NULL))
          { /* Leaf cell, plot: */
            affirm((LOCH(*n) == NULL) && (HICH(*n) == NULL), "ch");
            pswr_rectangle(ps, 
              LO(cell[0]), HI(cell[0]), 
              LO(cell[1]), HI(cell[1]), 
              TRUE, TRUE
            );
          }
        else
          { /* Non-leaf cell, recurse: */
            interval_t save = cell[r % d];
            affirm((LOCH(*n) != NULL) && (HICH(*n) != NULL), "ch");
            cell[r % DIM] = interval_split(&save, 0);
            do_plot(r+1, LOCH(*n), cell);
            cell[r % DIM] = interval_split(&save, 1);
            do_plot(r+1, HICH(*n), cell);
            cell[r % DIM] = save;
          }
      }

    double *rgb = &(o->fill.c[0]);
    pswr_set_fill_color(ps, rgb[0], rgb[1], rgb[2]);
    do_plot(0, sub, rootCell);
  }

 void plot_2d_bezier_array
  ( PSStream *ps,
    PlotOptions *o,
    interval_t cell[],
    dg_degree_t g, 
    double mag
  )
  { double cR = (o->color ? 1.000 : 0.400);
    double cG = (o->color ? 0.000 : 0.400);
    double cB = (o->color ? 0.000 : 0.400);
    pswr_set_fill_color(ps, cR,cG,cB);
    int ix, iy;
    for (ix = 0; ix <= g; ix++)
      { double tx = ((double)ix)/((double)g);
        double x = (1.0-tx)*LO(cell[0]) + tx*HI(cell[0]);
        for (iy = 0; iy <= g; iy++) 
          { double ty = ((double)iy)/((double)g);
            double y = (1.0-ty)*LO(cell[1]) + ty*HI(cell[1]);
            pswr_dot(ps, x, y, mag*0.30, TRUE, TRUE);
          }
      }
  }

void plot_2d_bezier_linkages
  ( PSStream *ps,
    PlotOptions *o,
    interval_t cell[], 
    dg_cont_t c,
    dg_degree_t g
  )
  { double frac = 1.0/((double)g);
    double xlo = LO(cell[0]), xhi = HI(cell[0]);
    double ylo = LO(cell[1]), yhi = HI(cell[1]);
    
    /* Vertical linkages: */
    { int ix;
      double dy = frac*(yhi - ylo);
      for (ix = 0; ix <= g; ix++)
        { double tx = ((double)ix)/((double)g);
          double x = (1.0-tx)*xlo + tx*xhi;
          pswr_segment(ps, x, ylo-dy, x, ylo+c*dy);
          pswr_segment(ps, x, yhi-c*dy, x, yhi+dy);
        }
    }
    
    /* Horizontal linkages: */
    { int iy;
      double dx = frac*(xhi - xlo);
      for (iy = 0; iy <= g; iy++) 
        { double ty = ((double)iy)/((double)g);
          double y = (1.0-ty)*ylo + ty*yhi;
          pswr_segment(ps, xlo-dx, y, xlo+c*dx, y);
          pswr_segment(ps, xhi-c*dx, y, xhi+dx, y);
        }
    }
  }

void plot_2d_locus_star
  ( PSStream *ps,
    PlotOptions *o, 
    dg_locus_t E,
    dg_tree_node_t *G,
    interval_t rootCell[],
    double hstep
  )
  {
    mdg_dim_t dE = dg_locus_dimension(DIM, E);
    mdg_cell_index_t kE[4];
    int nkE = ipow(2, DIM - dE);
    dg_locus_star(DIM, E, kE);
    
    interval_t cell[DIM]; /* Cell of {K(E)} that is being painted */
    ix_pos_t kpos;          /* Linearized index of {cell} in {K(E)}. */
            
    for (kpos = 0; kpos < nkE; kpos++)
      { dg_cell_box_root_relative(DIM, kE[kpos], cell);
        box_box_map(DIM, cell, rootCell, cell);
        hatch_rectangle(ps, cell, hstep);
      }
  }

void plot_2d_locus
  ( PSStream *ps,
    PlotOptions *o, 
    dg_locus_t E,
    dg_tree_node_t *G,
    interval_t rootCell[]
  )
  {
    mdg_dim_t dE = dg_locus_dimension(DIM, E);
    interval_t box[DIM]; /* Box of {E}. */
    dg_locus_box_root_relative(DIM, E, box);
    box_box_map(DIM, box, rootCell, box);
    switch (dE)
      { 
        case 0:
          pswr_set_pen(ps,  0.0,0.0,0.0, 0.5,  0.0, 0.0);
          pswr_set_fill_color(ps,  0.0, 0.0, 0.0);
          pswr_dot
            ( ps, 
              LO(box[0]), LO(box[1]),
              1.0,
              TRUE, TRUE      
            );
          break;
        case 1:
          pswr_set_pen(ps,  0.0,0.0,0.0, 0.5,  0.0, 0.0);
          pswr_segment
            ( ps, 
              LO(box[0]), HI(box[0]),
              LO(box[1]), HI(box[1])
            );
          break;
        case 2:
          pswr_set_pen(ps,  0.0,0.0,0.0, 0.15,  0.0, 0.0);
          double *rgb = &(o->fill.c[0]);
          double dim = 0.75;
          pswr_set_fill_color(ps, dim*rgb[0], dim*rgb[1], dim*rgb[2]);
          pswr_rectangle
            ( ps, 
              LO(box[0]), HI(box[0]), 
              LO(box[1]), HI(box[1]), 
              TRUE, TRUE
            );
          break;
        default:
          affirm(FALSE, "bad locus dimension");
      }
  }
  
void plot_2d_tent
  ( PSStream *ps,
    PlotOptions *o, 
    udg_pulse_kind_t pkind[],
    udg_cont_t c[],
    udg_degree_t g[],
    dg_tent_t t,
    dg_tree_node_t *G,
    interval_t rootCell[], 
    mdg_rank_t rank,
    mdg_rank_t plotDepth,
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *cR, double *cG, double *cB
  )
  { mdg_rank_t tent_rank = dg_cell_rank(t.ck/*ell*/);
    
    /* Get the grid sizes {gsz} at rank {tent_rank}: */
    mdg_grid_size_t gsz[DIM];
    dg_grid_size(DIM, tent_rank, gsz);
    
    /* Get maximum pulse index {pix[i]} for each axis: */
    udg_pulse_mother_index_t mix[DIM];
    { int ax; 
      for (ax = 0; ax < DIM; ax++)
       { mix[ax] = dg_pulse_max_index(pkind[ax], c[ax], g[ax]); }
    }
    
    /* Get the pulse index {pix[i]} for each axis: */
    udg_pulse_mother_index_t pix[DIM];
    udg_tent_pulse_indices(DIM, mix, t.tix, pix);
    
    /* Get the anchor cell's position: */
    mdg_grid_pos_t pos[DIM];
    dg_cell_position(DIM, t.ck/*ell*/, pos);
    
    /* Get the tent's support size {psz[i]} for each axis: */
    mdg_grid_size_t psz[DIM];
    dg_tent_mother_supp_size(DIM, pkind, c, g, pix, psz);
    
    auto void t_eval(double x[], double fx[]);
      /* Evaluates the tent function {t} at the point {x[0..DIM-1]}
        (in absolute coordinates), stores the result
        in {fx[0]}. */
    
    auto void do_plot(mdg_rank_t r, dg_tree_node_t *n, mdg_cell_index_t k);
/*    
    fprintf(stderr, "plotting tent ");
    dg_tent_print(stderr, DIM, pkind, c, g, t);
    fprintf(stderr, "\n");
    */
    do_plot(rank, G, 1);

    void t_eval(double x[], double fx[])
      { double z[DIM];
        box_point_unmap(DIM, x, rootCell, z);
        dg_tent_eval(DIM, pkind, c, g, pix, gsz, NULL, pos, z, 0, fx);
      }

    void do_plot(mdg_rank_t r, dg_tree_node_t *n, mdg_cell_index_t k)
      { fprintf(stderr, "{");
        if (n == NULL)
          { /* Omitted cell (?), do nothing */  }
        else if (! dg_pack_intersects_cell(DIM, k, t.ck/*ell*/, psz))
          { /* Skip cell */ }
        else if (r == tent_rank)
          { interval_t CB[DIM];
            dg_cell_print(stderr, DIM, k);
            /* Get cell's absolute box {cell}: */
            dg_cell_box_root_relative(DIM, k, CB);
            box_box_map(DIM, CB, rootCell, CB);
            /* Plot the tent function within {cell}: */
            plot_2d_function(ps, o, 
              t_eval, 
              n, CB, r, 
              plotDepth, 
              fStart, fStep, kMin, kMax, cR, cG, cB
            );
          }
        else if ((LOCH(*n) == NULL) || (HICH(*n) == NULL))
          { /* Leaf cell, do nothing: */
            affirm((LOCH(*n) == NULL) == (HICH(*n) == NULL), "ch");
          }
        else
          { /* Non-leaf cell, recurse: */
            affirm((LOCH(*n) != NULL) && (HICH(*n) != NULL), "ch");
            do_plot(r+1, LOCH(*n), 2*k);
            do_plot(r+1, HICH(*n), 2*k+1);
          }
        fprintf(stderr, "}");
      }
  }
  

void plot_2d_spline
  ( PSStream *ps,
    PlotOptions *o, 
    udg_pulse_family_t fam[],
    dg_cont_t c,
    dg_degree_t g,
    dg_tent_vec_t tv,
    double_vec_t coef,
    dg_tree_node_t *G,
    interval_t rootCell[], 
    mdg_rank_t rank,
    mdg_rank_t plotDepth,
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *cR, double *cG, double *cB
  )
  { 
    int d = 2;
    
    interval_t cell[d]; /* Cell of {K(E)} that is being painted */
            
    /* Should do something smarter... */
    
    auto void spln(double x[], double fx[]);
    /* Evaluates the spline {SUM { coef[i]*tv[i] }} at the point {x[0..DIM-1]}
      (in absolute coordinates), stores the result in {fx[0..DIM-1]}. */
    
    void spln(double x[], double fx[])
      { double z[d], s = 0.0;
        box_point_unmap(d, x, rootCell, z);
        int i;
        int nel = sizeof(tv)/sizeof(dg_tent_vec_t);
        for (i = 0; i < /*tv.*/nel; i++)
          { dg_tent_t t = tv.e[i]; 
            double c = coef.e[i];
            udg_pulse_t p[d];
            dg_tent_unpack(d, fam, &t, p);
            double tx;
            dg_tent_eval(d, fam, p, R, x, ord, f);
            s += c * tx;
          }
        fx[0] = s;
      }

    auto void do_plot(mdg_rank_t r, dg_tree_node_t *n, mdg_cell_index_t k);
    
    void do_plot(mdg_rank_t r, dg_tree_node_t *n, mdg_cell_index_t k)
      { 
        fprintf(stderr, "{");
        if (n == NULL)
          { /* Omitted cell (?), do nothing */  }
        else if ((LOCH(*n) == NULL) || (HICH(*n) == NULL))
          { /* Leaf cell, plot spline: */
            affirm((LOCH(*n) == NULL) == (HICH(*n) == NULL), "ch");
            mdg_cell_print(stderr, DIM, k);
            /* Get cell's absolute box {cell}: */
            mdg_cell_box_root_relative(DIM, k, cell);
            box_box_map(DIM, cell, rootCell, cell);
            /* Plot the tent function within {cell}: */
            plot_2d_function(ps, o, 
              spln, 
              n, cell, r, 
              plotDepth, 
              fStart, fStep, kMin, kMax, cR, cG, cB
            );
          }
        else
          { /* Non-leaf cell, recurse: */
            affirm((LOCH(*n) != NULL) && (HICH(*n) != NULL), "ch");
            do_plot(r+1, LOCH(*n), 2*k);
            do_plot(r+1, HICH(*n), 2*k+1);
          }
        fprintf(stderr, "}");
      }
    fprintf(stderr, "plotting spline\n");
    do_plot(rank, G, 1);
  }

//#warning necessario
void plot_2d_function
  ( PSStream *ps,
    PlotOptions *o, 
    Field func,
    dg_tree_node_t *G,
    interval_t rootCell[], 
    mdg_rank_t rank,
    mdg_rank_t plotDepth,
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *cR, double *cG, double *cB
  )
  {
    auto void do_plot(mdg_rank_t r, dg_tree_node_t *n, interval_t cell[]);
    
    void do_plot(mdg_rank_t r, dg_tree_node_t *n, interval_t cell[])
      { fprintf(stderr, "{");
        if (n == NULL)
          { /* Do nothing */ }
        else if ((LOCH(*n) == NULL) || (HICH(*n) == NULL))
          { /* Leaf cell, plot: */
            affirm((LOCH(*n) == NULL) == (HICH(*n) == NULL), "ch");
            fprintf(stderr, "�");
            pswr_set_pen(ps, 0,0,0,  0.10,  0.0, 0.0);
            plot_2d_function_in_cell(ps, 
              o, cell, r, func, 
              plotDepth, 
              fStart, fStep, kMin, kMax,
              cR, cG, cB
            );
            pswr_set_pen(ps, 0,0,0,  0.15,  0.0, 0.0);
            pswr_rectangle(ps, 
              LO(cell[0]), HI(cell[0]), 
              LO(cell[1]), HI(cell[1]),
              FALSE, TRUE
            );
          }
        else
          { /* Non-leaf cell, recurse: */
            interval_t save = cell[r % DIM];
            affirm((LOCH(*n) != NULL) && (HICH(*n) != NULL), "ch");
            cell[r % DIM] = interval_split(&save, 0);
            do_plot(r+1, LOCH(*n), cell);
            cell[r % DIM] = interval_split(&save, 1);
            do_plot(r+1, HICH(*n), cell);
            cell[r % DIM] = save;
          }
        fprintf(stderr, "}");
      }

    do_plot(rank, G, rootCell);
  }

// void plot_2d_mother_vertex_tent
//   ( PSStream *ps,
//     PlotOptions *o, 
//     udg_pulse_kind_t pkind,
//     dg_cont_t c,
//     dg_degree_t g,
//     dg_tent_mother_index_t tix,    /* Index of tent in set. */
//     mdg_grid_size_t gsz[],    /* X- and Y-domain periods (number of grid intervals). */
//     interval_t fr[],    /* X-, Y-, and function ranges of interest. */
//     mdg_rank_t plotDepth,
//     double fStart,  /* Synchronize isolines with this level. */
//     double fStep,   /* Isoline spacing. */
//     int kMin,       /* Minimum isoline index. */
//     int kMax,       /* Maximum isoline index. */
//     double *cR, double *cG, double *cB
//   )
//   { int M = g - c;
//             
//     auto void t_eval(double x[], double fx[]);
//     /* Evaluates the vertex mother tent function number {tix} at the point {x[0..DIM-1]}
//       (in absolute coordinates), stores the result in {fx[0]}. */
//     
//     void t_eval(double x[], double fx[])
//       { double v = 1.0;
//         int i;
//         for (i = 0; i < DIM; i++) 
//           { udg_pulse_mother_index_t ixi = tix % M; tix /= M;
//             v *= dg_pulse_mother_sum(pkind, c, g, ixi, gsz[i], x[i]);
//           }
//         fx[0] = v;
//       }
//     
//     /* Convert {fr[0..1]} to cell indices (possibly negative): */
//     int lox = (int)(1.0001 * LO(fr[0]));
//     int hix = (int)(1.0001 * HI(fr[0]));
//     int loy = (int)(1.0001 * LO(fr[1]));
//     int hiy = (int)(1.0001 * HI(fr[1]));
//     int ix, iy;
//     
//     for (ix = lox; ix < hix; ix++)
//       for (iy = loy; iy < hiy; iy++)
//         { 
//           interval_t cell[DIM] = 
//             { (interval_t){{ (double)ix, (double)(ix+1) }},
//               (interval_t){{ (double)iy, (double)(iy+1) }}
//             };
//           pswr_set_pen(ps, 0,0,0,  0.10,  0.0, 0.0);
//           plot_2d_function_in_cell(ps, 
//             o, cell, 0, t_eval, 
//             plotDepth, 
//             fStart, fStep, kMin, kMax,
//             cR, cG, cB
//           );
//           pswr_set_pen(ps, 0,0,0,  0.15,  0.0, 0.0);
//           pswr_rectangle(ps, 
//             LO(cell[0]), HI(cell[0]), 
//             LO(cell[1]), HI(cell[1]),
//             FALSE, TRUE
//           );
//         }
//   }
//  
void plot_2d_function_in_cell
  ( PSStream *ps,
    PlotOptions *o, 
    interval_t cell[], /* Cell where to plot. */
    mdg_rank_t rank,       /* Cell rank in grid. */
    Field eval,         /* Function to plot. */
    int plotDepth,      /* Subdivision depth for plotting. */
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *cR, double *cG, double *cB
  )
  {
    auto void do_plot
      ( mdg_rank_t r, 
        double c00[],
        double c01[],
        double c10[],
        double c11[]
      );
    
    void do_plot
      ( mdg_rank_t r, 
        double c00[],
        double c01[],
        double c10[],
        double c11[]
      )
      { if (r >= plotDepth) 
          { /* Small enough, plot: */
            dg_bands_in_quadrilateral
              ( ps,
                c00, c01, c10, c11,
                DIM+1,
                fStart, fStep, kMin, kMax,
                cR, cG, cB
              );
            dg_isolines_in_quadrilateral
              ( ps,
                c00, c01, c10, c11,
                DIM+1,
                fStart, fStep, kMin, kMax
              );
          }
        else
          { /* Still too big, recurse: */
            double m0[DIM+1], m1[DIM+1];
            if (r % DIM == 0)
              { m0[0] = (c00[0] + c10[0])/2;
                m0[1] = (c00[1] + c10[1])/2;
                eval(m0, &(m0[2]));
                m1[0] = (c01[0] + c11[0])/2;
                m1[1] = (c01[1] + c11[1])/2;
                eval(m1, &(m1[2]));
                do_plot(r+1, c00, c01, m0, m1);
                do_plot(r+1, m0, m1, c10, c11);
              }
            else
              { m0[0] = (c00[0] + c01[0])/2;
                m0[1] = (c00[1] + c01[1])/2;
                eval(m0, &(m0[2]));
                m1[0] = (c10[0] + c11[0])/2;
                m1[1] = (c10[1] + c11[1])/2;
                eval(m1, &(m1[2]));
                do_plot(r+1, c00, m0, c10, m1);
                do_plot(r+1, m0, c01, m1, c11);
              }
          }
      }

    /* Compute cell corners, perturbed slightly inwards: */
    double c00[DIM+1], c01[DIM+1], c10[DIM+1], c11[DIM+1];
    { double lo[DIM], hi[DIM];
      int j;
      for (j = 0; j < DIM; j++)
        { lo[j] = LO(cell[j]);
          hi[j] = HI(cell[j]);
          double eps = 0.0001*(hi[j] - lo[j]);
          lo[j] += eps;
          hi[j] -= eps;
        }
      c00[0] = c01[0] = lo[0]; 
      c10[0] = c11[0] = hi[0];
      c00[1] = c10[1] = lo[1]; 
      c01[1] = c11[1] = hi[1];
    }
    /* Evaluate function at corners: */
    eval(c00, &(c00[DIM]));
    eval(c01, &(c01[DIM]));
    eval(c10, &(c10[DIM]));
    eval(c11, &(c11[DIM]));
    do_plot(rank, c00, c01, c10, c11);
  }
  
//#warning necessario
double compute_2d_plot_depth
  ( interval_t rootCell[], /* The root cell (client units). */
    double scale,           /* Plot scale (mm per client unit). */
    double meshSize         /* Nominal plot step (mm). */
  )
  { double dx = HI(rootCell[0]) - LO(rootCell[0]);
    double dy = HI(rootCell[1]) - LO(rootCell[1]);
    double diam = scale*hypot(dx, dy);
    int dp = DIM*(int)ceil(log(diam/meshSize)/log(2.0));
    if (dp < 0) { dp = 0; }
    fprintf(stderr, "compute_2d_plot_depth: dp = %d\n", dp);
    return dp;
  }