addpath('C:\Users\labadmin\Desktop\k-Wave')
clear all;
clc;
%experimental properties
water_speed = 1480;
gly_speed = 1560;
gt = 1 %glycerin thickness in.
gt = gt*2.54*1e-2; %glycerin thickness m
% grid properties
pts_per_lam = 8;
c0_min = water_speed;
f_max = 500000;
x_size = .4;
dx = c0_min/(pts_per_lam*f_max);
dy = dx;
Nx = round(x_size/dx);
Ny = round(Nx/2);
% check prime factors
checkFactors(Ny - 40,Ny + 40)
checkFactors(Nx - 20, Nx + 20)
close all
%%
%set final grid properties
Nx = 1024;
Ny = 512;
%%
% create the grid and time array
kgrid = makeGrid(Nx,dx,Ny,dy);
CFL = .2;
c0_max = gly_speed;
dt = CFL*dx/c0_max;
kgrid.t_array = 0:dt:120e-6;
%medium properties
medium.sound_speed = water_speed*ones(Nx,Ny); % [m/s] Water
medium.density = 997*ones(Nx,Ny); % [kg m^-3] Water
%medium.BonA = 5*ones(Nx,Ny); %assume linear at first
%medium.alpha_coeff = .0022*ones(Nx,Ny); %assume lossless at first
%medium.alpha_power = 1.05*ones(Nx,Ny); %assume lossless at first
medium.alpha_mode = 'no_dispersion';
load('Elem_X')
load('Elem_Z')
coords = [3 7 13 20 28 85 88 92 98 105]; %elem #
%shift in x to fit in the kgrid
x_shift = min(Elem_X(coords)) - min(min(kgrid.x)) -.01;
y_shift = min(Elem_Z(coords)) - min(min(kgrid.y)) -.01;
[grid_data,order_id,re_order_id] = cart2grid(kgrid,[Elem_X(coords)-x_shift Elem_Z(coords)-y_shift]');
[val id]= max(Elem_Z(coords));
[rows cols] = find(grid_data);
gly_start = round(cols(order_id(id)) + .013/dx);
gt = round(gt/dy);
% medium.sound_speed(28:727,gly_start:gly_start+gt) = gly_speed; % [m/s] Glycerin
% medium.density(28:727,gly_start:gly_start+gt) = 1034.2; % [kg/m^3] % Glycerin
medium.sound_speed_ref = c0_min; %ensures phase error is bounded
%check dt
kmax = pi/dx;
dt < 2/(medium.sound_speed_ref*kmax)*asin(medium.sound_speed_ref/c0_max)
%%
% build transducer
diam = floor(.02/dx);
if rem(diam,2) == 0
diam = diam + 1;
end
[rows,cols] = find(grid_data);
bm = makeMultiArc([Nx, Ny],[rows cols], 9999, diam, [409 377]); %geometric focus array in gridpoints
imagesc(bm)
axis equal
source.p_mask = bm;
source_freq = 0.5e6; % [Hz]
source_mag = -10; % [Pa]
tc = gauspuls('cutoff',source_freq,0.6,[],-50);
time = -tc : kgrid.dt : kgrid.t_array(end) ;
pulse = gauspuls(time,source_freq,0.6);
source.p = filterTimeSeries(kgrid, medium, source_mag*pulse);
sensor.mask = zeros(Nx,Ny);
sensor.mask(190,205) = 1;
% Define data types to be stored in sensor_data
sensor.record = {'p', 'p_final','p_max_all','p_min_all'}; %[ p is acoustic pressure at each sensor, p_final is acoustic field across grid]
%Define Record start index
sensor.record_start_index = 1; %Start recording at kgrid.t_array( x )
%% Run the Sim
input_args = {'DataCast', 'single','PMLInside',false}; %put PML outside of grid
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,input_args{:});