%establishing taylor vortex
R = 0.2;
N = 64;
T = 2;
j = (-N/2):((N/2)-1):N;
x = 2*j/N;
y = 2*j/N;

[xn,yn] = meshgrid(x,y);

r = sqrt(xn.^2 + yn.^2);
wxy = (2.-(xn.^2 + yn.^2)./(R.^2)).*exp(0.5.*(1.-(xn.^2 + yn.^2)./(R.^2)))./R;

%fft of wxy
wmn = fftshift(fft2(wxy,N,N)/(N.^2));

%ODDball wavenumber zero
wmn(1,:) = 0;
wmn(64,:) = 0;
wmn(:,1) = 0;
wmn(:,64) = 0;

%establishing function phi_mn and fft of u,v
m = j./2;
n = j./2;
[p,q] = meshgrid(m,n);

Re = 200;
t0 = linspace(0,0.1,64);

%theoritical omega value
w_st = (2./(1+((2.*t0)./(Re.*R.^2))).^2).*(1-(r./R)./2.*(1+((2.*t0)./(Re.*R.^2)))).*exp(0.5.*(1-(r./R)./(1+((2.*t0)./(Re.*R.^2)))));


w_st(1,:) = 0;
w_st(64,:) = 0;
w_st(:,1) = 0;
w_st(:,64) = 0;

phi_mn = wmn./(((2.*pi())).^2).*(p.^2 + q.^2);
phi_mn(33,33) = 0;
umn = 1i.*(2.*pi().*q).*phi_mn;
vmn = -1i.*(2.*pi().*p).*phi_mn;

%convolution and dealiasing
M = N*3/2;
upad = zeros(M,M);
vpad = zeros(M,M);
wpad = zeros(M,M);
indpad = (17:80);
upad(indpad,indpad) = umn;
vpad(indpad,indpad) = vmn;
wpad(indpad,indpad) = wmn;
uij2 = ifft(ifftshift(upad).*M.^2);
vij2 = ifft(ifftshift(vpad).*M.^2);
wij2 = ifft(ifftshift(wpad).*M.^2);

%oddball wavenumber zero
uw1 = uij2.*wij2;
uw2 = fftshift(fft2(uw1,M,M)./(M.^2));
uw2(1:16,:) = [];
uw2(65:80,:) = [];
uw2(:,1:16) = [];
uw2(:,65:80) = [];
%
uw2(1,:) = 0;
uw2(64,:) = 0;
uw2(:,1) = 0;
uw2(:,64) = 0;

vw1 = vij2.*wij2;
vw2 = fftshift(fft2(vw1,M,M)./(M.^2));
vw2(1:16,:) = [];
vw2(65:80,:) = [];
vw2(:,1:16) = [];
vw2(:,65:80) = [];
%
vw2(1,:) = 0;
vw2(64,:) = 0;
vw2(:,1) = 0;
vw2(:,64) = 0;

w2 = wij2.*wij2;
wmn2 = fftshift(fft2(w2,M,M)./(M.^2));
wmn2(1:16,:) = [];
wmn2(65:80,:) = [];
wmn2(:,1:16) = [];
wmn2(:,65:80) = [];
%
wmn2(1,:) = 0;
wmn2(64,:) = 0;
wmn2(:,1) = 0;
wmn2(:,64) = 0;


%reshaping matrices to column matrix for ODE45
pN = reshape(p,[4096,1]);
qN = reshape(q,[4096,1]);
wxyN = reshape(wxy,[4096,1]);
wmnN = reshape(wmn2,[4096,1]);
uwN = reshape(uw2,[4096,1]);
vwN = reshape(vw2,[4096,1]);

%differential eqn
wdotN = -1i.*2.*pi().*pN.*(uwN)-1i.*2.*pi().*qN.*(vwN)-1./Re.*4.*(pi().^2).*(pN.^2 + qN.^2);
[t,w] = ode45(@(t,w) wdotN,[0,1],wxyN);

%reshape back to NxN matrix
wf = reshape(w(1,:),[64,64]);

%plot resulting in omega and theoriticak/analytical omega
hold on
plot(x(33:64),wf(33,33:64));
plot(x(33:64),w_st(33,33:64));