% Plot the soliton and its Fourier transform

format long;
beta = 1;

L = 25;
N = 200;     % N must be even here
h = L/N;
absc = -L/2:h:(L/2-h);
        % Invariant: length(absc) == N
[X,Y,Z] = meshgrid(absc,absc,absc);

xi = 1:N;
xifreq0 = floor(N/2+1);
xi = xi - xifreq0;
[XI1,XI2,XI3] = meshgrid(xi);
absXIsquared = XI1.^2 + XI2.^2 + XI3.^2;

[soliton,Petv_consts,EL_tests1,EL_tests2] = Petviashvili(beta, N, L);
s_x = squeeze(soliton((N/2+1):end,(N/2+1),(N/2+1)));
h = L/N;
x = 0:h:(L/2-h);

F_soliton = abs(fftshift(fftn(soliton)));
s_xi = squeeze(F_soliton((N/2+1):end,(N/2+1),(N/2+1)));
s_xi = (L/N)^3 * s_xi;
H = 1/L;
xi = 2*pi*(0:H:((N/L)/2-H));

figure(1); clf;
semilogy(x,s_x,'k-');

figure(2); clf;
semilogy(xi,s_xi,'k-');

figure(3); clf;
semilogy(1-Petv_consts,'k.'); hold on;
semilogy(EL_tests1,'k-.'); hold on;
semilogy(EL_tests2,'k-');
