/*******************************************************************************
+
+  LEDA 3.5
+
+  _spring.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
#include <LEDA/graph_alg.h>

#include <math.h>
#include <LEDA/array2.h>


#define FREPULSE(d) ((k2 > d) ? kk/d : 0)

static float log_2(int x)
{ float l = 0;
  while (x)
  { l++;
    x >>= 1;
   }
  return l/2;
}


void SPRING_EMBEDDING(const graph& G, node_array<double>& xpos, 
                                      node_array<double>& ypos,
                                      double xleft, double xright, 
                                      double ybottom, double ytop,
                                      int iterations)
{ list<node> L;
  SPRING_EMBEDDING(G,L,xpos,ypos,xleft,xright,ybottom,ytop,iterations); }


void SPRING_EMBEDDING(const graph& G, const list<node>& fixed_nodes,
                                      node_array<double>& xpos, 
                                      node_array<double>& ypos,
                                      double xleft, double xright, 
                                      double ybottom, double ytop,
                                      int iterations)
{

  if (xleft >= xright || ybottom >= ytop)
      error_handler(1,"SPRING_EMBDDING: illegal bounds.");

  node_array<list_item>	lit(G); 
  node_array<bool>  fixed(G,false);

  node u,v;
  edge e;

  forall(v,fixed_nodes) fixed[v] = true;

  int c_f = 1;

  double width  = xright - xleft;
  double height = ytop - ybottom; 

  double tx_null = width/50;
  double ty_null = height/50;
  double tx = tx_null;
  double ty = ty_null;

  double k = sqrt(width*height / G.number_of_nodes()) / 2;
  double k2 = 2*k;
  double kk = k*k;

  int ki = int(k);

  if (ki == 0) ki = 1;

  //build  matrix of node lists

  int xA = int(width / ki + 1);
  int yA = int(height / ki + 1);

  array2<list<node> > A(-1,xA,-1,yA);

  forall_nodes(v,G)
  { int i = int((xpos[v] - xleft)   / ki);
    int j = int((ypos[v] - ybottom) / ki);
    if (i >= xA || i < 0) error_handler(1,"spring: node out of range");
    if (j >= yA || j < 0) error_handler(1,"spring: node out of range");
    lit[v] = A(i,j).push(v);
  }


  while (c_f < iterations)
  {
    node_array<double>	xdisp(G,0); 
    node_array<double>	ydisp(G,0); 

   // repulsive forces

   forall_nodes(v,G) 
   { int i = int((xpos[v] - xleft)   / ki);
     int j = int((ypos[v] - ybottom) / ki);

     double xv = xpos[v];
     double yv = ypos[v];

     for(int m = -1; m <= 1; m++)
      for(int n = -1; n <= 1; n++)
       forall(u,A(i+m,j+n))
       { if(u == v) continue;
         double xdist = xv - xpos[u];
         double ydist = yv - ypos[u];
         double dist = sqrt(xdist * xdist + ydist * ydist);
         if (dist < 1e-3) dist = 1e-3;
         xdisp[v] += FREPULSE(dist) * xdist / dist;
         ydisp[v] += FREPULSE(dist) * ydist / dist;
        }
      xdisp[v] *=  (double(rand_int(750,1250))/1000.0);
      ydisp[v] *=  (double(rand_int(750,1250))/1000.0);
     }

 

   // attractive forces

   forall_edges(e,G)
   { node u = G.source(e);
     node v = G.target(e);
     double xdist=xpos[v]-xpos[u];
     double ydist=ypos[v]-ypos[u];
     double dist=sqrt(xdist*xdist+ydist*ydist);

     float f = (G.degree(u)+G.degree(v))/6.0;

     dist /= f;

     xdisp[v]-=xdist*dist/k;
     ydisp[v]-=ydist*dist/k;
     xdisp[u]+=xdist*dist/k;
     ydisp[u]+=ydist*dist/k;
    }


   // preventions

   forall_nodes(v,G)
   { 
     if (fixed[v]) continue;

     int i0 = int((xpos[v] - xleft)/ki);
     int j0 = int((ypos[v] - ybottom)/ki);

     double xd= xdisp[v];
     double yd= ydisp[v];
     double dist = sqrt(xd*xd+yd*yd);

     if (dist < 1) dist = 1;

     xd = tx*xd/dist;
     yd = ty*yd/dist;

     double xp = xpos[v] + xd;
     double yp = ypos[v] + yd;

     int i = i0;
     int j = j0;
     
     if (xp > xleft && xp < xright) 
     { xpos[v] = xp;
       i = int((xp - xleft)/ki);
      }

     if (yp > ybottom && yp < ytop) 
     { ypos[v] = yp;
       j = int((yp - ybottom)/ki);
      }
 
     if (i != i0 || j != j0)
     { if (lit[v] == nil) error_handler(1,"delete nil item");
       A(i0,j0).del_item(lit[v]);
       lit[v] = A(i,j).push(v);
      }
 
    }

   tx = tx_null / log_2(c_f); 
   ty = ty_null / log_2(c_f);

   c_f++;
  }

}
