function [ amp, phi, varargout ] = calib( adcp, ship, varargin)
%--------------------------------------------------------------------------
% SYNOPSIS
%
%   [ amp, phi [,time] ] = calib(adcp, ship, file, title, year_base)
%
% DESCRIPTION
%
%   Calculate calibration factor for velocities in ADCP and SHIP. These are
%   matrices with the following columns
%
%      1. time in decimal days
%      2. eastward velocity (or similar)
%      3. northward velocity (or similar)
%
%   Data are interpolated on the smallest common time interval. The longest
%   data set is interpolated to times of the shortest. A phase and
%   amplitude factor are determined for each profile. A first estimate is
%   calculated by taking the mean value of these. A second estimate is
%   determined by rejecting all values outside the expected range
%
%      val = mean +/- sdev
%
%   determined from the first estimate. If FILE is specified, all data are
%   written to file FILE.cal. If TITLE is specified, data are plotted and
%   saved into FILE.ps. Parameter YEAR_BASE is used for axis-label
%   annotation. 
%
%--------------------------------------------------------------------------

% Defaults:

file      = '';
title     = 'UNSPECIFIED';
year_base = -1;

% Arguments:

if nargin > 2,
  file = varargin{1};
end
if nargin > 3,
  title = varargin{2};
end
if nargin > 4,
  year_base = varargin{3};
end

% Find a common time interval:
t_min = max( min(adcp(:,1)), min(ship(:,1)));
t_max = min( max(adcp(:,1)), max(ship(:,1)));
i_adcp = find((adcp(:,1) >= t_min) & (adcp(:,1) <= t_max));
i_ship = find((ship(:,1) >= t_min) & (ship(:,1) <= t_max));
adcp = adcp(i_adcp,:);
ship = ship(i_ship,:);

% Interpolate the set with most data to the other one times
if length(i_adcp) < length(i_ship),
  ship = interpolate(ship, adcp);
  time = adcp(:,1);
else
  adcp = interpolate(adcp, ship);
  time = ship(:,1);
end

% Calculate complex velocities
a_ship = abs(ship(:,2) + i*ship(:,3));
p_ship = unwrap(angle(ship(:,2) + i*ship(:,3)));
a_adcp = abs(adcp(:,2) + i*adcp(:,3));
p_adcp = unwrap(angle(adcp(:,2) + i*adcp(:,3)));

% Calculate correction:
amp   = a_ship ./ a_adcp;
phi   = p_ship - p_adcp;
imask = ~isnan(amp) .* ~isnan(phi);
igood = find(imask);
amp   = amp(igood);
phi   = phi(igood);
time  = time(igood);

% First estimate:
a_mn = mean(amp);
a_sd = std(amp);
p_mn = mean(phi);
p_sd = std(phi);

% Second estimate:
a_ii   = find(abs(amp - a_mn) < a_sd);
amp_mn = mean(amp(a_ii));
amp_sd = std(amp(a_ii));
p_ii   = find(abs(phi - p_mn) < p_sd);
phi_mn = mean(phi(p_ii));
phi_sd = std(phi(p_ii));

% Save results:
if length(file),
  save_cal([ time, amp, phi ], ...
      [ a_mn  , a_sd  , length(amp)  ; ...
    p_mn  , p_sd  , length(phi)  ; ...
    amp_mn, amp_sd, length(a_ii) ; ...
    phi_mn, phi_sd, length(p_ii) ], ...
      file, title, year_base);
end

% Plot results
plot_cal([ time, amp, phi ], ...
  [ a_mn  , a_sd  , length(amp)  ; ...
    p_mn  , p_sd  , length(phi)  ; ...
    amp_mn, amp_sd, length(a_ii) ; ...
    phi_mn, phi_sd, length(p_ii) ], ...
  file, title, year_base);

% Arguments:
if nargout > 2,
  varargout{1} = time;
end

return

%--------------------------------------------------------------------------

function new = interpolate(old, ref);

new(:,1) = ref(:,1);
new(:,2) = interp1(old(:,1), old(:,2), ref(:,1));
new(:,3) = interp1(old(:,1), old(:,3), ref(:,1));

return;

%--------------------------------------------------------------------------

function save_cal(data, stat, file, title, year_base);

fid = fopen([ file '.cal'], 'a');

% Stamp:
fprintf(fid, '\n%% ADCP CALIBRATION: %s\n', date);

% Header:
fprintf(fid, '%% File : %s', file);
fprintf(fid, '%% Title: %s', title);
fprintf(fid, '%% Year : %4d', year_base);
fprintf(fid, '%% First Estimation:\n');
fprintf(fid, '%% 1. amp = %6.4f +/- %6.4f (n=%-d)\n', stat(1,:));
fprintf(fid, '%% 1. phi = %6.4f +/- %6.4f (n=%-d)\n', stat(2,:));
fprintf(fid, '%% Second Estimation:\n');
fprintf(fid, '%% 2. amp = %6.4f +/- %6.4f (n=%-d)\n', stat(3,:));
fprintf(fid, '%% 2. phi = %6.4f +/- %6.4f (n=%-d)\n', stat(4,:));
writeasc(fid, '  ', '1.0E38', '%12.6f', data);
fclose(fid);

return;

%--------------------------------------------------------------------------

function plot_cal(data, stat, file, title, year_base);

% Time extent
t_min = min(data(:,1));
t_max = max(data(:,1));

% Figure
figure;
orient tall;

% Amplitude data
subplot(221);
hold on;
plot(data(:,1), data(:,2), '.');
plot([ t_min, t_max], [ stat(1,1), stat(1,1) ], '-');
plot([ t_min, t_max], [ stat(1,1)+stat(1,2), stat(1,1)+stat(1,2) ], '--');
plot([ t_min, t_max], [ stat(1,1)-stat(1,2), stat(1,1)-stat(1,2) ], '--');
grid on;
if year_base > 0,
  xlabel(sprintf('Time (decimal days in %4d)', year_base));
else
  xlabel('Time (decimal days)');
end
ylabel('Amplitude');
g = get(gca, 'title');
set(g, 'string', sprintf('Pass 1: %-6.4f +/- %-6.4f (n=%-d)', stat(1,:)));


% Amplitude histogram
subplot(222);
hist(data(:,2), 20);
grid on;
xlabel(sprintf('Amplitude (n=%-d)', stat(1,3)));
ylabel('Occurences');
g = get(gca, 'title');
set(g, 'string', sprintf('Pass 2: %-6.4f +/- %-6.4f (n=%-d)', stat(3,:)));

% Phase data
subplot(223);
hold on;
plot(data(:,1), data(:,3), '.');
plot([ t_min, t_max], [ stat(2,1), stat(2,1) ], '-');
plot([ t_min, t_max], [ stat(2,1)+stat(2,2), stat(2,1)+stat(2,2) ], '--');
plot([ t_min, t_max], [ stat(2,1)-stat(2,2), stat(2,1)-stat(2,2) ], '--');
grid on;
if year_base > 0,
  xlabel(sprintf('Time (decimal days in %4d)', year_base));
else
  xlabel('Time (decimal days)');
end
ylabel('Phase');
g = get(gca, 'title');
set(g, 'string', sprintf('Pass 1: %-6.4f +/- %-6.4f (n=%-d)', stat(2,:)));

% Phase histogram
subplot(224);
hist(data(:,3), 20);
grid on;
xlabel(sprintf('Phase (n=%-d)', stat(2,3)));
ylabel('Occurences');
g = get(gca, 'title');
set(g, 'string', sprintf('Pass 2: %-6.4f +/- %-6.4f (n=%-d)', stat(4,:)));

% Usual stuf
if length(title),
  plotstuf('-title', title);
end

% Save the plot:
if length(file),
  eval([ 'print -dps -append ' file '.ps' ]);
end

return;
 

