/* Last edited on 2009-08-24 14:11:52 by stolfi */

void bz_plot_2D_fill_face_shade
  ( PSStream *ps,
    double c00[], 
    double c01[],
    double c10[], 
    double c11[],
    bz_patch_rdim_t n,
    int plotDepth,
    double fMin,
    double fMax
  )
  {
    int i0, i1, j;
    int steps = ipow(2, plotDepth);
    double p[bz_patch_MAX_RDIM], q[bz_patch_MAX_RDIM], r[bz_patch_MAX_RDIM], s[bz_patch_MAX_RDIM];
    for (i0 = 1; i0 <= steps; i0++)
      { for (i1 = 0; i1 <= steps; i1++)
          { double u0 = ((double)i0)/((double)steps);
            double u1 = ((double)i1)/((double)steps);
            double v0 = ((double)i0-1)/((double)steps);
            double v1 = ((double)i1)/((double)steps);
            for (j = 0; j < n; j++)
              { p[j] = dg_bilinear(u0,u1,c00[j],c01[j],c10[j],c11[j]);
                r[j] = dg_bilinear(v0,v1,c00[j],c01[j],c10[j],c11[j]);

              }
            if (i1 > 0) { dg_shade_quadrilateral(ps, p, q, r, s, n, fMin, fMax); }
            for (j = 0; j < n; j++) { q[j] = p[j]; s[j] = r[j]; }
          }
      }
  }

void bz_plot_2D_fill_face_bands
  ( PSStream *ps,
    double c00[], 
    double c01[],
    double c10[], 
    double c11[],
    bz_patch_rdim_t n,
    int plotDepth,
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *R, double *G, double *B
  )
  {
    demand(n == 3, "can't plot bands with more than one function value"); 
    int i0, i1, j;
    int steps = ipow(2, plotDepth);
    double p[bz_patch_MAX_RDIM], q[bz_patch_MAX_RDIM], r[bz_patch_MAX_RDIM], s[bz_patch_MAX_RDIM];
    for (i0 = 1; i0 <= steps; i0++)
      { for (i1 = 0; i1 <= steps; i1++)
          { double u0 = ((double)i0)/((double)steps);
            double u1 = ((double)i1)/((double)steps);
            double v0 = ((double)i0-1)/((double)steps);
            double v1 = ((double)i1)/((double)steps);
            for (j = 0; j < n; j++)
              { p[j] = dg_bilinear(u0,u1,c00[j],c01[j],c10[j],c11[j]);
                r[j] = dg_bilinear(v0,v1,c00[j],c01[j],c10[j],c11[j]);

              }
            if (i1 > 0) 
              { dg_bands_in_quadrilateral
                  (ps, p, q, r, s, n, fStart, fStep, kMin, kMax, R, G, B);
              }
            for (j = 0; j < n; j++) { q[j] = p[j]; s[j] = r[j]; }
          }
      }
  }

void bz_plot_2D_fill_face_isolines
  ( PSStream *ps,
    double c00[], 
    double c01[],
    double c10[], 
    double c11[],
    bz_patch_rdim_t n,
    int plotDepth,
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax        /* Maximum isoline index. */
  )
  {
    demand(n == 3, "can't plot isolines with more than one function value"); 
    int i0, i1, j;
    int steps = ipow(2, plotDepth);
    double p[bz_patch_MAX_RDIM], q[bz_patch_MAX_RDIM], r[bz_patch_MAX_RDIM], s[bz_patch_MAX_RDIM];
    for (i0 = 1; i0 <= steps; i0++)
      { for (i1 = 0; i1 <= steps; i1++)
          { double u0 = ((double)i0)/((double)steps);
            double u1 = ((double)i1)/((double)steps);
            double v0 = ((double)i0-1)/((double)steps);
            double v1 = ((double)i1)/((double)steps);
            for (j = 0; j < n; j++)
              { p[j] = dg_bilinear(u0,u1,c00[j],c01[j],c10[j],c11[j]);
                r[j] = dg_bilinear(v0,v1,c00[j],c01[j],c10[j],c11[j]);

              }
            if (i1 > 0) 
              { dg_isolines_in_quadrilateral(ps, p, q, r, s, n, fStart, fStep, kMin, kMax); }
            for (j = 0; j < n; j++) { q[j] = p[j]; s[j] = r[j]; }
          }
      }
  }