%	Transformation from circle to ellipse
%	z=2*lambda*cosh(omega)
%	lambda=(1/2)*(a^2-b^2)^(1/2)
%	a and b real, a > b
%	omega=r+i.*theta
%	Domain is between theta = 0 and theta = 2*pi for r => ro
%	ro = (1/2)*log((a+b)/(a-b))
clear all
in=input('(x/a)^2+(y/b)^2=1, Input value of [a b], a > b: ');
a=in(1);
b=in(2);
if a<b
   error('a > b necessary for ececution')
end
clf reset
ro = (1/2)*log((a+b)/(a-b));
c=(1/2)*(a+b);
lambda=(1/2)*(a^2-b^2)^(1/2);
rmax=acosh(3*a/(2*lambda));
[th,r]=meshgrid((0:10:360)*pi/(180),ro:0.05*rmax:rmax);
omega=r+i.*th;
%z=2*lambda.*cosh(omega);
%x=real(z);
%y=imag(z);
x=2*lambda.*cosh(r).*cos(th);
y=2*lambda.*sinh(r).*sin(th);
z=x+i.*y;
contour(x,y,r,20,'r')
hold on
contour(x,y,th,39,'k')
plot(x(1,:),y(1,:),'g')
plot(x(:,1),y(:,1),'k')
axis image
xlabel('x')
ylabel('y')
title(['Confocal Ellipses and Conjugate Hyperbolae, a = '...
      ,num2str(a),' b = ',num2str(b)])
