Code covered by the BSD License  

Highlights from
High intensity focused ultrasound simulator

High intensity focused ultrasound simulator

by

 

Simulates high intensity focused ultrasound beams and heating effects in layered media

BHT_plots(z,r,t,t_peak,Tpeak,Tmax_mat,Dmat,T0)
%% Authored by Joshua Soneson 2007, updated 2010
function[] = BHT_plots(z,r,t,t_peak,Tpeak,Tmax_mat,Dmat,T0)
% produces the following plots:
%	temperature contours when peak temp is reached
%	cumulative thermal dose contours on log scale
%	peak temp vs. time 

% cumulative thermal dose near focus:
x = max(max(Dmat));
P = floor(log10(x/2.4));
figure
axes('FontSize',18)
r2 = [-flipud(r(1:length(r)-1));r];
Dmat2 = [flipud(Dmat(1:length(r)-1,:));Dmat];
[C,H]=contour(z,r2,Dmat2,[.1 .1]);
n=length(C(1,:));
fill(C(1,2:n),C(2,2:n),'w');
xlim = get(gca,'XLim');
ylim = get(gca,'YLim');
[C,H]=contour(z,r2,Dmat2,[240 240],'r');
n=length(C(1,:));
fill(C(1,2:n),C(2,2:n),'r')
%hold on
%m=0.24;
%for p=1:P
%  m = [m,m*10^p];
%end
%contour(z,r,Dmat,m,'k','LineWidth',2);
axis([xlim(1),xlim(2),ylim(1),ylim(2)])
xlabel('z (cm)')
ylabel('r (cm)')
title('Thermal Dose (CEM43C)')
grid on

% temperature distribution near focus:
figure
peak_temp = max(Tmax_mat(1,:));
peak_index = find(Tmax_mat(1,:)==peak_temp);
peak_distance = 0.01*round(100*z(peak_index));
axes('FontSize',18)
Tmax_mat2 = [flipud(Tmax_mat(1:length(r)-1,:));Tmax_mat];
[C,H]=contour(z,r2,Tmax_mat2+T0,12);
%[C,H]=contour(z,r2,Tmax_mat+T0,[39.8 39.8]);
n=length(C(1,:));
if(n<2) else
  C(1,1) = C(1,2);
  C(2,1) = C(2,2);
  for nn=1:n-1
    if(abs(C(1,nn+1)-C(1,nn))>1) C(:,nn+1) = C(:,nn); end  
  end
  fill(C(1,:),C(2,:),'w');
end
xlim = get(gca,'XLim');
ylim = get(gca,'YLim');
%v = 40:5:100;
[C,H]=contour(z,r2,Tmax_mat2+T0,12,'LineWidth',2);
axis([xlim(1),xlim(2),ylim(1),ylim(2)])
%caxis([40 100])
colorbar;
xlabel('z (cm)')
ylabel('r (cm)')
title(['Temp (C), @ t = ',num2str(t_peak),' s, z = ',num2str(peak_distance),' cm'])
grid on

% peak temp vs. time:
figure
axes('FontSize',18)
plot(t,Tpeak+T0,'r','LineWidth',2)
ylim = get(gca,'YLim');
axis([0,t(length(t)),T0,ylim(2)])
xlabel('t (sec)')
ylabel('Temp (C)')
grid on

Contact us