% This produces a movie in the avi format showing the solution of 
% u_t + 6u^2u_x + u_xxx - (b(x,t)u)_x = 0 on [-pi,pi].
% The code is a simple modification of Trefethen's example p27.m
% The initial data is given by a double soliton with parameters
% a1, a2 (centers) and c1, c2 (velocities).
% The arguments of are time t, with parameters given in a vector form
% [a1,a2,c1,c2], a function bb = b(x,t), and log_2 of the number of 
% x grid points; the defauls is 8 but it may not work for all data.
% The last input variable is either 0 or 1: for 0 it shows 
% the solution of the PDE only and for 1 it shows the effective
% dynamics simultaneously with the solution. The default value is 1.
%
% The code effdyn.m has to be in the same directory of option 1 is
% used.
% 
% Here is an example
% 
% Bmovie(@(x,t) 40*cos(2*x+3-10^2*t)+30*sin(x+1+10^2*t),0.05,[-1,-2,6,10],8)
%
% The avi file is returned under the name bmovie.avi and can then be 
% changed or converted.


function bmovie(fun,tmax,ZZ,N,del)

     if (nargin < 4)

     N = 8;

     end


     if (nargin < 5)

     del = 1;

     end

N=2^N;

% Set up grid and two-soliton initial data:

dt = .4/N^2; x = (2*pi/N)*(-N/2:N/2-1)';

drawnow, set(gcf,'renderer','zbuffer')

u=sqrt(2)*qu(x,ZZ);

v = fft(u); k = [0:N/2-1 0 -N/2+1:-1]'; ik3 = 1i*k.^3;

% Solve PDE and plot results:
nplt = 100; nmax = round(tmax/dt);
udata = u; tdata = 0; 
for n = 1:nmax

t = n*dt; g = -i*dt*k; 

bb=fun(x,t);

E = exp(dt*ik3/2); E2 = E.^2;
a = g.*fft(real( ifft(     v    ) ).^3 - ifft( v ).*bb);
b = g.*fft(real( ifft(E.*(v+a/2)) ).^3 - ifft(E.*(v+a/2)).*bb); 
c = g.*fft(real( ifft(E.*v + b/2) ).^3 - ifft(E.*v + b/2).*bb);     
d = g.*fft(real( ifft(E2.*v+E.*c) ).^3 -  ifft(E2.*v+E.*c).*bb);
v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
if mod(n,nplt) == 0 
u = real(ifft(v)); %waitbar(n/nmax)
     udata = [udata u]; tdata = [tdata t];
    end
  end

udata = udata/sqrt(2);

     [K,P] = size(tdata);

T=tdata;

z1=0;
z2=0;

for j=1:P

z1=max(z1,max(udata(:,j)));
z2=min(z1,min(udata(:,j)));

end



if (del == 1) 

[T,Y]=effdyn(fun,T,ZZ,0);


w1=0;
w2=0;

for j=1:P

w1=max(w1,max(2*Y(j,[3:4])));
w2=min(w2,min(2*Y(j,[3:4])));

end


for j=1:P
   plot(x,udata(:,j),'k',x,qu(x,Y(j,:)),'LineWidth',2);
     xlim([-pi,pi])
     ylim([z2-1,z1+1])
     pause(0.02)
     fr(j)=getframe;
end

else

for j=1:P
    plot(x,udata(:,j),'k','LineWidth',2);
     xlim([-pi,pi])
     ylim([z2-1,z1+1])
     pause(0.02)
     fr(j)=getframe;
end

end

movie2avi( fr, 'bmovie.avi', 'fps', 5, 'quality', 100 ) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function qu = qu(x,ZZ)

     k1 = ZZ(3)*(x - ZZ(1));
     k2 = ZZ(4)*(x - ZZ(2));

     ga1 =  exp(-k1);
     ga2 = - exp(-k2);

     m11=(1+ga1.^2)./(2*ZZ(3));
     m22=(1+ga2.^2)./(2*ZZ(4));
     m12=(1+ga1.*ga2)./(ZZ(3)+ZZ(4));


     M=m11.*m22 - m12.^2;
     M1=ga1.*(m12 - m22) + ga2.*(m12-m11);
     qu = M1./M;

