% --------------------------------------------
%
% speImp.m
%
% solves shortpulse equation as in
% Sakovich paper
%
% PDE:      u_{xt} = u + c*(u^3)_{tt} 
%
%
% Initial:  u(x=0,t) = sakovich soliton
%           (found from analytical formula
%            using iteration)
% B.C:      vanishing for |x|->infinity
%
% written by T.Schaefer (June 2009)
%
%---------------------------------------------

function [x,u,uStart,uAna]=speImp()

% Set up parameters

xMin =  -256;            % range in x
xMax =   256;

nMin = -2048;            % grid
nMax = -nMin-1;
nPoints = nMax-nMin+1;

dx = (xMax-xMin)/(nPoints);
x = linspace(xMin,xMax-dx,nPoints);
dom = 2*pi/(xMax-xMin);

ep = 0.2;

dt= 25.625*ep/410;

tStart = 0.0;

stepNum = 410; 
itNum = 3;

tEnd = dt*stepNum;

% Set up random vector

% PDE parameters

c = 1/6;

% Set up vector for FFT

om  = dom*[0:nPoints/2,-nPoints/2+1:-1]; % with zero frequency
omC = dom*[1:nPoints/2,-nPoints/2+1:-1]; % without zero frequency

% Initial condition

u = sakAna(x,0,0.3);
uStart = u;

% Analytical Solutions

uAna = sakAna(x,tEnd,0.3);

% Create output of initial condition

uStart = u;             

% Start Integration, we work in Fourier space

u = fft(u);

for stepCtr=1:stepNum

    uIterate = u;
    
    for itCtr=1:itNum
        
    tCurr=stepCtr*dt;
        
    % Linear step
        
    help = 0.5*(uIterate+u);

    % We have to exclude the zero-frequency mode

    helpC = help(2:nPoints);
    helpC = -i./omC.*helpC;
    helpL = [0,helpC];
        
    % Nonlinear step
        
    helpNL = ifft(help);
    helpNL = helpNL.^3;
    helpNL = fft(helpNL);
    helpNL = i*om.*helpNL;

    uIterate = u + dt*(helpL+c*helpNL);
    
    end
    
    u = uIterate;
    
     %Filter in Fourier domain
    
     %kCut = 30;
     %kCut2 = 1000;
     %uOld(1:kCut)=0;
     %uOld(nPoints-kCut:nPoints)=0;
     %uOld(nPoints/2-kCut2:nPoints/2+kCut2)=0;
     
     %u(1:kCut)=0;
     %u(nPoints-kCut:nPoints)=0;
     %u(nPoints/2-kCut2:nPoints/2+kCut2)=0;
end         

u = real(ifft(u));


% function for Sakovich soliton

function u = sakAna(x,t,m)

% computes the analytical solution (Sakovich-soliton)
% at time t for a vector x and a soliton parameter m
% using iteration

n = sqrt(1-m^2);

% Iteration for y

itNum = 100;
y = x;

for k=1:itNum
    B = m*sin(2*n*(y-t))-n*sinh(2*m*(y+t));
    C = m^2*sin(n*(y-t)).^2+n^2*cosh(m*(y+t)).^2;
    y = x - 2*n*m*B./C;
end

A = m*sin(n*(y-t)).*sinh(m*(y+t))+n*cos(n*(y-t)).*cosh(m*(y+t));
C = m^2*sin(n*(y-t)).^2+n^2*cosh(m*(y+t)).^2;
u = 4*m*n*A./C;
