/* Routines for doing DLA */
/* Absolute derivative of function is dfunction */

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

static complex cx_1={1.0,0.0};
static complex cx_1pi={1.0,1.0};
static complex cx_1mi={1.0,-1.0};

/* Mobius map from Sigma to R = right half plane */
/* m(z) = (z-1)/(z+1) */
complex m(z)
	complex z;
{
	return(div(sub(z,cx_1),add(z,cx_1)));
}

double dm(z)
	complex z;
{
	double d;
	d=(z.x+1)*(z.x+1) + z.y*z.y;
	return(2/d);
}

/* Mobius map from R to Sigma */
/* mi(z) = (1+z)/(1-z) */
complex mi(z)
	complex z;
{
	return(div(add(z,cx_1),sub(cx_1,z)));
}

double dmi(z)
	complex z;
{
	double d;
	d=(z.x-1)*(z.x-1) + z.y*z.y;
	return(2/d);
}

/* Basic slit mapping, Sqrt[z^2+1], R->R */
complex slit(z)
	complex z;
{
	complex w;
	w = add(cx_1,mult(z,z));
	if(z.y > 0)   w=contsqrt(w,cx_1pi);
	else	      w=contsqrt(w,cx_1mi);
	return(w);
}

double dslit(z)
	complex z;
{
	return(cabs(div(z,slit(z))));
}

/* Basic bend mapping, z + Sqrt[z^2+1], R->R */ 
complex bend(z)
	complex z;
{	
	return(add(z,slit(z)));
}

double dbend(z)
	complex z;
{	
	return(cabs(add(cx_1,div(z,slit(z)))));
}

/* Map Sigma -> Sigma - lune centered at t*/
complex f(z,t,a)
	complex z, t;
	double a;
{
	complex w;
	double r, s;

	r = 2*(1+a)/(a*(2+a));
	s = 1/(r+sqrt(r*r+1));
	w = m(div(z,t));
	w.x *= r; w.y *= r;
	w = bend(w);
	w.x *= s; w.y *= s;
	w = mult(t,mi(w));
	return(w);
}

double df(z,t,a)
	complex z, t;
	double a;
{
	complex w;
	double d, r, s;

	r = 2*(1+a)/(a*(2+a));
	s = 1/(r+sqrt(r*r+1));
	w = m(div(z,t));
	d = dm(div(z,t));
	w.x *= r; w.y *= r;
	d *= r*dbend(w);
	w = bend(w);
	w.x *= s; w.y *= s;
	d *= s*dmi(w);
	return(d);
}

/* Capacity -- radius of mapping */
double capf(a)
	double a;
{
	double c, r, s;

	r=2*(1+a)/(a*(2+a));
        c=sqrt(1+r*r)/r;
	return(c);
}

/* Argument of critical point */
double cp(a)
	double a;
{
	complex z;
	z.x = 1+a;
	z.y = -1.0;
	z = mult(z,z);
	z = mult(z,z);
	z.x = -z.x;
	z.y = -z.y;
	return(arg(z));
}

/* Argument of critical value */
double cv(a)
	double a;
{
	complex z;
	double r, s;

	r=2*(1+a)/(a*(2+a));
	s = 1/(r+sqrt(r*r+1));
	z.x = 0.0;
	z.y = s;
	z = mi(z);
	return(arg(z));
}
