Hello all,
I've seen that some people asked before me for this topic but I don't find the correct solution that works for me.
I'm trying to obtain a relation among radius of radiating surface and parametric difference non linear beam generation.
I've two main question about code that I show below.
- I obtain a pressure of total beam pattern greater than I putted like a entry in my code. I've read in this forum that probably this can be solved using dirichlet option in source.p_mode but do not work in my case (Probably I miss something..). In any case I solved this problem finding a relationship between dx & dy and Nx & Ny.
How can I solved this problem ?
- With this relation that I found the maximum frequency that I can simulate is ~ 300kHz that is enough for difference frequency non linear beam but I would like observe also 2·f2 beam.
How can I increase the frequency range validity of the simulation without losing good amplitude of pressure emitted (in order to no obtain greater amplitude that I putted like amplitude of bifrequency signal).
I put my code below.
Thanks in advance for your help.
Cheers
***************************************************************************************
Amp = [100 1000 10000 100000 500000 750000 1000000 2000000 3000000 4000000 5000000];
for i = 1:11
USE_STATISTICS = false;
f_max = 170000;
c0_min = 1482;
points_per_wavelength = 3.5;
x_size = 0.25;
y_size = 0.5;
dx = (c0_min/(points_per_wavelength*f_max));
Nx = round(x_size/dx);
dy = (c0_min/(points_per_wavelength*f_max));
Ny = round(y_size/dy);
kgrid = makeGrid(Nx, dx, Ny, dy);
[kgrid.t_array, dt] = makeTime(kgrid, c0_min)
medium.sound_speed = 1482*ones(Nx, Ny); % [m/s]
medium.density = 1000*ones(Nx, Ny); % [kg/m^3]
att1 = 2.17e-3; %Attenuation dB/cm·MHz2
medium.alpha_coeff=att1*ones(Nx, Ny);
medium.alpha_power=2;
medium.BonA = 4.96;
radio_dist = 0.0375;
r = radio_dist/dx;
em_pos_y = round(0.01/dy);
source.p_mask = zeros(Nx, Ny);
source.p_mask(round((Nx/2)-r):round((Nx/2)+r),em_pos_y) = 1;
f1 = 120000;
f2 = 170000;
Mag_Pa = Amp(i);
signal = Mag_Pa.*(0.5.*sin(2*pi*f1.*kgrid.t_array) + 0.5.*sin(2*pi*f2.*kgrid.t_array));
source.p = signal;
source.p_mode = 'dirichlet';
source.p = signal(1:375);
source.p = filterTimeSeries(kgrid, medium, source.p);
%% Define a Cartesian sensor mask Directividad
sensor.mask = [1, 1, Nx, Ny].';
sensor.record = {'p','p_max'};
input_args = { 'DisplayMask', source.p_mask, 'RecordMovie', true, 'MovieType', 'image', 'MovieName', 'Directividad_Piston_2D', 'PMLInside', false, 'PlotScale', [-Mag_Pa-0.2*Mag_Pa, Mag_Pa+0.2*Mag_Pa]};
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,input_args{:});
% Post-Procesado
[freq, amp_spect] = spect(sensor_data.p, 1/dt, 'Dim', 3);
%Compute and plot the fundamental beam
[f1_value, f1_index] = findClosest(freq, f1);
beam_pattern_f1 = amp_spect(:, :, f1_index);
figure;
subplot(2,2,1); imagesc((kgrid.y_vec+max(kgrid.y_vec)).*100, (kgrid.x_vec+max(kgrid.x_vec)).*100, sensor_data.p_max);
colormap(jet(1024));
colorbar
ylabel('x-position [cm]');
xlabel('y-position [cm]');
title('Total Beam Pattern');
axis image;
%Compute and plot the difference beam
diff_freq = f2-f1;
[f2_value, f2_index] = findClosest(freq, diff_freq);
beam_pattern_f2 = amp_spect(:, :, f2_index);
subplot(2,2,2);imagesc((kgrid.y_vec+max(kgrid.y_vec)).*100, (kgrid.x_vec+max(kgrid.x_vec)).*100, beam_pattern_f2, [0 0.025*Mag_Pa]);
colormap(jet(1024));
colorbar
ylabel('x-position [cm]');
xlabel('y-position [cm]');
axis image;
title('f_2-f_1 Beam Pattern');
%Compute and plot the difference beam
diff2_freq = 2*f1;
[f3_value, f3_index] = findClosest(freq, diff2_freq);
beam_pattern_f3 = amp_spect(:, :, f3_index);
subplot(2,2,3);imagesc((kgrid.y_vec+max(kgrid.y_vec)).*100, (kgrid.x_vec+max(kgrid.x_vec)).*100, beam_pattern_f3);
colormap(jet(1024));
colorbar
ylabel('x-position [cm]');
xlabel('y-position [cm]');
axis image;
title('2·f1 Beam Pattern');
%Compute and plot the difference beam
f1_f2_freq = f1+f2;
[f4_value, f4_index] = findClosest(freq, f1_f2_freq);
beam_pattern_f4 = amp_spect(:, :, f4_index);
subplot(2,2,4);imagesc((kgrid.y_vec+max(kgrid.y_vec)).*100, (kgrid.x_vec+max(kgrid.x_vec)).*100, beam_pattern_f4);
colormap(jet(1024));
colorbar
ylabel('x-position [cm]');
xlabel('y-position [cm]');
axis image;
title('f1 + f2 Beam Pattern');
end