clear all
n=6;
m=2^n;
h=2*pi/m;
x=(0:h:2*pi-h)';
X=repmat(x,1,m);
Y=repmat(x',m,1);
k=[0:m/2-1 -m/2:-1];
 % Dx und Dy im Fourierraum

 %% easy test 
u0=1;v0=0;
w0=real(ifft2(Dx.*fft2(v0)-Dy.*fft2(u0)));
w=NavierStokes(w0,0.003,10,400,10);
M=plotNS(w);

 %% Kelvin-Helmholtz Instabilitaet
u0=0.5.*(1+tanh(10.* ... 
   (1-4./(2.*pi).*abs(pi-Y)))).* ... 
   (1+0.5.*sin(2.*X));
v0=0;
w0=real(ifft2(Dx.*fft2(v0)-Dy.*fft2(u0)));
w=NavierStokes(w0,0.001,50,1000,10);
M=plotNS(w);

 %% Wirbel-Dipole
psi0=0.01;
 % Wirbel1
xv=2.*pi./4;
yv=2.*pi./2+0.05;
lv=0.4.*sqrt(2).*min([xv,yv,2.*pi-xv,2.*pi-yv]);
u01=-2.*(Y-yv)./(lv.^2).*psi0.*... 
    exp(-((X-xv).^2+(Y-yv).^2)./lv.^2);
v01=2.*(X-xv)./(lv.^2).*psi0.*... 
    exp(-((X-xv).^2+(Y-yv).^2)./lv.^2);
 % Wirbel2
xv=2.*pi./4;
yv=2.*pi./2-0.05;
lv=0.4.*sqrt(2).*min([xv,yv,2.*pi-xv,2.*pi-yv]);
u02=2.*(Y-yv)./(lv.^2).*psi0.*... 
    exp(-((X-xv).^2+(Y-yv).^2)./lv.^2);
v02=-2.*(X-xv)./(lv.^2).*psi0.*... 
    exp(-((X-xv).^2+(Y-yv).^2)./lv.^2);
 % Dipole
u0=u01+u02;
v0=v01+v02;

w0=real(ifft2(Dx.*fft2(v0)-Dy.*fft2(u0)));
w=NavierStokes(w0,0.001,2000,500,10);
M=plotNS(w);
