format compact
format short e
clear

t_0 = 0; t_f = 5; M = 10; N = 32; 
dt = (t_f-t_0)/M;
NPrint = M;
% choose focussing (K=1) or defocussing (K=-1):
K = -1;
% NOTE: the signs of K are different here and in the code than in the class
% notes.  I didn't have a chance to go through and make everything 
% self-consistent...

dx = 2*pi/N;
x = 0:dx:2*pi;
x = x';
A = 4;
B = 3; 
C = 0;
u0 = A*exp(i*(B*x-t_0*(-K*A^2+B^2)+C));
% here's the exact solution for comparison purposes
uf = A*exp(i*(B*x-t_f*(-K*A^2+B^2)+C));

% A comment on using matlab and complex numbers...  First of all, don't use
% "i" in your for loops.  Or, if you don't want to think about that, then
% whenever you want to refer to sqrt(-1) refer to "1i".
%
% also, note that (1+1i)' = (1-1i).  So transposing a matrix also involves
% taking a complex conjugate.  On a practical level, this means that
% exp(i*x') and exp(i*x)' are different even though they have the same
% dimensions!
[u1,x,t1,kk,amp1,M1,H1] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u2,x,t2,kk,amp2,M2,H2] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u3,x,t3,kk,amp3,M3,H3] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u4,x,t4,kk,amp4,M4,H4] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u5,x,t5,kk,amp5,M5,H5] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u6,x,t6,kk,amp6,M6,H6] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u7,x,t7,kk,amp7,M7,H7] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u8,x,t8,kk,amp8,M8,H8] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u9,x,t9,kk,amp9,M9,H9] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u10,x,t10,kk,amp10,M10,H10] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u11,x,t11,kk,amp11,M11,H11] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);

% to watch movies of instabilities
figure(1), clf
for j=0:20, plot(kk,amp2(:,1+j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp3(:,1+2*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp4(:,1+4*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp5(:,1+8*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp6(:,1+16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp7(:,1+2*16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp8(:,1+4*16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp9(:,1+8*16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp10(:,1+16*16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp11(:,1+32*16*j)), figure(1), pause(1), end

% compare the computed solution to the exact solution for evidence
% that the schemes are convergent.
err(1) = max(abs(u1(:,length(t1))-uf));
err(2) = max(abs(u2(:,length(t2))-uf));
err(3) = max(abs(u3(:,length(t3))-uf));
err(4) = max(abs(u4(:,length(t4))-uf));
err(5) = max(abs(u5(:,length(t5))-uf));
err(6) = max(abs(u6(:,length(t6))-uf));
err(7) = max(abs(u7(:,length(t7))-uf));
err(8) = max(abs(u8(:,length(t8))-uf));
err(9) = max(abs(u9(:,length(t9))-uf));
err(10) = max(abs(u10(:,length(t10))-uf));
err(11) = max(abs(u11(:,length(t11))-uf))


% for initial data
%    u0 = A*exp(i*(B*x+0*(-K*A^2+B^2)+C));
% make plots of the instability (defocussing case)
% figure(1)
% clf
% plot(kk,amp5(:,1),'o-');
% hold on
% plot(kk,amp5(:,11),'o-');
% plot(kk,amp5(:,21),'o-');
% plot(kk,amp5(:,31),'o-');
% print -dps split_schroedinger_unstable_spec.ps
% 
% figure(2)
% clf
% plot(kk,amp6(:,1),'o-');
% hold on
% plot(kk,amp6(:,21),'o-');
% plot(kk,amp6(:,41),'o-');
% plot(kk,amp6(:,61),'o-');
% print -dps split_schroedinger_stable_spec.ps

% repeat the experiment using the focussing Schrodinger equation
t_0 = 0; t_f = 5; M = 10; N = 32; 
dt = (t_f-t_0)/M;
NPrint = M;
% choose focussing (K=1) or defocussing (K=-1):
K = 1;
% NOTE: the signs of K are different here and in the code than in the class
% notes.  I didn't have a chance to go through and make everything 
% self-consistent...
dx = 2*pi/N;
x = 0:dx:2*pi;
x = x';
A = 4;
B = 3; 
C = 0;
u0 = A*exp(i*(B*x-t_0*(-K*A^2+B^2)+C));
% here's the exact solution for comparison purposes.  Because of the 
% Benjamin-Fier instability, if I'm going to look for convergence, I need
% to look at a time before the BF instability kicks in.  So I'll look
% at time t=.25 .  I've chosen this time by looking at the power spectrum
% of u11 (see below) and things are okay at this time.
uf = A*exp(i*(B*x-(0.25)*(-K*A^2+B^2)+C));

[u1,x,t1,kk,amp1,M1,H1] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u2,x,t2,kk,amp2,M2,H2] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u3,x,t3,kk,amp3,M3,H3] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u4,x,t4,kk,amp4,M4,H4] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u5,x,t5,kk,amp5,M5,H5] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u6,x,t6,kk,amp6,M6,H6] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u7,x,t7,kk,amp7,M7,H7] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u8,x,t8,kk,amp8,M8,H8] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u9,x,t9,kk,amp9,M9,H9] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u10,x,t10,kk,amp10,M10,H10] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);
dt = dt/2;
NPrint = 2*NPrint;
[u11,x,t11,kk,amp11,M11,H11] = split_schroedinger1(u0,K,t_0,t_f,dt,NPrint,N);

% to watch movies of instabilities
for j=0:20, plot(kk,amp2(:,1+j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp3(:,1+2*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp4(:,1+4*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp5(:,1+8*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp6(:,1+16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp7(:,1+2*16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp8(:,1+4*16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp9(:,1+8*16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp10(:,1+16*16*j)), figure(1), pause(1), end
for j=0:20, plot(kk,amp11(:,1+32*16*j)), figure(1), pause(1), end

% compare the computed solution to the exact solution for evidence
% that the schemes are convergent.
%
% I want to compare at time t = 0.25.  The u1 solution doesn't get
% this time but the subsequent ones do.
% err(1) = max(abs(u1(:,length(t1))-uf));
clear err
i1 = find( t2 == 1/4);
err(1) = max(abs(u2(:,i1)-uf));
i1 = find( t3 == 1/4);
err(2) = max(abs(u3(:,i1)-uf));
i1 = find( t4 == 1/4);
err(3) = max(abs(u4(:,i1)-uf));
i1 = find( t5 == 1/4);
err(4) = max(abs(u5(:,i1)-uf));
i1 = find( t6 == 1/4);
err(5) = max(abs(u6(:,i1)-uf));
i1 = find( t7 == 1/4);
err(6) = max(abs(u7(:,i1)-uf));
i1 = find( t8 == 1/4);
err(7) = max(abs(u8(:,i1)-uf));
i1 = find( t9 == 1/4);
err(8) = max(abs(u9(:,i1)-uf));
i1 = find( t10 == 1/4);
err(9) = max(abs(u10(:,i1)-uf));
i1 = find( t11 == 1/4);
err(10) = max(abs(u11(:,i1)-uf))

% for initial data
%    u0 = A*exp(i*(B*x+0*(-K*A^2+B^2)+C));
% make plots of the Benjamin-Feir instability (focussing case)
% figure(1)
% clf
% plot(kk,amp5(:,1),'o-');
% hold on
% plot(kk,amp5(:,11),'o-');
% plot(kk,amp5(:,21),'o-');
% plot(kk,amp5(:,31),'o-');
% print -dps split_schroedinger_focussing_unstable_spec.ps
%
% figure(1)
% clf
% plot(kk,amp6(:,1),'o-');
% hold on
% plot(kk,amp6(:,21),'o-');
% plot(kk,amp6(:,41),'o-');
% plot(kk,amp6(:,61),'o-');
% print -dps split_schroedinger_focussing_stable_spec.ps
% 
% figure(2)
% clf
% hold on
% plot(kk,amp6(:,81),'o-');
% plot(kk,amp7(:,161),'-');
% plot(kk,amp8(:,321),'-');
% plot(kk,amp9(:,641),'o-');
% print -dps split_schroedinger_ben_feir_spec.ps


% Difference between Particle Number at final time and initial time
abs(M1(length(t1))-M1(1))
abs(M2(length(t2))-M2(1))
abs(M3(length(t3))-M3(1))
abs(M4(length(t4))-M4(1))
abs(M5(length(t5))-M5(1))
abs(M6(length(t6))-M6(1))
abs(M7(length(t7))-M7(1))

% Difference between Hamiltonian at final time and initial time
abs(H1(length(t1))-H1(1))
abs(H2(length(t2))-H2(1))
abs(H3(length(t3))-H3(1))
abs(H4(length(t4))-H4(1))
abs(H5(length(t5))-H5(1))
abs(H6(length(t6))-H6(1))
abs(H7(length(t7))-H7(1))
abs(H8(length(t8))-H8(1))
abs(H9(length(t9))-H9(1))
abs(H10(length(t10))-H10(1))
abs(H11(length(t11))-H11(1))

% the ratio of those differences...
abs(H1(length(t1))-H1(1))/abs(H2(length(t2))-H2(1))
abs(H2(length(t2))-H2(1))/abs(H3(length(t3))-H3(1))
abs(H3(length(t3))-H3(1))/abs(H4(length(t4))-H4(1))
abs(H4(length(t4))-H4(1))/abs(H5(length(t5))-H5(1))
abs(H5(length(t5))-H5(1))/abs(H6(length(t6))-H6(1))
abs(H6(length(t6))-H6(1))/abs(H7(length(t7))-H7(1))
abs(H7(length(t7))-H7(1))/abs(H8(length(t8))-H8(1))
abs(H8(length(t8))-H8(1))/abs(H9(length(t9))-H9(1))
abs(H9(length(t9))-H9(1))/abs(H10(length(t10))-H10(1))
abs(H10(length(t10))-H10(1))/abs(H11(length(t11))-H11(1))

