function adcpcal1(config)
% adcpcal1(config)
% called by adcpcal (the stub copied to the cal/watertrk directory)

% all this used to be in the function called "adcpcal" in cal/watertrk.
% - moved most of it here in this function
% JH 2001/01/09
% EF 2001/10/24 debugged, along with the stub.

%passed from adcpcal
clip_ph   = config.clip_ph  ;
clip_amp  = config.clip_amp ;
clip_var  = config.clip_var ;
clip_dt   = config.clip_dt  ;

%these were the optional arguments from adcpcal
if (isfield(config, 'trange')),   trange =  config.trange ; end
if (isfield(config, 'dvrange')),  dvrange = config.dvrange ; end

if isstr(config.cal_),
   cal = read_asc(config.cal_, 11);
else
   cal = config.cal_;
end

if (isempty(cal))
   error('no data in watertrack')
end



[m, n] = size(cal);
j_ph = n;               % column index of phase
j_amp = n-1;            % column index of amplitude
j_t = 2;                % column index of time
j_du = 5;               % column index of delta-u
j_dv = 6;               % column index of delta-v
j_dt = n-2;             % column index of time shift
j_mvar = n-3;           % column index of variance
j_var = n-4;           % column index of variance

% Fri  12-15-1989  changed from min var to variance

if  (isempty(trange))
   trange(1) = min(cal(:, j_t));
   trange(2) = max(cal(:, j_t));
end
if isempty(dvrange)
   dvrange = [-100 100 -100 100]; % default: liberal ranges in m/s
end


i0 = min(find(cal(:, j_t) >= trange(1)));
i1 = max(find(cal(:, j_t) <= trange(2)));
range = i0:i1;
mr = length(range);

calr = cal(range, :);

% --> Find the raw medians for use in editing:
rmed_ph  = median(calr(:, j_ph));
rmed_amp = median(calr(:, j_amp));
rmed_dt  = median(calr(:, j_dt));

% --> Edit by finding the row index range for which amp and phase
%     are good.
one_col = ones(mr, 1);

a_good = (  ((calr(:, j_amp) - rmed_amp) < ( one_col * clip_amp)) ...
            & ((calr(:, j_amp) - rmed_amp) > (-one_col * clip_amp)));

p_good = (  ((calr(:, j_ph) - rmed_ph)   < ( one_col * clip_ph)) ...
            & ((calr(:, j_ph) - rmed_ph)   > (-one_col * clip_ph)));

dt_good = (  ((calr(:, j_dt) - rmed_dt)   < ( one_col * clip_dt)) ...
             & ((calr(:, j_dt) - rmed_dt)   > (-one_col * clip_dt)));

var_good = calr(:, j_mvar)  <  one_col * clip_var; %min var for WECOMA

dv_good = calr(:,j_du) > dvrange(1) & calr(:,j_du) < dvrange(2) ...
          & calr(:,j_dv) > dvrange(3) & calr(:,j_dv) < dvrange(4);

igood = find( p_good & a_good & dt_good & var_good & dv_good);

adcpcal2(calr, config.label, igood, j_t, j_ph, j_amp, j_dt);

% --> Statistics for edited amp and phase
emed_ph = median(calr(igood, j_ph));
emean_ph = mean(calr(igood, j_ph));
estd_ph = std(calr(igood, j_ph));
emed_amp = median(calr(igood, j_amp));
emean_amp = mean(calr(igood, j_amp));
estd_amp = std(calr(igood, j_amp));
emed_dt = median(calr(igood, j_dt));
emean_dt = mean(calr(igood, j_dt));
estd_dt = std(calr(igood, j_dt));
emed_mvar = median(calr(igood, j_mvar));
emean_mvar = mean(calr(igood, j_mvar));
estd_mvar = std(calr(igood, j_mvar));
emed_var = median(calr(igood, j_var));
emean_var = mean(calr(igood, j_var));
estd_var = std(calr(igood, j_var));
emed_du = median(calr(igood, j_du));
emean_du = mean(calr(igood, j_du));
estd_du = std(calr(igood, j_du));
emed_dv = median(calr(igood, j_dv));
emean_dv = mean(calr(igood, j_dv));
estd_dv = std(calr(igood, j_dv));
n_edit = length(igood);

tmean = mean(calr(igood, j_t));
p = polyfit(calr(igood, j_t)-tmean, calr(igood, j_ph), 1);
pa = polyfit(calr(igood, j_t)-tmean, calr(igood, j_amp), 1);


FID = fopen('adcpcal.out','a');

fprintf(FID, [config.label '\n   Time range %6.2f to %6.2f'], ...
        trange(1), trange(2));
fprintf(FID, ['\n   Calculation done at ' pctime ]);
fprintf(FID, '\n   delta-u min = %6.2f, max = %6.2f;', ...
        dvrange(1), dvrange(2));
fprintf(FID, '\n   delta-v min = %6.2f, max = %6.2f',    ...
        dvrange(3), dvrange(4));
fprintf(FID, '\n   clip_amp = %4.2f,  clip_ph = %4.1f', ...
        clip_amp, clip_ph);
fprintf(FID, '\n   clip_dt = %4.0f,  clip_var = %5.3f', ...
        clip_dt, clip_var);
fprintf(FID, '\nNumber of edited points: %3.0f out of %3.0f', ...
        n_edit, mr)
fprintf(FID, '\n   amp   = %6.4f  + %6.4f (t - %5.1f)', ...
        pa(2), pa(1), tmean);
fprintf(FID, '\n   phase = %6.2f  + %6.4f (t - %5.1f)', ...
        p(2), p(1), tmean);
fprintf(FID, '\n            median     mean      std');
fprintf(FID, '\namplitude  %7.4f  %7.4f  %7.4f',      ...
        emed_amp, emean_amp, estd_amp);
fprintf(FID, '\nphase      %7.4f  %7.4f  %7.4f',      ...
        emed_ph, emean_ph, estd_ph);
fprintf(FID, '\nnav - pc   %7.4f  %7.4f  %7.4f',      ...
        emed_dt, emean_dt, estd_dt);
fprintf(FID, '\nvar        %7.4f  %7.4f  %7.4f',  ...
        emed_var, emean_var, estd_var);
fprintf(FID, '\nmin var    %7.4f  %7.4f  %7.4f',  ...
        emed_mvar, emean_mvar, estd_mvar);
fprintf(FID, '\ndelta-u    %7.4f  %7.4f  %7.4f',  ...
        emed_du, emean_du, estd_du);
fprintf(FID, '\ndelta-v    %7.4f  %7.4f  %7.4f\n\n',  ...
        emed_dv, emean_dv, estd_dv);



fclose(FID);

