clear
%load Mar28_data.mat

N = 128*4;
dx = 2*pi/N;
x1 = 0:dx:2*pi-dx;
u0 = cos(x1);

t_0 = 0;
t_f = 1;
M = 1250;
D = 0;
[u1,x1,t1,kk1,amp1] = burgers_if(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:25, 
    jj = j*50+1;
    t1(jj)
    subplot(2,1,1)
    plot(x1,u1(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 0');
    subplot(2,1,2)
    plot(kk1,amp1(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end
t_0 = 0;
t_f = 1;
M = 2*1250;
D = 0;
[u2,x1,t2,kk1,amp2] = burgers_if(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:25, 
    jj = j*2*50+1;
    t2(jj)
    subplot(2,1,1)
    plot(x1,u2(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 0');
    subplot(2,1,2)
    plot(kk1,amp2(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end
t_0 = 0;
t_f = 1;
M = 4*1250;
D = 0;
[u3,x1,t3,kk1,amp3] = burgers_if(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:25, 
    jj = j*4*50+1;
    t3(jj)
    subplot(2,1,1)
    plot(x1,u3(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 0');
    subplot(2,1,2)
    plot(kk1,amp3(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end
u1_if = u1;
amp1_if = amp1;
u2_if = u2;
amp2_if = amp2;
u3_if = u3;
amp3_if = amp3;

t_0 = 0;
t_f = 1;
M = 1250;
D = 0;
[u1,x1,t1,kk1,amp1] = burgers_ab(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:25, 
    jj = j*50+1;
    t1(jj)
    subplot(2,1,1)
    plot(x1,u1(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 0');
    subplot(2,1,2)
    plot(kk1,amp1(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end
t_0 = 0;
t_f = 1;
M = 2*1250;
D = 0;
[u2,x1,t2,kk1,amp2] = burgers_ab(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:25, 
    jj = j*2*50+1;
    t2(jj)
    subplot(2,1,1)
    plot(x1,u2(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 0');
    subplot(2,1,2)
    plot(kk1,amp2(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end
t_0 = 0;
t_f = 1;
M = 4*1250;
D = 0;
[u3,x1,t3,kk1,amp3] = burgers_ab(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:25, 
    jj = j*4*50+1;
    t3(jj)
    subplot(2,1,1)
    plot(x1,u3(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 0');
    subplot(2,1,2)
    plot(kk1,amp3(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end
u1_ab = u1;
amp1_ab = amp1;
u2_ab = u2;
amp2_ab = amp2;
u3_ab = u3;
amp3_ab = amp3;

figure(1), clf, hold on
plot(kk1,amp1_if(:,101),'LineWidth',3)
plot(kk1,amp1_if(:,501),'LineWidth',3)
plot(kk1,amp1_if(:,751),'LineWidth',3)
axis([0,N/2,-17,2])
% t1(101),t1(501),t1(751)
% solution at times 0.08, 0.40, and 0.60
% print -dps Burgers_no_diffusion_EE_coarse.ps

figure(2), clf, hold on
plot(kk1,amp2_if(:,201),'LineWidth',3)
plot(kk1,amp2_if(:,1001),'LineWidth',3)
plot(kk1,amp2_if(:,1501),'LineWidth',3)
axis([0,N/2,-17,2])
% t2(201),t2(1001),t2(1501)
% solution at times 0.08, 0.40, and 0.60
% print -dps Burgers_no_diffusion_EE_medium.ps


figure(3), clf, hold on
plot(kk1,amp3_if(:,401),'LineWidth',3)
plot(kk1,amp3_if(:,2001),'LineWidth',3)
plot(kk1,amp3_if(:,3001),'LineWidth',3)
axis([0,N/2,-17,2])
% t3(401),t3(2001),t3(3001)
% solution at times 0.08, 0.40, and 0.60
% print -dps Burgers_no_diffusion_EE_fine.ps


figure(4), clf, hold on
plot(kk1,amp1_ab(:,101),'LineWidth',3)
plot(kk1,amp1_ab(:,501),'LineWidth',3)
plot(kk1,amp1_ab(:,751),'LineWidth',3)
axis([0,N/2,-17,2])


N = 128;
dx = 2*pi/N;
x1 = 0:dx:2*pi-dx;
u0 = cos(x1);

t_0 = 0;
t_f = 2;
M = 5000;

D = 2;
[u1,x1,t,kk1,amp1] = burgers_if(u0,D,t_0,t_f,M,N);
for j=0:40, 
    jj = j*125+1;
    subplot(2,1,1)
    plot(x1,u1(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 2');
    subplot(2,1,2)
    plot(kk1,amp1(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end

D = 1;
[u2,x1,t,kk1,amp2] = burgers_if(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:40, 
    jj = j*125+1;
    subplot(2,1,1)
    plot(x1,u2(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 1');
    subplot(2,1,2)
    plot(kk1,amp2(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end

D = 1/2;
[u3,x1,t,kk1,amp3] = burgers_if(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:40, 
    jj = j*125+1;
    subplot(2,1,1)
    plot(x1,u3(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 1/2');
    subplot(2,1,2)
    plot(kk1,amp3(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end

D = 1/4;
[u4,x1,t,kk1,amp4] = burgers_if(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:40, 
    jj = j*125+1;
    subplot(2,1,1)
    plot(x1,u4(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 1/4');
    subplot(2,1,2)
    plot(kk1,amp4(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end

D = 1/8;
[u5,x1,t,kk1,amp5] = burgers_if(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:40, 
    jj = j*125+1;
    subplot(2,1,1)
    plot(x1,u5(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 1/8');
    subplot(2,1,2)
    plot(kk1,amp5(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end

% need to point-double!
N = 2*N;
dx = 2*pi/N;
x2 = 0:dx:2*pi-dx;
u0 = cos(x2);
D = 1/16;
[u6,x2,t,kk2,amp6] = burgers_if(u0,D,t_0,t_f,M,N);
figure(1)
for j=0:40, 
    jj = j*125+1;
    subplot(2,1,1)
    plot(x2,u6(:,jj)), axis([0,2*pi,-1.1,1.1]), 
    title('D = 1/16');
    subplot(2,1,2)
    plot(kk2,amp6(:,jj)), axis([0,N/2,-18,2])
    figure(1), pause(1) 
end

% need to point-double!
N = 2*N;
dx = 2*pi/N;
x3 = 0:dx:2*pi-dx;
u0 = cos(x3);
D = 1/32;
[u7,x3,t,kk3,amp7] = burgers_if(u0,D,t_0,t_f,M,N);


% time t=1
t(2501)
figure(1), clf, hold on
plot(x1,u1(:,2501),'LineWidth',3)
%plot(x1,u2(:,2501),'LineWidth',3)
plot(x1,u3(:,2501),'LineWidth',3)
plot(x1,u4(:,2501),'LineWidth',3)
plot(x1,u5(:,2501),'LineWidth',3)
plot(x2,u6(:,2501),'LineWidth',3)
plot(x3,u7(:,2501),'LineWidth',3)
plot(x1,u3_ab(:,length(t3)),'--','LineWidth',3);
axis([0,2*pi,-1,1])
% print -dps Burgers_with_diffusion_shock_time.ps
figure(1)

figure(2), clf, hold on
plot(kk1,amp1(:,2501),'LineWidth',3)
%plot(kk1,amp2(:,2501))
plot(kk1,amp3(:,2501),'LineWidth',3)
plot(kk1,amp4(:,2501),'LineWidth',3)
plot(kk1,amp5(:,2501),'LineWidth',3)
plot(kk2,amp6(:,2501),'LineWidth',3)
plot(kk3,amp7(:,2501),'LineWidth',3)
plot(kk1,amp3_ab(:,length(t3)),'--','LineWidth',3);
axis([0,N/2,-17,2])
% print -dps Burgers_with_diffusion_PS_shock_time.ps
figure(2)



