function adcp2geo(theInputFile, theOutputFile, theStart, theEnd)

% adcp2geo -- Convert ADCP beam data to earth coordinates.
%  adcp2geo('theInputFile', 'theOutputFile', theStart, theEnd)
%   processes the ADCP "vel" variables in 'theInputFile' and
%   writes the results to the "vel" variables in 'theOutputFile'.
%   If the filenames are not given, a dialog is invoked to get
%   them.  The output file should be an exact copy of the input
%   file, suitably renamed.  This function does not create the
%   output file automatically.  If theStart and/or theEnd are
%   not given, they default to the starting and ending ensembles,
%   respectively.
 
% Copyright (C) 1998 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 23-Sep-1998 13:51:40.
% Updated    24-Sep-1998 10:48:03.

if nargin < 1, help(mfilename), end

if nargin < 1
	disp(' ## InputFile...')
	f = netcdf('nowrite');
else
	f = netcdf(theInputFile, 'nowrite');
end
if isempty(f), return, end

if nargin < 2
	disp(' ## OutputFile...')
	g = netcdf('write');
else
	g = netcdf(theOutputFile, 'write');
end
if isempty(g), close(f), return, end

disp([' ### ' mfilename])
disp([' ## Input:  ' name(f)])
disp([' ## Output: ' name(g)])

bin = size(f('bin'));
ensemble = size(f('ensemble'));

if nargin < 3, theStart = 1; end
if nargin < 4, theEnd = ensemble; end

if ischar(theStart), theStart = eval(theStart); end
if ischar(theEnd), theEnd = eval(theEnd); end

theStart = max(theStart, 1);
theEnd = min(theEnd, ensemble);

for i = 1:4
	vel{i} = f{['vel' int2str(i)]};   % Input ADCP velocities.
	cur{i} = g{['vel' int2str(i)]};   % Output currents and error.
end

% Defaults.

theElevations = [-70 -70 -70 -70];
theAzimuths = [270 90 0 180];

% #187 ADCP.

if 0
	theElevations = [-69.4300 -69.5200 -69.7300 -69.7700];
	theAzimuths = [270.3800 90.3700 359.6300 179.6200];
end

theHeading = f{'Hdg'};
thePitch = f{'Ptch'};
theRoll = f{'Roll'};
theOrientation = lower(f.orientation(:));
if isempty(theOrientation)
	theOrientation = 'up';
end
theBlankingDistance = 0;

theFillValue = fillval(vel{i});

lasterr('')

p = zeros(4, bin);
q = zeros(4, bin);

disp([' ## Processing ensembles ' int2str(theStart) ':' int2str(theEnd) '.'])

for i = theStart:theEnd
	if (0), disp([' ### ' int2str(i)]); end
	try
		drawnow
		count = i - theStart;
		remaining = theEnd - i + 1;
		if count > 1
			dt = 0.5 * dt + 0.5 * toc;
			tic
		elseif count == 1
			dt = toc;
			tic
		elseif count == 0
			tic
		end
		if count == 10  | (count > 10 & rem(remaining, 100) == 0)
			t = remaining * dt;
			dd = dhms(t/86400);
			disp([' ## Doing: #' int2str(i) ...
					'; Remaining: ' int2str(remaining) ...
					' ensembles; ' dd])
		end
		for k = 1:4
			p(k, :) = vel{k}(i, :);
		end
		p(p == theFillValue) = nan;
		if (1)
			q = bm2geo(p.', theElevations, theAzimuths, ...
				theHeading(i), thePitch(i), theRoll(i), ...
				theOrientation, theBlankingDistance).';   % Note transposes.
		else   % Quicker-and-dirtier.
			q(1, :) = [1 -1 0 0] * p;
			q(2, :) = [0 0 1 -1] * p;
			q(3, :) = [1 1 1 1] * p;
			q(4, :) = [1 1 -1 -1] * p;
			h = theHeading(i) * pi / 180;
			c = cos(h);
			s = sin(h);
			r = [c -s; s c];
			q(1:2, :) = r * q(1:2, :);
		end
		q(isnan(q)) = theFillValue;
		for k = 1:4
			cur{k}(i, :) = q(k, :);
		end
	catch
		disp([' ## Aborted: ' lasterr])
		disp([' ## While processing record #' int2str(i)])
		break
	end
end

g.transform(:) = 'EART';

close(f)
close(g)

try
	beep
catch
end
