%
%  2D swift-hoh with ETD -- set-up -- djm, 26 jan 2007
%

global kg lg xg yg Nx Ny epsi Cr Cq u uF Ny32 Nx32 ixpad iypad
warning off MATLAB:log:logOfZero

%  set-up choice (0 = continue from final values)
choice = 1;

%  run parameters
switch choice
    case{0,1}
        epsi = 0.2;  Cr = epsi^2;  Cq = 1.0;
        Nx = 64;  Ny = Nx;  Lx = sqrt(2)*2*pi;  Ly = (Ny/Nx)*Lx;
        dt = 0.2;  Nt = 10;  Nc = 50;
end

%  coordinates
ds = Lx/Nx;  xx = 0:ds:Lx-ds;  yy = 0:ds:Ly-ds;
[xg,yg] = meshgrid(xx,yy);

%  initial condition
switch choice
    case{0}
    case{1}
        disp('IV:  pitchfork bifurcation test');
        u = 0.3*sin(2*pi*xg/Lx).*sin(2*pi*yg/Ly);
        u = u./(1 + 0.3*sin(10*pi*xg/Lx).*sin(10*pi*yg/Ly));
        time = 0;  
        
        load deaColors;  cmap = red2blu;
        mkmovie = 0;  mov_dir = 'mov_dir';  fnum = 1;
end

%  spectral IVs
uF = fft2(u);  uF(Ny/2+1,:) = 0;  uF(:,Nx/2+1) = 0;

%  small spectral derivatives
kk = fftshift([0 -Nx/2+1:Nx/2-1])*(2*pi/Lx);
ll = fftshift([0 -Ny/2+1:Ny/2-1])*(2*pi/Ly);
klmax = max([max(kk(:)) max(ll(:))]);

[kg,lg] = meshgrid(kk,ll);

%  zero-padding for pv
Nx32=Nx*3/2;  Ny32=Ny*3/2;
ixpad = [1:Nx/2 Nx+1:Nx32];  iypad = [1:Ny/2 Ny+1:Ny32];

%  exponential TD set-up
lam = Cr - (Cq^2 - (kg.^2 + lg.^2)).^2;

c0 = exp(lam*dt);  ch = exp(lam*dt/2);

%  complex means
Mc = 32;  r = exp(i*2*pi*(reshape((1:Mc),1,1,Mc)-.5)/Mc);
LR = dt*lam(:,:,ones(Mc,1)) + r(ones(Ny,1),ones(Nx,1),:);

%  ab3
c1ab = dt*real(mean( (exp(LR) - 3 + (3*exp(LR) - 5)./(2*LR) + ...
    (exp(LR) - 1)./(LR.*LR))./LR ,3));
c2ab = dt*real(mean( (3 + (4 - 2*exp(LR))./LR + ...
    (2 - 2*exp(LR))./(LR.*LR))./LR,3));
c3ab = dt*real(mean( (-1 + (exp(LR) - 3)./(2*LR) + ...
    (exp(LR) - 1)./(LR.*LR))./LR ,3));

%  rk3
c1rk = dt*real(mean(           (exp(LR/2) - 1)./LR            ,3));
c2rk = dt*real(mean(           (exp(LR  ) - 1)./LR            ,3));
c3rk = dt*real(mean(  (-4-LR+exp(LR).*(4-3*LR+LR.^2))./LR.^3  ,3));
c4rk = dt*real(mean(    (4+2*LR+exp(LR).*(-4+2*LR))  ./LR.^3  ,3));
c5rk = dt*real(mean(  (-4-3*LR-LR.^2+exp(LR).*(4-LR))./LR.^3  ,3));

clear rhstemp MC r LR

%  time stepping
SwiftH_ab3