#include "defines.h"
#include <math.h>
/* Routines to construct the voronoi diagram */

extern struct VEC2 V2_sum();
extern struct VEC2 V2_sub();
extern struct VEC2 V2_times();
extern double V2_cprod();
extern struct VEC2 V2_cross();
extern double V2_dot();
extern double V2_magn();


circle_center(a,b,c,f)
    /* computes the center of the circle passing through A, B & C. */
    /* fills in the area and vc fields of f */
struct VEC2 a,b,c;
TRI_PTR f;
{
    struct VEC2 ac,bc;
    register double ac2,bc2,det;

    ac = V2_sub(a,c); ac2=V2_dotq(ac,ac);
    bc = V2_sub(b,c); bc2=V2_dotq(bc,bc);
    Area(f) = V2_cprodq(ac,bc);

    ac = V2_times(bc2,ac);
    bc = V2_times(ac2,bc);

    Corner(f) =  V2_sum(c,V2_times(-.5/Area(f),V2_cross(V2_sub(ac,bc))));

}

#define cntr(e) circle_center(origv(e), destv(e), destv(onext(e)),left(e))

find_vor_corners()
{
    EDGE_PTR e;
    SITE_PTR v;
    struct VEC2 center,edgev,inf;
    double vol;
    double min_len,tot_int_len,sol_vol,temp_len;

    for_tris(e)
	new_corner(sym(e));
	cntr(e);
    end_tris(e)

    /* Do outside edges to infinity */
    for_rring_a(One_bndry_edge,e)
	fflush(stdout);
	center = leftv(e); /* should be cntr(e) from above */
	edgev = V2_sub(destv(e), origv(e));
	info(e).del_length = V2_magnq(edgev);
	info(e).vor_length = -1.;
	info(e).coupling = 0.;
	inf = V2_sum(center,V2_times(1/info(e).del_length,V2_cross(edgev)));
	/* gives a length 1 vector in the right direction towards infinity */
	new_inf_corner(e);
	rightv(e) = inf;
    end_rring_a(One_bndry_edge,e);

    for_sites(v)
	Vol(v) = 0.0;
    end_sites(v)

    min_len = HUGE;
    tot_int_len = 0;
    for_int_edges(e)
	edgev = V2_sub(destv(e), origv(e));
	if ((info(e).del_length = V2_magnq(edgev)) < min_len)
	    min_len=info(e).del_length;
	edgev = V2_sub(leftv(e), rightv(e));
	temp_len = V2_magnq(edgev);
	if (!n_wulff)
	    info(e).vor_length = temp_len;
	else
	{
	    register double temp;
	    register int w;
	    info(e).vor_length = 0;
	    for(w=0;w<n_wulff;w++)
		if ((temp=V2_dotq(edgev,wulff[w])) > info(e).vor_length)
		{
		    info(e).vor_length = temp;
		    info(e).which_wulff = w;
		}
	}
	if (Type(orig(e)) != Type(dest(e)))
	    tot_int_len += info(e).vor_length;
	info(e).coupling = temp_len/info(e).del_length *
		(HtCond(dest(e))+HtCond(orig(e)))/2.0;
	vol = info(e).del_length * temp_len / 4.0;
	Vol(orig(e)) += vol; Vol(dest(e)) += vol;
    end_int_edges(e)
    /* Now volume is right for all sites, except bndry sites, whose
    volume should be infinite.  We have the finite 'piece' of it. */

    printf("Minimum site separation %f\n",min_len);
    printf("Total interface energy (length) %f\n",tot_int_len);
    sol_vol = 0;
    for_sites(v)
	if (Type(v))
	    sol_vol += Vol(v);
    end_sites(v)
    printf("Total volume of solid %f\n",sol_vol);
}