function [press_field]=fieldII_ciratten(f0, fs, c, attenfreq, x_real, y_real, z_real)
%FIELDII_CIRATTEN Calculate the pressure generated by FIELD II program 
%with attenuation after correction. Refer:
%http://www.es.oersted.dtu.dk/staff/jaj/field
%   Usage: [press_field]=fieldII_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 rectangular
set_field('use_rectangles',1);
set_field('use_triangles',0);
% set the sampling frequency and the speed of sound
set_field('c',c);
set_field('fs',fs);

atten = attenfreq*f0/1e6;   % db/cm  db/cm/MHz
%% set the attenuation
set_field('att', atten*100);  %% db/m
set_field('freq_att',atten*0.1151/f0);
set_field('att_f0', f0);
set_field('use_att',1);

wavelength = c/f0  %% changed

% Set initial parameters
R_coef = 4; % Radius of transducer
R = R_coef*wavelength;
ele_size=R/128; % Size of mathematical elements
% Define the transducer
Th = xdc_piston (R,ele_size);

% points = [x_real, y_real, z_real];

%%########### This part for the impulse response+
% field_h = zeros(size(t,2), size(y_real,2), size(z_vector,2));
tic
for ix = 1:size(x_real,2)
    for iy = 1:size(y_real,2)
        for iz = 1:size(z_real,2)
            
            points = [x_real(ix)  y_real(iy)  z_real(iz)];
            [h,start_time] = calc_h(Th,points); 
        
            t = start_time:1/fs:start_time+(length(h)-1)/fs;
            t = t';          
            F_impulse = sum(h.*exp(-i*2*pi*f0*t))*1/fs;   %% 1/fs is the interval of integrand  
            press_field(ix,iy,iz) = F_impulse;    
        end
    end
end
toc

press_field = press_field*fs;
% press_field = press_field/max(max(max(abs(press_field))));
% figure,plot(t*1e6, hh/max(hh))