http://www.k-wave.org/forum/topic/setting-incident-intensity-of-acoustic-plane-wave Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:31:45 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php 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/ <p>Hi Jon,</p> <p>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). </p> <p>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>&#39;PlotScale&#39;, [-source_mag, source_mag]</code> to your optional inputs, you should see the same display (in the linear case) irrespective of your chosen intensity.</p> <p>Hope that helps,</p> <p>Brad. </p> 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/ <p>Hi Brad,</p> <p>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?</p> <p>Thanks for your help.</p> <p>Sincerely,</p> <p>Jon</p> <p>function planewave(distance)</p> <p>% distance defines lesion-to-detector distance (cm)</p> <p>dx = 1e-4; % grid point spacing in the x direction [meters]<br /> dy = 1e-4; % grid point spacing in the y direction [meters]<br /> Nx = 587; % number of pixels in x-direction<br /> Ny = 500; % number of pixels in y-direction<br /> kgrid = makeGrid(Nx, dx, Ny, dy); % make computational grid<br /> middle = 294;<br /> x0 = middle;<br /> frequency = 3.81;</p> <p>% define properties of the lesion<br /> y0 = round(500 - (distance*(1/100)*(1/dy)));<br /> medium.sound_speed = 1550*makeDisc(Nx,Ny,x0,y0,6); % [m/s]<br /> medium.density = 1040*makeDisc(Nx,Ny,x0,y0,6); % [kg/m^3]<br /> medium.alpha_coeff = .570*(frequency^1.3)*makeDisc(Nx,Ny,x0,y0,6); % [dB/(MHz^y*cm)]<br /> medium.alpha_power = 0; % value of power y<br /> lesionmask = makeDisc(Nx,Ny,x0,y0,6);</p> <p>medium.sound_speed(medium.sound_speed==0) = 1550;<br /> medium.density(medium.density==0) = 1015;<br /> medium.alpha_coeff(medium.alpha_coeff==0) = ((0.25*0.158*(frequency^1.7))...<br /> + (0.75*0.870*(frequency^1.5)));<br /> medium.alpha_power = 0; % value of power y</p> <p>% define sensor<br /> sensor.mask = zeros(Nx,Ny);<br /> sensor.mask(:,500) = 1;<br /> sensor.record = {'p'};</p> <p>%define source<br /> source.p_mask = zeros(Nx,Ny);<br /> source.p_mask(:,1) = 1;<br /> source_intensity = 1e-4; % [W/m^2]<br /> rho0 = 1015; % [kg/m^3]<br /> c0 = 1550; % [m/s]<br /> source_mag = sqrt(2*source_intensity*rho0*c0); % [Pa]<br /> source_freq = frequency*(1e6); % [Hz]</p> <p>% calculate the time length of one acoustic period<br /> T = 1/source_freq;</p> <p>% create the time array using a nice number of points<br /> points_per_period = 10;<br /> dt = T/points_per_period;<br /> t_end = (10/source_freq) + (3.4e-5);<br /> kgrid.t_array = 0:dt:t_end;</p> <p>% define a time varying sinusoidal source</p> <p>source.p = source_mag*sin(2*pi*source_freq.*kgrid.t_array);</p> <p>% smooth startup using a cosine ramp<br /> ramp_length = 40;<br /> ramp = (-cos((0:(ramp_length-1))*pi/ramp_length ) + 1)/2;<br /> source.p(1:length(ramp)) = source.p(1:length(ramp)) .* ramp;</p> <p>source.p = filterTimeSeries(kgrid, medium, source.p, 'PlotSpectrums', false, 'PPW', 2);</p> <p>num_cycles = 1;<br /> sensor.record_start_index = kgrid.Nt - round(num_cycles * T / dt) + 1;</p> <p>medium.sound_speed = smooth(kgrid, medium.sound_speed, true);<br /> medium.density = smooth(kgrid, medium.density, true);</p> <p>input_args = {'DisplayMask', (source.p_mask + sensor.mask + lesionmask), 'DataCast', 'single', 'PMLSize', [0, 20], 'PMLAlpha', [0, 2], 'PMLInside', false};</p> <p>sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});</p> <p>pressure = sensor_data.p;<br /> pressure = pressure.^2;<br /> pressure = transpose(pressure);<br /> pressure = mean(pressure);<br /> Iavg = pressure./(rho0*c0);<br /> xaxis = -293*0.1:0.1:293*0.1;<br /> plot(xaxis, Iavg)<br /> xlabel('Lateral Position (mm)')<br /> ylabel('Acoustic Intensity (W/m^2)')<br /> hold on<br /> end </p>