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{:});