/*
 * Copyright (c) 1992, The Geometry Center
 *                     University of Minnesota 
 *                     1300 South Second Street
 *                     Minneapolis, MN  55454
 *
 * email address: software@geom.umn.edu
 *
 * This software is copyrighted as noted above.  It may be freely copied,
 * modified, and redistributed, provided that the copyright notice is
 * preserved on all copies.
 *
 * There is no warranty or other guarantee of fitness for this software,
 * it is provided solely "as is".  Bug reports or fixes may be sent
 * to the authors, who may or may not act on them as they desire.
 *
 * You may not include this software in a program or other software product
 * without supplying the source, or without informing the end-user that the
 * source is available for no extra charge.
 *
 * If you modify this software, you should include a notice giving the
 * name of the person performing the modification, the date of modification,
 * and the reason for such modification.
 *
 *     The National Science and Technology Research Center for
 *      Computation and Visualization of Geometric Structures
 */

#include "link.h"
#include <stdio.h>
#include <math.h>

#define MAX_CXS          100     /* max # canonical rotations */
#define	EPSILON		.2
#define BASE_IT		25	/* iterate dynamics base_it + it_rate*nx**2 (nx = # of xings) */
#define IT_RATE		.1
#define PI	3.141592653589


/*
#define DEBUG
*/

#define DISTANCE(x1, x2)	\
	sqrt ((x1->x - x2->x)*(x1->x - x2->x) + (x1->y - x2->y)*(x1->y - x2->y))

    double angle_table[][4] = {
	PI, 0, PI/2, -PI/2,	/* orientation = RIGHT */
	PI, 0, -PI/2, PI/2	/* orientation = LEFT */
	};

    extern int tofromtable[][2];

    short orient_table[4][4] = {
	UNORIENT, UNORIENT, RIGHT, LEFT,
	UNORIENT, UNORIENT, LEFT, RIGHT,
	LEFT, RIGHT, UNORIENT, UNORIENT,
	RIGHT, LEFT, UNORIENT, UNORIENT
	};

void getXY(crossing *xing, int which, int num, double *x, double *y, region *bigreg)
{
  double x1, y1, x2, y2, p, q, hey;
  crossing *xing2;
  int dir;
  if (xing->edgeRegions[which][num] != bigreg)
    {
      *x = xing->edgeRegions[which][num]->x;
      *y = xing->edgeRegions[which][num]->y;
      return;
    }
  
  dir = (which == UP) ? N0 : N1;
  xing2 = xing->nhbr[dir];
  x1 = xing->x; y1 = xing->y;
  x2 = xing2->x; y2 = xing2->y;
  p = (x1 + x2)/2;
  q = (y1 + y2)/2;
  hey = sqrt(p*p + q*q);
  *x = (p*3)/hey; *y = (q*3)/hey;
}

void linkdyn(link *linkptr, short tangle)

{
    crossing *xingptr;
    region *rptr, *find_canonical_region(), *reg;
    double min, min2, angle, angle0, tstore[MAXX][6], rstore[MAXX][2];
    int i,j,k,m,n,p, edgeout, edgein, offset,offsets[MAX_CXS], which_offset=0,numx;
    short redraw=0;
    double angles[4], r, dx, dy, sum, angle_N0, x, y;
    
    if (!tangle) rptr = find_canonical_region(linkptr);

/* DRAW_NEW:*/

    /* initialize coordinates of crossing to origin */
    for (xingptr = linkptr->xinglist, i=0; i<linkptr->nxings; ++i, xingptr = NEXTXING(xingptr))
      {
	xingptr->x = xingptr->y = xingptr->phase = 0.0;
	xingptr->fixed = 0;
	xingptr->edgeX[0] = xingptr->edgeY[0]
	  = xingptr->edgeX[1] = xingptr->edgeY[1] = 0;
      }
    /* Ditto for regions */
    if (!tangle)
      {
	for (i=0, reg = linkptr->rlist[0]; i < linkptr->nregions;
	     i++, reg = linkptr->rlist[i])
	  reg->x = reg->y = 0.0;
	/* construct base polygon on unit circle */
	/* determine which direction (clockwise or cclockwise) to position the vertices */
	edgeout = rptr->dlist[0];
	edgein = rptr->dlist[rptr->nv-1];
	edgein = tofromtable[rptr->vlist[rptr->nv-1]->ud[edgein]][edgein & 1];
	if (orient_table[edgeout][edgein] !=  rptr->vlist[0]->orient)
	  for (xingptr = linkptr->xinglist, i=0; i<linkptr->nxings; ++i, xingptr = NEXTXING(xingptr))
	    xingptr->orient = (xingptr->orient == RIGHT) ? LEFT : RIGHT;
	
	/* position the crossings of this region on the unit circle */
	angle0 = 2*PI/rptr->nv;
	if (!redraw) offset=find_canonical_xing(rptr,offsets,&numx);
	for (p=0; p<rptr->nv; ++p)
	  {
	    i = (p+offset) % rptr->nv;
	    angle = angle0*p;
	    rptr->vlist[i]->x = cos(angle);
	    rptr->vlist[i]->y = sin(angle);
	    if (rptr->dlist[i] == N0)
	      rptr->vlist[i]->phase = angle + (PI/4.0);
	    else if (rptr->dlist[i] == P0)	
	      rptr->vlist[i]->phase = angle + 5*(PI/4.0);
	    else if (rptr->dlist[i] == N1)  
	      rptr->vlist[i]->phase = angle + 
		(rptr->vlist[i]->orient == RIGHT  ? 3*(PI/4.0) : -PI/4.0);
	    else if (rptr->dlist[i] == P1)  
	      rptr->vlist[i]->phase = angle + 
		(rptr->vlist[i]->orient == RIGHT  ? -(PI/4.0) : 3*PI/4.0);
	    else printf("linkdyn: Error in region reconstruction\n");
	    rptr->vlist[i]->fixed = 1;
	  }
      }
    else      /* pin down four corners (if tangle) */
      {
	double angles2[4][2];

	/* initialize angles of nw, ne, sw, and se */
	angles2[0][0] = 3*(PI/4.0); angles2[0][1]= -(PI/4.0);    /* angles2[i][0] gives angle of unit circle */
	angles2[1][0] = PI/4.0;     angles2[1][1] = -3*(PI/4.0); /* to place end.  angles2[i][1] gives phase  */
	angles2[2][0] = -3*(PI/4.0);angles2[2][1] =  PI/4.0;     /* assuming strand is directed inward      */
	angles2[3][0] = -(PI/4.0);  angles2[3][1] = 3*(PI/4.0);

	for (xingptr = linkptr->xinglist, i=0; i<4; xingptr = NEXTXING(xingptr), i++)
	  {
	    xingptr->x = cos(angles2[i][0]);
	    xingptr->y = sin(angles2[i][0]);
	    xingptr->phase = (angles2[i][1] +
			      ((xingptr->orient == TERMINUS) ? PI/2.0 : 0));
	    xingptr->fixed = 1;
	  }
      }

    /* iterate a fixed # of times, proportional to # of crossings */
    for (min = 0.0, n = 0; n < BASE_IT + IT_RATE*linkptr->nxings*linkptr->nxings; ++n)
      {
	int denominator[MAXX][6];
	int rdenom[MAXX][2];
	crossing *xptr, *prevXing;
	int weight;
	int which, dir;

	/* Update all the crossings */
	for (xingptr = linkptr->xinglist, i=0; i<linkptr->nxings; ++i, xingptr = NEXTXING(xingptr))
	  {
	    denominator[i][0] = denominator[i][1] = denominator[i][2]
	      = denominator[i][3] = denominator[i][4] = denominator[i][5] = 0;
	    tstore[i][0] = tstore[i][1] = tstore[i][2] = tstore[i][3]
	      = tstore[i][4] = tstore[i][5] = 0;

	    /* First, deal with the crossing's edge neighbors */

	    tstore[i][0] += xingptr->edgeX[UP];
	    tstore[i][1] += xingptr->edgeY[UP];

	    tstore[i][0] += xingptr->edgeX[DOWN];
	    tstore[i][1] += xingptr->edgeY[DOWN];

	    xptr = xingptr->nhbr[P0];
	    which = (xptr->nhbr[N0] == xingptr) ? UP : DOWN;
	    tstore[i][0] += xptr->edgeX[which];
	    tstore[i][1] += xptr->edgeY[which];

	    xptr = xingptr->nhbr[P1];
	    which = (xptr->nhbr[N0] == xingptr) ? UP : DOWN;
	    tstore[i][0] += xptr->edgeX[which];
	    tstore[i][1] += xptr->edgeY[which];

	    denominator[i][0] = denominator[i][1] = 4;

	    /* And the edge's crossing neighbors */

	    tstore[i][2] += xingptr->x;   /* Up edge */
	    tstore[i][3] += xingptr->y;
	    tstore[i][2] += xingptr->nhbr[N0]->x;
	    tstore[i][3] += xingptr->nhbr[N0]->y;
	    tstore[i][4] += xingptr->x;   /* Down edge */
	    tstore[i][5] += xingptr->y;
	    tstore[i][4] += xingptr->nhbr[N1]->x;
	    tstore[i][5] += xingptr->nhbr[N1]->y;
	    denominator[i][2] += 2;
	    denominator[i][3] += 2;
	    denominator[i][4] += 2;
	    denominator[i][5] += 2;
	    /* Next, deal with the crossing's region neighbors */
	    if (!tangle)
	      {
		for (j=0; j<4; j++)
		  {
		    if (rptr != xingptr->regions[j])
		      {
			denominator[i][0] += 1;
			denominator[i][1] += 1;
			
			tstore[i][0] += xingptr->regions[j]->x;
			tstore[i][1] += xingptr->regions[j]->y;
		      }
		  }
	      }

	    /* And the edge's region neighbors */
	    if (!tangle)
	      {
		getXY(xingptr, UP, 0, &x, &y, rptr);    /* Up edge */
		tstore[i][2] += x;        
		tstore[i][3] += y;
		getXY(xingptr, UP, 1, &x, &y, rptr);	    
		tstore[i][2] += x;
		tstore[i][3] += y;
		getXY(xingptr, DOWN, 0, &x, &y, rptr);  /* Down edge */
		tstore[i][4] += x;
		tstore[i][5] += y;
		getXY(xingptr, DOWN, 1, &x, &y, rptr);
		tstore[i][4] += x;
		tstore[i][5] += y;
		denominator[i][2] += 2;
		denominator[i][3] += 2;
		denominator[i][4] += 2;
		denominator[i][5] += 2;
	      }
	  }
	if (!tangle)
	  {
	    /* Update all the regions */
	    for (i=0, reg = linkptr->rlist[0]; i < linkptr->nregions;
		 i++, reg = linkptr->rlist[i])
	      {
		if (reg != rptr)
		  {
		    rdenom[i][0] = rdenom[i][1] = 0;
		    rstore[i][0] = rstore[i][1] = 0;
		    
		    /* Deal with crossing neighbors */
		    
		    for (j=0, xingptr = reg->vlist[0]; j < reg->nv;
			 j++, xingptr = reg->vlist[j])
		      {
			if (tangle && xingptr->n < 4) weight = 5; else weight = 1;
			
			rdenom[i][0] += weight;
			rdenom[i][1] += weight;
			
			rstore[i][0] += (xingptr->x) * weight;
			rstore[i][1] += (xingptr->y) * weight;
		      }
		    
		    /* Deal with edge neighbors */
		    /* (this part sucks rotten tomatoes) */
		    
		    for (prevXing = NULL, j=0,
			 xingptr = reg->vlist[0], dir = reg->dlist[0];
			 j < reg->nv; j++, xingptr = reg->vlist[j], dir = reg->dlist[j])
		      {
			if (prevXing)
			  {
			    /* Previous direction was a wrong-way, so take care
			       of previous edge now */
			    which = (xingptr->nhbr[N0] == prevXing) ? UP : DOWN;
			    rstore[i][0] += xingptr->edgeX[which];
			    rstore[i][1] += xingptr->edgeY[which];
			    rdenom[i][0]++; rdenom[i][1]++;
			  }
			/* Try to take care of current edge, or set flag saying
			   we'll do it later */
			switch(dir)
			  {
			  case N0:
			    rstore[i][0] += xingptr->edgeX[UP];
			    rstore[i][1] += xingptr->edgeY[UP];
			    rdenom[i][0]++; rdenom[i][1]++;
			    prevXing = NULL;
			    break;
			  case N1:
			    rstore[i][0] += xingptr->edgeX[DOWN];
			    rstore[i][1] += xingptr->edgeY[DOWN];
			    rdenom[i][0]++; rdenom[i][1]++;
			    prevXing = NULL;
			    break;
			  default:
			    prevXing = xingptr;
			    break;
			  }
		      }
		    if (prevXing)            /* Must do final edge from first point */
		      {
			xingptr = reg->vlist[0];
			which = (xingptr->nhbr[N0] == prevXing) ? UP : DOWN;
			rstore[0][0] += xingptr->edgeX[which];
			rstore[0][1] += xingptr->edgeY[which];
			rdenom[0][0]++; rdenom[0][1]++;
		      }
		  }
	      }
	  }
	/* record the results */
	for (xingptr = linkptr->xinglist, i=0; i<linkptr->nxings; ++i, xingptr = NEXTXING(xingptr))
	  {
	    if ( !(xingptr->fixed))	
	      {
		xingptr->x = (denominator[i][0] ? tstore[i][0]/denominator[i][0] : 0);
		xingptr->y = (denominator[i][1] ? tstore[i][1]/denominator[i][1] : 0);
	      }
	    xingptr->edgeX[UP] = (denominator[i][2]?tstore[i][2]/denominator[i][2]:0);
	    xingptr->edgeY[UP] = (denominator[i][3]?tstore[i][3]/denominator[i][3]:0);
	    xingptr->edgeX[DOWN]=(denominator[i][4]?tstore[i][4]/denominator[i][4]:0);
	    xingptr->edgeY[DOWN]=(denominator[i][5]?tstore[i][5]/denominator[i][5]:0);
	  }
	if (!tangle)
	  {
	    for (i=0, reg = linkptr->rlist[0]; i < linkptr->nregions;
		 i++, reg = linkptr->rlist[i])
	      {
		if (reg != rptr)
		  {
		    reg->x = (rdenom[i][0] ? rstore[i][0] / rdenom[i][0] : 0);
		    reg->y = (rdenom[i][1] ? rstore[i][1] / rdenom[i][1] : 0);
		  }
	      }
	  }
      }

        /* split apart crossings close in y */
        for (n = 0; n < BASE_IT + IT_RATE*linkptr->nxings*linkptr->nxings; ++n)
	  {
	    float *bigger, *smaller;
	    for (xingptr = linkptr->xinglist, i=0; i<linkptr->nxings; ++i, xingptr = NEXTXING(xingptr))
	      for (j=0; j<4; j++)
		if (xingptr->nhbr[j] != xingptr)
		  {
		    if (xingptr->nhbr[j]->y > xingptr->y)
		      {
			bigger = &(xingptr->nhbr[j]->y);
			smaller = &(xingptr->y);
		      }
		    else
		      {
			bigger = &(xingptr->y);
			smaller = &(xingptr->nhbr[j]->y);
		      }
		    if (*bigger-*smaller < 0.05)
		      *bigger += 0.07;
		  }
	  }
		

	/* compute minimum distance between adjacent nodes */
	for (min = 1.0, xingptr = linkptr->xinglist, i=0; 
	    i<linkptr->nxings; ++i, xingptr = NEXTXING(xingptr))
	      {
	      for (min2 = 1.0, j=0; j<4; ++j)
		{
		if (xingptr != xingptr->nhbr[j])
		    {
		    r = DISTANCE(xingptr, xingptr->nhbr[j]);
		    if (r == 0)
			{
			crossing *xx[2];
			double lx[2], ly[2], dx, dy;
			for (m = 0, k=0; k<4; ++k)	
			    if (xingptr->nhbr[j] != xingptr->nhbr[k])	
				{
				lx[m] = xingptr->nhbr[k]->x;	
				ly[m] = xingptr->nhbr[k]->y;	
				xx[m] = xingptr->nhbr[k];
				m++;
			 	if (m>1) break;
				}
			dx = lx[1] - lx[0];
			dy = ly[1] - ly[0];
			/* horrible kludge to fix special case problem */
			if ( !(linkptr->nxings %2))
			    {
			    xingptr->x = lx[0] + dx/2 + dy/4;
			    xingptr->y = ly[0] + dy/2 - dx/4;
			    xingptr->nhbr[j]->x = lx[0] + dx/2 - dy/4;
			    xingptr->nhbr[j]->y = ly[0] + dy/2 + dx/4;
			    }
	  		else
			    {
			    xingptr->x = lx[0] + dx/2 - dy/4;
			    xingptr->y = ly[0] + dy/2 + dx/4;
			    xingptr->nhbr[j]->x = lx[0] + dx/2 + dy/4;
			    xingptr->nhbr[j]->y = ly[0] + dy/2 - dx/4;
			    }
			r = DISTANCE(xingptr, xingptr->nhbr[j]);
			}
			
		    if (r < min2)  min2 = r;
		    }
		}
	      xingptr->sep = min2;
#ifdef DEBUG
	fprintf(stderr,"separation of crossing %d is %f f\n", xingptr->n, xingptr->sep);
#endif
	      if (min2 < min)	min = min2;
	    }

#ifdef DEBUG
	fprintf(stderr,"min distance on iteration %d is %f\n", n, min);
#endif
	linkptr->min_sep = min;
	
	/* now compute phases of the nodes */
	for (xingptr = linkptr->xinglist, i=0; i<linkptr->nxings; ++i, xingptr = NEXTXING(xingptr))
	  {
	  if ( !(xingptr->fixed))
	    {
	    /* compute angles for each edge */
	    for (j=0; j<4; ++j)
		{
		dx = xingptr->nhbr[j]->x - xingptr->x;
		dy = xingptr->nhbr[j]->y - xingptr->y;
		angles[j] = atan2(dy, dx);
		angles[j] -= angle_table[xingptr->orient][j];
		}
	    angle_N0 = angles[N0];
	    for (sum = 0, j=0; j<4; ++j) 	angles[j] -= angle_N0;
	    for (sum = 0, j=0; j<4; ++j)
		{
		while (angles[j] < PI)	angles[j] += 2*PI;
		while (angles[j] > PI)	angles[j] -= 2*PI;
		sum += angles[j];
		}
	    xingptr->phase = (sum)/4.0 + angle_N0;
	    }
#ifdef DEBUG
	printf("phase, crossing %d: %f\n",xingptr->n, xingptr->phase);
#endif
	}
    redraw = 0;
    
    /* draw all rots for this particular canreg */
/*    if (draw_all_rots && which_offset < numx)
      {
      linkdraw(linkptr);
	offset = offsets[++which_offset];
	redraw=1;
	goto DRAW_NEW;
      }*/
    
    /* do next canreg */
    which_offset = 0;
/*    if (draw_all_regs && which_canreg < ncan)
      {
	linkdraw(linkptr);
	rptr = canregs[++which_canreg];
	goto DRAW_NEW;
      }*/
    ncan = 0;
  }






    

/*
Unsuccessful attempt to let points on unit circle move on unit circle
		if (xingptr->fixed)
		    {
		    d = sqrt(tstore[i][0]*tstore[i][0] + tstore[i][1]*tstore[i][1]);
		    if (d > 0)	scale = 4.0/d;
		    else scale = 1.0;
		    xingptr->x = scale * tstore[i][0]/4;
		    xingptr->y = scale * tstore[i][1]/4;
		    }
		else
		    {
		    xingptr->x = tstore[i][0]/4;
		    xingptr->y = tstore[i][1]/4;
		    }
*/
