function w0=InitialW(n,p);
% INITIALW constructs initial vorticity profiles 
%   w0=InitialW(n,p); constructs an initial vorticity profile 
%   w0 for the Navier Stokes solver with 2^n points in the x 
%   and y direction. The maximum velocities are scaled to be one.
%   According to p the following profiles are chosen:
%
%   default: simple harmonic profile for testing purposes
%   p(1)=1 : random vorticity with gaussian amplitude about k0=1/2*kmax or p(2)
%            and width p(3)
%   p(1)=2 : flow past a cylinder profile with gaussian distribution
%            amplitude p(2) and decay p(3) and p(4) random noise 

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=repmat(i*k',1,m);                        % derivative operators
Dy=repmat(i*k,m,1);
DS=Dx.^2+Dy.^2;                             % Laplacian in Fourier
DS(1,1)=1; DD=1./DS; DD(1,1)=0;             % Laplace solver, constant=0
DS(1,1)=0;

if p(1)==1,
  kx=repmat(k',1,m);
  ky=repmat(k,m,1);
  if length(p)==1, p(2)=m/4; p(3)=1; end;
  wh=exp(-p(3)*(p(2)-sqrt(kx.^2+ky.^2)).^2).*exp(i*2*pi*rand(m)); 
  wh(1,1)=0;
  u0=real(ifft2(-Dy.*wh.*DD));              % corresponding velocity field
  v0=real(ifft2(Dx.*wh.*DD));               % for scaling purposes
  ma=max(max(abs(u0)));
  u0=u0/ma; v0=v0/ma;
  w0=real(ifft2(Dx.*fft2(v0)-Dy.*fft2(u0)));
elseif p(1)==2,
  if length(p)==1, p(2)=1;p(3)=1;p(4)=0.01; end;
%  v0=zeros(m,m);
  v0=p(4)*sin(X);                           % small disturbance
  u0=p(2)*exp(-p(3)*(Y-pi).^2);
  u0=u0/max(max(abs(u0)));                  % scale max to one
  uh=fft2(u0); 
  uh(1,1)=0;                                % mean zero
  u0=real(ifft2(uh)); 
  vh=fft2(v0);
  wh=Dx.*vh-Dy.*uh;
  w0=real(ifft2(wh));
else
  u0=sin(2*Y);
  v0=cos(X);
  wh=Dx.*fft2(v0)-Dy.*fft2(u0);
  w0=real(ifft2(wh));
end;
PlotAll(w0,0);



