% ジューコスフキー翼のまわりの流線を描くサンプルプログラム
% 翼の内外の判定の部分を改訂したもの
%-------------------------------------------------------

% パラメターの設定

u=1
a=1
b=0.8
alpha=0
beta=10

% alpha, beta をラジアンに直す

alpha=alpha*pi/180
beta=beta*pi/180

% Γの計算

gamma=-4*pi*u*a*sin(alpha+beta)

% 座標の設定

x=-3:0.025:3;
y=-1:0.025:1.5;

% 速度ポテンシャルの計算

for k=1:length(x)
  for l=1:length(y)
    z=x(k)+i*y(l);
    zz=(z+sqrt(z*z-4*b*b))/2;
    zeta=zz-b+a*exp(-i*beta); 
    if(abs(zeta)< a)
      zz=(z-sqrt(z*z-4*b*b))/2;
      zeta=zz-b+a*exp(-i*beta); 
    end
    if(abs(zeta)>=a)    
      f(k,l)=u*(zeta*exp(-i*alpha)+a*a/zeta*exp(i*alpha)-i*gamma/(2*pi)*log(zeta));
    else
      f(k,l)=0;
    end
  end
end

psi=imag(f);

% 等高線の表示

gset size ratio (2.5/6)

contour(psi,[-1.0:0.05:1.7],x,y)