function [press_zemanek] = zemanek_ciratten(f0, fs, c, attenfreq, x_real, y_real, z_real)
%ZEMANEK_CIRATTEN Calculate the pressure generated by the Zemanek method
%with attenuation. Refer: J. Zemanek, “Beam behavior within the nearfield
%of a vibrating piston,” J. Acoust. Soc. Am., vol. 49, no. 1, pp. 181–191,
%1971.
%   Usage: [press_zemanek] = zemanek_ciratten(f0, fs, c, x_real, y_real, z_real)
%   f0 - The central frequency of the pulse. Unit: Hz
%   fs - The sampling frequency. Unit: Hz
%   c  - The speed of sound. Unit: m/s
%   attenfreq - The frequency dependent attenuation. Unit: db/cm/MHz
%   x_real - The x coordinate in row. Unit: meter
%   y_real - The y coordinate in row. Unit: meter
%   z_real - The z coordinate in row. Unit: meter

%   Copyright 2006 MSU
%   $Vision: 1.0 $  $Date: 2006/06/15 $

wavelength = c/f0;
radius = 4*wavelength;
wavenum = 2*pi/wavelength;

atten = attenfreq*f0/1e6;   % db/cm  db/cm/MHz
atten = atten*0.1151;  % np/cm

n = 2*round(4*radius/wavelength);   %% subdivision in radial direction
m = 2*round(2*n*pi);           %% subdivision in angle direction

%% use a simple loop to calculate the core part of the Zemanek point source
%% method.
for ix = 1:size(x_real,2)
    for iy = 1:size(y_real,2)
        for iz = 1:size(z_real,2)
            
            r = sqrt((x_real(ix)).^2 + (y_real(iy)).^2);
            z = z_real(iz);
            
            X = r;
            Z = z;
            
            % delta_sigma = 1/(n*radius/wavelength);  %% changed
            delta_sigma = radius/n;
            delta_phi  = 2*pi/m;
            
            R = sqrt(Z.^2 + X.^2);
            theta = atan(X/Z);
            
            press_tmp = 0;
            
            for p=1:m
                for q=1:n
                    sigma_q = delta_sigma * (q-0.5);
                    phi_p = delta_phi * (p-0.5);
                    Rpq = sqrt(R.^2 + sigma_q.^2 - 2*R.*sigma_q.*sin(theta).*cos(phi_p));
                    delta_Sq = sigma_q.*delta_sigma.*delta_phi;
                    wavenew = 2*pi/wavelength - j*atten*100;   %% update the wavenum here to incorporate the attenuation
                    press_tmp = press_tmp + delta_Sq/Rpq*exp(-j*wavenew* Rpq);
                end
            end
            press_zemanek(ix,iy,iz) = press_tmp/(2*pi);
        end
    end
end