Code covered by the BSD License

# High intensity focused ultrasound simulator

### Josh Soneson (view profile)

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
```