function echo_info = publish_echogram(param,mdata,lay)
% echo_info = publish_echogram(param,mdata,lay)
%
% param: structure controlling the plotting of the echogram
%  .fig_hand: figure handle to use for plot
%  .num_x_tics: number of xtics across the bottom (the xtics give
%     the geographic info), usually 4 for figure copy and 6 for
%     full page print out
%  .depth: 1x2 double vector representing the limits of the depth
%     axis to plots (-500 to 3500 typical, deep antarctica -500 to 4000,
%     but for custom plots -100 to 2000 or whatever highlights the data
%     best
%  .elev_comp: if true, the echogram is elevation compensated so that
%     the plotted echogram is as if the aircraft flew at a constant
%     elevation
%  .er_ice: the permittivity used to create the depth axis, default is 
%     loaded from physical_constants.
% echo_info: structure with handles to plots/axes
%  .ah_echo_time: time axis handle
%  .ah_echo: depth axis handle
%  .h_surf: surface plot handle
%  .h_bot: bottom plot handle
%
% Example: see run_publish_echogram
%
% Author: John Paden
%
% See also: publish_map, run_publish_echogram

physical_constants;
if ~isfield(param,'er_ice')
  param.er_ice = er_ice;
end

% =======================================================================
% Create Surface and Bottom Variables
% =======================================================================
lay.Bottom  = lay.layerData{2}.value{2}.data;
lay.Surface = lay.layerData{1}.value{2}.data;
lay.Thickness = lay.Bottom-lay.Surface;
neg_idxs = find(lay.Thickness < 0 & isfinite(lay.Thickness));
if ~isempty(neg_idxs)
  warning('  Negative thickness detected');
  lay.Bottom(neg_idxs) = lay.Surface(neg_idxs);
end

% Interpolate layer onto GPS times of echogram
warning('off','MATLAB:interp1:NaNinY');
lay.Surface = interp1(lay.GPS_time,lay.Surface,mdata.GPS_time,'linear','extrap');
lay.Bottom = interp1(lay.GPS_time,lay.Bottom,mdata.GPS_time,'linear','extrap');
warning('on','MATLAB:interp1:NaNinY');

% Create fast-time correction vector
%   The layer file may contain a different elevation profile
%   than the data file. The surface and bottom need to be adjusted
%   to account for these differences.
elev_interp    = interp1(lay.GPS_time,lay.Elevation,mdata.GPS_time,'linear','extrap');
fast_time_correction = (mdata.Elevation - elev_interp)/(c/2);

% =======================================================================
% Elevation compensation
% =======================================================================

if param.elev_comp == 1
  max_elev = max(mdata.Elevation);
  dRange = max_elev - mdata.Elevation;
  dt = mdata.Time(2)-mdata.Time(1);
  dBins = round(dRange / (c/2) / dt);
  dtime = dRange/(c/2);
  zero_pad_len = max(abs(dBins));
  mdata.Data = cat(1,mdata.Data,zeros(zero_pad_len,size(mdata.Data,2)));
  mdata.Time = mdata.Time(1) + (mdata.Time(2)-mdata.Time(1)) * (0:1:size(mdata.Data,1)-1);
  
  for rline = 1:size(mdata.Data,2)   
    mdata.Data(:,rline) = interp1(mdata.Time, mdata.Data(:,rline), mdata.Time - dtime(rline), 'linear',0);
    mdata.Elevation(rline) = mdata.Elevation(rline) + dRange(rline);
    lay.Surface(rline) = lay.Surface(rline) + dtime(rline);
    lay.Bottom(rline) = lay.Bottom(rline) + dtime(rline);
  end
elseif param.elev_comp == 2
  % Compensate for surface variations
  % Choose filtering method for surface
  if 0
    surf_filt = medfilt1(lay.Surface,101);
  elseif 1
    if length(lay.Surface) >= 100
      [Bfilt,Afilt] = butter(2,0.02);
      surf_filt = filtfilt(Bfilt,Afilt,lay.Surface);
      tmp = polyval(polyfit(1:51,reshape(lay.Surface(1:51),[1 51]),2),1:51);
      surf_filt(1:50) = tmp(1:50) - tmp(51) + surf_filt(51);
      tmp = polyval(polyfit(1:51,reshape(lay.Surface(end-50:end),[1 51]),2),1:51);
      surf_filt(end-49:end) = tmp(2:51) - tmp(1) + surf_filt(end-50);
    else
      surf_filt = medfilt1(lay.Surface,round(length(lay.Surface)/20)*2+1);
    end
  elseif 0
    surf_filt = polyval(polyfit(1:length(lay.Surface),lay.Surface,4),1:length(lay.Surface));
  else
    surf_filt = lay.Surface;
  end
  lay.Surface = surf_filt;
  mdata.Elevation = mean(mdata.Elevation) + (surf_filt - mean(surf_filt)) * c/2;  
  
  max_elev = max(mdata.Elevation);
  dRange = max_elev - mdata.Elevation;
  dt = mdata.Time(2)-mdata.Time(1);
  dBins = round(dRange / (c/2) / dt);
  dtime = dRange/(c/2);
  zero_pad_len = max(abs(dBins));
  mdata.Data = cat(1,mdata.Data,zeros(zero_pad_len,size(mdata.Data,2)));
  mdata.Time = mdata.Time(1) + (mdata.Time(2)-mdata.Time(1)) * (0:1:size(mdata.Data,1)-1);
  
  for rline = 1:size(mdata.Data,2)   
    mdata.Data(:,rline) = interp1(mdata.Time, mdata.Data(:,rline), mdata.Time - dtime(rline), 'linear',0);
    mdata.Elevation(rline) = mdata.Elevation(rline) + dRange(rline);
    lay.Surface(rline) = lay.Surface(rline) + dtime(rline);
    lay.Bottom(rline) = lay.Bottom(rline) + dtime(rline);
  end
end

if isfield(param,'detrend') && ~isempty(param.detrend.mode) && param.detrend.mode == 1
  surf_bins = interp1(mdata.Time,1:length(mdata.Time),lay.Surface);
  % mdata.Data = mdata.Data ./ repmat(mean(mdata.Data,2),[1 size(mdata.Data,2)]);
  % mdata.Data = hanning(size(mdata.Data,1))*ones(.*mdata.Data
  % mdata.Data = mdata.Data .* repmat(detwinhanning(size(mdata.Data,1)),[size(mdata.Data,1) 1]);
  
%   wdep = 4000;
%   wind = [hanning(wdep) ; 0*ones(size(mdata.Data,1)-wdep,1)];
%   wind  = wind*ones(1,size(mdata.Data,2));
%   surf_bins = round(surf_bins);
%   shift = (surf_bins-80+wdep);
%   for ii=1:size(wind,2)
%     wind_tmp = flipud(wind(:,ii));
%     wind_tmp = circshift(wind_tmp,shift(ii));
%     wind(:,ii) = wind_tmp;
%   end
%   mdata.Data = mdata.Data.*wind+0.0001;
end

% Create depth axis
good_surface_vals = lay.Surface + fast_time_correction;
good_surface_vals = good_surface_vals(isfinite(good_surface_vals));
mean_surface_time = mean(good_surface_vals);
if ~isfinite(mean_surface_time)
  mean_surface_time = 0;
end

Depth = (mdata.Time-mean_surface_time)*c/2/sqrt(param.er_ice);

DBottom = lay.Bottom - mean_surface_time + fast_time_correction;
DBottom = DBottom*c/2/sqrt(param.er_ice);
DSurface = lay.Surface - mean_surface_time + fast_time_correction;
DSurface = DSurface*c/2/sqrt(param.er_ice);
Surface_Depth = DSurface;

% Limit depths according to input param.depth
depth_range = eval(param.depth);
depth_good_idxs = find(Depth >= depth_range(1) & Depth <= depth_range(end));

% Create echogram plot
figure(param.fig_hand); clf;

ah_echo_time = axes; % 2-way travel time axis
ah_echo = axes; % Depth axis with data
plot(ah_echo_time,zeros(size(mdata.Time(depth_good_idxs))),mdata.Time(depth_good_idxs)*1e6,'');
axis(ah_echo_time,'tight')
echogram_vals = lp(mdata.Data(depth_good_idxs,:));
imagesc([],Depth(depth_good_idxs),echogram_vals);  
set(ah_echo_time,'XTickLabel','');
set(ah_echo_time,'YAxisLocation','Right');
set(ah_echo_time,'YDir','Reverse');
xtl = create_standard_x_labels(mdata.Latitude,mdata.Longitude,mdata.Elevation,param.num_x_tics);
add_x_labels(ah_echo,xtl,{'distance','latitude','longitude'})
ylabel(ah_echo,sprintf('depth, e_r = %.2f (m)', param.er_ice));
ylabel(ah_echo_time,'2-way travel time (us)','VerticalAlignment','bottom','Rotation',270)
if isfield(param,'colormap') && ~isempty(param.colormap)
  colormap(param.colormap);
else
  colormap(1-gray);
end
if isfield(param,'caxis') && ~isempty(param.caxis)
  echogram_vals = sort(echogram_vals(round(linspace(1,numel(echogram_vals),1000))));
  caxis(echogram_vals(1+round(param.caxis*999)));
end

echo_title = title(ah_echo,sprintf('Data Frame ID: %s', param.frm_id),'Interpreter','none');

% Plot surface and bottom layer
hold on
echo_info.h_surf = plot(ah_echo,DSurface,'--m');
echo_info.h_bot = plot(ah_echo,DBottom,'--r');
hold off

echo_info.ah_echo_time = ah_echo_time;
echo_info.ah_echo = ah_echo;

return;

