function vvg = WaveGPU(N, Nsteps, hTopAxes, hBottomAxes, hStartStopBtn, ...
    hIterationText, hElapsedTimeText, plotStep, ticT)
%% Solving 2nd Order Wave Equation Using Spectral Methods
% This example solves a 2nd order wave equation: utt = uxx + uyy, with u =
% 0 on the boundaries. It uses a 2nd order central finite difference in
% time and a Chebysehv spectral method in space (using FFT).
%
% The code has been modified from an example in Spectral Methods in MATLAB
% by Trefethen, Lloyd N.
%
% Copyright 2015 The MathWorks, Inc.

if nargin == 2
    guiMode = false;
else
    guiMode = true;
    on_state = get(hStartStopBtn,'Max');
end

% Points in X and Y
x = cos(pi*(0:N)/N); % using Chebyshev points

% Send x to the GPU
x = gpuArray(x);

y = x';

% Calculating time step
dt = 6/N^2;

% Setting up grid
[xx,yy] = meshgrid(x,y);

% Calculate initial values
vv = exp(-40*((xx-.4).^2 + yy.^2));
vvold = vv;

ii = 2:N;
index1 = 1i*[0:N-1 0 1-N:-1];
index2 = -[0:N 1-N:-1].^2;

% Sending data to the GPU
dt = gpuArray(dt);
index1 = gpuArray(index1);
index2 = gpuArray(index2);

% Weights used for spectral differentiation via FFT
W1T = repmat(index1,N-1,1);
W2T = repmat(index2,N-1,1);
W3T = repmat(index1.',1,N-1);
W4T = repmat(index2.',1,N-1);

WuxxT1 = repmat((1./(1-x(ii).^2)),N-1,1);
WuxxT2 = repmat(x(ii)./(1-x(ii).^2).^(3/2),N-1,1);
WuyyT1 = repmat(1./(1-y(ii).^2),1,N-1);
WuyyT2 = repmat(y(ii)./(1-y(ii).^2).^(3/2),1,N-1);

% Start time-stepping
n = 1;

while n <= Nsteps
    
    V = [vv(ii,:) vv(ii,N:-1:2)];
    U = real(fft(V.')).';
    
    W1test = (U.*W1T).';
    W2test = (U.*W2T).';
    W1 = (real(ifft(W1test))).';
    W2 = (real(ifft(W2test))).';
    
    % Calculating 2nd derivative in x
    uxx(ii,ii) = W2(:,ii).* WuxxT1 - W1(:,ii).*WuxxT2;
    uxx([1,N+1],[1,N+1]) = 0;
    
    V = [vv(:,ii); vv((N:-1:2),ii)];
    U = real(fft(V));
    
    W1 = real(ifft(U.*W3T));
    W2 = real(ifft(U.*W4T));
    
    % Calculating 2nd derivative in y
    uyy(ii,ii) = W2(ii,:).* WuyyT1 - W1(ii,:).*WuyyT2;
    uyy([1,N+1],[1,N+1]) = 0;
    
    % Computing new value using 2nd order central finite difference in
    % time
    vvnew = 2*vv - vvold + dt*dt*(uxx+uyy);
    vvold = vv; 
    vv = vvnew;
    
    if guiMode
        if ~isequal(get(hStartStopBtn, 'Value'), on_state)
            break
        end
        if n == 1 || mod(n, plotStep) == 0
            updateUI(gather(xx), gather(yy), gather(vv), gather(vvold), n, N, hTopAxes, ...
                hBottomAxes, hIterationText, hElapsedTimeText, ticT);
        end
    end
    n = n + 1;
end

% Gather vvg back from GPU memory when done
vvg = gather(vv);
