% This script solves u_{tt} = mu*u_{xx} + f_x(t,u)
% by a pseudospectral method.

clear all;

%==========================================================================
%
% In "run_twoway.m" set
%
%   t0, tf, dtplot, mu
%   x0, x1, ord, pltord
%   epsilon, ltrunc
%   initial condition w0
%
%==========================================================================

%% Set time parameters
%%
t0 = 0;             % initial time
tf = 1.867;             % final time
mu = 4;             % bifurcation parameter

dtplot = (tf-t0)/101;       % time between solution plots
tspan = t0:dtplot:tf;

%% Set space parameters
%%
x0 = 0;             % left endpoint
x1 = 2*pi;          % right endpoint
ord = 13;            % n = 2^ord  space grid points
pltord = 9;         % use 2^pltord points for surface plots

xlength = x1-x0;        
nord = 2^ord;
dx = xlength/nord;
x = x0+dx*(0:(nord-1));
lx = length(x);
k = make_k(lx);

%% Set viscosity parameters
%%
epsilon = dx;       % coefficient of numerical viscosity
ltrunc = 1;         % viscosity applied to wavenumbers k > ltrunc
% ltrunc = round(sqrt(lx));

%% Define initial conditions  w0 = [u(x,0),u_t(x,0)]
%%
w0 = wzero(x,k,lx,mu);

%% Solve the equation
%%
%    options = odeset('AbsTol', 1.0e-6, 'RelTol', 1.0e-6);
%    [T, u] = ode45 (@F, [t0 tf], u0, options);
%[t,w] = ode15s(@h,tspan,w0,[],k,lx,ltrunc,epsilon,mu);
[t,w] = ode45(@h,tspan,w0,[],k,lx,ltrunc,epsilon,mu);

%% Plot the solution
%%
clf
u = w(:,1:lx);
 
%% Surface plot
%%  
    
    figure(1);
    
    lt = length(t);
    
    pltord = min(pltord,ord);
    dxplot = 2^(ord-pltord);
    
    nxplot = 1:dxplot:lx;
    lxp = length(nxplot);
    xplot = x(nxplot);
    uplot = u(:,nxplot);
    
    waterfall(xplot,t,uplot)
    view(10,70)
%    umax = max(max(u));
%    umin = min(min(u));   
%    axis([0 2*pi t0 tf umin umax])
    axis tight
    xlabel('x'), ylabel('t'), zlabel('u')
    grid off

    figure(2);
    
    meshc(xplot,t,uplot)
    view(10,70)
    axis tight
    xlabel('x'), ylabel('t'), zlabel('u')
    grid off
   

nmovie = 1;

if (nmovie==1)
    figure(3);
    hold off;
    
    lt = length(t); 
    
    umax = max(max(u));
    umin = min(min(u));
    textx = x1 -(x1-x0)/5;
    textu = umax - (umax-umin)/10;
 
    for j = 1:lt
      umplot = u(j,:);
      plot(x,umplot);
      title(strcat('time =',num2str(t(j))));
      axis([x0 x1 umin umax]);
      text(textx,textu,strcat('time =',num2str(t(j))));   
      xlabel('x'); ylabel('u');
      mvplot(j) = getframe;
    end
    
movie2avi(mvplot,'hiz.avi','fps',4,'quality',100,'compression','None');
movie(mvplot,1);
    
end