function [twiddle, dda, dem_vel] = get_atbiasp(data, varargin)
% [twiddle, dda, dem_vel] = get_atbiasp(data,varargin)
% detrend amplitude; then look at the product of 
%           "dda"                  "dem_vel"
%      d/dz (demeaned amp) *  (foreward velocity)
% 'data' comes from [data, config] = agetmat('a_'); 
% varargin is name, value pairs:
% name          value    what
% boxlen        15       boxfilter length       
% noise_tol     2        noise floor is the deepest value + this
% pfc           pfc.cfg.bigtarget_amp     25      from aflagit_setup
%               pfc.cfg.bigtarget_offset   5
%               pfc.cfg.bigtarget_avg     2 
%

fprintf('getting along-track bias parameter\n')
tic
cfg.boxlen            = 15;
cfg.noise_tol         = 2;     %tolerance for finding the noise floor
cfg.pfc.cfg.bigtarget_amp = 25;
cfg.pfc.cfg.bigtarget_offset = 5;
cfg.pfc.cfg.bigtarget_avg = 2;

cfg = fillstruct(cfg, varargin, 0);


pfc = cfg.pfc;  %profile flags configuration

if (~isstruct(data))
   error('must pass in data and config structures from agetmat')
end

%take out a linear fit to log10(amp) for amp greater than bottom amp+count
[bm1, b] = bigtarget(data.amp, pfc.cfg.bigtarget_amp, ...
                     pfc.cfg.bigtarget_offset, ...
                     pfc.cfg.bigtarget_avg);
okalldata_mask =  ~bottommask(bm1, b, 0.85, 1); %zero out the bottom; 1 for data


%need to find the value at the bottom of each profile
%- first propagate deepest value downward: do this by  appending a layer of zeros 
%     at the bottom and finding the location where diff==1; 
%     then propagate the value from there downward.
tempmask = okalldata_mask;
tempmask = [tempmask; 0*tempmask(1,:)];  %1 for data; 0 for bottom 
bottomi = find(diff(tempmask) == -1); %diff (N+1) - Nth value

noisefloorM = tempmask(1:end-1,:).*data.amp;     %ones where the data are; 0 else
    noisefloor = noisefloorM(bottomi);             %should be 1xNprofs
    noisefloorM = ones(size(data.bins(:))) * noisefloor(:)';
floordatai = find(data.amp < noisefloorM + cfg.noise_tol);
okalldata_mask(floordatai) = 0;

%now take out a linear fit to log10(amp profile above noise floor)
dem_amp = zeros(size(data.amp));
for profi = 1:length(dem_amp(1,:))
   mask = okalldata_mask(:,profi);

   if (length(find(mask==1)) > 2) 
      xx = data.bins(mask)';
      yy = log10(data.amp(mask,profi));
      coef = polyfit(xx, yy, 1);
      trend = polyval(coef, xx);
   else
      trend = 0;
   end
   res = 0*data.amp(:,profi); 
   res(mask) = data.amp(mask, profi) - 10.^(trend);
   dem_amp(:,profi) =  res - boxfilt(res, cfg.boxlen);
end


dem_vel = zeros(size(data.fabs_sm));
for profi = 1:length(data.fabs_sm(1,:))
   dem_vel(:,profi) = data.fabs_sm(:,profi) - ...
       boxfilt(data.fabs_sm(:,profi), cfg.boxlen);
end

dda =  diff(dem_amp);
twiddle = sign(dda).*dda.*(table1([data.bins', dem_vel] , mid(data.bins')));
twiddle = [twiddle; NaN*twiddle(1,:)];


%clg
%for profi = 1:length(dda(1,:))
%   h=plot([0 0]',[-50 0]','color',.5*[1 1 1], 'linewidth',2);
%   hold on
%   h=plot(dda(:,profi), -mid(data.bins),'k');lw(h,2)
%   h=plot(dem_vel(:,profi), -data.bins,'g');lw(h,2)
%   xl(-15,15)
%   paused
%   clg
%end


toc