function [MSL,Dstd] = RDsurface(rawdata,ADCP_offset,progPath,DepthFile);

%function [MSL,Dstd] = RDsurface(progPath,rawdata,DepthFile,ADCP_offset);
%This function will run the RDI surface program to obtain an approximate
%water depth.  If file names are not given, they will be requested.
%
%INPUTS:
%	progPath = the full path for the RDI surface program on your computer
%	rawdata = ADCP Binary data file
%	DepthFile = the output data file name
%	ADCP_offset = the height of the ADCP above the bottom in meters!
%						if not give will be pulled from netcdf file.
%
%Outputs:
%	MSL = the mean sea level
%	Dstd = the standard of deviation of the depth which should 
%			 be an approximate tidal fluctuation

% Written by Jessica M. Cote
% for the U.S. Geological Survey
% Coastal and Marine Geology Program
% Woods Hole, MA
% http://woodshole.er.usgs.gov/
% Please report bugs to jcote@usgs.gov

%version
%updated 10-Aug-1999 16:20:42


if nargin < 1, rawdata='';, end
if nargin < 2, ADCP_offset = '';, end
if nargin < 3, progPath='';,end
if nargin < 4, DepthFile ='';, end

if isempty(rawdata), rawdata ='*';, end
if isempty(DepthFile), DepthFile = '*';, end
if isempty(ADCP_offset), ADCP_offset = NaN;, end

%where are we currently
s=pwd;

%Get binary ADCP data file
if any(rawdata == '*')
   [dataFile, dataPath] = uigetfile(rawdata, 'Select Binary ADCP File:');
   if ~any(dataFile), return, end
   if dataPath(end) ~= filesep, dataPath(end+1) = filesep; end
   rawdata = [dataPath dataFile];
else
   [dataPath, dataname, ext]=fileparts(rawdata);
   dataFile=[dataname ext];
end

%even if it is given we need to make sure

%get the ADCP trasnducer offset out of *.cdf file
if isnan(ADCP_offset)
   [theFile, thePath] = uigetfile('*.cdf', 'Select Netcdf ADCP File to find ADCP offset:');
	if ~any(theFile), return, end
	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	cdfFile = [thePath theFile];
   f=netcdf(cdfFile);
   ADCP_offset=f{'D'}.xducer_offset_from_bottom(:);
   close(f)
end

%where is the program located
if isempty(progPath)
prompt={'In put Path:'};
title='Where is RDI surface program?';
lineNo=1;
DefAns={['drive' filesep 'folder' filesep]};
dlgresult=inputdlg(prompt,title,lineNo,DefAns);
   progPath=char(dlgresult{1});   
   
end

if progPath(end) ~= filesep, progPath(end+1) = filesep; end

   
%Create output file for depths
if any(DepthFile == '*')
   [thePath,theFile,ext] = fileparts(rawdata);
  DepthFile = fullfile(Path, [theFile '.dat']);
end

[Dpath, Dname, dext]=fileparts(DepthFile);
Dfile=([Dname dext]);

%SURFACE InFile OutFile [Ensembles Skip Bin1 LastBin]
% Copy the ADCP file to a the RDI directory.
cpfile = fullfile(progPath,dataFile)

if isunix
		eval(['!cp ' rawdata ' ' cpfile])
	elseif any(findstr(lower(computer), 'pcwin')) | isVMS
		eval(['!copy ' rawdata ' ' cpfile])
	elseif any(findstr(lower(computer), 'mac')) & ...
			exist('aduplicate') == 2
		feval('aduplicate', rawdata, cpfile)
	else
		fcopy(rawdata, cpfile)
      
   end
   
%this all has to happen in the same directory as the surface program
	try
		eval(['cd ' progPath])
	catch
   	prompt={'Correct Path:'};
		title='Path incorrect or missing end filesep?';
		lineNo=1;
		DefAns={progPath};
		dlgresult=inputdlg(prompt,title,lineNo,DefAns);
   	progPath=char(dlgresult{1});   
      if progPath(end) ~= filesep, progPath(end+1) = filesep; end
      eval(['cd ' progPath]);
   end
   
disp('Running RDI surface program')   
eval(['!surface ' dataFile ' ' Dfile])
delete(cpfile)
   
datfile = mfriend(Dfile);
surfdat=load(datfile);

%move some stuff around
%if ~isequal(Dpath,progPath)
	if isunix
		eval(['!mv ' Dfile ' ' DepthFile])
	elseif any(findstr(lower(computer), 'pcwin')) | isVMS
		eval(['!move ' Dfile ' ' DepthFile])
	elseif any(findstr(lower(computer), 'mac')) & ...
			exist('aduplicate') == 2
      feval('aduplicate', Dfile, DepthFile)
      delete(Dfile)
	else
		fcopy(Dfile,DepthFile)
      delete(Dfile)
   end
%end

%go back to where we were before
cd(s)

%these are 2 ways of finding the good data,
%the first is by the data quality given by rdi
[m,n]=size(surfdat);
qual=find(surfdat(:,7)==4);

%the second is just taking everything greater than 0,
%which looks as though it should be the same.
Dm=surfdat(:,6);
%idgood=find(Dm>0);
%Dm=Dm(idgood);
Dm=Dm(qual);

Davg=mean(Dm);
Dmin=min(Dm);
Dmax=max(Dm);
Dstd=std(Dm);

%are we really getting only the good data?
p=figure;, plot(Dm,'r.'), hold

%Depth_constraints.Are_these_depths_reasonable={'radiobutton',0}
%Depth_constraints.Minimum_depth_in_meters={num2str(Dmin)};
%Depth_constraints.Maximum_depth={num2str(Dmax)};
%Depth_constraints=uigetinfo(Depth_constraints);

prompt={'Minimum depth in meters:',...
   'Maximum depth in meters:','Are these depths reasonable?'};
title='Check surface program output';
lineNo=1;
DefAns={num2str(Dmin),num2str(Dmax),'N'};
dlgresult=inputdlg(prompt,title,lineNo,DefAns);
   
   if ~isequal(Dmin,str2num(dlgresult{1})) | ~isequal(Dmax,str2num(dlgresult{2})) 
   Dmin=str2num(dlgresult{1});
   Dmax=str2num(dlgresult{2});
	idgood=find(Dm > Dmin & Dm < Dmax);
	Davg=mean(Dm(idgood));
	Dstd=std(Dm(idgood));
	clf, plot(Dm(idgood),'r.');,hold;
	Dm=Dm(idgood);
   
   elseif upper(char(dlgresult{3}))=='N';       
      prompt={'Minimum depth allowed:',...
   'Maximum depth allowed:'};
   title='Input valid range for depth';
   lineNo=1;
   DefAns={num2str(Dmin),num2str(Dmax)};
   dlgresult=inputdlg(prompt,title,lineNo,DefAns);
   
   Dmin=str2num(dlgresult{1});
   Dmax=str2num(dlgresult{2});
	idgood=find(Dm > Dmin & Dm < Dmax);
	Davg=mean(Dm(idgood));
	Dstd=std(Dm(idgood));
	clf, plot(Dm(idgood),'r.');
	Dm=Dm(idgood);
   
   
end

MSL=Davg+ADCP_offset;
disp(['ADCP measured ' num2str(MSL) ' m from surface to the seabed (mean sea level)']);
disp(['The tidal range is approximately ' num2str(Dstd) ' m']);
disp([' '])

pause(3)
close(p);
