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

TDNL(u,U,X,K,J,c,cutoff,Ppos,Pneg,I_td)
%% Authored by Joshua Soneson 2007, updated 2010
function[u,U,Ppos,Pneg,I_td] = TDNL(u,U,X,K,J,c,cutoff,Ppos,Pneg,I_td)
% converts spectrum to one cycle of the time-domain waveform and 
% integrates the invicid Burger's equation using upwind/downwind 
% method with periodic boundary conditions. TDNL stands for Time
% Domain NonLinear.

% set peak and trough values to zero; in case they are not assigned later
if(K==1)	% linear case - do nothing
else		% nonlinear case - enter loop
  for j=J:-1:1
    % execute nonlinear step only if amplitude is above cutoff;
    % row j=1 is always computed so plots look nice and smooth!
    I_td(j) = 0;
    if(sqrt(u(j,1)*u(j,1) + u(J+j,1)*u(J+j,1))>cutoff | j==1)
      % convert from sin/cos representation to complex exponential 	
      U(2:K+1) = u(j,:) - i*u(j+J,:);		
      U(2*K:-1:K+2) = u(j,1:K-1) + i*u(j+J,1:K-1); 
      U(1) = 0;
      % transform to time domain:
      U = K*real(ifft(U));
      I_td(j) = trapz(U.^2);
      % determine how many steps necessary for CFL<1 (CFL<0.8 to be safe).
      P = ceil(c*max(U)/0.8);
      % Nonlinear integration (upwind/downwind) algorithm 
      for p=1:P
        for k=1:2*K		
          if(U(k)<0)
            if(k==1)
              X(k) = U(k) + c*(U(1)*U(1) - U(2*K)*U(2*K))/P/2;
            else
              X(k) = U(k) + c*(U(k)*U(k) - U(k-1)*U(k-1))/P/2;
            end
          else
            if(k==2*K)
              X(k) = U(k) + c*(U(1)*U(1) - U(k)*U(k))/P/2;
            else
              X(k) = U(k) + c*(U(k+1)*U(k+1) - U(k)*U(k))/P/2;
            end
          end
        end
        U = X;
      end
      I_td(j) = I_td(j) - trapz(X.^2); 
      % transform back to frequency domain:
      Ppos(j) = max(X);
      Pneg(j) = min(X);
      X = fft(X)/K;
      % convert back to sin/cos representation:
      u(j,:) = real(X(2:K+1));
      u(j+J,:) = -imag(X(2:K+1));
    end
  end
end

Contact us