/* See figtools.h */
/* Last edited on 2009-01-06 03:26:38 by stolfi */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include <dg_grid.h>
#include <dg_tree.h>
#include <dg_pulse.h>
#include <bz_basic.h>
#include <bz_patch.h>
#include <dg_plot.h>
#include <affirm.h>
#include <pswr.h>

#include <figtools.h>

#define mm (72.27/25.4)
  /* One millimeter in points. */

/* INTERNAL PROTOTYPES */

/* IMPLEMENTATIONS */

PSStream *ps = NULL;
  /* "The" Postscript file for single-document plotting. */

double hsz_old = -1.0, vsz_old = 1.0;
  /* Size of last figure (pt) for single_document plotting. */

void start_document(PlotOptions *o)
  { if (o->eps)
      { ps = pswr_new_stream(txtcat(o->outName, "-"), NULL, TRUE,  NULL,  NULL,         TRUE,  640.0, 480.0); }
    else
      { ps = pswr_new_stream(txtcat(o->outName, "-"), NULL, FALSE, "doc", o->paperSize, FALSE, 0.0, 0.0); }
  }

void start_figure
  ( PlotOptions *o,
    char *figTag,
    char *title,
    interval_t bbox[2],    /* Client plot window. */
    double scale            /* Plot scale (mm per client unit). */
  )
  { 
    fprintf(stderr, "=== %s =========================\n", figTag);
    
    double xsz = HI(bbox[0]) - LO(bbox[0]);
    double ysz = HI(bbox[1]) - LO(bbox[1]);

    /* Figure size in pt: */
    double hsz = scale * xsz * mm;
    double vsz = scale * ysz * mm;

    /* Figure margin width in pt: */
    double hmg = o->margin[0] * mm;
    double vmg = o->margin[1] * mm;

    affirm(ps != NULL, "no output stream");
    if (o->eps)
      { double htot = hsz + 2*hmg;
        double vtot = vsz + 2*vmg;
        pswr_set_canvas_size(ps, htot, vtot);
        pswr_new_canvas(ps, figTag);
        pswr_set_canvas_layout(ps, hsz, vsz, FALSE, hmg, vmg, 0, 1, 1);
        pswr_new_picture
          ( ps,
            LO(bbox[0]), HI(bbox[0]),
            LO(bbox[1]), HI(bbox[1])
          );
      }
    else
      { /* Check whether figure size changed: */
        if ((hsz != hsz_old) || (vsz != vsz_old))
          { /* Start new canvas: */
            pswr_set_canvas_layout(ps, hsz, vsz, FALSE, hmg, vmg, 2, 0, 0);
            hsz_old = hsz; vsz_old = vsz;
          }
        /* Start a new figure in the same canvas: */
        pswr_new_picture
          ( ps,
            LO(bbox[0]), HI(bbox[0]),
            LO(bbox[1]), HI(bbox[1])
          );
        pswr_set_pen(ps, 1,0,0, 0.25, 0,0);
        pswr_add_caption(ps, figTag, 0.0);
        if ((title != NULL) && ((*title) != '\000') && (strcmp(title,figTag) != 0))
          { pswr_add_caption(ps, title, 0.0); }
     }
   
    pswr_set_pen(ps, 0,0,0, 0.15, 0,0);
  }
  
void finish_figure(PlotOptions *o)
  {  }

void finish_document(PlotOptions *o)
  { if (ps != NULL)
      { pswr_close_stream(ps); ps = NULL; }
  }

void hatch_rectangle(PSStream *ps, interval_t B[], double hstep)
  {
    int i;
    double dx = HI(B[0]) - LO(B[0]);
    double dy = HI(B[1]) - LO(B[1]);
    double nx = -dy, ny = dx;
    int nh = (int)ceil(hypot(dx,dy)/hstep);
    fprintf(stderr, 
      "hatch hstep = %.5f d = (%.5f, %.5f) nh = %d\n",
      hstep, dx, dy, nh
    );
    if (nh < 3) { nh = 3; }
    for (i = 1; i < nh; i++)
      { double t = ((double) i)/((double)nh);
        double xm = LO(B[0]) + t*dx;
        double ym = LO(B[1]) + t*dy;
        double lim;
        /* Find {ta} st. {(xm,ym) - ta*(nx,ny)} is on low-right border */
        double ta = 4.0;
        lim = xm - HI(B[0]); if (ta*nx < lim) { ta = lim/nx; }
        lim = ym - LO(B[1]); if (ta*ny > lim) { ta = lim/ny; }
        /* Find {tb} st. {(xm,ym) + ta*(nx,ny)} is on up-left border */
        double tb = 4.0;
        lim = LO(B[0]) - xm; if (tb*nx < lim) { tb = lim/nx; }
        lim = HI(B[1]) - ym; if (tb*ny > lim) { tb = lim/ny; }
        /* Plot segment: */
        pswr_segment(ps, 
            xm - ta*nx, ym - ta*ny,
            xm + tb*nx, ym + tb*ny
          );
      }
  }

/* Matrix for the isometric projection: */
#define XPR0 (-0.7071067811865475)
#define YPR0 (-0.4082482904638630)
#define ZPR0 (+0.5773502691896258)

#define XPR1 (+0.7071067811865475)
#define YPR1 (-0.4082482904638630)
#define ZPR1 (+0.5773502691896258)

#define XPR2 (00.0000000000000000)
#define YPR2 (+0.8164965809277260)
#define ZPR2 (+0.5773502691896258)


void persp_project
  ( double x, double y, double z, 
    double *xp, double *yp, double *zp,
    double vdist
  )
  {
    (*xp) = XPR0*x + XPR1*y + XPR2*z;
    (*yp) = YPR0*x + YPR1*y + YPR2*z;
    (*zp) = ZPR0*x + ZPR1*y + ZPR2*z;
    if (vdist > 0.0) 
      { double mag = vdist/(vdist - (*zp));
        (*xp) *= mag; (*yp) *= mag; (*zp) *= mag;
      }
  }

  
void print_1d_cell(FILE *wr, interval_t cell)
  { fprintf(wr, 
      "[%+8.5f _ %+8.5f]\n",
      LO(cell), HI(cell)
    );
  }
  
void print_2d_cell(FILE *wr, interval_t cell[])
  { fprintf(wr, 
      "[%+8.5f _ %+8.5f]×[%+8.5f _ %+8.5f]\n",
      LO(cell[0]), HI(cell[0]), 
      LO(cell[1]), HI(cell[1])
    );
  }
  
void print_3d_cell(FILE *wr, interval_t cell[])
  { fprintf(wr, 
      "[%+8.5f _ %+8.5f]×[%+8.5f _ %+8.5f]×[%+8.5f _ %+8.5f]",
      LO(cell[0]), HI(cell[0]), 
      LO(cell[1]), HI(cell[1]), 
      LO(cell[2]), HI(cell[2])
    );
  }

bool_t check_pulse_index(dg_pulse_family_t *fam, dg_pulse_mother_index_t pix)
  { if (pix >= fam->nmp)
      { fprintf(stderr, "** mother %d invalid for family ", pix);
        dg_pulse_family_print(stderr, fam);
        fprintf(stderr, "\n");
        return FALSE;
      }
    else 
      { return TRUE; }
  }

interval_t get_pulse_range
  ( dg_pulse_kind_t pkind, /* Kind of mother spline. */
    dg_cont_t c,           /* Continuity class of spline. */
    dg_degree_t g,         /* Degree of spline. */
    dg_pulse_mother_index_t pix,  /* Puse index. */
    dg_grid_size_t gsz      /* Period of spline (#intervals). */
  )
  { interval_t fr = (interval_t){{ +INF, -INF }};
    dg_pulse_family_t fam = dg_pulse_family(pkind, c, g);
    int m = dg_pulse_supp_count(&fam, pix, gsz);
    dg_pulse_t pu = dg_pulse(&fam, pix, gsz, 0);
    double bz[g+1];
    int ix;
    for (ix = 0; ix < m; ix++)
      { dg_pulse_to_bezier(&fam, &pu, ix, bz);
        bz_index_t k;
        for (k = 0; k <= g; k++)
          { double bzk = bz[k];
            if (bzk < LO(fr)) { LO(fr) = bzk; }
            if (bzk > HI(fr)) { HI(fr) = bzk; }
          }
      }
    return fr;
  }

interval_t get_common_pulse_range
  ( dg_pulse_kind_t pkind,  /* Pulse kind. */
    dg_cont_t c,            /* Continuity class of spline. */
    dg_degree_t gmax,       /* Max degree to consider. */
    dg_grid_size_t gsz       /* Num of cells in grid. */
  )
  { interval_t fr = (interval_t){{ +INF, -INF }};
    dg_degree_t g;
    for (g = 0; g <= gmax; g++)
      { dg_pulse_family_t fam = dg_pulse_family(pkind, c, g);
        dg_pulse_mother_index_t pix;
        for (pix = 0; pix < fam.nmp; pix++)
          { interval_t frp = get_pulse_range(pkind, c, g, pix, gsz);
            box_join(2, &fr, &frp, &fr);
          }
      }
    return fr;
  } 

interval_t multiply_ranges(interval_t *x, interval_t *y)
  {
    interval_t z;
    if (LO(*x) >= 0)
      {	if (LO(*y) >= 0)
	  { LO(z) = LO(*x) * LO(*y);
	    HI(z) = HI(*x) * HI(*y);
	  }
	else if (HI(*y) <= 0)
	  { LO(z) = HI(*x) * LO(*y);
	    HI(z) = LO(*x) * HI(*y);
	  }
	else
	  { LO(z) = HI(*x) * LO(*y);
	    HI(z) = HI(*x) * HI(*y);
	  }
      }
    else if (HI(*x) <= 0)
      {	if (LO(*y) >= 0)
	  { LO(z) = LO(*x) * HI(*y);
	    HI(z) = HI(*x) * LO(*y);
	  }
	else if (HI(*y) <= 0)
	  { LO(z) = HI(*x) * HI(*y);
	    HI(z) = LO(*x) * LO(*y);
	  }
	else
	  { LO(z) = LO(*x) * HI(*y);
	    HI(z) = LO(*x) * LO(*y);
	  }
      }
    else
      { if (LO(*y) >= 0)
	  { LO(z) = LO(*x) * HI(*y);
	    HI(z) = HI(*x) * HI(*y);
	  }
	else if (HI(*y) <= 0)
	  { LO(z) = HI(*x) * LO(*y);
	    HI(z) = LO(*x) * LO(*y);
	  }
	else
	  { { double z1 = LO(*x) * HI(*y);
              double z2 = HI(*x) * LO(*y);
              LO(z) = (z1 < z2 ? z1 : z2);
            }
	    { double z1 = LO(*x) * LO(*y);
              double z2 = HI(*x) * HI(*y);
              HI(z) = (z1 > z2 ? z1 : z2);
            }
	  }
      }
    return z;
  }