/* Last edited on 2002-04-27 16:45:06 by stolfi */
/* Test of libdygrid routines. */

#include <stdio.h>
#include <math.h>
#include <dygrid.h>
#include <isolines.h>
#include <js.h>
#include <ps.h>

/*
  A coordinate axis. */
  
void dyg_begin_figure(
  FILE *f,
  dyg_dimensions* dim
);
void dyg_end_figure(FILE *f);
void dyg_plot_2d_tree(
  FILE *f, dyg_header *h, 
  dyg_dimensions* dim,
  axis long_axis,
  int deep_max
);

void dyg_plot_2d_subtree(
  FILE *f, dyg_node *p, 
  dyg_dimensions* dim,
  axis long_axis,
  int deep,
  int deep_max
);

void dyg_draw_point_region(
  FILE* f, dyg_header *h, 
  dyg_dimensions* dim,
  int deep, double x, double y
);

double length(double dx, double dy, double dz);
bool toobig(dyg_dimensions* dim);

int dummy_data_adjust(dyg_node* node);

int main(int argc, char **argv)
{
  dyg_header *h; 
  dyg_node_list* list = 0;

  FILE *f = open_write("gridplot.eps");
//  iso_triangle_list *triangles, *tritemp;
  dyg_dimensions dim;//, d;
  iso_coords p;

  axis long_axis = x_axis;
  int deep_max = argc > 1 ? atoi(argv[1]) : 14;
//  double i;
  int param;

  double Step = argc > 4 ? atof(argv[4]) : 0.0;

  p.axis[X] = argc > 2 ? atof(argv[2]) : 0.0;
  p.axis[Y] = argc > 3 ? atof(argv[3]) : 0.0;
  p.axis[Z] = 0;


  dim.min[X] = 0.0;
  dim.max[X] = 1.0;
  dim.min[Y] = 0.0;
  dim.max[Y] = sqrt(2.0)/2;
  dim.min[Z] = 0.0;
  dim.max[Z] = 0.0;


  fprintf(stderr,"\n*\n*\n*\n");
  
  h = dyg_new_tree(2, &dim, x_axis);

  dyg_shatter_node(h->root, deep_max, toobig, dummy_data_adjust, &param, &dim, long_axis, false);
  dyg_begin_figure(f, &dim);
  /*
  list = dyg_find_adjacents_cells(h->root, &dim, h->inicial_axis, 1, deep_max, p.axis);
  if( list == NULL )
    fprintf(stderr, "\nList == NULL ");
  for(aux = list; aux != NULL; aux = aux->next)
    {
      dyg_get_cell_dimension(aux->node->index, &dim, h->inicial_axis, &d);      

      triangles = iso_split_cell_in_triangles(h->root, &dim);
      tritemp = 0;
      i = 0.1;
      for( ; triangles != NULL; triangles = triangles->next)
	{
	  if( tritemp )
	    {
	      free(tritemp->triangle);
	      free(tritemp);
	    }
	       
	  ps_fill_triangle(f,triangles->triangle->P.axis[X], triangles->triangle->P.axis[Y], triangles->triangle->Q.axis[X], 
                           triangles->triangle->Q.axis[Y], triangles->triangle->R.axis[X], triangles->triangle->R.axis[Y], 
			   i * i, 1 - i, i );	  
	  
 	  ps_draw_segment(f, triangles->triangle->P.axis[X], triangles->triangle->P.axis[Y], triangles->triangle->Q.axis[X], 
			  triangles->triangle->Q.axis[Y]);
	  ps_draw_segment(f, triangles->triangle->R.axis[X], triangles->triangle->R.axis[Y], triangles->triangle->Q.axis[X], 
			  triangles->triangle->Q.axis[Y]);
	  ps_draw_segment(f, triangles->triangle->P.axis[X], triangles->triangle->P.axis[Y], triangles->triangle->R.axis[X], 
			  triangles->triangle->R.axis[Y]);
	  i +=  0.1;

	  tritemp = triangles;
	}

      if( tritemp )
	{
	  free(tritemp->triangle);
	  free(tritemp);
	}
      //       ps_fill_and_draw_rectangle(f, d.min[X], d.max[X], d.min[Y], d.max[Y], 255,0,0);
    }

  dyg_flush_list(list);
  */

  ps_fill_dot(f, p.axis[X], p.axis[Y], 0.5, 0,0,0);

  fprintf(stderr, "\n\nIso Start" );
  iso_isolines( f, h, 0, Step, 0, true);
  dyg_plot_2d_tree(f, h, &dim, long_axis, MAX_INT);


  dyg_end_figure(f);

  dyg_free_tree(h);
  free(h);
  fprintf(stderr, "\n\nOK so far...\n");
  return(0);
}

double length(double dx, double dy, double dz)
{
  return sqrt(dx*dx + dy*dy + dz*dz);
}

bool toobig(dyg_dimensions* dim)
/*
  An arbitrary splitting criterion. */
{
  double xctr = (dim->min[X] + dim->max[X])/2.0;
  double yctr = (dim->min[Y] + dim->max[Y])/2.0;
  double r = length((dim->max[X] - dim->min[X])/2.0, (dim->max[Y] - dim->min[Y])/2.0, 0.0);
  double d1 = length((xctr - 0.2)/9.0, (yctr - 0.3)/9.0, 0.003);
  double d2 = length((xctr - 0.6)/5.0, (yctr - 0.5)/5.0, 0.003);
  double d3 = length((yctr - (0.9 + 0.4*sin(2.0*xctr)))/4.0, 0.0, 0.003);
  double maxr = 3.0/(1.0/d1 + 1.0/d2 + 1.0/d3);
  return (r > maxr);
}

void dyg_begin_figure(
  FILE *f,
  dyg_dimensions* dim
)
/*
  Prepares an Encapsulated Postscript file for drawing a dyadic grid. */
{
  double dx = dim->max[X] - dim->min[X];
  double ex = 0.1*dx;
  double toth = (dx + 2*ex)*450.0;
  
  double dy = dim->max[Y] - dim->min[Y];
  double ey = 0.1*dy;
  double totv = (dy + 2*ey)*450.0;
  
  ps_begin_figure(f, 
    dim->min[X]-ex, dim->max[X]+ex, dim->min[Y]-ey, dim->max[Y]+ey, 
    0.0, toth, 0.0, totv,
    5,7
  );
}

void dyg_end_figure(FILE *f)
/*
  Terminates an Encapsulated Postscript drawing of a dyadic grid. */
{
  ps_end_figure(f);
}
  
void dyg_plot_2d_tree(
  FILE *f, dyg_header *h, 
  dyg_dimensions* dim,
  axis long_axis,
  int deep_max
)
{
  affirm(h->d == 2, "wrong tree dimension");
  ps_set_pen(f, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0);
  ps_draw_rectangle(f, dim->min[X], dim->max[X], dim->min[Y], dim->max[Y]);
  dyg_plot_2d_subtree(f, h->root, dim, long_axis, 1, deep_max);
}
  
void dyg_plot_2d_subtree(
  FILE *f, dyg_node *p, 
  dyg_dimensions* dim,
  axis long_axis,
  int deep,
  int deep_max
)
{
  dyg_dimensions d;

  d.min[X] = dim->min[X];
  d.max[X] = dim->max[X];
  d.min[Y] = dim->min[Y];
  d.max[Y] = dim->max[Y];

  if ((p != NULL) && ((p->c0 != NULL) || (p->c1 != NULL)) && deep < deep_max)
    { 
      double xmid = (dim->max[X] + dim->min[X])/2.0;
      double ymid = (dim->max[Y] + dim->min[Y])/2.0;

      if( p->virtual )
	ps_set_pen(f, 0.0, 0.0, 1.0, 0.1, 0.0, 0.0);
      else
	ps_set_pen(f, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0);

      ps_fill_dot(f, xmid, ymid, 0.5, 0, 0, 0);


      if (long_axis == x_axis)
      { 
          ps_draw_segment(f, xmid, dim->min[Y], xmid, dim->max[Y]);
	  d.max[X] = xmid;
          dyg_plot_2d_subtree(f, p->c0, &d, y_axis, deep + 1, deep_max);

	  d.max[X] = dim->max[X];
	  d.min[X] = xmid;
          dyg_plot_2d_subtree(f, p->c1, &d, y_axis, deep + 1, deep_max);
      }
      else
      { 
          ps_draw_segment(f, dim->min[X], ymid, dim->max[X], ymid);
	  d.max[Y] = ymid;
          dyg_plot_2d_subtree(f, p->c0, &d, x_axis, deep + 1, deep_max);

	  d.max[Y] = dim->max[Y];
	  d.min[Y] = ymid;

          dyg_plot_2d_subtree(f, p->c1, &d, x_axis, deep + 1, deep_max);
      }
    }
}

void dyg_draw_point_region(
  FILE* f, dyg_header *h, 
  dyg_dimensions* dim,
  int deep, double x, double y
)
{
  axis naxis;
  double aux;
  int d = 1;
 
  dyg_node* p =  h->root;
  naxis = h->inicial_axis;
  while( d < deep && p->c0)
  {
    if( naxis == x_axis )
      {
       aux = (dim->max[X] + dim->min[X])/2.0;
       if( x >= aux )
	 {
           dim->min[X] = aux;
           p = p->c1;
	 }
       else if( x < aux )
	 {
            dim->max[X] = aux;
	    p = p->c0;
	 }
      }
    else
      {
       aux = (dim->max[Y] + dim->min[Y])/2.0;
       if( y >= aux )
	 {
           dim->min[Y] = aux;
	   p = p->c1;
	 }
       else if( y < aux )
	 {
	   dim->max[Y] = aux;        
	   p = p->c0;
	 }
      }
    d++;
    naxis = next_axis(naxis);
  }

  ps_set_pen(f, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0);
  ps_fill_and_draw_rectangle(f, dim->min[X], dim->max[X], dim->min[Y], dim->max[Y], 0,0,255);
}  

int dummy_data_adjust(dyg_node* node)
{
  int ret = 0;

  if( node->data == NULL )
    node->data = malloc(sizeof(double));
  else
  {
    ret = sizeof(*((double*)node->data)) / sizeof(double);
    node->data = realloc(node->data, (ret + 1) * sizeof(double));
  }
  
  ret = (int)(((double*)(node->data))[ret] = node->index);
  return ret;
}