k-Wave User Forum » User Favorites: MGCain

k-Wave User Forum » User Favorites: MGCain mgcain Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 09:15:30 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php MGCain on "Simulate ultrasound pressure on piezoelectric inclusions in medium" http://www.k-wave.org/forum/topic/simulate-ultrasound-pressure-on-piezoelectric-inclusions-in-medium#post-7582 Mon, 08 Jun 2020 14:12:57 +0000 MGCain 7582@http://www.k-wave.org/forum/ Thanks for this Brad. I will explore simulation with k-wave to determine pressure fields and model the piezo response analytically. I can also use COMSOL later on if necessary. Thanks again. Bradley Treeby on "Simulate ultrasound pressure on piezoelectric inclusions in medium" http://www.k-wave.org/forum/topic/simulate-ultrasound-pressure-on-piezoelectric-inclusions-in-medium#post-7571 Sat, 06 Jun 2020 18:43:11 +0000 Bradley Treeby 7571@http://www.k-wave.org/forum/ Hi MGCain, If you're interested in small inclusions and coupling with piezoelectricity, k-Wave is probably not the right tool. A finite-element model that allows you to compute coupled physics might be the way to go. Brad. MGCain on "Simulate ultrasound pressure on piezoelectric inclusions in medium" http://www.k-wave.org/forum/topic/simulate-ultrasound-pressure-on-piezoelectric-inclusions-in-medium#post-7537 Mon, 01 Jun 2020 14:58:05 +0000 MGCain 7537@http://www.k-wave.org/forum/ Hi, I am interested in using k-wave to sinmulate an ultrasound pressure wave on dispersed micron-sized inclusions of piezoelectric material inside a polymer matrix. In the simplest formalism I could launch a plane wave at one edge and measure pressure at all points inside my 3D volume. From pressure I can deduce the piezo response (I expect!). Any advice please on the specific approach/ example set I could start to explore? Thanks. rainbowwhale on "Simulation of a flat reflector" http://www.k-wave.org/forum/topic/simulation-of-a-flat-reflector#post-7374 Tue, 07 Apr 2020 03:24:36 +0000 rainbowwhale 7374@http://www.k-wave.org/forum/ CFL number should be calculated using highest speed of sound. By using highest speed of sound, current CFL number is about 4.13. I tried using dt as 0.05/(PPP*F0) which has CFL number of about 0.2 and it worked. Ryuji on "Simulation of a flat reflector" http://www.k-wave.org/forum/topic/simulation-of-a-flat-reflector#post-7363 Sat, 04 Apr 2020 02:07:08 +0000 Ryuji 7363@http://www.k-wave.org/forum/ Hi, thanks for providing the great library. I would like to simulate a small 40kHz-ultrasound transducer with a flat reflector in 3D. When I set the medium heterogeneous (air and acrylic), I cannot see any propagation even in the air. Is this because of the large change in impedance? I would appreciate it if you could tell me how to make my code work. Thanks advance! Best regards, Ryuji % ========================================================================= % sorce code % ========================================================================= clear all; % ========================================================================= % DEFINE LITERALS % ========================================================================= % select which code to run % 1: MATLAB CPU code % 2: MATLAB GPU code % 3: C++ code MODEL = 2; % scale parameter (changes the resolution of the simulation) SC = 1; % grid parameters PML_SIZE = 10; % size of the perfectly matched layer [grid points] Nx = 128*SC - 2*PML_SIZE; % number of grid points in the x direction Ny = 64*SC - 2*PML_SIZE; % number of grid points in the y direction Nz = 64*SC - 2*PML_SIZE; % number of grid points in the z direction % sampling parameters PPW = 10*SC; % points per wavelength (this defines the grid size) PPP = 10*SC; % points per period (this defines the CFL) T_END = 1e-3; % total compute time [s] (this must be long enough to reach steady state) RECORD_PERIODS = 2; % number of periods to record % medium parameters C0 = 346; % sound speed [m/s] RHO0 = 1.2; % density [kg/m^3] CR = 1.43e3; % sound speed in the reflector [m/s] RHOR = 1.18e3; % density of the reflector [kg/m^3] % source parameters F0 = 40e3; % source frequency [Hz] SOURCE_RADIUS = 5e-3; % piston radius [m] SOURCE_MAG = 0.5e6; % source pressure [Pa] % ========================================================================= % CREATE GRID % ========================================================================= % calculate the grid spacing based on the PPW and F0 dx = C0 / (PPW * F0); % [m] dt = 1 / (PPP*F0); % [s] SOURCE_RADIUS = SOURCE_RADIUS / dx; % calculate the CFL disp(['CFL = ' num2str(C0*dt/dx)]); % create the computational grid kgrid = makeGrid(Nx, dx, Ny, dx, Nz, dx); % create the time array kgrid.t_array = 0:dt:T_END; % assign medium properties % medium.sound_speed = C0; % medium.density = RHO0; medium.sound_speed = C0 * ones(Nx, Ny, Nz); medium.density = RHO0 * ones(Nx, Ny, Nz); % assign reflector properties medium.sound_speed(Nx, :, :) = CR; medium.density(Nx, :, :) = RHOR; % assign the reference sound speed to the background medium % medium.sound_speed_ref = C0; % ========================================================================= % CREATE SOURCE % ========================================================================= % create piston source mask piston = makeDisc(Ny, Nz, Ny/2, Nz/2, SOURCE_RADIUS); source.p_mask = zeros(Nx, Ny, Nz); source.p_mask(1, :, :) = piston; % create time varying source source.p = SOURCE_MAG * sin(2*pi*F0*kgrid.t_array); % filter source with up-ramp ramp_length = round((2*pi/F0)/dt); ramp = (-cos( (0:(ramp_length-1))*pi/ramp_length ) + 1)/2; source.p(1:ramp_length) = ramp .* source.p(1:ramp_length); % ========================================================================= % CREATE SENSOR MASK % ========================================================================= % set sensor mask to record central axis, not including the source point % sensor.mask = zeros(Nx, Ny, Nz); % sensor.mask(2:end, Ny/2, Nz/2) = 1; sensor.mask = [2, 1, Nz/2, Nx, Ny, Nz/2].'; % record the maximum pressure sensor.record = {'p', 'p_rms'}; % average only the final few periods when the field is in steady state sensor.record_start_index = kgrid.Nt - RECORD_PERIODS*PPP + 1; % ========================================================================= % RUN k-WAVE SIMULATION % ========================================================================= % run code switch MODEL case 1 % MATLAB CPU sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ... 'DataCast', 'single', 'PMLInside', false, ... 'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG); case 2 % MATLAB GPU sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ... 'DataCast', 'gpuArray-single', 'PMLInside', false, ... 'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG); case 3 % C++ sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, 'PMLInside', false); otherwise end % ========================================================================= % ANALYTICAL SOLUTION % ========================================================================= % calculate the wavenumber k = 2*pi*F0./C0; % define radius and axis a = SOURCE_RADIUS*dx; x = (kgrid.x_vec(2:end, :) - min(kgrid.x_vec(:))); % calculate the analytical solution for a piston in an infinite baffle % for comparison (Eq 5-7.3 in Pierce) r_a = sqrt(x.^2 + a^2); p_ref = SOURCE_MAG*abs(-2i*exp(1i*k*(x+r_a)/2).*sin((k*r_a - k*x)/2)); % ========================================================================= % VISUALISATION % ========================================================================= % plot the pressure along the focal axis of the piston figure; plot(sensor_data.p(end,:), 'bx'); img = gather(sensor_data.p_rms); image(img, 'CDataMapping', 'scaled') daspect([1 1 1]) colormap default colorbar Bradley Treeby on "Planewave simulation" http://www.k-wave.org/forum/topic/planewave-simulation#post-5481 Fri, 29 Apr 2016 10:47:03 +0000 Bradley Treeby 5481@http://www.k-wave.org/forum/ Hi Rajesh, From memory you can set both of these values to inf to use an unfocused transducer. Alternatively, you could also create the source and sensor masks manually (i.e., set <code>source.p_mask</code> and <code>sensor.mask</code> rather than using the <code>kWaveTransducer</code> class), and define a single source time series to drive all of the elements. Brad. Bradley Treeby on "Particle acoustophoresis at "larger" scales using ultrasonic waves" http://www.k-wave.org/forum/topic/particle-acoustophoresis-at-larger-scales-using-ultrasonic-waves#post-5478 Fri, 29 Apr 2016 10:35:37 +0000 Bradley Treeby 5478@http://www.k-wave.org/forum/ Hi Andres, The main challenge with any grid based method (including k-Wave) is the accuracy with which you represent the heterogeneous structure of the material on the discrete grid. In some cases, the actual structure isn't that important, just the statistical effect it has on the propagating wave field. k-Wave doesn't have a problem with this type of geometry per-se, but you might want to perform a convergence study. That is, check to what degree the simulation output depends on the grid discretisation. If you want another model to compare to, you could also try <a href="http://www.simsonic.fr/">SimSonic</a>. I've got a wrapper function that lets you call SimSonic2D with k-Wave inputs if that is helpful. Regarding elastic vs fluid simulation, I don't have much experience with the situation you are interested in, but if the media can support shear, then an elastic simulation is probably necessary. In any case, it should be easy enough to run the simulation with the k-Wave fluid and elastic codes and compare the results. Hope that helps, Brad rajeshlvvl on "Planewave simulation" http://www.k-wave.org/forum/topic/planewave-simulation#post-5472 Mon, 18 Apr 2016 04:32:17 +0000 rajeshlvvl 5472@http://www.k-wave.org/forum/ Hi Brad, We aretrying to simulate planewave propagation in a medium using 3D model. The way we are trying to generate a planewave is by giving very high values for azimuthal focus and elevation focus as below. Is the approach right? Or is there any other way to generate them? Thanks -Rajesh transducer.focus_distance = 3000000e-3; % focus distance [m] transducer.elevation_focus_distance = 3000000e-3;% focus distance in the elevation plane [m] transducer.steering_angle = 0; % steering angle [degrees] transducer.steering_angle_max = 0; % maximum steering angle [degrees] % apodization transducer.transmit_apodization = 'Rectangular'; transducer.receive_apodization = 'Rectangular'; transducer.active_elements = ones(transducer.number_elements, 1); transducer.input_signal = input_signal; AndresMM on "Particle acoustophoresis at "larger" scales using ultrasonic waves" http://www.k-wave.org/forum/topic/particle-acoustophoresis-at-larger-scales-using-ultrasonic-waves#post-5457 Fri, 08 Apr 2016 09:23:12 +0000 AndresMM 5457@http://www.k-wave.org/forum/ Hi everyone, I am trying to model the acoustophoretic effect of ultrasound excitation on a column of water containing a particle bed with sizes in the range of 20-90 micrometers. To do this I want to find the pressure and speed time-averaged effects (P_rms and U_rms) and calculate the Gorkov potential, its gradient and obtain the acoustic radiation forces field. My idea is to afterwards couple this force field to an external CFD simulation and see how different sized particles will behave under different flow conditions. Something similar has already been done for microfluidic setups (<a href="http://pubs.rsc.org/en/content/articlelanding/2012/lc/c2lc40733g#" rel="nofollow">http://pubs.rsc.org/en/content/articlelanding/2012/lc/c2lc40733g#</a>!divAbstract) or (<a href="http://pubs.rsc.org/en/content/articlelanding/2015/lc/c5lc00866b#" rel="nofollow">http://pubs.rsc.org/en/content/articlelanding/2015/lc/c5lc00866b#</a>!divAbstract) (However they used a FEM simulation, COMSOL multi-physics). I want to use kwave for the steps of finding P_rms and U_rms for my geometry. However I have some doubts regarding the modeling of the particle bed. What I have done so far is create a porous bed on the bottom of the column. My 2D-mask looks like this <a href="http://imgur.com/a/vpy1B" rel="nofollow">http://imgur.com/a/vpy1B</a> , where the transducer is the top small line and is in direct contact with the main medium (water), the vertical lines are the tank walls (PMMA), and the bottom is a steel bottom (which acts as a wave reflector). - How well can KWave manage such a "complex" geometry including holes and porosities? - Am I right in believing that an elastic simulation is strictly necessary to find the resulting P_rms and U_rms of the geometry? - Would you have any advice regarding the problem in general? Cheers, Andrés Bradley Treeby on "Simulation of ultrasound propagating from soft tissue to air" http://www.k-wave.org/forum/topic/simulation-of-ultrasound-propagating-from-soft-tissue-to-air#post-5450 Tue, 05 Apr 2016 08:55:34 +0000 Bradley Treeby 5450@http://www.k-wave.org/forum/ Hi Jingba, First, it's worth pointing out that the situation you are simulating has an impedance mismatch on the order of 1e6, which makes it challenging to deal with computationally. The NaNs are because the code is unstable with the grid parameters you have used. There are a few problems with your code: (1) The definition of dx should use the minimum sound speed in the medium, which in your case will be air (not 1500 m/s). You can still use the sound speed in tissue to define your value of <code>t_end</code>. (2) The CFL needs to be reduced for the code to be stable. Something like the following should work <pre><code>% create the time array t_end = (Nx*dx)*1.5/c0_tissue; % [s] CFL = 0.05; kgrid.t_array = makeTime(kgrid, c0_tissue, CFL, t_end);</code></pre> (3) You should be aware of the frequency content of your source, and the maximum frequency you have chosen for the grid. If you plot the frequency spectrum of <code>transducer_TX.input_signal</code>, you will see that it has energy above your chosen maximum frequency. You could either increase the maximum frequency, reduce the source frequency, or remove the unsupported frequencies using <code>filterTimeSeries</code>. Hope that helps, Brad Bowen Jing on "Simulation of ultrasound propagating from soft tissue to air" http://www.k-wave.org/forum/topic/simulation-of-ultrasound-propagating-from-soft-tissue-to-air#post-5431 Fri, 18 Mar 2016 08:36:15 +0000 Bowen Jing 5431@http://www.k-wave.org/forum/ Hi, I am trying to model the reflection and refraction of ultrasound at a interface between soft tissue and air. The sound speed and density of air are 340 m/s and 1.29 Kg/m3 respectively. the size of the interface is much larger than the wavelength of ultrasound source. ultrasound frequency is on the order of 10 MHz. However, the ultrasound echo received by the linear transducer is a matrix of NaN (not a number). The Matlab code is shown below. I don't know why it didn't work out. <pre><code>clc clear all close all % simulation settings DATA_CAST = 'double'; % set to 'single' or 'gpuArray-single' to speed up computations %% ======================================================================== = % DEFINE THE K-WAVE GRID % ======================================================================== = % compute grid spacing based on a desired x_size and f_max x_size = 50e-3; % [m] desired domain size on x direction (depth...) c0_min = 1500; % [m/s] f_max = 15e6; % [Hz] maximum frequency of ultrasound interested points_per_wavelength = 2; dx = c0_min/(points_per_wavelength*f_max); % [m] dy = dx; % [m] dz = dx; % [m] pml_x_size = 20; % [grid points] pml_y_size = 20; % [grid points] pml_z_size = 20; % [grid points] Nx_at_smallest_prime_factor = 256; Ny_at_smallest_prime_factor = 128; Nz_at_smallest_prime_factor = 128; % set number of grid points not including the PML Nx = Nx_at_smallest_prime_factor - 2*pml_x_size; % [grid points] Ny = Ny_at_smallest_prime_factor - 2*pml_y_size; % [grid points] Nz = Nz_at_smallest_prime_factor - 2*pml_z_size; % [grid points] % create the k-space grid kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz); % create the time array t_end = (Nx*dx)*1.5/c0_min; % [s] kgrid.t_array = makeTime(kgrid, c0_min, [], t_end); %% ======================================================================== = % DEFINE THE MEDIUM PARAMETERS % ======================================================================== = % define the properties of the air c0_air = 340; % [m/s] rho0_air = 1.29; % [kg/m^3] % medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] % medium.alpha_power = 1.5; % medium.BonA = 6; medium.sound_speed = c0_air*ones(Nx, Ny, Nz); medium.density = rho0_air*ones(Nx, Ny, Nz); % define the properties of the soft tissue c0_tissue = 1540; % [m/s] rho0_tissue = 1000; % [kg/m^3] std = 0.005; speed_map_tissue = c0_tissue*(ones(Nx, Ny, Nz) + std*randn(Nx, Ny, Nz)); density_map_tissue = rho0_tissue*(ones(Nx, Ny, Nz) + std*randn(Nx, Ny, Nz)); % define the region of tissue region_tissue = zeros(Nx, Ny, Nz); radius = Nx/2 -10; for i = 1:Nz region_tissue(:,:,i) = makeDisc(Nx, Ny, 1, Ny/2, radius); end % define the properties of tissue medium.sound_speed(region_tissue == 1) = speed_map_tissue(region_tissue == 1); medium.density(region_tissue == 1) = density_map_tissue(region_tissue == 1); % show the sound speed map of compound medium, the air and tissue figure imagesc(medium.sound_speed(:,:,44)) %% ======================================================================== = % DEFINE THE INPUT SIGNAL % ======================================================================== = % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = 12e6; % [Hz] tone_burst_cycles = 2; % create the input signal using toneBurst input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,... 'Envelope', 'Gaussian', 'Plot',true); % scale the source magnitude by the source_strength divided by the % impedance (the source is assigned to the particle velocity) input_signal = (source_strength./(c0_tissue*rho0_tissue)).*input_signal; %% ======================================================================== = % DEFINE THE TX TRANSDUCER % ======================================================================== = % the linear array transducer for TX transducer_TX.number_elements = 32; % total number of transducer elements transducer_TX.element_width = 2 ; % width of each element [grid points] transducer_TX.element_length = 20 ; % length of each element [grid points] transducer_TX.element_spacing = 0; % spacing (kerf width) between the elements [grid points] transducer_TX.radius = inf; % radius of curvature of the transducer [m] transducer_width_TX = transducer_TX.number_elements*transducer_TX.element_width ... + (transducer_TX.number_elements - 1)*transducer_TX.element_spacing; % transducer_TX.position: X = 1 transducer_TX.position = round([1, Ny/2 - transducer_width_TX/2, Nz/2 - transducer_TX.element_length/2]); % properties used to derive the beamforming delays transducer_TX.sound_speed = 1540; % sound speed [m/s] transducer_TX.focus_distance = inf; % focus distance [m] transducer_TX.elevation_focus_distance = 20e-3;% focus distance in the elevation plane [m] transducer_TX.steering_angle = 0; % steering angle [degrees] % apodization transducer_TX.transmit_apodization = 'Rectangular'; transducer_TX.receive_apodization = 'Rectangular'; % use all the element to TX transducer_TX.active_elements = ones(transducer_TX.number_elements, 1); % append input signal used to drive the transducer transducer_TX.input_signal = input_signal; % create the transducer using the defined settings transducer_TX_class = makeTransducer(kgrid, transducer_TX); transducer_TX_class.properties transducer_TX_class.plot %% ======================================================================== = % RUN THE SIMULATION % ======================================================================== = % set the input settings input_args = {'PMLInside', false, 'PMLSize', [pml_x_size, pml_y_size, pml_z_size],... 'DataCast', DATA_CAST}; [sensor_data] = kspaceFirstOrder3D(kgrid, medium, transducer_TX_class, transducer_TX_class, input_args{:}); element_data = transducer_TX_class.combine_sensor_data(sensor_data); save element_data.mat element_data</code></pre> radel on "Piezoelectric propagation" http://www.k-wave.org/forum/topic/piezoelectric-propagation#post-5327 Thu, 26 Nov 2015 03:35:19 +0000 radel 5327@http://www.k-wave.org/forum/ Thanks Brad. For the time being I will work with the elastic models, approximating a piezoelectric material. Bradley Treeby on "Piezoelectric propagation" http://www.k-wave.org/forum/topic/piezoelectric-propagation#post-5323 Tue, 24 Nov 2015 15:24:42 +0000 Bradley Treeby 5323@http://www.k-wave.org/forum/ Hi Radel, k-Wave can model the propagation of elastic waves within solid materials (using <code>pstdElastic2D</code> and <code>pstdElastic3D</code>). However, it doesn't include a model of piezoelectricity if that's what you are interested in simulating. Brad. radel on "Piezoelectric propagation" http://www.k-wave.org/forum/topic/piezoelectric-propagation#post-5318 Tue, 24 Nov 2015 03:23:56 +0000 radel 5318@http://www.k-wave.org/forum/ Hi All, Is it possible to simulate acoustic wave propagation in a piezoelectric material (anisotropic-Lithium niobate) using K-wave? Thanks in advance! Regards, Radel Bradley Treeby on "Setting Incident Intensity of Acoustic Plane Wave" http://www.k-wave.org/forum/topic/setting-incident-intensity-of-acoustic-plane-wave#post-4568 Tue, 10 Jun 2014 11:12:40 +0000 Bradley Treeby 4568@http://www.k-wave.org/forum/ Hi Jon, The behaviour you see is expected - what you are seeing is the underlying behaviour of the so-called band-limited interpolant. This is what is actually propagated by the simulation (see section III in <a href="http://www.k-wave.org/papers/2011-Treeby-JASA.pdf">this paper</a> or Sec. 2.8 in the k-Wave manual if you want some more explanation). You only see it for the high-intensity simulation because the default plot scale is [-1, 1], while your pressure magnitude is ~1e5. If you add <code>'PlotScale', [-source_mag, source_mag]</code> to your optional inputs, you should see the same display (in the linear case) irrespective of your chosen intensity. Hope that helps, Brad. jrosenfield on "Setting Incident Intensity of Acoustic Plane Wave" http://www.k-wave.org/forum/topic/setting-incident-intensity-of-acoustic-plane-wave#post-4565 Mon, 09 Jun 2014 19:37:38 +0000 jrosenfield 4565@http://www.k-wave.org/forum/ Hi Brad, You probably recall that you were helping me some time ago to simulate monochromatic plane waves propagating through simple phantoms. I would like to set my time-averaged incident intensity to be 1 W/cm^2 (1e4 W/m^2) in my current simulation. I am including the code for your reference. When I set the incident intensity to be very low (1e-4 W/m^2), the simulation visually appears to run correctly. However, if I set the incident intensity to be the desired 1e4 W/m^2, the simulation does not seem to run correctly. In particular, the movie of the 2D simulation appears to show that the pressure field propagates the entire length of the phantom in one time step, and drastic changes in pressure can be observed. You only have to look at the first few time steps to see this. Could you please run the attached code and let me know what might be causing the apparent instability in the higher-intensity simulation? Thanks for your help. Sincerely, Jon function planewave(distance) % distance defines lesion-to-detector distance (cm) dx = 1e-4; % grid point spacing in the x direction [meters] dy = 1e-4; % grid point spacing in the y direction [meters] Nx = 587; % number of pixels in x-direction Ny = 500; % number of pixels in y-direction kgrid = makeGrid(Nx, dx, Ny, dy); % make computational grid middle = 294; x0 = middle; frequency = 3.81; % define properties of the lesion y0 = round(500 - (distance*(1/100)*(1/dy))); medium.sound_speed = 1550*makeDisc(Nx,Ny,x0,y0,6); % [m/s] medium.density = 1040*makeDisc(Nx,Ny,x0,y0,6); % [kg/m^3] medium.alpha_coeff = .570*(frequency^1.3)*makeDisc(Nx,Ny,x0,y0,6); % [dB/(MHz^y*cm)] medium.alpha_power = 0; % value of power y lesionmask = makeDisc(Nx,Ny,x0,y0,6); medium.sound_speed(medium.sound_speed==0) = 1550; medium.density(medium.density==0) = 1015; medium.alpha_coeff(medium.alpha_coeff==0) = ((0.25*0.158*(frequency^1.7))... + (0.75*0.870*(frequency^1.5))); medium.alpha_power = 0; % value of power y % define sensor sensor.mask = zeros(Nx,Ny); sensor.mask(:,500) = 1; sensor.record = {'p'}; %define source source.p_mask = zeros(Nx,Ny); source.p_mask(:,1) = 1; source_intensity = 1e-4; % [W/m^2] rho0 = 1015; % [kg/m^3] c0 = 1550; % [m/s] source_mag = sqrt(2*source_intensity*rho0*c0); % [Pa] source_freq = frequency*(1e6); % [Hz] % calculate the time length of one acoustic period T = 1/source_freq; % create the time array using a nice number of points points_per_period = 10; dt = T/points_per_period; t_end = (10/source_freq) + (3.4e-5); kgrid.t_array = 0:dt:t_end; % define a time varying sinusoidal source source.p = source_mag*sin(2*pi*source_freq.*kgrid.t_array); % smooth startup using a cosine ramp ramp_length = 40; ramp = (-cos((0:(ramp_length-1))*pi/ramp_length ) + 1)/2; source.p(1:length(ramp)) = source.p(1:length(ramp)) .* ramp; source.p = filterTimeSeries(kgrid, medium, source.p, 'PlotSpectrums', false, 'PPW', 2); num_cycles = 1; sensor.record_start_index = kgrid.Nt - round(num_cycles * T / dt) + 1; medium.sound_speed = smooth(kgrid, medium.sound_speed, true); medium.density = smooth(kgrid, medium.density, true); input_args = {'DisplayMask', (source.p_mask + sensor.mask + lesionmask), 'DataCast', 'single', 'PMLSize', [0, 20], 'PMLAlpha', [0, 2], 'PMLInside', false}; sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); pressure = sensor_data.p; pressure = pressure.^2; pressure = transpose(pressure); pressure = mean(pressure); Iavg = pressure./(rho0*c0); xaxis = -293*0.1:0.1:293*0.1; plot(xaxis, Iavg) xlabel('Lateral Position (mm)') ylabel('Acoustic Intensity (W/m^2)') hold on end