
/* We consider SL_2(C) adjoined orientation-reversing 
   maps.  The latter are distinguished by det=-1
   and represent f(z) = (aw+b)/(cw+d) where w=z-bar.
   We compute the action on circles.

   Lines a represented by circles (c,r) with r < 0 and
   with c = closest point to the origin, and with
   r = angle in revolutions [-1,0].
*/

#include "cx.h"
#include "circles.h"
#include "lim.h"
#include <math.h>

#define INF 1e20

complex det(x)
	matrix x;
{
	return(sub(mult(x.a,x.d),mult(x.b,x.c)));
}

matrix make_sl2(x)
	matrix x;
{
	complex d;

	d = recip(cx_sqrt(det(x)));
	x.a = mult(x.a,d);
	x.b = mult(x.b,d);
	x.c = mult(x.c,d);
	x.d = mult(x.d,d);
	return(x);
}

matrix conjmat(x)
	matrix x;
{
	x.a = cx_conj(x.a);
	x.b = cx_conj(x.b);
	x.c = cx_conj(x.c);
	x.d = cx_conj(x.d);
	return(x);
}

matrix inverse(x)
	matrix x;
{
	matrix z;
	complex d;
	d = det(x);
	z.a = x.d;
	z.b.x = -x.b.x;
	z.b.y = -x.b.y;
	z.c.x = -x.c.x;
	z.c.y = -x.c.y;
	z.d = x.a;
	if(d.x > 0) return(z);
	z.a.x = -z.a.x;
	z.b.x = -z.b.x;
	z.c.x = -z.c.x;
	z.d.x = -z.d.x;
	return(z);
}

matrix product(x,y)
	matrix x,y;
{
	matrix z;
	complex d;

	d = det(x);
	if(d.x < 0) y = conjmat(y);
	z.a = add(mult(x.a,y.a),mult(x.b,y.c));
	z.b = add(mult(x.a,y.b),mult(x.b,y.d));
	z.c = add(mult(x.c,y.a),mult(x.d,y.c));
	z.d = add(mult(x.c,y.b),mult(x.d,y.d));
	return(z);
}

matrix circle_to_matrix(circ)
	circle circ;
{
	matrix z;
	if(circ.r <= 0) return(line_to_matrix(circ));
	z.a.x = circ.c.x/circ.r; 
	z.a.y = circ.c.y/circ.r;
	z.b.x = circ.r-(circ.c.x*circ.c.x+circ.c.y*circ.c.y)/circ.r;
	z.b.y = 0;
	z.c.x = 1/circ.r;
	z.c.y = 0;
	z.d.x = -circ.c.x/circ.r;
	z.d.y =  circ.c.y/circ.r;
	return(z);
}

matrix line_to_matrix(circ)
	circle circ;
{
	matrix z;
	complex unit;
	static complex i = {0.0,1.0};

	unit.x = cos(-PI*circ.r);
	unit.y = sin(-PI*circ.r);
	z.a = mult(i,unit);
	z.b = mult(i,sub(  mult(cx_conj(unit),circ.c),
			   mult(unit,cx_conj(circ.c))   ));
	z.c.x = 0.0; z.c.y = 0.0; 
	z.d = mult(i,cx_conj(unit));
	return(z);
}

circle matrix_to_circle(x)
	matrix x;
{
	circle circ;

	if(x.c.x < LINEFUZZ && x.c.x > -LINEFUZZ) return(matrix_to_line(x));
	circ.c = divide(x.a,x.c);
	circ.r = 1/x.c.x;
	if(circ.r < 0) circ.r = -circ.r;
	return(circ);
}

circle matrix_to_line(x)
	matrix x;
{
	circle circ;

	circ.r =   -arg(x.a)/PI - 1.5;
	circ.c.x = -x.b.x/2; 
	circ.c.y = -x.b.y/2; 
	circ.c = mult(circ.c,x.a);
	return(circ);
}

circle image_circle(x,c)
	matrix x;
	circle c;
{
	matrix xi, y;

	y = circle_to_matrix(c);
	xi = inverse(x); 
	y = product(x,product(y,xi));
	c = matrix_to_circle(y);
	return(c);
}

/* Vector subroutines */
vector vadd(v,w)
	vector v, w;
{
	vector s;
	s.x = v.x + w.x; s.y = v.y + w.y; s.z = v.z + w.z;
	return(s);
}

vector vsub(v,w)
	vector v, w;
{
	vector s;
	s.x = v.x - w.x; s.y = v.y - w.y; s.z = v.z - w.z;
	return(s);
}

double vabs(v)
	vector v;
{
	return(sqrt(v.x*v.x + v.y*v.y + v.z*v.z));
}

/* Line subroutines */
double cot(a)
	double a;
{
	return(tan(PI/2-a));
}
	
/* Return x so line contains (x,y) */
double meetx(circ,y)
	circle circ;
	double y;
{
	double rots, x;

	rots = -circ.r;
	if(rots == 0.0 || rots == 1.0) return(INF);
	x = (y-circ.c.y)*cot(PI*rots) + circ.c.x; 
	return(x);
}

/* Return y so line contains (x,y) */
double meety(circ,x)
	circle circ;
	double x;
{
	double rots, y;

	rots = -circ.r;
	if(rots == 0.5) return(INF);
	y = (x-circ.c.x)*tan(PI*rots) + circ.c.y; 
	return(y);
}

/* Sphere subroutines */
vector complex_to_s2(z)
	complex z;
{
	vector v;
	double r;

	r = z.x*z.x + z.y*z.y;
	v.x = 2*z.x/(1+r);
	v.y = 2*z.y/(1+r);
	v.z = (1-r)/(1+r);
	return(v);
}

/* Convert circle in plane to circle on S^2 */
/* Method is to find a pair of opposite points on circle */
/* Return norm vector; put radius in R^3 in r*/
vector circle_to_s2(circ,r)
	circle circ;
	double *r;
{
	double ac, s;
	complex w1, w2, c, cn;
	vector v, v1, v2;

	if(circ.r <= 0) return(line_to_s2(circ,r));

	c = circ.c; 
	s = circ.r;
	ac = cx_abs(c);
/* Circle centered at origin */
	if(ac == 0.0)
	{	v.x = v.y = 0; v.z = 1;
		if(s>1.0) v.z = -1;
		*r = 2*s/(1+s*s);
		return(v);
	}
/* Else find points closest and farthest from origin */
	cn.x = c.x/ac; cn.y = c.y/ac;
	w1.x = c.x+s*cn.x; w1.y = c.y+s*cn.y; 
	w2.x = c.x-s*cn.x; w2.y = c.y-s*cn.y; 

	v1 = complex_to_s2(w1);
	v2 = complex_to_s2(w2);
	v  = vadd(v1,v2);
	*r = vabs(vsub(v1,v2))/2;
/* Take care if we have a great circle */
	if(vabs(v) < VEPS) {
		s = sqrt(v1.x*v1.x+v1.y*v1.y);
		v.x = -v1.x * v1.z / s;
		v.y = -v1.y * v1.z / s;
		v.z = s;
	}
	return(v);
}

/* Convert line in plane to circle on S^2 */
vector line_to_s2(circ,r)
	circle circ;
	double *r;
{
	vector p, v;
	static vector north_pole = {0.0,0.0,-1.0};
	double norm;

	p = complex_to_s2(circ.c);
	if(p.x == 0.0 && p.y == 0.0) 
	{	v.y = cos(-PI*circ.r); 
		v.x = -sin(-PI*circ.r);
		v.z = 0.0;
		*r  = 1.0;
	}
	else
	{	v  = vadd(p,north_pole);
		*r = vabs(vsub(p,north_pole))/2;
	}
	return(v);
}
