%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Make sound field snapshots and movies using the
% DREAM toolbox.
%
%
% Fredrik Lingvall : 2005-10-25
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if findstr(computer,'linux') && exist('OCTAVE_VERSION')
  disp('***************************************************************');
  disp('*');
  disp('* Your are running Octave on a Linux platform!');
  disp('* You need MPlayer/MEncoder and Gnuplot 4.x ');
  disp('* to run this script!');
  disp('*');
  disp('***************************************************************');
  pause(4);

  % Detect the number of CPUs/Cores on the system.
  [dummy,tmp_str]= system('cat /proc/cpuinfo | grep model | grep name');
  n_cpus = size(strfind(tmp_str,'model'),2);
  fprintf('\n*** Detected a %d cpu system ***\n\n',n_cpus);
else
  n_cpus = 1;
end

% Paths to eps, png, and jpegs.
eps_path = '../eps/';
png_path = 'png/';


eps_file = '50_m20';

% Switches
COMPUTE_SIR = 1;
GENERATE_EPS = 0;
MAKE_MOVIE = 1;
AUTO_SCALE = 0;


% 25 MHz
Ts = 0.04;

if COMPUTE_SIR

  % Soundspeed in water (in [m/s]).
  cp = 1500;

  % Length of A-scans.
  K = 1000;

  %
  % Define the array.
  %

  % Spacing between array elements (array pitch).
  arr_samp = 1;                         % [mm]

  % Element size.
  a = 0.9;                              % [mm]
  b = 1.0;                              % [mm]

  % Transmit.
  lt = 8;
  xt =  arr_samp*(-lt+0.5:lt-0.5);
  trans = [xt; a*ones(1,2*lt);];

  fprintf('Array aperture = %1.1f [mm]\n',xt(length(xt))-xt(1));


  % Number of transducer elements.
  L = size(xt,2);

  % Time vector.
  us = (0:Ts:Ts*40000);

  % Max width of ROI.
  x_max = 25;

  idx_z_min = 1;
  idx_z_max = 918;

  % Spatial sampling resolution.
  spat_samp_x = 0.3; % Horizontal resolution [mm].
  spat_samp_z = 0.3; % Vertical resoltion [mm].

  %
  % Horizontal grid.
  %
  x = -x_max:spat_samp_x:x_max;
  N = length(x);

  %
  % Vertical grid.
  %
  z_min = us(idx_z_min)*cp/1e3;
  z_max = us(idx_z_max)*cp/1e3;
  z = (z_min:spat_samp_z:z_max);
  M = length(z);

  %
  % Observation grid matrix.
  %
  [X,Z] = meshgrid(x,z);
  Y = zeros(size(X));
  xo = X(:);
  yo = Y(:);
  zo = Z(:);
  Ro = [xo yo zo];

  % Number of array elements.
  L = size(trans,2);

  % Number of observation points
  MN = size(Ro,1);

  fprintf('z_min = %1.1f [mm]\n',z_min);
  fprintf('z_max = %1.1f [mm]\n',z_max);

  fprintf('KL x MN = %d x %d \n',K*L,M*N);


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % Electrical impulse response.
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  he_nt=50;
  t = (0:((he_nt-1)))*Ts;

  f0 = 3;                             % Center frequency [MHz].
  fprintf('\nf = %1.2f [MHz]\n',f0);
  lambda = cp/f0/1e3; % [mm].
  fprintf('lambda = %1.2f [mm]\n',lambda);

  t0 = 0.55;                 % Time delay to max amplitude [mus].
  a_n = 20;                        % Envelop parameter.
  h_e = -exp(-a_n.*(t-t0).^2).*cos(2.*pi.*f0.*t);

  % Number of samples to max amplitude of pulse.
  [mx,p_delay] = max(h_e);

  f_tmp = abs(freqz(h_e,1,1024));
  h_e = h_e/max(f_tmp); % Unity gain at center freq.
  f_e = abs(freqz(h_e,1,1024));
  f = (0:1023)/1045/Ts/2;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % DREAM parameters.
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % Descretization parameters.
  dx = 0.05; 				% [mm]
  dy = 0.05; 				% [mm]
  dt = Ts; 				% [us]
  nt = 1300;   				% Length of single-path spatial impulse response vector.
  s_par = [dx dy dt nt];


  fprintf('\nDREAM length nt = %d\n\n',nt);

  % Material parameters.
  v     = 1.0; 				% Normal velocity.
  alfa  = 0; 				% Absorbtion (dB/cm Hz).
  m_par = [v cp alfa];

  % Element size for transmit [mm].
  geom_par = [a b];

  % Grid function (position vectors of the transmit elements).
  gx = trans(1,:);
  gx = gx(:);
  gy = zeros(length(gx),1);
  gz = zeros(length(gx),1);
  G = [gx gy gz];

  %
  % User defined focusing.
  %

  foc_dep = 50;
  theta = -20;
  foc_met = 'ud';

  xp = foc_dep*tan(theta*pi/180)*ones(L,1);
  zp = foc_dep*ones(L,1);

  r_max = 0;
  for l=1:L
    tmp_max = norm([xp(l) zp(l)]' - [gx(l) 0]');
    if tmp_max > r_max
      r_max = tmp_max;
    end
  end

  focal = zeros(L,1);
  for n=1:L
    focal(n) = (r_max - norm([xp(n) zp(n)]' - [gx(n) 0]'))*1000/cp;
  end

  fprintf('Array focused at  = %1.1f [mm]\n',foc_dep);
  fprintf('Array steered = %1.1f [deg]\n',theta);

  % Beam steering
  steer_met = 'off'; % Done above by foc_met = 'ud' above.

  phi = 0; 				% Angle in y-direction.
  steer_par = [theta phi];              % Not used here.

  apod_met = 'off';

  % Apodization
  apod_met = 'off';
  apod = ones(length(gx),1); 		% Not used here.
  win_par = 1; 				% Not used here.


  z_min = min(Ro(:,3));                 % Depth to closest observation point.
  delay = z_min/cp*1000;                % Single path delay.

  mx_td = max(focal);
  fprintf('----------------------------------------------------------------------\n');
  fprintf(' Max steer delay = %1.2f [us] (= %1.2f [mm], = %d [samples])\n',mx_td, mx_td*1000/cp,...
          round(mx_td/Ts));
  fprintf('----------------------------------------------------------------------\n');

  %
  % Compute transmit SIRs.
  %

  eval('tic');

  if n_cpus == 1
    % Serial.
    H_t = dream_arr_rect(Ro,geom_par,G,s_par,delay,m_par,foc_met,...
                         focal,steer_met,steer_par,apod_met,apod,win_par,'stop');
  else
    % Parallel.
    H_t = dream_arr_rect_p(Ro,geom_par,G,s_par,delay,m_par,foc_met,...
                           focal,steer_met,steer_par,apod_met,apod,win_par,n_cpus,'stop');
  end

  t = toc;
  if t > 60
    fprintf('elapsed time: dream_arr_rect = %1.1f [min]\n',t/60);
  else
    fprintf('elapsed time: dream_arr_rect = %1.1f [s]\n',t);
  end

  %
  % Compute pressure response (i.e., convolve with electro-acoustical impulse response).
  %

  eval('tic');

  %P_t = conv_p(H_t,h_e,n_cpus);
  P_t = fftconv_p(H_t,h_e,n_cpus); % If you have fftw 3.x on your system.

  t = toc;
  if t > 60
    fprintf('elapsed time: fftconv_p = %1.1f [min]\n',t/60);
  else
    fprintf('elapsed time: fftconv_p = %1.1f [s]\n',t);
  end
  fprintf('----------------------------------------------------------------------\n');
end

l_z = length(z);

if exist('OCTAVE_VERSION')
  figure(1);
  clf
else
  fig=figure(1);
  set(fig,'DoubleBuffer','on');
  set(gca,'xlim',[min(x) max(x)],'ylim',[z(1) z(l_z)],...
          'NextPlot','replace','Visible','off')

end

kk = 0;
for k=1:3:1100
  kk = kk+1;

  fprintf('Snapshot taken at %1.1f [us]\n', delay+k*Ts)

  ds = Ts*cp/1000;
  u_samp = ceil(spat_samp_z/ds); % n times undersampled.

  if u_samp > 1
    % Take mean of u_samp samples.
    if (round(k-u_samp/2) >=1)
      o = sum(abs(hilbert(P_t(round(k-u_samp/2:k+u_samp/2),:))))'/u_samp;
    else
      o = sum(abs(hilbert(P_t(round(k:k+u_samp/2),:))))'/u_samp;
    end
    O = reshape(o,M,N);
  else
    o = P_t(k,:);
    O = reshape(o,M,N);
  end

  mx = max(abs(o));

  if AUTO_SCALE
    O = 64*abs(O)/mx;
    Im = abs(O);
  else
    g = 0.7;
    Im = g*abs(O);
  end

  if exist('OCTAVE_VERSION')
    image(x,z,Im);
    axis('ij')
    xlabel('x [mm]');
    ylabel('z [mm]');
    legend('off');
  else
    h = image(x,z,Im);
    set(gca,'FontSize',16);
    xlabel('x [mm]','FontSize',20);
    ylabel('z [mm]','FontSize',20);
    set(gca,'YTick', 0:10:80)
    drawnow
  end
  grid('off');

  t_str = sprintf('Wavefield Snapshot at %1.1f [us]', delay+k*Ts);
  title(t_str);

  if MAKE_MOVIE

    if exist('OCTAVE_VERSION')
      print(sprintf ('%s/snapshot_%.5d.png', png_path, k));

      %if k<10
      %  eval(['print("' png_path 'snapshot_000' num2str(k) '.png", "-dpng");']);
      %elseif k<100
      %  eval(['print("' png_path 'snapshot_00' num2str(k) '.png", "-dpng");']);
      %elseif k<1000
      %  eval(['print("' png_path 'snapshot_0' num2str(k) '.png", "-dpng");']);
      %else
      %  eval(['print("' png_path 'snapshot_' num2str(k) '.png", "-dpng");']);
      %end
    else % Matlab
      % Add a frame to Matlab movie.
      drawnow;
      snapshots(kk) = getframe;
    end
  end

end % for

if MAKE_MOVIE

  if findstr(computer,'linux') && exist('OCTAVE_VERSION')
    system(['mencoder -quiet mf://png/*.png -mf w=640:h=480:fps=' ...
            '25:type=png -ovc lavc -lavcopts vcodec=wmv1 -oac copy -' ...
            'o snapshots.avi']);
  else
    % Create an AVI movie from MATLAB movie.
    movie2avi(snapshots,'snapshots.avi');
  end

end

t = toc;
if t > 60
  fprintf('elapsed time = %1.1f [min]\n',t/60);
else
  fprintf('elapsed time = %1.1f [s]\n',t);
end

if GENERATE_EPS
  t_str = sprintf('%1.1f', delay+k*Ts);
  if exist('OCTAVE_VERSION')
    eval(['print("' eps_path 'snapshot_' eps_file '_'  t_str '.eps", "-depsc");']);
  else
    eval(['print -deps ' eps_path 'snapshot_' eps_file '_' t_str '.eps']);
  end
end

%
% Electrical impulse response.
%

figure(2);
clf
plot(us(1:length(h_e)),h_e/max(h_e));
grid('on')
axis([0 max(us(length(h_e))) -1.2 1.2])

title('Pulse');

if exist('OCTAVE_VERSION')
  legend('off');
  xlabel('t [us]');
  ylabel('Normalized Amplitude');
else
  xlabel('t [\mus]','FontSize',20);
  ylabel('Normalized Amplitude','FontSize',20);
  set(gca,'FontSize',16);
end

if GENERATE_EPS
  if exist('OCTAVE_VERSION')
    eval(['print("' eps_path 'imp_resp.eps", "-depsc");']);
  else
    eval(['print -depsc ' eps_path 'imp_resp.eps']);
  end
end

figure(3);
clf
plot(f,f_e);
grid('on')
axis([0 max(f) 0 1.2])

if exist('OCTAVE_VERSION')
  legend('off');
  xlabel('f [MHz]');
  ylabel('Normalized Amplitude');
else
  xlabel('f [MHz]','FontSize',20);
  ylabel('Normalized Amplitude','FontSize',20);
  set(gca,'FontSize',16);
end

title('Pulse Spectrum');

if GENERATE_EPS
  if exist('OCTAVE_VERSION')
    eval(['print("' eps_path 'freq_imp_resp.eps", "-depsc");']);
  else
    eval(['print -depsc ' eps_path 'freq_imp_resp.eps']);
  end
end