%  w02wwave.m -- djm -- 20 jan 2006
%  - modified considerably for scattering experiments
%
%Solving the wave equation
%     u_tt = c^2*u_xx 
% on a infinite domain,
%   with u(x,0) = f(x),
%        u'(x,0) = 0,
%        and c is discontinous function
%
% Numerical method:
%     Centered difference in time, and centered difference in space
%
% Wilson Au
% http://www.math.sfu.ca/~wilsona/
% Jan 2006
%
% 
%   Parameters:        line26 - line31
%   Initial condition: line47 - line49
%   Plotting Setup:    line70 - line72
%

tic;  choice = 2;
%close all

% Parameters
N = 8*512;                                   % Number of point
L1 = -16*pi;  L2 = 16*pi;                  % Domain (i.e. [L1 L2])

% dx
dx = (L2-L1)/(N-1);                        % dx
x = L1:dx:L2;                              % spatial grid

% Wavespeed
switch choice
    case{1}
        c_m = 1;                                   % c-
        c_p = 2;                                   % c+
        jump  = 0;                                 % location of discontinuity
        jump2 = pi;                                % location of discontinuity

        C = (c_m-c_p)*(sign(x-jump) == -1)+c_p;    % Wavespeed function
        A_T = 2*c_p/(c_m+c_p);                     % Amplitude of transmitted wave
        A_R = (c_p-c_m)/(c_m+c_p);                 % Amplitude of reflected wave
    case{2}
        c_m = 1;                                   % c background
        c_p = 2;                                   % c max
        jump  = 0;                                 % location of discontinuity
        jump2 = pi;                                % location of discontinuity
        C = c_m./(((x>jump)&(x<jump2)).*sin(x) *(c_m/c_p - 1) + 1); 
        A_T = 1.0;  A_R = 1.0;
    case{3}
        c_m = 1;                                   % c background
        c_p = 2;                                   % c max
        jump  = 0;                                 % location of discontinuity
        jump2 = pi;                                % location of discontinuity
        C = c_m./((sech(2*(x-pi/2)).^2) *(c_m/c_p - 1) + 1);
        A_T = 1.0;  A_R = 1.0;
end

% dt
dt = sqrt(0.99*dx^2/max(C)^2);             % dt
tmax = 30;                                 % End time
numIter = ceil(tmax/dt);                   % Number of iterations

% CFL condition: c^2*dt^2/dx^2 < 1
s = C.^2*dt^2/dx^2;

% Initial condition
switch choice
    case{1}
        center = -20;                              % Center of the peak
        A = 1.5;                                   % Amplitude of the peak
        one_sided = 1;                             % one_sided = 1 for 1-sided wave
        % one_sided = 0 for 2-sided wave
        two_sided = 1-one_sided;
        u0 = A*sech(0.5*(x-center));

        if one_sided == 1
            % "One sided" wave
            u1 = A*sech(0.5*(x-center-c_m*dt)); Amp = max(u1);
        elseif two_sided == 1
            % "Two sided" wave
            u1 = 0.5*s.*(u0([2:N N])+u0([1 1:end-1])) + (1-s).*u0;
        end
    case{2,3}
%        omega = 14*pi/8;  pc = -24;  pw = 8;
%        omega = 13/6;
        pc  = -4*pi;  pw = 3*pi;
        ue  = (1+  cos(1*(x-pc)))/2.*(abs(x-pc      )<=pw);
        uep = ( -1*sin(1*(x-pc)))/2.*(abs(x-pc      )<=pw);

        ue  = max(ue,abs((x-pc))<=pw-pi);
        uep = uep.*(abs(x-pc)>pw-pi);
        
        u0  = ue.*cos(omega*(x-pc+pw)+pi/2);
        u0p = -C .* (uep.*cos(omega*(x-pc+pw)+pi/2) ...
              - omega*ue.*sin(omega*(x-pc+pw)+pi/2));
        
        u1  = u0 + dt*u0p + (0.5*(dt/dx)^2) * (C.^2).*[0 diff(u0,2) 0];
end
Amp = max(u1)/2;
  
% Initialize
old_u = u0; cur_u = u1;
new_u = u0;
u_data = [cur_u]; t_data = [0];

% Plotting Setup
p1 = 40;                       % Number of animation frame (See Figure 2)
p2 = 20;                       % Number of waterfall line (See Figure 4)
p3 = 5;                        % Number of profile plot (See Figure 3)
numPlot = ceil(numIter/p1);    % Number of iteration between each output
numFrame1 = ceil(numIter/p2);  % Number of iteration between each waterfall
numFrame2 = ceil(p2/p3);       % Number of iteration between each profile
maxA = Amp*max(A_R,A_T); minA = Amp*min(A_R,A_T);
top = min(min(u0),minA); bottom = max(max(u0),maxA);
Y1 = min(min(u0),minA)-0.1*abs(top-bottom);
Y2 = max(max(u0),maxA)+0.1*abs(top-bottom);

% Plotting
%figure(1), clf;
% plot(x,C,'k-');
% text((jump+L1)/2,(c_p+c_m)/2,...
%     ['\bf c_- = ',num2str(c_m)],...
%     'HorizontalAlignment','center',...
%     'FontSize',20);
% text((L2+jump)/2,(c_p+c_m)/2,...
%     ['\bf c_+ = ',num2str(c_p)],...
%     'HorizontalAlignment','center',...
%     'FontSize',20);
% xlabel(['\bf x']);
% title(['\bf Discontinuous Wavespeed']);
% axis([L1 L2 min(c_m,c_p)-0.1*abs(c_p-c_m)-((c_p-c_m)==0) ...
%     max(c_m,c_p)+0.1*abs(c_p-c_m)+((c_p-c_m)==0)])

figure(2), clf;
plot(x,u0,'r-',[jump  jump ],[-1 6],'g:',[jump2 jump2],[-1 6],'g:');
text((jump+L1)/2,5.5,...
    ['\bf c_- = ',num2str(c_m)],...
    'HorizontalAlignment','center',...
    'FontSize',20);
text((L2+jump)/2,5.5,...
    ['\bf c_+ = ',num2str(c_p)],...
    'HorizontalAlignment','center',...
    'FontSize',20);
xlabel(['\bf x']); ylabel(['\bf u(x,t)']);
title(['\bf Initial Condition']);
axis([L1 L2 Y1 Y2])

t1 = toc;
%pause
tic;

scnsize = get(0,'ScreenSize');
if ispc == 1
    h = waitbar(0,'\bf please wait...', ...
        'Position',[scnsize(3)/10 scnsize(4)/10 270 56.25]);
else
    h = waitbar(0,'\bf please wait...', ...
        'Position',[scnsize(3)/10 scnsize(4)/10 347 56.25]);
end

% Main Loop
for k = 1:numIter
    % Calculate
    s = C.^2*dt^2/dx^2;
    new_u = s.*(cur_u([2:N N])+cur_u([1 1:end-1])) + 2*(1-s).*cur_u - old_u;


    if mod(k,numPlot) == 0 | k == numIter
        % Plot
        figure(2);
%        plot(x,new_u,'r-',x,u0,'k-',[jump jump],[-1 6],'g:', ...
        plot(x,new_u,'r-', ...
            [jump jump],[-1 6],'g:', ...
            [jump2 jump2],[-1 6],'g:', ...
            [L1 L2],Amp*[A_T A_T],'b:', ...
            [L1 L2],Amp*[A_R A_R],'b:');
        text((jump+L1)/2,5.5,...
             ['\bf c_- = ',num2str(c_m)],...
             'HorizontalAlignment','center',...
             'FontSize',20);
        text((L2+jump)/2,5.5,...
             ['\bf c_+ = ',num2str(c_p)],...
             'HorizontalAlignment','center',...
             'FontSize',20);
        xlabel(['\bf x']); ylabel(['\bf u(x,t)']);
        title(['\bf t = ',num2str(k*dt)]);
        axis([L1 L2 Y1 Y2])
        drawnow;
        
    end

    % Store data
    if mod(k,numFrame1) == 0
        u_data = [new_u; u_data];
        t_data = [k*dt; t_data];   
    end
    
    % Update
    old_u = cur_u;
    cur_u = new_u;
    
    waitbar(k/numIter,h);
end
close(h);

% figure(3), clf;
% hold on
% [m,n] = size(u_data);
% plot(x,u_data(m,:),'k-');
% plot([jump jump],[Y1 Y2],'g:','LineWidth',2);
% for k = 1:numFrame2:m-numFrame2
%     plot(x,u_data(k,:),'r-');
% end
% text((jump+L1)/2,Y2-0.1*(Y2-Y1),...
%      ['\bf c_- = ',num2str(c_m)],...
%      'HorizontalAlignment','center',... 
%      'FontSize',20);
% text((L2+jump)/2,Y2-0.1*(Y2-Y1),...
%      ['\bf c_+ = ',num2str(c_p)],...
%      'HorizontalAlignment','center',...
%      'FontSize',20);
% xlabel('\bf x'), ylabel('\bf u(x,t)');        
% title( ['\bf Solution profiles'] );
% axis([L1 L2 Y1 Y2])
% grid on;
% hold off

% figure(4), clf;
% waterfall(x,t_data,u_data), colormap(1e-6*[1 1 1]); view(-20,25)      
% xlabel('\bf x'), ylabel('\bf t'), zlabel('\bf u(x,t)');
% title(['\bf Profile']);
% axis([L1 L2 0 tmax Y1 Y2]) 
% grid off;
% set(gca,'ytick',[0 tmax]);
% pbaspect([1 1 .13])

t2 = toc;
disp(['Elapsed time is ',num2str(t1+t2),' seconds.']);

backscatter0 = sum((   u0(3:N/2+1)-   u0(1:N/2-1)).^2)/(dx);
backscatter  = sum((new_u(3:N/2+1)-new_u(1:N/2-1)).^2)/(dx);
format long e;  backscatter/backscatter0
data = [data ; omega backscatter/backscatter0 max(abs(new_u(21*N/64:25*N/64)))]
format short