function [press_zemanek,press_fnm,press_field,press_field_be] = test_fun()
%TEST_FUN gives the basic usage of the three codes for the FNM, FieldII and
%point source superposition methods.

% basic parameters
f0 = 2.0e6;         % Excitation frequency. Unit: Hz
fs = 1024e6;        % The sampling frequency. Unit: Hz
c = 1500;           % speed of sound in water. Unit: m/s

attenfreq = 0.5;    % attenuation. db/cm/MHz

wavelength = c/f0;  % The wavelength of the wave
radius = 4*wavelength;  % The radius of the circular piston

% for the on axis situation
x_vector = 0;
y_vector = 0;
z_vector = 0.1:0.025:4; %% on asix


% % for the plane that parallel to the piston face
% x_vector = -1.5:0.075:1.5;
% y_vector = -1.5:0.075:1.5;
% z_vector = 0.5; %% YX plane

x_unit = radius;
y_unit = radius;
z_unit = radius.^2/wavelength;

x_real = x_vector*x_unit;
y_real = y_vector*y_unit;
z_real = z_vector*z_unit;

% For the FNM
[press_fnm]=FNM_ciratten(f0, fs, c, attenfreq, x_real, y_real, z_real);
press_fnm = squeeze(press_fnm);

% For the point source method (zemanek)
[press_zemanek] = zemanek_ciratten(f0, fs, c, attenfreq, x_real, y_real, z_real);
press_zemanek = squeeze(press_zemanek);

% For the Field II program after correction
[press_field]=fieldII_ciratten(f0, fs, c, attenfreq, x_real, y_real, z_real);
press_field = squeeze(press_field);

% If you want to compute the pressure field for the Field II before
% correction, just change the code I have listed in my website in file
% "fieldII_ciratten.m".


% Plot the figures
if (size(press_fnm,1)==1|size(press_fnm,2)==1)
    figure,plot(z_vector, abs(squeeze(press_fnm-press_field))/max(max(abs(squeeze(press_fnm)))),'r-'), hold on
    plot(z_vector, abs(squeeze(press_fnm-press_zemanek))/max(max(abs(squeeze(press_zemanek)))),'g:'), hold off
    axis([0 4 0 0.5])
    xlabel('z (unit of a^2/\lambda)');
    ylabel('normalized error');
    legend('Field II after correction','zemanek method');
else
    figure,mesh(x_vector,y_vector, abs(squeeze(press_fnm-press_field))/max(max(abs(squeeze(press_fnm))))) 
    xlabel('x (unit of a)')
    ylabel('y (unit of a)')
    zlabel('normalized error');
    title('Field II after correction')
    
    figure,mesh(x_vector,y_vector, abs(squeeze(press_fnm-press_zemanek))/max(max(abs(squeeze(press_fnm)))))
    xlabel('x (unit of a)');
    ylabel('y (unit of a)');
    zlabel('normalized error');
    title('Zemanek method');
end