diff --git a/README.md b/README.md
index 7d2c196..8deabd4 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,6 @@
# MATLAB Discrete Trigonometric Transform Library
+[](https://www.mathworks.com/matlabcentral/fileexchange/75071-matlab-discrete-trigonometric-transform-library)
+
## Author
This library is written by Bradley Treeby, University College London. Contact: b.treeby@ucl.ac.uk
@@ -31,4 +33,17 @@ This program is free software: you can redistribute it and/or modify it under th
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
-You should have received a copy of the GNU General Public License along with this program. If not, see .
\ No newline at end of file
+You should have received a copy of the GNU General Public License along with this program. If not, see .
+
+## Change Log
+
+* v1.1 (21 April 2020):
+ * Fixed bug in `gradientDtt1D` (missing `numDim` function)
+ * Added `align_output` option to `gradientDtt1D`
+ * Updated `example_wave_eq_pstd_1D` to enforce correct boundary position
+ * Added `example_wave_eq_pstd_1D_non_reflecting` to simulate general boundary conditions
+ * Added `example_wsws_gradient` to illustrate grid shifting with and without `align_output`
+* Added `example_wave_eq_pstd_2D_neumann` and `example_wave_eq_pstd_2D_dirichlet`
+
+* v1.0 (17 April 2020):
+ * Initial release
\ No newline at end of file
diff --git a/examples/example_gradient_calculation.m b/examples/example_gradient_calculation.m
index c7a0ec8..66b5439 100644
--- a/examples/example_gradient_calculation.m
+++ b/examples/example_gradient_calculation.m
@@ -1,9 +1,9 @@
% DESCRIPTION:
% This example script creates different periodic sequences of known
% symmetry, and then numerically calculates the gradients using
-% gradientDTT1D. The numerical gradient is then compared with the known
+% gradientDtt1D. The numerical gradient is then compared with the known
% analytical gradient. A test sequence is created for all supported DTT
-% symmetries. The grid shifting functionality of gradientDTT1D can also
+% symmetries. The grid shifting functionality of gradientDtt1D can also
% be tested by changing the value of shift.
%
% ABOUT:
@@ -13,7 +13,7 @@
%
% Copyright (C) 2013-2020 Bradley Treeby
%
-% See also dtt1D, gradientDTT1D
+% See also dtt1D, gradientDtt1D
% This program is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
@@ -68,7 +68,7 @@
end
% compute gradient
- dfdx_num = gradientDTT1D(f, dx, test_case, shift);
+ dfdx_num = gradientDtt1D(f, dx, test_case, shift);
dfdx_analytical = -sin(x_out);
% plot functions
@@ -114,7 +114,7 @@
end
% compute gradient
- dfdx_num = gradientDTT1D(f, dx, test_case, shift);
+ dfdx_num = gradientDtt1D(f, dx, test_case, shift);
dfdx_analytical = -sin(x_out);
% plot functions
@@ -160,7 +160,7 @@
end
% compute gradient
- dfdx_num = gradientDTT1D(f, dx, test_case, shift);
+ dfdx_num = gradientDtt1D(f, dx, test_case, shift);
dfdx_analytical = -sin(x_out);
% plot functions
@@ -206,7 +206,7 @@
end
% compute gradient
- dfdx_num = gradientDTT1D(f, dx, test_case, shift);
+ dfdx_num = gradientDtt1D(f, dx, test_case, shift);
dfdx_analytical = -sin(x_out);
% plot functions
@@ -252,7 +252,7 @@
end
% compute gradient
- dfdx_num = gradientDTT1D(f, dx, test_case, shift);
+ dfdx_num = gradientDtt1D(f, dx, test_case, shift);
dfdx_analytical = cos(x_out);
% plot functions
@@ -298,7 +298,7 @@
end
% compute gradient
- dfdx_num = gradientDTT1D(f, dx, test_case, shift);
+ dfdx_num = gradientDtt1D(f, dx, test_case, shift);
dfdx_analytical = cos(x_out);
% plot functions
@@ -344,7 +344,7 @@
end
% compute gradient
- dfdx_num = gradientDTT1D(f, dx, test_case, shift);
+ dfdx_num = gradientDtt1D(f, dx, test_case, shift);
dfdx_analytical = cos(x_out);
% plot functions
@@ -390,7 +390,7 @@
end
% compute gradient
- dfdx_num = gradientDTT1D(f, dx, test_case, shift);
+ dfdx_num = gradientDtt1D(f, dx, test_case, shift);
dfdx_analytical = cos(x_out);
% plot functions
diff --git a/examples/example_transform_wsws_sequence.m b/examples/example_transform_wsws_sequence.m
index 5735e12..827b82e 100644
--- a/examples/example_transform_wsws_sequence.m
+++ b/examples/example_transform_wsws_sequence.m
@@ -16,9 +16,9 @@
% ABOUT:
% author - Bradley Treeby
% date - 31 March 2020
-% last update - 2 April 2020
+% last update - 20 April 2020
%
-% See also dtt1D, gradientDTT1D
+% See also dtt1D, gradientDtt1D
% This program is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
diff --git a/examples/example_wave_eq_pstd.m b/examples/example_wave_eq_pstd_1D.m
similarity index 73%
rename from examples/example_wave_eq_pstd.m
rename to examples/example_wave_eq_pstd_1D.m
index e2beabd..55ef445 100644
--- a/examples/example_wave_eq_pstd.m
+++ b/examples/example_wave_eq_pstd_1D.m
@@ -1,16 +1,22 @@
% DESCRIPTION:
% This example script solves the 1D wave equation (written as two
% coupled first-order equations) using a DTT-based PSTD method subject
-% to different combinations of Dirichlet and Neumann boundary
-% conditions at each end of the domain. The spatial gradients are
-% calculated using discrete cosine transforms (DCTs) and discrete sine
-% transforms (DST) chosen to enfore the required boundary condition.
-% Time integration is performed using a first-order accurate foward
-% difference scheme. A classical Fourier PSTD method is also used with
-% periodic bounday conditions for comparison. Further details are given
-% in [1].
+% to four different combinations of Dirichlet and Neumann boundary
+% conditions at each end of the domain.
%
-% This script reproduces Figure 5 of [1].
+% The spatial gradients are calculated using discrete cosine transforms
+% (DCTs) and discrete sine transforms (DST) chosen to enfore the
+% required boundary condition. Time integration is performed using a
+% first-order accurate foward difference scheme. A classical Fourier
+% PSTD method is also used with periodic bounday conditions for
+% comparison.
+%
+% For the DTT-based simulations, the position of the boundary is
+% assumed to be at the first and last grid points. To account for the
+% different implied symmetries for the different transform types, the
+% calculations are performed on grids of different sizes.
+%
+% Further details are given in [1].
%
% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral
% time-domain (PSTD) methods for the wave equation: Realising boundary
@@ -19,11 +25,11 @@
% ABOUT:
% author - Bradley Treeby
% date - 31 March 2020
-% last update - 1 April 2020
+% last update - 20 April 2020
%
% Copyright (C) 2020 Bradley Treeby
%
-% See also dtt1D, gradientDTT1D
+% See also dtt1D, gradientDtt1D
% This program is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
@@ -44,11 +50,12 @@
% DEFINE LITERALS
% =========================================================================
-% shift variables used by gradientDTT1D
-shift_pos = 1;
-shift_neg = 2;
+% shift and align variables used by gradientDtt1D (shift = 1 and shift = 2
+% return the same output if align_output = false)
+shift_output = 1;
+align_output = false;
-% transform types used by gradientDTT1D
+% transform types used by gradientDtt1D
DCT1 = 1; % WSWS
DCT2 = 2; % HSHS
DCT3 = 3; % WSWA
@@ -118,6 +125,7 @@
% select figure
figure(plot_fig);
+ hold on;
% preallocate matrix to store data for waterfall plot
p_waterfall = zeros(waterfall_snapshots, Nx);
@@ -132,57 +140,86 @@
% Neumann-Neumann / WSWS
T1 = DCT1;
T2 = DST2;
+
+ % set indices for representative sample
+ ind1 = 1;
+ ind2 = Nx;
+
+ % set title
waterfall_title = 'Neumann-Neumann (WSWS)';
+ title(waterfall_title);
case 2
% Neumann-Dirichlet / WSWA
T1 = DCT3;
T2 = DST4;
+
+ % set indices for representative sample
+ ind1 = 1;
+ ind2 = Nx - 1;
+
+ % set title
waterfall_title = 'Neumann-Dirichlet (WSWA)';
+ title(waterfall_title);
case 3
% Dirichlet-Neumann / WAWS
T1 = DST3;
T2 = DCT4;
+
+ % set indices for representative sample
+ ind1 = 2;
+ ind2 = Nx;
+
+ % set title
waterfall_title = 'Dirichlet-Neumann (WAWS)';
+ title(waterfall_title);
case 4
- % Neumann-Neumann / WSWS
+ % Dirichlet-Dirichlet / WAWA
T1 = DST1;
T2 = DCT2;
+
+ % set indices for representative sample
+ ind1 = 2;
+ ind2 = Nx - 1;
+
+ % set title
waterfall_title = 'Dirichlet-Dirichlet (WAWA)';
+ title(waterfall_title);
end
- % assign initial condition for the pressure
- p = p0;
+ % assign initial condition for p, just selecting representative sample
+ p = p0(ind1:ind2);
% run the model backwards for dt/2 to calculate the initial condition
% for u, to take into account the time staggering between p and u
- u = u0 + (dt / 2) / rho0 * gradientDTT1D(p0, dx, T1, shift_pos);
+ u = u0(1:end - 1) + (dt / 2) / rho0 * gradientDtt1D(p, dx, T1, shift_output, align_output);
% calculate pressure in a loop
for time_ind = 1:Nt
% update u
- u = u - dt / rho0 * gradientDTT1D(p, dx, T1, shift_pos);
+ u = u - dt / rho0 * gradientDtt1D(p, dx, T1, shift_output, align_output);
% update p
- p = p - dt * rho0 * c0^2 * gradientDTT1D(u, dx, T2, shift_neg);
+ p = p - dt * rho0 * c0^2 * gradientDtt1D(u, dx, T2, shift_output, align_output);
% plot pressure field
if ~rem(time_ind, plot_freq)
- plot(x, p, 'k-');
+ cla;
+ plot(x(ind1:ind2), p, 'k-');
set(gca, 'YLim', [-1, 1]);
drawnow;
end
% save the pressure field
if ~rem(time_ind, Nt / waterfall_snapshots)
- p_waterfall(waterfall_ind, :) = p;
+ p_waterfall(waterfall_ind, ind1:ind2) = p;
waterfall_ind = waterfall_ind + 1;
end
@@ -209,6 +246,10 @@
% select figure
figure(plot_fig);
+% set title
+waterfall_title = 'Periodic';
+title(waterfall_title);
+
% preallocate matrix to store data for waterfall plot
p_waterfall = zeros(waterfall_snapshots, Nx);
@@ -242,6 +283,7 @@
% plot pressure field
if ~rem(time_ind, plot_freq)
+ cla;
plot(x, p, 'k-');
set(gca, 'YLim', [-1, 1]);
drawnow;
@@ -264,7 +306,7 @@
ylabel('time');
zlabel('pressure');
xlabel('position');
-title('Periodic');
+title(waterfall_title);
grid off
% close animation figure
diff --git a/examples/example_wave_eq_pstd_1D_non_reflecting.m b/examples/example_wave_eq_pstd_1D_non_reflecting.m
new file mode 100644
index 0000000..d94b433
--- /dev/null
+++ b/examples/example_wave_eq_pstd_1D_non_reflecting.m
@@ -0,0 +1,284 @@
+% DESCRIPTION:
+% This example script solves the 1D wave equation (written as two
+% coupled first-order equations) using a DTT-based PSTD method subject
+% to four different combinations of Dirichlet and Neumann boundary
+% conditions at each end of the domain.
+%
+% The solutions for different boundary conditions are then summed to
+% give different reflection coefficients at each end of the domain,
+% including non-reflecting boundaries. This builds on
+% example_wave_eq_pstd.m.
+%
+% Further details are given in [1].
+%
+% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral
+% time-domain (PSTD) methods for the wave equation: Realising boundary
+% conditions with discrete sine and cosine transforms", 2020.
+%
+% ABOUT:
+% author - Bradley Treeby
+% date - 19 April 2020
+% last update - 20 April 2020
+%
+% Copyright (C) 2020 Bradley Treeby
+%
+% See also dtt1D, gradientDtt1D
+
+% This program is free software: you can redistribute it and/or modify it
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your
+% option) any later version.
+%
+% This program is distributed in the hope that it will be useful, but
+% WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+% Public License for more details.
+%
+% You should have received a copy of the GNU General Public License along
+% with this program. If not, see .
+
+clearvars;
+
+%% ========================================================================
+% DEFINE LITERALS
+% =========================================================================
+
+% shift and align variables used by gradientDtt1D (shift = 1 and shift = 2
+% return the same output if align_output = false)
+shift_output = 1;
+align_output = false;
+
+% transform types used by gradientDtt1D
+DCT1 = 1; % WSWS
+DCT2 = 2; % HSHS
+DCT3 = 3; % WSWA
+DCT4 = 4; % HSHA
+DST1 = 5; % WAWA
+DST2 = 6; % HAHA
+DST3 = 7; % WAWS
+DST4 = 8; % HAHS
+
+% number of time snapshots to include in waterfall plot
+waterfall_snapshots = 80;
+
+% number of boundary conditions
+bc_num = 4;
+
+%% ========================================================================
+% DEFINE SIMULATION SETTINGS
+% =========================================================================
+
+% set the grid size
+Nx = 256; % grid size [m]
+dx = 1/Nx; % grid spacing [m]
+
+% set the medium properties
+c0 = 1500; % sound speed [m/s]
+rho0 = 1000; % density [kg/m^3]
+
+% set CFL and number of time steps
+CFL = 0.2;
+Nt = 2720;
+
+%% ========================================================================
+% DEFINE INITIAL CONDITIONS
+% =========================================================================
+
+% spatial grid
+x = (0:Nx - 1) * dx;
+
+% properties of Gaussian initial condition
+width = Nx * dx / 28;
+offset = Nx * dx / 3;
+
+% define initial pressure distribution on regular grid as a Gaussian
+p0 = exp(-((x - offset) / width).^2);
+
+%% ========================================================================
+% RUN SIMULATIONS USING DTT-BASED PSTD METHOD
+% =========================================================================
+
+% calculate the time step
+dt = CFL * dx / c0;
+
+% normalised plot axes
+x_ax = (0:Nx - 1) / (Nx - 1);
+t_ax = (0:waterfall_snapshots - 1) / (waterfall_snapshots - 1);
+
+% preallocate matrix to store simulation data
+p_four_bc = zeros(bc_num, Nx, waterfall_snapshots);
+
+% loop over boundary conditions
+for bc_ind = 1:bc_num
+
+ % initialise index for storing waterfall data
+ waterfall_ind = 1;
+
+ % define which transforms to use
+ switch bc_ind
+ case 1
+
+ % Neumann-Neumann / WSWS
+ T1 = DCT1;
+ T2 = DST2;
+
+ % set indices for representative sample
+ ind1 = 1;
+ ind2 = Nx;
+
+ case 2
+
+ % Neumann-Dirichlet / WSWA
+ T1 = DCT3;
+ T2 = DST4;
+
+ % set indices for representative sample
+ ind1 = 1;
+ ind2 = Nx - 1;
+
+ case 3
+
+ % Dirichlet-Neumann / WAWS
+ T1 = DST3;
+ T2 = DCT4;
+
+ % set indices for representative sample
+ ind1 = 2;
+ ind2 = Nx;
+
+ case 4
+
+ % Dirichlet-Dirichlet / WAWA
+ T1 = DST1;
+ T2 = DCT2;
+
+ % set indices for representative sample
+ ind1 = 2;
+ ind2 = Nx - 1;
+
+ end
+
+ % assign initial condition for p, just selecting representative sample
+ p = p0(ind1:ind2);
+
+ % run the model backwards for dt/2 to calculate the initial condition
+ % for u, to take into account the time staggering between p and u, in
+ % this case setting u0 = 0
+ u = (dt / 2) / rho0 * gradientDtt1D(p, dx, T1, shift_output, align_output);
+
+ % calculate pressure in a loop
+ for time_ind = 1:Nt
+
+ % update u
+ u = u - dt / rho0 * gradientDtt1D(p, dx, T1, shift_output, align_output);
+
+ % update p
+ p = p - dt * rho0 * c0^2 * gradientDtt1D(u, dx, T2, shift_output, align_output);
+
+ % save the pressure field
+ if ~rem(time_ind, floor(Nt / waterfall_snapshots))
+ p_four_bc(bc_ind, ind1:ind2, waterfall_ind) = p;
+ waterfall_ind = waterfall_ind + 1;
+ end
+
+ end
+
+end
+
+%% ========================================================================
+% COMBINE SIMULATIONS WITH DIFFERENT BC
+% =========================================================================
+
+% form matrix of the different possible boundary conditions, where +1
+% corresponds to a positive image source, and -1 to a negative image
+% source, thus the four columns correspond to Neumann-Neumann,
+% Neumann-Dirichlet, Dirichlet-Neumann, and Dirichlet-Dirichlet
+M = [1 1 -1 -1;
+ 1 1 1 1;
+ 1 -1 1 -1;
+ 1 -1 -1 1];
+
+ % open a new figure window
+plot_fig = figure;
+
+% loop over different reflection coefficients
+for ind = 1:2
+
+ % set the desired reflection coefficient
+ switch ind
+ case 1
+
+ % non-reflecting boundaries
+ R_l = 0;
+ R_r = 0;
+
+ case 2
+
+ % partially reflecting boundary
+ R_l = 0;
+ R_r = 0.5;
+
+ end
+
+ % form r vector (this corresponds to the image source amplitudes)
+ r = [R_l, 1, R_r, R_r * R_l].';
+
+ % solve for the weights for each of the pre-calculated fields
+ w = 1/4 * M.' * r %#ok
+
+ % preallocate matrix
+ p_waterfall = zeros(Nx, waterfall_snapshots);
+
+ % loop over time (waterfall) index
+ for waterfall_ind = 1:waterfall_snapshots
+
+ % form the pressure field by summing weighted fields
+ p_waterfall(:, waterfall_ind) = w.' * p_four_bc(:, :, waterfall_ind);
+
+ % plot the fields and the summation
+ figure(plot_fig);
+
+ subplot(5, 1, 1);
+ plot(x_ax, p_four_bc(1, :, waterfall_ind), 'k-');
+ set(gca, 'YLim', [-1, 1]);
+ title('Neumann-Neumann (WSWS)');
+
+ subplot(5, 1, 2);
+ plot(x_ax, p_four_bc(2, :, waterfall_ind), 'k-');
+ set(gca, 'YLim', [-1, 1]);
+ title('Neumann-Dirichlet (WSWA)');
+
+ subplot(5, 1, 3);
+ plot(x_ax, p_four_bc(3, :, waterfall_ind), 'k-');
+ set(gca, 'YLim', [-1, 1]);
+ title('Dirichlet-Neumann (WAWS)');
+
+ subplot(5, 1, 4);
+ plot(x_ax, p_four_bc(4, :, waterfall_ind), 'k-');
+ set(gca, 'YLim', [-1, 1]);
+ title('Dirichlet-Dirichlet (WAWA)');
+
+ subplot(5, 1, 5);
+ plot(x_ax, p_waterfall(:, waterfall_ind), 'k-');
+ set(gca, 'YLim', [-1, 1]);
+ title(['Sum r = [' num2str(r(1)) ', ' num2str(r(2)) ', ' num2str(r(3)) ', ' num2str(r(4)) ']']);
+ drawnow;
+
+ end
+
+ % waterfall plot of evolution of field
+ figure;
+ waterfall(x_ax, t_ax, p_waterfall.');
+ view(10, 70);
+ colormap([0, 0, 0]);
+ set(gca, 'ZLim', [-1, 1], 'FontSize', 12);
+ ylabel('time');
+ zlabel('pressure');
+ xlabel('position');
+ title(['r = [' num2str(r(1)) ', ' num2str(r(2)) ', ' num2str(r(3)) ', ' num2str(r(4)) ']']);
+ grid off;
+
+end
+
+% close animation figure
+close(plot_fig);
\ No newline at end of file
diff --git a/examples/example_wave_eq_pstd_2D_dirichlet.m b/examples/example_wave_eq_pstd_2D_dirichlet.m
new file mode 100644
index 0000000..30d7565
--- /dev/null
+++ b/examples/example_wave_eq_pstd_2D_dirichlet.m
@@ -0,0 +1,217 @@
+% DESCRIPTION:
+% This example script solves the 2D wave equation (written as two
+% coupled first-order equations) using a DTT-based PSTD method subject
+% to Dirichlet boundary conditions on each side of the domain.
+%
+% The spatial gradients are calculated using discrete cosine transforms
+% (DCTs) and discrete sine transforms (DST) chosen to enfore the
+% required boundary condition. Time integration is performed using a
+% first-order accurate foward difference scheme.
+%
+% Due to the symmetry of the DTTs used (WAWA), the position of the
+% boundary is one point outside the first and last grid points in the
+% domain.
+%
+% Further details are given in [1].
+%
+% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral
+% time-domain (PSTD) methods for the wave equation: Realising boundary
+% conditions with discrete sine and cosine transforms", 2020.
+%
+% ABOUT:
+% author - Bradley Treeby
+% date - 21 April 2020
+% last update - 21 April 2020
+%
+% Copyright (C) 2020 Bradley Treeby
+%
+% See also dtt1D, gradientDtt1D
+
+% This program is free software: you can redistribute it and/or modify it
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your
+% option) any later version.
+%
+% This program is distributed in the hope that it will be useful, but
+% WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+% Public License for more details.
+%
+% You should have received a copy of the GNU General Public License along
+% with this program. If not, see .
+
+clearvars;
+
+% =========================================================================
+% DEFINE LITERALS
+% =========================================================================
+
+% transform types used by dtt1D
+DCT1 = 1; % WSWS
+DCT2 = 2; % HSHS
+DCT3 = 3; % WSWA
+DCT4 = 4; % HSHA
+DST1 = 5; % WAWA
+DST2 = 6; % HAHA
+DST3 = 7; % WAWS
+DST4 = 8; % HAHS
+
+% plot frequency for snapshots of the field (time steps)
+plot_freq = 250;
+
+% =========================================================================
+% DEFINE SIMULATION SETTINGS
+% =========================================================================
+
+% set the grid size (assuming grid is square)
+Nx = 256; % grid size [m]
+dx = 1/Nx; % grid spacing [m]
+
+% set the medium properties
+c0 = 1500; % sound speed [m/s]
+rho0 = 1000; % density [kg/m^3]
+
+% set CFL and number of time steps
+CFL = 0.2;
+Nt = 1000;
+
+% =========================================================================
+% DEFINE INITIAL CONDITIONS
+% =========================================================================
+
+% spatial grid
+x = (0:Nx - 1) * dx;
+
+% properties of Gaussian initial condition
+width = Nx * dx / 28;
+offset = Nx * dx / 3;
+
+% define initial pressure distribution on regular grid as a Gaussian
+p = 0.5 * exp(-((x - offset) / width).^2);
+p = p.' * p;
+
+% define initial velocity to be zero
+ux = 0;
+uy = 0;
+
+% =========================================================================
+% DEFINE DTT VARIABLES
+% =========================================================================
+
+% compute the implied period of the input function
+M = 2 .* (Nx + 1);
+
+% define wavenumber vectors for WAWA / DST-I
+k_wawa = 2 .* pi .* (1:1:(M/2 - 1)).' ./ (M .* dx);
+
+% define wavenumber vectors for HSHS / DCT-II
+k_hshs = 2 .* pi .* (0:1:(M/2 - 1)).' ./ (M .* dx);
+
+% =========================================================================
+% RUN SIMULATIONS USING DTT-BASED PSTD METHOD
+% =========================================================================
+
+% calculate the time step
+dt = CFL * dx / c0;
+
+% create custom colour map
+cm = [bone(256); flipud(hot(256))];
+
+% calculate pressure in a loop
+for time_ind = 1:Nt
+
+ % ---------------
+ % calculate dpdx
+ % ---------------
+
+ % forward transform over x
+ dpdx = dtt1D(p, DST1, 1);
+
+ % spectral derivative
+ dpdx = bsxfun(@times, k_wawa, dpdx);
+
+ % WAWA -> HSHS, prepend zero
+ dpdx = [zeros(1, Nx); dpdx]; %#ok
+
+ % inverse transform over x, C2^-1 = C3
+ dpdx = dtt1D(dpdx, DCT3, 1) ./ M;
+
+ % ---------------
+ % calculate dpdy
+ % ---------------
+
+ % forward transform over y
+ dpdy = dtt1D(p, DST1, 2);
+
+ % spectral derivative
+ dpdy = bsxfun(@times, k_wawa.', dpdy);
+
+ % WAWA -> HSHS, prepend zero
+ dpdy = [zeros(Nx, 1), dpdy]; %#ok
+
+ % inverse transform over y, C2^-1 = C3
+ dpdy = dtt1D(dpdy, DCT3, 2) ./ M;
+
+ % ---------------
+ % update u
+ % ---------------
+
+ % update components of particle velocity
+ ux = ux - dt / rho0 * dpdx;
+ uy = uy - dt / rho0 * dpdy;
+
+ % ---------------
+ % calculate duxdx
+ % ---------------
+
+ % forward transform over x
+ duxdx = dtt1D(ux, DCT2, 1);
+
+ % spectral derivative
+ duxdx = bsxfun(@times, -k_hshs, duxdx);
+
+ % HSHS -> WAWA, remove left endpoint
+ duxdx = duxdx(2:end, :);
+
+ % inverse transform over x, S1^-1 = S1
+ duxdx = dtt1D(duxdx, DST1, 1) ./ M;
+
+ % ---------------
+ % calculate duydy
+ % ---------------
+
+ % forward transform over y
+ duydy = dtt1D(uy, DCT2, 2);
+
+ % spectral derivative
+ duydy = bsxfun(@times, -k_hshs.', duydy);
+
+ % HSHS -> WAWA, remove left endpoint
+ duydy = duydy(:, 2:end);
+
+ % inverse transform over y, S1^-1 = S1
+ duydy = dtt1D(duydy, DST1, 2) ./ M;
+
+ % ---------------
+ % update p
+ % ---------------
+
+ % update p
+ p = p - dt * rho0 * c0^2 * (duxdx + duydy);
+
+ % ---------------
+ % plot
+ % ---------------
+
+ % plot snapshots of pressure field
+ if (time_ind == 1) || ~rem(time_ind, plot_freq)
+ figure;
+ imagesc(x, x, p, 0.075 * [-1, 1]);
+ axis image;
+ colormap(cm);
+ box on;
+ set(gca, 'XTick', 0:0.25:1, 'YTick', 0:0.25:1, 'XTickLabel', {}, 'YTickLabel', {});
+ drawnow;
+ end
+
+end
\ No newline at end of file
diff --git a/examples/example_wave_eq_pstd_2D_neumann.m b/examples/example_wave_eq_pstd_2D_neumann.m
new file mode 100644
index 0000000..b6f987f
--- /dev/null
+++ b/examples/example_wave_eq_pstd_2D_neumann.m
@@ -0,0 +1,216 @@
+% DESCRIPTION:
+% This example script solves the 2D wave equation (written as two
+% coupled first-order equations) using a DTT-based PSTD method subject
+% to Neumann boundary conditions on each side of the domain.
+%
+% The spatial gradients are calculated using discrete cosine transforms
+% (DCTs) and discrete sine transforms (DST) chosen to enfore the
+% required boundary condition. Time integration is performed using a
+% first-order accurate foward difference scheme.
+%
+% Due to the symmetry of the DTTs used (WSWS), the position of the
+% boundary is at the first and last grid points in the domain.
+%
+% Further details are given in [1].
+%
+% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral
+% time-domain (PSTD) methods for the wave equation: Realising boundary
+% conditions with discrete sine and cosine transforms", 2020.
+%
+% ABOUT:
+% author - Bradley Treeby
+% date - 20 April 2020
+% last update - 21 April 2020
+%
+% Copyright (C) 2020 Bradley Treeby
+%
+% See also dtt1D, gradientDtt1D
+
+% This program is free software: you can redistribute it and/or modify it
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your
+% option) any later version.
+%
+% This program is distributed in the hope that it will be useful, but
+% WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+% Public License for more details.
+%
+% You should have received a copy of the GNU General Public License along
+% with this program. If not, see .
+
+clearvars;
+
+% =========================================================================
+% DEFINE LITERALS
+% =========================================================================
+
+% transform types used by dtt1D
+DCT1 = 1; % WSWS
+DCT2 = 2; % HSHS
+DCT3 = 3; % WSWA
+DCT4 = 4; % HSHA
+DST1 = 5; % WAWA
+DST2 = 6; % HAHA
+DST3 = 7; % WAWS
+DST4 = 8; % HAHS
+
+% plot frequency for snapshots of the field (time steps)
+plot_freq = 250;
+
+% =========================================================================
+% DEFINE SIMULATION SETTINGS
+% =========================================================================
+
+% set the grid size (assuming grid is square)
+Nx = 256; % grid size [m]
+dx = 1/Nx; % grid spacing [m]
+
+% set the medium properties
+c0 = 1500; % sound speed [m/s]
+rho0 = 1000; % density [kg/m^3]
+
+% set CFL and number of time steps
+CFL = 0.2;
+Nt = 1000;
+
+% =========================================================================
+% DEFINE INITIAL CONDITIONS
+% =========================================================================
+
+% spatial grid
+x = (0:Nx - 1) * dx;
+
+% properties of Gaussian initial condition
+width = Nx * dx / 28;
+offset = Nx * dx / 3;
+
+% define initial pressure distribution on regular grid as a Gaussian
+p = 0.5 * exp(-((x - offset) / width).^2);
+p = p.' * p;
+
+% define initial velocity to be zero
+ux = 0;
+uy = 0;
+
+% =========================================================================
+% DEFINE DTT VARIABLES
+% =========================================================================
+
+% compute the implied period of the input function
+M = 2 .* (Nx - 1);
+
+% define wavenumber vectors for WSWS / DCT-I
+k_wsws = 2 .* pi .* (0:1:M/2).' ./ (M .* dx);
+
+% define wavenumber vectors for HAHA / DST-II
+k_haha = 2 .* pi .* (1:1:M/2).' ./ (M .* dx);
+
+% =========================================================================
+% RUN SIMULATIONS USING DTT-BASED PSTD METHOD
+% =========================================================================
+
+% calculate the time step
+dt = CFL * dx / c0;
+
+% create custom colour map
+cm = [bone(256); flipud(hot(256))];
+
+% calculate pressure in a loop
+for time_ind = 1:Nt
+
+ % ---------------
+ % calculate dpdx
+ % ---------------
+
+ % forward transform over x
+ dpdx = dtt1D(p, DCT1, 1);
+
+ % spectral derivative
+ dpdx = bsxfun(@times, -k_wsws, dpdx);
+
+ % WSWS -> HAHA, remove left endpoint
+ dpdx = dpdx(2:end, :);
+
+ % inverse transform over x, S2^-1 = S3
+ dpdx = dtt1D(dpdx, DST3, 1) ./ M;
+
+ % ---------------
+ % calculate dpdy
+ % ---------------
+
+ % forward transform over y
+ dpdy = dtt1D(p, DCT1, 2);
+
+ % spectral derivative
+ dpdy = bsxfun(@times, -k_wsws.', dpdy);
+
+ % WSWS -> HAHA, remove left endpoint
+ dpdy = dpdy(:, 2:end);
+
+ % inverse transform over y, S2^-1 = S3
+ dpdy = dtt1D(dpdy, DST3, 2) ./ M;
+
+ % ---------------
+ % update u
+ % ---------------
+
+ % update components of particle velocity
+ ux = ux - dt / rho0 * dpdx;
+ uy = uy - dt / rho0 * dpdy;
+
+ % ---------------
+ % calculate duxdx
+ % ---------------
+
+ % forward transform over x
+ duxdx = dtt1D(ux, DST2, 1);
+
+ % spectral derivative
+ duxdx = bsxfun(@times, k_haha, duxdx);
+
+ % HAHA -> WSWS, prepend zero
+ duxdx = [zeros(1, Nx); duxdx]; %#ok
+
+ % inverse transform over x, C1^-1 = C1
+ duxdx = dtt1D(duxdx, DCT1, 1) ./ M;
+
+ % ---------------
+ % calculate duydy
+ % ---------------
+
+ % forward transform over y
+ duydy = dtt1D(uy, DST2, 2);
+
+ % spectral derivative
+ duydy = bsxfun(@times, k_haha.', duydy);
+
+ % HAHA -> WSWS, prepend zero
+ duydy = [zeros(Nx, 1), duydy]; %#ok
+
+ % inverse transform over y, C1^-1 = C1
+ duydy = dtt1D(duydy, DCT1, 2) ./ M;
+
+ % ---------------
+ % update p
+ % ---------------
+
+ % update p
+ p = p - dt * rho0 * c0^2 * (duxdx + duydy);
+
+ % ---------------
+ % plot
+ % ---------------
+
+ % plot snapshots of pressure field
+ if (time_ind == 1) || ~rem(time_ind, plot_freq)
+ figure;
+ imagesc(x, x, p, 0.075 * [-1, 1]);
+ axis image;
+ colormap(cm);
+ box on;
+ set(gca, 'XTick', 0:0.25:1, 'YTick', 0:0.25:1, 'XTickLabel', {}, 'YTickLabel', {});
+ drawnow;
+ end
+
+end
\ No newline at end of file
diff --git a/examples/example_wsws_gradient.m b/examples/example_wsws_gradient.m
new file mode 100644
index 0000000..e7bf125
--- /dev/null
+++ b/examples/example_wsws_gradient.m
@@ -0,0 +1,113 @@
+% DESCRIPTION:
+% This example script shows a WSWS gradient calculated on regular and
+% staggered grids, both with and without alignment.
+%
+% Further details are given in [1].
+%
+% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral
+% time-domain (PSTD) methods for the wave equation: Realising boundary
+% conditions with discrete sine and cosine transforms", 2020.
+%
+% ABOUT:
+% author - Bradley Treeby
+% date - 20 April 2020
+% last update - 20 April 2020
+%
+% Copyright (C) 2020-2020 Bradley Treeby
+%
+% See also dtt1D, gradientDtt1D
+
+% This program is free software: you can redistribute it and/or modify it
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your
+% option) any later version.
+%
+% This program is distributed in the hope that it will be useful, but
+% WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+% Public License for more details.
+%
+% You should have received a copy of the GNU General Public License along
+% with this program. If not, see .
+
+clearvars;
+
+% create function with WSWS symmetry
+Nx = 10;
+dx = pi/(Nx-1);
+x_in = 0:dx:pi;
+f = cos(x_in);
+
+% =========================================================================
+% GRADIENT WITH ALIGNMENT
+% =========================================================================
+
+% define x vectors
+x_out_def = x_in;
+x_out_pos = x_in + dx/2;
+x_out_neg = x_in - dx/2;
+
+% compute gradient using DCT-I (WSWS symmetry)
+dfdx_def = gradientDtt1D(f, dx, 1, 0);
+dfdx_pos = gradientDtt1D(f, dx, 1, 1);
+dfdx_neg = gradientDtt1D(f, dx, 1, 2);
+
+% analytical derivative
+x_an = -dx:dx/5:2*pi - dx;
+dfdx_an = -sin(x_an);
+
+% plot and compare
+figure;
+plot(...
+ x_an, dfdx_an, 'k-', ...
+ x_out_def, dfdx_def, 'r.', ...
+ x_out_pos, dfdx_pos, 'bo', ...
+ x_out_neg, dfdx_neg, 'g.');
+legend(...
+ 'Analytical Derivative', ...
+ 'DTT shift = 0', ...
+ 'DTT shift = 1', ...
+ 'DTT shift = 2', ...
+ 'Location', 'Best');
+xlabel('x');
+ylabel('f''(x)');
+title('DCT-I WSWS (align output = true)');
+set(gca, 'XLim', [-dx, 2*pi - dx]);
+grid on;
+
+% =========================================================================
+% GRADIENT WITHOUT ALIGNMENT
+% =========================================================================
+
+% define x vectors
+x_out_def = x_in(2:end-1);
+x_out_pos = x_in(1:end-1) + dx/2;
+x_out_neg = x_in(1:end-1) + dx/2;
+
+% compute gradient using DCT-I (WSWS symmetry)
+dfdx_def = gradientDtt1D(f, dx, 1, 0, false);
+dfdx_pos = gradientDtt1D(f, dx, 1, 1, false);
+dfdx_neg = gradientDtt1D(f, dx, 1, 2, false);
+
+% analytical derivative
+x_an = -dx:dx/5:2*pi - dx;
+dfdx_an = -sin(x_an);
+
+% plot and compare
+figure;
+plot(...
+ x_an, dfdx_an, 'k-', ...
+ x_out_def, dfdx_def, 'r.', ...
+ x_out_pos, dfdx_pos, 'bo', ...
+ x_out_neg, dfdx_neg, 'g.');
+legend(...
+ 'Analytical Derivative', ...
+ 'DTT shift = 0', ...
+ 'DTT shift = 1', ...
+ 'DTT shift = 2', ...
+ 'Location', 'Best');
+xlabel('x');
+ylabel('f''(x)');
+title('DCT-I WSWS (align output = false)');
+set(gca, 'XLim', [-dx, 2*pi - dx]);
+grid on;
\ No newline at end of file
diff --git a/gradientDTT1D.m b/gradientDtt1D.m
similarity index 53%
rename from gradientDTT1D.m
rename to gradientDtt1D.m
index c06d2aa..49fa19e 100644
--- a/gradientDTT1D.m
+++ b/gradientDtt1D.m
@@ -1,8 +1,8 @@
-function deriv = gradientDTT1D(f, dx, dtt_type, shift)
+function deriv = gradientDtt1D(f, dx, dtt_type, shift, align_output)
%GRADIENTDTT1D Calculate gradient using discrete trigonometric transforms.
%
% DESCRIPTION:
-% gradientDTT1D computes the spectral gradient of a 1D input function
+% gradientDtt1D computes the spectral gradient of a 1D input function
% using discrete trigonometric transforms (DTTs). The DTTs used to
% transform the input to and from the frequency domain are chosen based
% on the symmetry of the input f, defined using dtt_type. The gradient
@@ -14,47 +14,55 @@
% need to be defined). This means that in some cases the basis
% functions weights are trimmed or appended after the forward
% transform. To ensure the output of gradientDTT1D is the same length
-% as the input f, additional values are appended or removed from the
-% calculated gradient after the inverse transform is calculated. These
-% values are known from the symmetry of the output sequence.
+% as the input f, additional values can be appended or removed from the
+% calculated gradient after the inverse transform is calculated by
+% setting the optional input pad to true (the default). These values
+% are known from the symmetry of the output sequence.
%
-% For additional details, see [1].
+% For additional details on gradient calculation using DTTs, see [1].
%
% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral
% time-domain (PSTD) methods for the wave equation: Realising boundary
% conditions with discrete sine and cosine transforms", 2020.
%
% INPUTS:
-% f - matrix or vector to find the gradient of
-% dx - grid point spacing
-% dtt_type - type of discrete trigonometric transform, this should
-% correspond to the assumed input symmetry of the input
-% function, where:
+% f - Vector to find the gradient of.
+% dx - Grid point spacing.
+% dtt_type - Type of discrete trigonometric transform. This should
+% correspond to the assumed input symmetry of the input
+% function, where:
%
-% 1: DCT-I WSWS
-% 2: DCT-II HSHS
-% 3: DCT-III WSWA
-% 4: DCT-IV HSHA
-% 5: DST-I WAWA
-% 6: DST-II HAHA
-% 7: DST-III WAWS
-% 8: DST-IV HAHS
+% 1: DCT-I WSWS
+% 2: DCT-II HSHS
+% 3: DCT-III WSWA
+% 4: DCT-IV HSHA
+% 5: DST-I WAWA
+% 6: DST-II HAHA
+% 7: DST-III WAWS
+% 8: DST-IV HAHS
%
% OPTIONAL INPUTS:
-% shift - integer controlling whether derivative is shifted to a
-% staggered grid (default = 0), where
+% shift - Integer controlling whether derivative is shifted to a
+% staggered grid (default = 0), where
%
-% 0: no shift
-% 1: shift by + dx/2
-% 2: shift by - dx/2
+% 0: no shift
+% 1: shift by + dx/2
+% 2: shift by - dx/2
+%
+% align_output - Boolean controlling whether the returned values are
+% padded and trimmed based on the implied symmetry so
+% the output is the same length as the input (default =
+% true). Note, if align_output is false, then the
+% gradient calculated for shift = 1 and shift = 2 will
+% be the same.
%
% OUTPUTS:
-% dfdx - gradient of the input function
+% dfdx - Gradient of the input function.
%
% ABOUT:
-% author - Bradley Treeby
-% date - 23 April 2013
-% last update - 2 April 2020
+% author - Bradley Treeby
+% date - 23 April 2013
+% last update - 20 April 2020
%
% Copyright (C) 2013-2020 Bradley Treeby
%
@@ -84,15 +92,22 @@
DST4 = 8; % HAHS
% check for shift input
-if nargin < 4
- shift = false;
+if (nargin < 4) || isempty(shift)
+ shift = 0;
end
-% make sure the input is 1D
-if numDim(f) > 1
- error('Input function must be 1D.');
+% check for pad input
+if (nargin < 5) || isempty(align_output)
+ align_output = true;
end
+% check inputs
+validateattributes(f, {'numeric'}, {'real', 'vector', 'nonsparse'}, 'gradientDTT1D', 'f', 1);
+validateattributes(dx, {'numeric'}, {'real', 'scalar', 'nonsparse'}, 'gradientDTT1D', 'dx', 2);
+validateattributes(dtt_type, {'numeric'}, {'integer', 'scalar', 'nonsparse', '>=', 1, '<=', 8}, 'gradientDTT1D', 'dtt_type', 3);
+validateattributes(shift, {'numeric'}, {'integer', 'scalar', 'nonsparse', '>=', 0, '<=', 2}, 'gradientDTT1D', 'shift', 4);
+validateattributes(align_output, {'logical'}, {'scalar'}, 'gradientDTT1D', 'align_output', 5);
+
% reshape to be a row vector
f = reshape(f, 1, []);
@@ -341,161 +356,163 @@
% add back in the implied values so the output is the same length as the
% input
-switch dtt_type
- case 1
-
- switch shift
- case 0
-
- % WSWS -> WAWA, add both endpoints
- deriv = [0, deriv, 0];
-
- case 1
-
- % WSWS -> HAHA, shift right, mirror right endpoint
- deriv = [deriv, -deriv(end)];
-
- case 2
-
- % WSWS -> HAHA, shift left, mirror left endpoint
- deriv = [-deriv(1), deriv];
-
- end
-
- case 2
-
- switch shift
- case 0
-
- % HSHS -> HAHA, no change
-
- case 1
-
- % HSHS -> WAWA, shift right, append zero
- deriv = [deriv, 0];
-
- case 2
-
- % HSHS -> WAWA, shift left, prepend zero
- deriv = [0, deriv];
-
- end
-
- case 3
-
- switch shift
- case 0
-
- % WSWA -> WAWS, prepend zero, remove right endpoint
- deriv = [0, deriv(1:end - 1)];
-
- case 1
-
- % WSWA -> HAHS, shift right, no change
-
- case 2
-
- % WSWA -> HAHS, shift left, mirror left endpoint, remove
- % right endpoint
- deriv = [-deriv(1), deriv(1:end - 1)];
-
- end
-
- case 4
-
- switch shift
- case 0
-
- % HSHA -> HAHS, no change
-
- case 1
-
- % HSHA -> WAWS, shift right, no change
-
- case 2
-
- % HSHA -> WAWS, shift left, prepend zero, remove right
- % endpoint
- deriv = [0, deriv(1:end - 1)];
-
- end
-
- case 5
-
- switch shift
- case 0
-
- % WAWA -> WSWS, remove both endpoints
- deriv = deriv(2:end - 1);
-
- case 1
-
- % WAWA -> HSHS, shift right, remove left endpoint
- deriv = deriv(2:end);
-
- case 2
-
- % WAWA -> HSHS, shift left, remove right endpoint
- deriv = deriv(1:end - 1);
-
- end
-
- case 6
-
- switch shift
- case 0
-
- % HAHA -> HSHS, no change
-
- case 1
-
- % HAHA -> WSWS, shift right, remove left endpoint
- deriv = deriv(2:end);
-
- case 2
-
- % HAHA -> WSWS, shift left, remove right endpoint
- deriv = deriv(1:end - 1);
-
- end
-
- case 7
-
- switch shift
- case 0
-
- % WAWS -> WSWA, remove left endpoint, append zero
- deriv = [deriv(2:end), 0];
-
- case 1
-
- % WAWS -> HSHA, shift right, remove left endpoint, mirror
- % right endpoint
- deriv = [deriv(2:end), -deriv(end)];
-
- case 2
-
- % WAWS -> HSHA, shift left, no change
-
- end
-
- case 8
-
- switch shift
- case 0
-
- % HAHS -> HSHA, no change
-
- case 1
-
- % HAHS -> WSWA, shift right, remove left endpoint, append
- % zero
- deriv = [deriv(2:end), 0];
-
- case 2
-
- % HAHS -> WSWA, shift left, no change
-
- end
-
+if align_output
+ switch dtt_type
+ case 1
+
+ switch shift
+ case 0
+
+ % WSWS -> WAWA, add both endpoints
+ deriv = [0, deriv, 0];
+
+ case 1
+
+ % WSWS -> HAHA, shift right, mirror right endpoint
+ deriv = [deriv, -deriv(end)];
+
+ case 2
+
+ % WSWS -> HAHA, shift left, mirror left endpoint
+ deriv = [-deriv(1), deriv];
+
+ end
+
+ case 2
+
+ switch shift
+ case 0
+
+ % HSHS -> HAHA, no change
+
+ case 1
+
+ % HSHS -> WAWA, shift right, append zero
+ deriv = [deriv, 0];
+
+ case 2
+
+ % HSHS -> WAWA, shift left, prepend zero
+ deriv = [0, deriv];
+
+ end
+
+ case 3
+
+ switch shift
+ case 0
+
+ % WSWA -> WAWS, prepend zero, remove right endpoint
+ deriv = [0, deriv(1:end - 1)];
+
+ case 1
+
+ % WSWA -> HAHS, shift right, no change
+
+ case 2
+
+ % WSWA -> HAHS, shift left, mirror left endpoint, remove
+ % right endpoint
+ deriv = [-deriv(1), deriv(1:end - 1)];
+
+ end
+
+ case 4
+
+ switch shift
+ case 0
+
+ % HSHA -> HAHS, no change
+
+ case 1
+
+ % HSHA -> WAWS, shift right, no change
+
+ case 2
+
+ % HSHA -> WAWS, shift left, prepend zero, remove right
+ % endpoint
+ deriv = [0, deriv(1:end - 1)];
+
+ end
+
+ case 5
+
+ switch shift
+ case 0
+
+ % WAWA -> WSWS, remove both endpoints
+ deriv = deriv(2:end - 1);
+
+ case 1
+
+ % WAWA -> HSHS, shift right, remove left endpoint
+ deriv = deriv(2:end);
+
+ case 2
+
+ % WAWA -> HSHS, shift left, remove right endpoint
+ deriv = deriv(1:end - 1);
+
+ end
+
+ case 6
+
+ switch shift
+ case 0
+
+ % HAHA -> HSHS, no change
+
+ case 1
+
+ % HAHA -> WSWS, shift right, remove left endpoint
+ deriv = deriv(2:end);
+
+ case 2
+
+ % HAHA -> WSWS, shift left, remove right endpoint
+ deriv = deriv(1:end - 1);
+
+ end
+
+ case 7
+
+ switch shift
+ case 0
+
+ % WAWS -> WSWA, remove left endpoint, append zero
+ deriv = [deriv(2:end), 0];
+
+ case 1
+
+ % WAWS -> HSHA, shift right, remove left endpoint, mirror
+ % right endpoint
+ deriv = [deriv(2:end), -deriv(end)];
+
+ case 2
+
+ % WAWS -> HSHA, shift left, no change
+
+ end
+
+ case 8
+
+ switch shift
+ case 0
+
+ % HAHS -> HSHA, no change
+
+ case 1
+
+ % HAHS -> WSWA, shift right, remove left endpoint, append
+ % zero
+ deriv = [deriv(2:end), 0];
+
+ case 2
+
+ % HAHS -> WSWA, shift left, no change
+
+ end
+
+ end
end
\ No newline at end of file