function [press_fnm]=FNM_ciratten(f0, fs, c, attenfreq, x_real, y_real, z_real)
%FNM_CIRATTEN Calculate the pressure generated by the fast nearfield method
%with attenuation. Refer:  R. J. McGough, T. V. Samulski, and J. F. Kelly,
%"An e_cient grid sectoring method for calculations of the nearfield pressure
%generated by a circular piston,” J. Acoust. Soc. Am., vol. 115, no. 5, 
%pp. 1942–1954, 2004.
%   Usage: [press_fnm]=FNM_ciratten(f0, fs, c, attenfreq, 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 $


% Use to calculate the Gauss abscissas and weights 
load gauss1to500.mat;  %% abscissas; weights
gauss_abscissas = abscissas';  %% transpose   %%% notice keep consistence with dll file
gauss_weight     = weights;
max_interval = pi;
min_interval = 0;
abscissas_all = (max_interval-min_interval)/2*gauss_abscissas + (max_interval+min_interval)/2;
weights_all = gauss_weight*(max_interval-min_interval)/2;
 

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

wavelength = c/f0  %% changed
R_coef = 4; % Radius of transducer
R = R_coef*wavelength;

tic
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);
            
            abscissas = abscissas_all(500,1:500);
            weights   = weights_all(500,1:500);
            
            tmp1 = (r*cos(abscissas)-R)./(r.^2+R.^2-2*R*r*cos(abscissas));
            wav_new = 2*pi/wavelength - j*atten*100;
            tmp2 = exp(-j*wav_new*sqrt(r.^2+R.^2-2*R*r*cos(abscissas)+z.^2)) - exp(-j*wav_new*z);
            
          
            press_fnm(ix,iy,iz) = R/(j*wav_new*pi)*sum(tmp1.*tmp2.*weights);             
        end
    end
end
toc

press_fnm = press_fnm;
% press_fnm = press_fnm/max(max(max(abs(press_fnm))));