%	Potential flow past ellipse
%	Velocity at infinity is (-U,-V) 
%	or -U-i*V=-magv*exp(-i*alpha)
%	w(z)=-magv*(a+b)*cosh(omega-ro+i*alpha)	
%	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
in1=input('(x/a)^2+(y/b)^2=1, Input value of [a b], a > b: ');
a=in1(1);
b=in1(2);
if a<b
   error('a > b necessary for ececution')
end
in2=input('-U-i*V=-magv*exp(-i*alpha), Input [magv alpha]: ');
magv=in2(1);
alpha=in2(2);
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(2*a/(2*lambda));
[th,r]=meshgrid((0:3:360)*pi/(180),ro:0.01*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)])
pause 
clf reset
w=z;  %allocate memory
w=-magv*(a+b).*(cosh(r-ro).*cos(th+alpha)...
   +i.*sinh(r-ro).*sin(th+alpha));
contour(x,y,real(w),21,'r')
hold on
contour(x,y,imag(w),21,'k')
plot(x(1,:),y(1,:),'g')
axis image
xlabel('x')
ylabel('y')
title('Potentials and Streamlines for Flow Past Ellipse')     
