(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialization Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)



Off[General::spell1];
<< Graphics`Graphics`
<<Graphics`Colors`
<<DiscreteMath`ComputationalGeometry`

one[lijst_] := Module[
  {i,l},
	l = Length[lijst];
  Table[lijst[[i,1]],{i,l}]
]

two[lijst_] := Module[
  {i},
  Table[lijst[[i,2]],{i,Length[lijst]}]
]

seriestopoint [function_] := Module[ 
		{s, point, i},
			
		s = Series[function,{t,0,11}];
		point = 0;
		For[ i=1, i<= 10,
			sc= SeriesCoefficient[s,i];
			If[sc != 0,
			  point = N[point +sc  10^(-i+1 )];
			];
			i++
		];
		point
]

	SetAttributes[seriestopoint,Listable];

showseries[ lijst_] := Module[
		{points, popltext, i},
		
		points = seriestopoint[lijst];
		popltext = Table[{i,points[[i]]+{.035,.035}},{i,1,Length[points]}];
		Show[Graphics[ {RGBColor[1,0,0],PointSize[0.015],Point/@points,
			RGBColor[0,0,0],
			Apply[Text,popltext,1],
			Line[{{-1.3,0},{1.3,0}}] ,
			Line[ {{0,1.3},{0,-1.3}}]       } ],
			AspectRatio -> Automatic,
			PlotRange -> {-1.3,1.3}]
]

showtrajectories[lijst_] := Module[
		{l, listh, h, pos},
		l = Length[lijst];
		listh = lijst /. t -> h;
		h = .59;
		pos =Table[{i,listh[[i]]+{.03,.03}},{i,l}] ;
		ParametricPlot[Evaluate[lijst], {t,0,1}, Prolog -> {Apply[Text,pos,1]}, 
      Axes -> False, ImageSize -> 500]
	]		

issmallerseries[f1_, f2_] := Module[
		(* Returns True iff f1 < f2 *)
		{ii, pow, c1, c2},
		
		If[ f1 == f2, Return[True] ];
		If[f1 == 0,  
       Return[ f1 < firstcoef[f2] ] ];
		If[ f2 == 0,
       Return[ firstcoef[f1] < f2 ] ] ;
		
		 ii = 1;
		 While[ True,
			  pow = Series[{f1,f2}, {t,0,ii}];
			  c1 = SeriesCoefficient[ pow[[1]],ii];
			  c2 = SeriesCoefficient[ pow[[2]],ii];
			  If[ c1 != c2,
				  Return[ c1 < c2]
			  ];
			  ii++;
			];
			  
	]

sortseriespairs [lijst_]:= Module[
		{modlist, firstst, tt, singlepos, next, i},
		
		modlist = Sort[lijst, issmallerseries[#1[[1]], #2[[1]]]&];
		firsts = one[modlist]; 
		tt = Table[
        Flatten[Position[firsts, firsts[[i]],1 ] ], {i,1,Length[firsts]} ];
		singlepos = {};
		next = {};
		For[ i=1, i <= Length[tt],
			If[tt[[i]] != next,
				next = tt[[i]];
				singlepos = Append[singlepos, next];
			];
		i++;
		];
  	Flatten[
      Table[ Sort[Part[modlist,singlepos[[i]]], 
          issmallerseries[#1[[2]], #2[[2]]]&], {i,1,Length[singlepos]}],1]
	
		]

lo[a_,b_] := Module[
		(* returns true iff a < b, w.r.t. the bigo ordening *)
		{},
		If[ a < 0 &&  b >= 0, Return[True]];
		If[ a==0 && b > 0, Return[True]];
		If[ a < 0 && b < a, Return[True]];
		If[ a > 0 && b  > 0 &&  b < a, Return[True]];
		
		False
		]

findsortedpos[list_, entry_] := Module[
		{notfound, l, i},
		notfound = True;
		l = Length[list];
		i=1;
		While[ notfound && i <= l,
			If[issmallerseries[entry,list[[i]]],
				notfound = False;
				Return[i];
			];
			i++;
		];
		i
		]

containso[ interval_, coef_] := Module[
		{a, b, c},
		(* returns true iff c in {a.b}^o w.r.t. the bigo ordening *)
		
		a = interval[[1]];
		b = interval[[2]];
		c = coef;
		
		(* Print[interval,"  ",coef];	*)
	   
    If[ lo[b,a], Return["Error, empty interval..." ]];
		
		(* first check if a is bigo-smaller than c *)
	  If[  lo[a,c] && lo[c,b],
			Return[True];
		];
		
		False
		]

sign[l_] := Module[
		{i, dd, cc, nn},
		dd = Det[Evaluate[Table[Join[{1},l[[i]]],{i,1,3}]]];
		If[PolynomialQ[dd,t],
		  cc=CoefficientList[dd,t];
			nn = Table[0,{i,Length[cc]}];
			If[ cc == nn,
				Return[0];
			];
	   ];
		dd
		]

ca[cl_] := Module[
        {m,i},
        m = Table[Append[cl[[i]],1],{i,1,3}];
        Det[m]
]
cd[cl_] := Module[
        {m,i},
         m = 
      Table[{Plus @@{cl[[i]][[1]]^2, cl[[i]][[2]]^2}, cl[[i]][[2]], 1},
                {i, 1, 3}];
        -Det[m]
]
ce[cl_] := Module[
        {m,i},
         m = 
      Table[{Plus @@{cl[[i]][[1]]^2, cl[[i]][[2]]^2}, cl[[i]][[1]], 1},
                {i, 1, 3}];
        Det[m]
]
cf[cl_] := Module[
        {m,i},
         m = Table[{Plus @@{cl[[i]][[1]]^2, cl[[i]][[2]]^2}, cl[[i]][[1]],
                 cl[[i]][[2]]}, {i, 1, 3}];
        -Det[m]
]


hasnontrivcenter[lijst_] := Module[
		(* Returns <True> i.a.o.i. the circle defined by <lijst> 
        has its center outside the origin *)
		{a,d,e},
		a = ca[lijst];
		d = cd[lijst];
		If[ Abs[Limit[d/a, t -> 0]] > 0,
			Return[True];
		];
		e = ce[lijst];
		If[Abs[Limit[e/a, t -> 0]] > 0,
			Return[True];
		];
		Return[False]
	]

firstcoef[function_] := Module[
		{depth=0,coef=0,s},
	  
		If[function == 0,
		  Return[0]
		];
		
		While[ coef == 0,
			depth++;
			s = Series[function, {t,0,depth}];
			coef = SeriesCoefficient[s,depth];
		];
		
		coef
	]

	SetAttributes[firstcoef,Listable];

firstcoefpoly[function_] := Module[
		(* Finds the coefficient belonging to the term of lowest
			  degre *)
		{depth=0,coef=0},
	  
		If[function == 0,
		  Return[0]
		];
		
		While[ coef == 0,
			depth++;
			coef= Coefficient[function,t,depth];
		];
		
		coef
	]

	SetAttributes[firstcoefpoly,Listable];

bigoplsign[function_] := Module[
		{depth=0,coef=0,s},
	  
		If[function == 0,
		  Return[0]
		];
		
		While[ coef == 0,
			depth++;
			s = Series[function, {t,0,depth}];
			coef = SeriesCoefficient[s,depth];
		];
		
		Sign[coef] depth
	]

SetAttributes[bigoplsign,Listable];

maxtree[lijst_] := Module[
		(* looks for the biggest series *)
	   
		{cand, zeros,  depth, pow,coef,max,possies},
		
		If[ Length[lijst] == 1,
			Return[lijst];
		];
		
		cand = lijst;
		
		(* First ckeck for zeros *)
	  zeros = Position[cand,0];
		cand =Delete[cand,zeros];
	
		depth = 0;
		While[ Length[cand] > 1 ,
		  depth++;
			pow =  Series[cand,{t,0,depth}];
			coef = Map[SeriesCoefficient[#1,depth]&,pow];
			max = Max[coef];
			possies = Flatten[Position[coef, max]];
			cand = Part[cand, possies];
	  ];
		
		cand = cand[[1]];
			
		If[ Length[zeros] > 0,
			If[ firstcoef [cand] < 0,
			  Return[0];
			];
		];
			
		cand
					
]


mintree[lijst_] := Module[
		(* looks for the smallest series *)
	   
		{cand, zeros, depth, pow,coef,min,	possies},
		
		If[ Length[lijst] == 1,
			Return[lijst];
		];
		
		cand = lijst;
		
		(* First ckeck for zeros *)
	  zeros = Position[cand,0];
		cand =Delete[cand,zeros];
		
		depth = 0;
		While[ Length[cand] > 1 ,
		  depth++;
			pow = Series[cand,{t,0,depth}];
			coef = Map[SeriesCoefficient[#1,depth]&,pow];
			min = Min[coef];
		  possies = Flatten[Position[coef, min]];
		 cand = Part[cand, possies];
	  ];
		
		cand = cand[[1]];
			
		If[ Length[zeros] > 0,
			If[ firstcoef [cand] > 0,
			  Return[0];
			];
		];
			
		cand
			
]

serieshull1D[lijst_]:= Module[
		{},
		{mintree[Union[lijst]],maxtree[Union[lijst]]}
]


	hullcand[lijst_] := Module[

  {l, i,li, extr,	oextr, cand, done},
  
	(* compute maxima in both x- and y-direction *)
   l = Length[lijst];
	 If[ l <= 3,
			Return[Range[l]];
		];
		
   extr = {serieshull1D[one[lijst]],serieshull1D[two[lijst]]};
   oextr = bigoplsign[extr];
	If[ Not[Sign[oextr] == {{-1,1},{-1,1}}],
				(* Print["One sided input..."]; *)
				Return[Range[l]];
	];
   (* Print[extr]; *)
   Print["O extremes: ", oextr]; 
   
   cand ={};
   For[ i=1, i <= l,

     done = False;
     li = lijst[[i]];

     (* 2. if not, check for O(max) or O(min). *)

     (* a. w.r.t x-coord. *)
     
      If[ Not[done ]&&
         
          Not[containso[oextr[[1]], bigoplsign[li[[1]] ] ] ],
         cand = Append[cand, i];
         done = True;
     ];

     (* b. w.r.t. y-coord. *)
     
      If[ Not[done]&&
         
          Not[ containso[oextr[[2]], bigoplsign[li[[2]] ] ] ],
         cand = Append[cand,i];
     ];

     i++
   ];

cand
]


positions[l1_,l2_] := Module[
		{i},
		Flatten[Table[Position[l1,l2[[i]] ],{i,1,Length[l2]}]]
	]

rc[lijst_] := Module[
		(* computes the rc of the normal of the line define by lijst[[1]] and \
lijst[[2]] *)
		{below, top, phi},
		below = lijst[[2,2]]-lijst[[1,2]];
		top = lijst[[1,1]]-lijst[[2,1]];
		If[ below == 0, 
			If[ top == 0, 
				 Return["0/0 expression encountered..."]
			];
			 Return[ Sign[firstcoef[top]] Infinity];
		];
   	phi = Limit[(top/below),t->0];
		 If[phi == ComplexInfinity,
			  phi = Infinity
			];
		phi
	]

anglerc[lijst_] := Module[
		{},
		ArcTan[rc[lijst]]
	]

upperhullseries[modlist_] := Module[ 
		(* Computes the upper hull of modlist, 
      that should be a sorted list of trajectories *)
	{lupper, n, i},
	lupper = {modlist[[1]],modlist[[2]]};
  n = Length[modlist];
  For[ i=3, i <= n,
    lupper = Append[ lupper, modlist[[i]] ];
    While[ 
        Length[lupper] > 2 &&
      
          Sign[firstcoef[sign[ Take[lupper, -3] ] ] ]> 0,
        lupper = Delete[lupper, Length[lupper] - 1];
        (* Print["upper hull: ", positions[modlist,lupper]];*)
    ];
    i++
  ];
lupper
]
		

lowerhullseries[modlist_] := Module[
		(* Computes the lower hull of modlist, 
      that should be a sorted list of trajectories *)
		{llower, i, n},
		n = Length[modlist];
		llower = { modlist[[n]],modlist[[n-1]]};
    For[i=n-2, i>=1,
      (* Print[i]; *)
      llower = Append[ llower, modlist[[i]] ];
      While[ 
        Length[llower] > 2 &&
        
          Sign[firstcoef[sign[Take[llower, -3] ] ] ] > 0,
          llower = Delete[llower,Length[llower] - 1];
          (* Print["lower hull: ",positions[modlist, llower]]; *)
      ];
      i--;
   ];
	llower
]

filteredges[list_, edgel_, hull_] := Module[
		(* leaves only those sites in edgel that are global w.r.t the edge vertices \
*)
		{l,filtlist, i,j},
		l = Length[hull];
		filtlist =  Table[{}, {i,l}];
		For[i=1, i<=l,
			For[j=1, j <= Length[edgel[[i]]],
				If[ hasnontrivcenter[ 
            list[[Append[hull[[{i,Mod[i+1,l] /. 0 -> l}]],
                  edgel[[i,j]] ] ]] ],
					filtlist[[i]] = Append[ filtlist[[i]], edgel[[i,j]] ]
					];
				j++;
				];
			i++;
			];
		filtlist
		]

zonetestpoly[a_,b_, lijst_] := Module[
		{perpab,bhat,ahat,inzone,i,l1,l2,s1,s2},
		perpab = {a[[2]]-b[[2]],b[[1]]-a[[1]]};
		bhat = a + perpab;
		ahat = b + perpab;
		inzone = {};
		For[i=1, i <= Length[lijst],
			l1 = {bhat,a,lijst[[i]]};
			l2 = {ahat,b, lijst[[i]]};
			s1 = Sign[firstcoefpoly[sign[l1]]];
			s2 = Sign[firstcoefpoly[sign[l2]]];
			If[s1 s2< 0,
				inzone = Append[inzone, i];
			];
			i++;
		];
		
		inzone
		]

pointlinedistance[u_,v_,w_] := Module[
		(* Let l be the line defined by the points u and v.
					The module returns d(w,l),where w is a third point *)
		{s,top},
		s= (v[[1]]-u[[1]])^2 + (v[[2]]-u[[2]])^2;
		top = (v[[2]]-u[[2]])(w[[1]]-u[[1]])-
				(v[[1]]-u[[1]])(w[[2]]-u[[2]]);
		(* Take abs value *)
		If[firstcoef[top]== -1,
			top = - top;
			];
		 top/Sqrt[s]
		]

permsort[list_,p_] := Module[
		(* sorts the elements in <list> according to the function <p>
				 returns as a second argument the ordered list and additionaly, 
      as the first argument, the inverse permutation *)
		{l,extlist,sortlist,
      perm},
		l = Length[list];
		extlist = Table[{i,list[[i]]},{i,l}];
		sortlist = Sort[extlist,p[#1[[2]],#2[[2]]]&];
		perm = Table[sortlist[[i,1]],{i,l}];
		sortlist = Table[sortlist[[i,2]],{i,l}];
		{perm,sortlist}
	]

replacebyseries[list_] := Module[
		(* replaces the functions in list> 
        by their power series expansions in such a way that all functions \
remain different and that differences in x- and y-
          coordinates are expressed in the power series expansions as well  *)
\
	
		{l, h1, lh1, i, j, k, seriesone, h2, lh2, seriestwo},
		
		l = Length[list];
		If[ l > Length[Union[list]],
			Return["Error, all sites should be different... "];
		];
		h1 = one[list];
		lh1 = Length[Union[h1]];
		i=1;
		seriesone := Series[h1,{t,0,i}];
		While[Length[Union[Normal[seriesone]]] < lh1,
	     i++;
		];
		seriesone = Normal[seriesone];

    h2 = two[list];
    lh2 = Length[Union[h2]];
     j=1;
		seriestwo := Series[h2,{t,0,j}];
     While[Length[Union[Normal[seriestwo]]] < lh2,
        j++;
     ];
		seriestwo = Normal[seriestwo];
		
	Table[{seriesone[[k]],seriestwo[[k]]},{k,1,l}]

]

bisectorangle[b_,e_,ip_] := Module[
		(* Returns 0<=alpha<2 Pi such that alpha represents the
			  directions of the outward bisector of b and e.  *)
	{alpha, u, signzero, 
      signu},
	alpha = ArcTan[rc[{b,e}]];
	u = {Cos[alpha],Sin[alpha]};
	signzero = Sign[firstcoef[sign[{b,e, ip}]]];
	signu = Sign[firstcoef[sign[{b,e,u}]]];
	If[ signzero == signu,
			alpha = alpha + Pi;
	];
	alpha
	]

insidepoint[abc_] := Module[
		(* Returns the vector 1/2 a+1/4 + 1/4 c = 
        A + 1/2(b-a) + 
          1/4(c-b) and this is a vector that is contained in the interior of \
the triangle abc  *)
	{},
	1/4(2 abc[[1]]+abc[[2]]+abc[[3]])
]

voronoilimittwo[lijst_,pos_, globlijst_] := Module[
		{worklist, phi, angles},
		worklist = lijst[[pos]];
		phi = bisectorangle[worklist[[1]],worklist[[2]],{0,0}];
		angles= {phi, phi + Pi};
		showlimithull[pos, angles, globlijst]
	]

voronoilimitmult[lijst_,pos_, globlijst_] := Module[
  {ip,angles, worklist, replijst, lua,l, i, b, e,  
   dirhull, uniqangles, newgloblijst, j,aa, ab},

	replijst= replacebyseries[lijst];
	(* Easy case: |pos| = 2. *)
	l = Length[pos];
	If[ l ==2 ,
			Return[ voronoilimittwo[ replijst, pos, globlijst] ];
	];
	
(* Compute all bisector angles. For this, compute some interior point <ip> *)

    worklist = replijst[[pos]];

If[ sign[worklist[[ {1,2,3} ]] ] == 0,
   Return["## Not yet implemented..."];
];
ip = Expand[insidepoint[worklist[[{1,2,3}]]  ]];

angles = {};	
For[i=1,  i <= l,
  b = worklist[[i]];
  e = worklist[[ Mod[i+1,l] /. 0 -> l ]];
  angles = Append[angles,bisectorangle[b,e, ip]];
  i++;
];

(* If all angles are different, we can start drawing *)

    lua = Length[Union[angles]];
If[ l== lua,
			Return[ showlimithull[ pos, angles, globlijst] ];
];
		
(* If not all angles are different, 
      the direction hull differs from the limit hull and we have to adapt 
  i. dirhull ii. angles iii. globlijst *)	

  dirhull = {};
  uniqangles = {};
  newgloblijst = Table[{},{i,1,lua}];
  j=0;
  For[ i=1, i <= l,
    aa = angles[[Mod[i-1,l] /. 0 -> l]];
    ab = angles[[i]];
    If[ aa != ab, 
      j++;
      dirhull = Append[dirhull, pos[[i]] ];
      uniqangles = Append[uniqangles,  angles[[i]] ];
      newgloblijst[[j]] = 
          Join[globlijst[[i]],newgloblijst[[j /. 0 -> lua]] ] ;
    ];
    (* else *)
    If[ aa == ab,
       newgloblijst[[j /. 0 -> lua ]] = 
            Join[newgloblijst[[j  /. 0 -> lua]], {pos[[i]]},
       globlijst[[i]] ];
    ] ;
    i++;
  ];

 showlimithull[dirhull, uniqangles, newgloblijst]
]

showlimithull[dirhull_, angles_,newgloblijst_] := Module[
		{i, l, epos, o, lines, fnumbers, lnumbers, mults},
		
		Print["direction hull ", dirhull];
		Print[angles];
		Print[newgloblijst];
		
		l = Length[dirhull];
	
		For[ i=1, i<= l,
       Print["edge: ", dirhull [[i]], "-",dirhull [[Mod[i+1,l] /. 0 -> l]],
         " bisector angle: ", angles[[i]], 
         " globals: ",newgloblijst[[i]] ];
      i++;
    ];
		
		epos = Append[dirhull ,dirhull [[1]]];
		o= {0,0};
		lines= Table[Line[{o,{Cos[angles[[i]]],Sin[angles[[i]] ] } }],{i,1,l}];
		fnumbers = Table[ Text[epos[[i]],{Cos[angles[[i]]+Pi/24], Sin[angles[[i]] 
         +Pi/24] }],{i,1,l}];
		lnumbers = Table[ Text[epos[[i+1]],{ Cos[angles[[i]] - Pi/24], 
           Sin[angles[[i]] - Pi/24] }],{i,1,l}];
		mults = Table[Text[1 +Length[newgloblijst[[i]]],1.1{  Cos[angles[[i]]], 
             Sin[angles[[i]]]}],{i,1,l}];
		
		Show[Graphics[   {Red,lines ,Black, fnumbers , lnumbers ,Blue, mults}], 
      AspectRatio -> Automatic]
	]

direction[point_] := Module[
  {x,y,phi},
  x = point[[1]];
  y = point[[2]];
  If[ x == 0,
    If[ y== 0,
      Return["Error, 0/0 expression encountered..."];
    ];
    Return[Infinity];
  ];
  phi = Limit[ y/x, t -> 0];
  If[ phi == ComplexInfinity,
    phi = Infinity;
  ];
  phi
]

squarenorm[lijst_] := Module[
  {},
  lijst[[1]]^2 + lijst[[2]]^2
]

handleglobaledge[b_,e_,cand_] := Module[
  {bn={0,0}, en, candn, l, phi, mat, dist, ps, permutation,
   globalsites, c, i, xc2, allxcs, j, q, r},

  If[ cand == {},
    Return[ {} ]
  ];

  (* translate the scene s.t. <b>=(0,0) *)
  l = Length[cand];
  en = e-b;
  candn = Table[cand[[i]] - b, {i,l}];

  (* compute the direction of the line segment be *)
 
  phi = direction[en];
  phi = ArcTan[phi];

  (* rotate the scene around <0> over an angle <phi>, where
     -<phi> is the angle between the positive x-axis and
     the line segment <be> *)
  
    mat = - { { Cos[phi], -Sin[phi] }, { Sin[phi], Cos[phi] } };
  en = mat.en;
  candn=Table[mat.candn[[i]],{i,l}];

  (* Get the distances of the candidates to the line <be>. *)
  
    dist =  Table[pointlinedistance[{0,0},en,candn[[i]] ], {i,1,l}];

  (* order candidates according to increasing distance *)
  
    ps = permsort[ dist, issmallerseries ];
  (* ps[[1]] is the order of the candidates and ps[[2]] are
     the distances it self *)
  permutation = ps[[1]];
  Print["Permutation ", permutation];
  dist=ps[[2]];
 (* initialise <global sites> *)
  (* 
      We use the square of the distance to 0! *)
  
    globalsites = { {{0,0},0}, {en, squarenorm[en]} };

  (* We can add the first candidate immediately! But this is
     not completely trivial *)

  c := candn[[ permutation[[i]] ]];
  q := globalsites[[j-1,1]];
  r := globalsites[[j,1]];


  (* Now add the candidates one by one *)
  For[i=1, i <= l,
    Print["i ",i];
    (* get xc^2 *)
    xc2 = squarenorm[c] - dist[[i]]^2;
    (* only continue if xcs doesn't occur yet *)
    
      allxcs = two[globalsites];
    If[ MemberQ[allxcs, xc2] == False,
      (* find proper interval *)
      j = findsortedpos[ allxcs, xc2];
      (* test angle with <q> *)
      q = globalsites[[j-1,1]];
      Print["<qc ", anglerc[{q,c}] ];
      If[ Abs[anglerc[{q,c}]] == Pi/2,
        (* test angle with <r> *)
        r = globalsites[[j,1]];
          Print["<rc ", anglerc[{r,c}]];
        If[ Abs[anglerc[{r,c}]] == Pi/2,
          (* global! so add <c> to <globalsites> at position <j> *)
          
            globalsites = Insert[ globalsites, {c, xc2}, j];
        ];
      ];
    ];
    i++;
  ];

  (* Now we want to have some more nice output *)
  
    newglobals = one[globalsites];
  (* we now the first and the last *)
  newglobals = Delete[newglobals,-1];
  newglobals = Delete[newglobals,1];
  pos = positions[candn,newglobals]

]
		

global3[list_] := Module[
		(* THIS IS THE PROGRAM THAT COMPUTES ALL GLOBAL SITES *)
  {lijst, modlist, 
      sortcand, unsortedcand, lijstcand, lhull, lupper, n, i, llower, jj, 
      alledges},
  
  (* We want all sites to be different *)
  
    If[ Length[list] > Length[Union[list]],
    Return["Error, all sites should be different... "];
  ];

  lijst = replacebyseries[list];
  
  modlist = lijst[[hullcand[lijst]]];
  unsortedcand = positions[lijst,modlist];
  lijstcand = lijst[[unsortedcand]];
  (* Sort the candidate sites *)
  modlist = sortseriespairs[modlist];
  sortcand =  positions[lijst,modlist];
  Print["sorted candidates: ", sortcand];

  (* Compute the upper hull *)
  lupper = upperhullseries[modlist]; 
  Print["\n","final upper hull ", positions[lijst,lupper]];

  (* Compute the lower hull *)
  llower = lowerhullseries[modlist];
  Print["final lower hull ", positions[lijst,llower]];

  (* join the upper and the lower hull *)
  
    jj = Join[lupper, Delete[llower,{{1},{Length[llower]}}]];
  lhull = Flatten[Table[Position[lijst,jj[[i]] ], {i,1,Length[jj]} ]];
  hl = Length[lhull];

  (* Now the zone test *)

  alledges = {};
  For[ i=1, i <= hl,
    b = lijst[[lhull[[i]]]];
    e = lijst[[lhull[[Mod[i+1,hl] /.  0 -> hl  ]]  ]];
    z = zonetestpoly[b, e, lijstcand];
    z = unsortedcand[[z]];
    alledges = Append[alledges, z];
    i++;
  ];
  
  Print["alledges ",alledges];

  (* Now we filter *)
  Print["\nfiltering..."];
  alledges = filteredges[lijst, alledges, lhull];
  Print["Candidates for the global hull: ", alledges];
  
  (* We apply the direction test on the left candidates *)

  For[ i=1, i <= hl,
    If[ Length[ alledges[[i]] ] > 1,
      b = lijst[[lhull[[i]]]];
      e = lijst[[lhull[[Mod[i+1,hl] /.  0 -> hl  ]]  ]];
      alledges[[i]] = 
            alledges[[i]][[handleglobaledge[b,e,lijst[[alledges[[i]] ]] ] ]];
    ];
    i++;
  ];

  Print["new global sites ", alledges];
 
  {lhull, alledges}	

]

Print["Finished!"];



serieshull[list_] := Module[
		(* THIS IS THE PROGRAM THAT COMPUTES ALL GLOBAL SITES *)
  {length, labels, 
      sortlist,  lhull, lupper, n, hl,i, llower, jj, alledges, b, e, z},
  
  (* We want all sites to be different *)
    length=Length[list];
      If[ length > Length[Union[list]],
          Return["Error, all sites should be different... "];
        ];
    labels=Range[length];
    
      sortlist = sortseriespairs[list];
    
      (* Compute the upper hull *)
      lupper = upperhullseries[sortlist]; 
      Print["\n","final upper hull ", positions[list,lupper]];
    
      (* Compute the lower hull *)
      llower = lowerhullseries[sortlist];
      Print["final lower hull ", positions[list,llower]];
    
      (* join the upper and the lower hull *)
      
    jj = Join[lupper, Delete[llower,{{1},{Length[llower]}}]];
      lhull = Flatten[Table[Position[list,jj[[i]] ], {i,1,Length[jj]} ]];
      hl = Length[lhull];
    
      (* Now the zone test *)
    
      alledges = {};
      For[ i=1, i <= hl,
          b = list[[lhull[[i]]]];
          e = list[[lhull[[Mod[i+1,hl] /.  0 -> hl  ]]  ]];
          z = zonetestpoly[b, e, list];
          alledges = Append[alledges, z];
          i++;
        ];
      
      Print["alledges ",alledges];
    
      (* Now we filter *)
      Print["\nfiltering..."];
      alledges = filteredges[list, alledges, lhull];
      Print["Candidates for the global hull: ", alledges];
      
      (* We apply the direction test on the left candidates *)
    
      For[ i=1, i <= hl,
          If[ Length[ alledges[[i]] ] > 1,
              b = list[[lhull[[i]]]];
              e = list[[lhull[[Mod[i+1,hl] /.  0 -> hl  ]]  ]];
              
        alledges[[i]] = 
          alledges[[i]][[handleglobaledge[b,e,list[[alledges[[i]] ]] ] ]];
            ];
          i++;
        ];
    
      Print["new global sites ", alledges];
     
      {lhull, alledges}	
    
    ]

allatonce[lijst_] := Module[
		{},
		gl =serieshull[lijst];
		voronoilimitmult[ lijst, gl[[1]],gl[[2]]]
	]