function hh = catfit(p)

if(length(p)<12)
	hh(1:5) = NaN;
	hh = hh';
	return;
end

p = princomp(p')'*p;

xo = p(1,:);
yo = -p(2,:);
x = [xo; yo];
n = length(xo);

a = 4;
c = -4;
tx = mean(xo);
ty = mean(yo);
theta = 0;
ts = 0:n-1;
lamb = 100;
nr = Inf;
it = 0;

while(it < 100)
	it = it +1;
	
	h = [a c tx ty theta ts(2:end)]';
	
	xc =  [tx + cos(theta)*(c + ts) - a*cosh((c + ts)/a)*sin(theta)
		ty + sin(theta)*(c + ts) + a*cosh((c + ts)/a)*cos(theta)];
	
	da = [(sinh((c + ts)/a)*sin(theta).*(c + ts))/a - cosh((c + ts)/a)*sin(theta)
		cosh((c + ts)/a)*cos(theta) - (sinh((c + ts)/a)*cos(theta).*(c + ts))/a];
 
 
	dc =	[cos(theta) - sinh((c + ts)/a)*sin(theta)
			 sin(theta) + sinh((c + ts)/a)*cos(theta)];
 
	dtx = repmat([1; 0], 1,n);
	dty = repmat([0;1], 1, n);
	dtheta = [- sin(theta)*(c + ts) - a*cosh((c + ts)/a)*cos(theta)
	cos(theta)*(c + ts) - a*cosh((c + ts)/a)*sin(theta)];
 
	J1 = [da(:) dc(:) dtx(:) dty(:) dtheta(:)];
	J2 = zeros(2*n, n-1);
	J2(3:(2*n+2):end) = cos(theta) - sinh((c + ts(2:end))/a)*sin(theta);
	J2(4:(2*n+2):end) = sin(theta) + sinh((c + ts(2:end))/a)*cos(theta);
	
	r = xc(:) - x(:);
	
	J = [J1 J2];
	
	nrc = norm(r);
	if(abs(nr - nrc)< 0.001)
		break;
	end
	nr = nrc;
	
	JJ = J'*J;
	h = h - inv(JJ + lamb*diag(diag(JJ)))*J'*r;
	lamb = lamb*0.8;
	
	xx = [cos(theta) -sin(theta); sin(theta) cos(theta)]'*xc;
	
	a = h(1);
	c = h(2);
	tx = h(3);
	ty = h(4);
	theta = h(5);
	ts(2:end) = h(6:end)';
end
hh = [a c tx ty theta]';






