% Script file adcpview.m
%
%       ADCP Data Quick View 
%	for Self Contained ADCP applications
%
%	Read data converted to netCDF format by rdi2cdf.m
%	This is data from Workhorse and Broadband ADCP's
%
%	Written for use on a PC-notebook with only 8 MB RAM,
%	this m-file will restrict the total amount of data
%	that can be read at once with the variable MAXDATA.
%
%	This file assumes that the following were recorded:
%	velocity, correlation, intensity and % good.  
%	If your data was recorded in earth coordinates:
%		Beam 1 = to the East
%		Beam 2 = to the North
%		Beam 3 = away from the instrument
%		Beam 4 = error velocity
%
%	To invole averaging, set a variable called AVERAGE in your
%	workspace to the number of ensembles to average.
%
%	TO EXECUTE THIS FILE, TYPE ADCPMENU TO CALL THE GUI MENU
%

%	Supporting m-files and their dependents:
%		rdi2cdf.m	RDI to netCDF format converter
%		adcpmenu.m	menu driver which calls this script
%		mexcdf.mex	netCDF access routines
%		mcdump.m	netCDF file header dumper
%		julian.m	convert from gregorian to julian time
%		gregoria.m  convert from julian to gregorian time
%		timeplt.m	plot data on Gregorian axis
%		bm2earth.m	transform ADCP data recorded in beam to earth coordinates,
%					accounting for ADCP pitch, roll, heading, sound speed, etc.

% Written by Marinna Martini
% for the U.S. Geological Survey
% Marine and Coastal Program
% Woods Hole Center, Woods Hole, MA
% http://woodshole.er.usgs.gov/
% Please report bugs to mmartini@usgs.gov
%
% gregoria, julian and timeplt routines written by Rich Signell
% netCDF attributes and variable names follow EPIC conventions
%
% information on netCDF may be obtained at
% http://www.unidata.ucar.edu
% it's free!!%
%
% 10/29/96	Version 0.0
%

echo off
MAXDATA=800000;
XFREQS=[76.8 153.6 307.2 614.4 1228.8 2400.0];
ANGLES = [15 20 30 0];
BEAMS=['1' '2' '3' '4'];
NBEAMS=4;
XDUCER_HEIGHT=2;	% MAB of ADCP transducers
% offset indicators
VELOCITY = 3;
CORRELATION = 4;
INTENSITY = 5;
GOOD = 6;
% data names
DNAMES=['vel';'cor';'AGC';'PGd'];
DIRECTS = ['East '; 'North'; 'Vert.'; 'Error'];
% ops flags
if exist('AVERAGE') ~= 1, AVERAGE = 1; end
if exist('FIRSTBST') ~=1, FIRSTBST = 1; end

% get the handle of the menu window to be able to hide it later
% it is the current figure because it called this script
menuptr=get(0,'currentfigure');

if ((cop == 1) | (exist('file') ~= 1))
	[file, pth]=uigetfile('*.cdf','Get ADCP netCDF Data file');
	path_file = [pth file];
	if file==0, break; end % user aborted
	disp(['New file: ' path_file]);
	% and get other pertinent data
	cdfid = mexcdf('OPEN',path_file,'NOWRITE');
	if cdfid < 0, error(['Failed to open ' path_file]); end
	% get some basic data
	xfreq = mexcdf('ATTGET',cdfid,'NC_GLOBAL','frequency'); 
	sensor_config = mexcdf('ATTGET',cdfid,'NC_GLOBAL','sensor_configuration');
	junk = mexcdf('ATTGET',cdfid,'NC_GLOBAL','orientation');
	if strcmp(junk,'DOWN'), beams_up=0; 
	else beams_up=1; end		
	beam_angle = mexcdf('ATTGET',cdfid,'NC_GLOBAL','beam_angle'); 
	[junk, nbins, rcode] = mexcdf('DIMINQ',cdfid,'bin');
	binsize = mexcdf('ATTGET',cdfid,'D','bin_size');	% meters
	blank = mexcdf('ATTGET',cdfid,'D','blanking_distance');		% meters
	junk = mexcdf('ATTGET',cdfid,'NC_GLOBAL','transform');	% coordinate transformation
	if strcmp(junk,'BEAM'), transform = 0;
	elseif strcmp(junk,'INST'), transform = 1;
	elseif strcmp(junk,'SHIP'), transform = 2;
	else transform = 3; end	
	heading_bias=mexcdf('ATTGET',cdfid,'NC_GLOBAL','heading_bias')./100;	% degrees
	junk = mexcdf('ATTGET',cdfid,'NC_GLOBAL','beam_pattern');	% coordinate transformation
	if strcmp(junk,'convex'), convex = 1;
	else convex = 0;
	end	
	[junk, nens, rcode] = mexcdf('DIMINQ',cdfid,'ensemble');
	XDUCER_HEIGHT =	mexcdf('ATTGET',cdfid,'D','xducer_offset_from_bottom');
	center_first_bin = mexcdf('ATTGET',cdfid,'D','center_first_bin'); 
	if XDUCER_HEIGHT == 0,
		disp('Transducer height above bottom = 0.')
		XDUCER_HEIGHT = input('Enter new value: ');
	end
	D=(0:nbins-1).*binsize+center_first_bin+XDUCER_HEIGHT;
	% ------ done getting pertinent data ------
	% do we need to limit loading?
	corner = [0 0];
	edges = [nens nbins];
	if nens*nbins/AVERAGE > MAXDATA,
		n=menu('Too much data!  Load partial:','Bins','Ensembles');
		set(menuptr,'visible','off');
		if n==1, 
			disp(sprintf('There are %d bins,\n',nbins))
			disp(sprintf('Only %d can be loaded into memory at one time\n',...
				floor(MAXDATA/nens)))
			while 1,
				corner(2)=input('Starting bin: ')-1;
				edges(2)=input('Number of bins: ');
				if edges(2)*nens > MAXDATA, 
					disp('Too many bins...')
				else 
					break;
				end
			end
		elseif n==2,
			disp(sprintf('There are %d ensembles,\n',nens))
			disp(sprintf('Only %d can be loaded into memory at one time\n',...
				floor(MAXDATA/edges(2))))
			while 1,
				corner(1)=input('Starting ensemble: ')-1;
				edges(1)=input('Number of ensembles: ');
				if edges(2)*edges(1) > MAXDATA,
					disp('Too many ensembles...')
				else
					break;
				end
			end	
		end
		set(menuptr,'visible','on');
	end
	nens = edges(1);
	if AVERAGE > 1, 
		nens=floor((nens-FIRSTBST)/AVERAGE); 
		corner(1) = corner(1)+FIRSTBST;
	end
	nbins = edges(2);
	mexcdf('CLOSE',cdfid);
	disp('Done assessing new file');
	disp(setstr(7));
	clear junk i j k buf
	if AVERAGE > 1,
		disp(sprintf('Averaging is set to %i',AVERAGE))
		disp(sprintf('starting at record %i',FIRSTBST))
	end
end
%
% --- view Fixed Leader and Header Data
%
if (cop == 2)

	% use mexddump to display header info
	% thanks to John Evans
	mcdump(path_file,'nc_global');
	%mcdump(path_file,'D');

end
%
% ------ Plot Variable Leader data
%
if cop == 3,
	n=menu('Select data to plot:','All','Temperature','Sound Speed','Orientation',...
		'Xmit Voltage & Current','DAC','VDD3 & VDD1','Main Battery Voltage');
	t=get(0,'Children');
	if length(t)~=1,
		% close the other windows to save memory
		for i=2:max(t),
			close(i)
		end
	end
	% get time information
	if exist('TIM') ~= 1, 
		if AVERAGE > 1,
			TIM = mcgetavg(path_file,'TIM',corner(1),edges(1),AVERAGE)';	
		else
			TIM = mcgetvar(path_file,'TIM',corner(1),edges(1))';	
		end
		% do some massage
		missing_ens=find(TIM < -99999);
		TIM(missing_ens)=[];
		nens=length(TIM);
	end
	if n==1 | n==2,
		% get temperature data
		if exist('Tx') ~= 1,
			if AVERAGE > 1,
				[Tx, Tx_name, Tx_units, Tx_sts] = mcgetavg(path_file,'Tx',corner(1),edges(1), AVERAGE);
				Tx=Tx'; Tx_sts=Tx_sts';
			else
				[Tx, Tx_name, Tx_units] = mcgetvar(path_file,'Tx',corner(1),edges(1));
			end
			% do some massage
			Tx(missing_ens)=[];
		end
		% plot 
		figure('name',Tx_name,'numbertitle','off');
	    timeplt(TIM,Tx);
    	ylabel([Tx_name ': ' Tx_units])
		title(file)
	end
	if n==1 | n==3,
		% get the speed of sound data
		if exist('sv') ~= 1,
			if AVERAGE > 1,
				[sv, sv_name, sv_units, sv_sts] = mcgetavg(path_file,'sv',corner(1),edges(1), AVERAGE);
				sv=sv';  sv_sts=sv_sts';
			else
				[sv, sv_name, sv_units] = mcgetvar(path_file,'sv',corner(1),edges(1));
			end
			% do some massage
			sv(missing_ens)=[];
		end
		% plot 
		figure('name',sv_name,'numbertitle','off');
	    timeplt(TIM,sv);
    	ylabel([sv_name ': ' sv_units])
		title(file)
	end
	if n==1 | n==4,
		% get orientation data
		if exist('Hdg') ~= 1 | exist('Ptch') ~= 1 | exist('Roll') ~= 1,
			if AVERAGE > 1,
				[Hdg, Hdg_name, Hdg_units, Hdg_sts] = mcgetavg(path_file,'Hdg',corner(1),edges(1),AVERAGE);
				[Ptch, Ptch_name, Ptch_units, Ptch_sts] = mcgetavg(path_file,'Ptch',corner(1),edges(1),AVERAGE);
				[Roll, Roll_name, Roll_units, Roll_sts] = mcgetavg(path_file,'Roll',corner(1),edges(1),AVERAGE);
				Hdg=Hdg'; Hdg_sts=Hdg_sts';
				Ptch=Ptch';	Ptch_sts=Ptch_sts';
				Roll=Roll'; Roll_sts=Roll_sts';
			else
				[Hdg, Hdg_name, Hdg_units] = mcgetvar(path_file,'Hdg',corner(1),edges(1));
				[Ptch, Ptch_name, Ptch_units] = mcgetvar(path_file,'Ptch',corner(1),edges(1));
				[Roll, Roll_name, Roll_units] = mcgetvar(path_file,'Roll',corner(1),edges(1));
			end
			% do some massage
			Hdg(missing_ens)=[];
			Ptch(missing_ens)=[];
			Roll(missing_ens)=[];
		end
		figure('name','Orientation','numbertitle','off');
	    timeplt(TIM,[Hdg Ptch Roll],[1 2 3],...
	    	str2mat(Hdg_name, Ptch_name, Roll_name));
			% ['Heading'; 'Pitch  '; 'Roll   ']);
		text('units','normalized','position',[1.05 0.01],...
			'string','all units degrees','rotation',-90)		
		title(file)
	end
	if n==1 | n==5,
		% get energy output data
		if exist('xmitc')~=1 | exist('xmitv')~=1,
			if AVERAGE > 1,
				[xmitc, xmitc_name, xmitc_units, xmitc_sts] = mcgetavg(path_file,'xmitc',corner(1),edges(1),AVERAGE);
				[xmitv, xmitv_name, xmitv_units, xmitv_sts] = mcgetavg(path_file,'xmitv',corner(1),edges(1),AVERAGE);
				xmitc=xmitc'; xmitc_sts=xmitc_sts';
				xmitv=xmitv'; xmitv_sts=xmitv_sts';
			else
				[xmitc, xmitc_name, xmitc_units] = mcgetvar(path_file,'xmitc',corner(1),edges(1));
				[xmitv, xmitv_name, xmitv_units] = mcgetvar(path_file,'xmitv',corner(1),edges(1));
			end
			% do some massage
			xmitc(missing_ens)=[];
			xmitv(missing_ens)=[];
		end
		% plot
		figure('name','ADCP Xmit Data','numbertitle','off');
	    timeplt(TIM,[xmitc xmitv],[1 2],str2mat(xmitc_name, xmitv_name));
		text('units','normalized','position',[1.05 0.7],...
			'string',xmitv_units,'rotation',-90)		
		text('units','normalized','position',[1.05 -0.5],...
			'string',xmitc_units,'rotation',-90)		
		title(file)
	end
	if n==1 | n==6,
		% get data
		if exist('dac')~=1,
			if AVERAGE > 1,
				[dac, dac_name, dac_units, dac_sts] = mcgetavg(path_file,'dac',corner(1),edges(1),AVERAGE);
			else
				[dac, dac_name, dac_units] = mcgetvar(path_file,'dac',corner(1),edges(1));
				dac=dac'; dac_sts=dac_sts';			
			end
			% do some massage
			dac(missing_ens)=[];
		end
		% plot 
		figure('name',dac_name,'numbertitle','off');
	    timeplt(TIM,dac);
    	ylabel([dac_name ': ' dac_units])
		title(file)
	end
	if n==1 | n==7,
		% get voltage data
		if exist('VDD1')~=1 | exist('VDD3')~=1,
			if AVERAGE > 1,
				[VDD1, VDD1_name, VDD1_units, VDD1_sts] = mcgetavg(path_file,'VDD1',corner(1),edges(1),AVERAGE);
				[VDD3, VDD3_name, VDD3_units, VDD3_sts] = mcgetavg(path_file,'VDD3',corner(1),edges(1),AVERAGE);
				VDD1=VDD1'; VDD1_sts=VDD1_sts';
				VDD3=VDD3'; VDD3_sts=VDD3_sts';
			else
				[VDD1, VDD1_name, VDD1_units] = mcgetvar(path_file,'VDD1',corner(1),edges(1));
				[VDD3, VDD3_name, VDD3_units] = mcgetvar(path_file,'VDD3',corner(1),edges(1));
			end
			% do some massage
			VDD1(missing_ens)=[];
			VDD3(missing_ens)=[];
		end
		figure('name','ADCP Voltage Data','numbertitle','off');
	    timeplt(TIM,[VDD1 VDD3],[1 2],str2mat(VDD1_name, VDD3_name));
		text('units','normalized','position',[1.05 0.7],...
			'string',VDD1_units,'rotation',-90)		
		text('units','normalized','position',[1.05 -0.5],...
			'string',VDD3_units,'rotation',-90)		
		title(file)
	end
	if n==1 | n==8,
		% get data
		if AVERAGE > 1,
			[VDC, VDC_name, VDC_units, VDC_sts] = mcgetavg(path_file,'VDC',corner(1),edges(1),AVERAGE);
			VDC=VDC'; VDC_sts=VDC_sts';
		else
			[VDC, VDC_name, VDC_units] = mcgetvar(path_file,'VDC',corner(1),edges(1));
		end
		% do some massage
		VDC(missing_ens)=[];
		% plot 
		figure('name','Main battery voltage: ','numbertitle','off');
	    timeplt(TIM,VDC);
    	ylabel(['Main battery voltage: ' VDC_units])
		title(file)
	end
end

%
% ------ Plot Beam related data
%
if cop == 4,
	buf=['   Velocity';'Correlation';'  Intensity';'     % Good'];
	k=menu('Select data to plot:','Velocity','Correlation','Intensity','% Good');
	% load the data
	% get time information
	if exist('TIM') ~= 1, 
		if AVERAGE > 1,
			TIM = mcgetavg(path_file,'TIM',corner(1),edges(1),AVERAGE)';	
		else
			TIM = mcgetvar(path_file,'TIM',corner,edges)';	
		end
		% do some massage
		missing_ens=find(TIM < -99999);
		TIM(missing_ens)=[];
		nens=length(TIM);
	end
	% build a single holder for all the beam data
	if exist(DNAMES(k,:)) ~= 1,
		eval([ DNAMES(k,:), '=zeros(NBEAMS*nbins,nens);' ]);
		if AVERAGE > 1, eval([ DNAMES(k,:), '_sts=zeros(NBEAMS*nbins,nens);' ]); end
		j=1;
		for i=1:nbins:NBEAMS*nbins,
			if AVERAGE > 1,
				eval(['[',DNAMES(k,:),'(i:i+nbins-1,:),',DNAMES(k,:),'_name,',...
					DNAMES(k,:),'_units, ',DNAMES(k,:),'_sts]=mcgetavg(path_file,''',...
					[DNAMES(k,:) int2str(j)],''',corner,edges,AVERAGE);']);
			else
				eval(['[',DNAMES(k,:),'(i:i+nbins-1,:),',DNAMES(k,:),'_name,',...
					DNAMES(k,:),'_units]=mcgetvar(path_file,''',...
					[DNAMES(k,:) int2str(j)],''',corner,edges);']);
			end
			j=j+1;
		end
		% remove bad data for matlab convenience
		eval([DNAMES(k,:),'(:,missing_ens)=[];']);
		eval(['flags=find(',DNAMES(k,:),'==-999999);']);
		eval([DNAMES(k,:),'(flags) = ones(size(flags)).*NaN;']);
	end
	clear flags
	% data returned as bin*beam by ensemble
	if exist(DNAMES(k,:))==1, 
		figure('name',buf(k,:),'numbertitle','off');
	    eval(['lptr_',DNAMES(k,1),'=timeplt(TIM,',DNAMES(k,:),'(1,:));']);
    	if k==1, 
    		if transform == 3,
	    		ylabel('mm/s, Earth Coordinates'),
	    	else
	    		ylabel('mm/s, Beam Coordinates'),
    		end		
    	end
		title([file ':  ' buf(k,:)]);
		eval(['tptr_',DNAMES(k,1),'=text(0.05,0.05,['' Beam 1, Bin '',int2str(corner(2)+1)],''units'',''normalized'');']);
		dtype=k;
		% put up controls
		% text labels
		uicontrol(gcf,'style','text','backgroundcolor','green',...
			'units','normalized','position',[0.6 0.95 0.04 0.04],...
			'string',sprintf('%d',corner(2)+1));
		uicontrol(gcf,'style','text','backgroundcolor','green',...
			'units','normalized','position',[0.96 0.95 0.04 0.04],...
			'string',sprintf('%d',corner(2)+edges(2)));
		eval(['uptrbin_',DNAMES(k,1),'=uicontrol(gcf,''style'',''text'',''backgroundcolor'',''green'',',...
			'''units'',''normalized'',''position'',[0.88 0.95 0.04 0.04],',...
			'''string'',sprintf(''%d'',corner(2)+1));'])
		uicontrol(gcf,'style','text','backgroundcolor','green',...
			'units','normalized','position',[0.1 0.95 0.04 0.04],...
			'string','1');
		uicontrol(gcf,'style','text','backgroundcolor','green',...
			'units','normalized','position',[0.46 0.95 0.04 0.04],...
			'string',sprintf('%d',NBEAMS));
		eval(['uptrbeam_',DNAMES(k,1),'=uicontrol(gcf,''style'',''text'',''backgroundcolor'',''green'',',...
			'''units'',''normalized'',''position'',[0.41 0.95 0.04 0.04],''string'',''1'');'])
		% sliders to select new values
		eval([ ' uptrbin1_' DNAMES(k,1) '=uicontrol(''Style'',''slider'',''Position'',[0.6 0.91 0.4 0.04], ' ...
			' ''Units'',''normalized'',''Max'',corner(2)+edges(2),''Min'',corner(2)+1,''Value'',corner(2)+1, ' ...
			' ''Userdata'',[uptrbin_' DNAMES(k,1) '],''CallBack'', ' ...
			' [''newbin=ceil(get(gco,''''value''''));'', ' ...
			' ''t=get(gco,''''userdata'''');'', ' ...
			' ''set(t,''''string'''',sprintf(''''%d'''',newbin));'', ' ...
			' ]); ' ...
			]);
		eval([ ' uptrbeam1_' DNAMES(k,1) '=uicontrol(''Style'',''slider'',''Position'',[0.1 0.91 0.4 0.04], ' ...
			' ''Units'',''normalized'',''Max'',NBEAMS,''Min'',1,''Value'',1, ' ...
			' ''Userdata'',[uptrbeam_' DNAMES(k,1) '],''CallBack'', ' ...
			' [''newbeam=ceil(get(gco,''''value''''));'', ' ...
			' ''t=get(gco,''''userdata'''');'', ' ...
			' ''set(t,''''string'''',sprintf(''''%d'''',newbeam));'', ' ...
			' ]); ' ...
			]);
		% buttons to activate selection
		eval([ ' uicontrol(gcf,''style'',''pushbutton'',''units'',''normalized'', ' ...
			' ''position'',[0.66 0.95 0.22 0.04],''string'',''Show New Bin>>'', ' ...
			' ''userdata'',[uptrbin1_' DNAMES(k,1) ', lptr_' DNAMES(k,1) ', corner(2)+edges(2), ' ...
			' uptrbeam1_' DNAMES(k,1) ', tptr_' DNAMES(k,1) '],''Callback'', ' ...
			' [''t=get(gco,''''userdata'''');'', ' ...	% get the pointer to the slider
			' ''newbin=ceil(get(t(1),''''value''''));'', ' ... % get the value of the bin slider
			' ''cbeam=ceil(get(t(4),''''value''''));'', ' ... % get the	value of the beam slider
			' ''set(t(2),''''ydata'''',' DNAMES(k,:) '((cbeam-1)*nbins+newbin,:));'', ' ...
			' ''set(t(5),''''string'''',['''' Beam '''',int2str(cbeam),'''' Bin '''',int2str(newbin)]);'', ' ...
			' ]); '...
			]);
		eval([ ' uicontrol(gcf,''style'',''pushbutton'',''units'',''normalized'', ' ...
			' ''position'',[0.16 0.95 0.25 0.04],''string'',''Show New Beam>>'', ' ...
			' ''userdata'',[uptrbin1_' DNAMES(k,1) ', lptr_' DNAMES(k,1) ', corner(2)+edges(2), ' ...
			' uptrbeam1_' DNAMES(k,1) ', tptr_' DNAMES(k,1) '],''Callback'', ' ...
			' [''t=get(gco,''''userdata'''');'', ' ...
			' ''cbin=ceil(get(t(1),''''value''''));'', ' ...
			' ''newbeam=ceil(get(t(4),''''value''''));'', ' ...
			' ''set(t(2),''''ydata'''',' DNAMES(k,:) '((newbeam-1)*nbins+cbin,:));'', ' ...
			' ''set(t(5),''''string'''',['''' Beam '''',int2str(newbeam),'''' Bin '''',int2str(cbin)]);'', ' ...
			' ]); '...
			]);
		eval([ 'clear lptr_' DNAMES(k,1) ' tptr_' DNAMES(k,1) ' dtype uptrbin1_' DNAMES(k,1) ...
			' uptrbin_' DNAMES(k,1) ' uptrbeam_' DNAMES(k,1) ' uptrbeam1_' DNAMES(k,1) ]);
	end
	clear i j k
end
%
% ------- plot on a 3D mesh surface
%
if cop==5,
    buf=['   Velocity';'Correlation';'  Intensity';'     % Good'];
    k=menu('Select data to plot:','Velocity','Correlation','Intensity','% Good');
    if (k==1) & (transform == 3),
	    b=menu('Select:',DIRECTS(1,:),DIRECTS(2,:),DIRECTS(3,:),DIRECTS(4,:));
	else
    	b=menu('Select the beam to plot:','1','2','3','4','avg of all');
    end
    if exist(DNAMES(k,:))==1, 
		figure('name',buf(k,:),'numbertitle','off');
		if b == 5,
			eval([DNAMES(k,:),'_beam_avg = zeros(nbins,nens);']);
			for i=1:NBEAMS,
				eval([DNAMES(k,:),'_beam_avg = ',DNAMES(k,:),...
					'_beam_avg+',DNAMES(k,:),'(((i-1)*nbins)+1:i*nbins,:)./NBEAMS;']);
			end
	   		eval(['pcolor((TIM-TIM(1)).*24,'...
				,'D,'...
				,DNAMES(k,:),'_beam_avg);']);
	
		else
	   		eval(['pcolor((TIM-TIM(1)).*24,'...
				,'D,'...
				,'reshape(',DNAMES(k,:),'(((b-1)*nbins)+1:b*nbins,:),nbins,nens));']);
		end
		tbuf=gregoria(TIM(1));
        xlabel(sprintf('Hours since %i/%i/%i %i:%i:%i',tbuf(2),tbuf(3),tbuf(1),tbuf(4),tbuf(5),tbuf(6)));
        ylabel('Meters Above Bottom');
    	cptr = colorbar('vert');
    	tptr = get(cptr,'xlabel');
    	if (k==1) & (transform == 3),
	        title([file,', ',buf(k,:),':  ',DIRECTS(b,:)])
	    else
	    	if (b == 5),
		        title([file,', ',buf(k,:),':  Avg of all beams'])
		    else
		        title([file,', ',buf(k,:),':  Beam ',int2str(b)])
		    end
	    end
    	if k==1, 
    		zlabel('mm/s'), 
	    	set(tptr,'string','mm/s');
		elseif k==2 | k==4, 
			zlabel('%');
	    	set(tptr,'string','%');
        end
	end
    clear cptr tptr b k buf t tbuf
end

%
% ------- close plot windows
%
if cop==6,
	t=get(0,'children');
	for i=2:max(t),
		close(i)
	end
end
%
% -------- statistics
%
if cop==7,
	% do statistics by bin and beam
	buf=['   Velocity';'Correlation';'  Intensity';'     % Good'];
	for k=2:4,
		if exist(DNAMES(k,:))==1,
			eval([DNAMES(k,:),'_mean=zeros(nbins,NBEAMS);']);
			eval([DNAMES(k,:),'_std=zeros(nbins,NBEAMS);']);
			for i=1:NBEAMS,
				for j=1:nbins,
					eval(['junk=',DNAMES(k,:),'((i-1)*nbins+j,:);']);
					flags = find(isnan(junk)==0);
			    	eval([DNAMES(k,:),'_mean(j,i) = mean(junk(flags));']);
					eval([DNAMES(k,:),'_std(j,i) = std(junk(flags));']);
				end
			end
			disp('')
			disp([buf(k,:),' Statistics']);
			disp('Bin   Means 1....4          STDs  1....4')
			frmt='%2.0f    %3.0f %3.0f %3.0f %3.0f      %3.0f %3.0f %3.0f %3.0f';
			for i=corner(2)+1:corner(2)+edges(2),
				eval(['disp(sprintf(frmt, i, ',DNAMES(k,:),'_mean(i,:), ',DNAMES(k,:),'_std(i,:)));']);
			end
			clear frmt
		else
			disp(['No ',buf(k,:),' data loaded'])
		end
	end
	clear frmt buf i j k junk
end

%
% ------ convert from beam to earth
%
if cop==8,
	if transform == 3,
		disp('Data already in Earth Coordinates');
	elseif transform == 0,
		% get time information
		if exist('TIM') ~= 1, 
			disp('Loading time data into MATLAB workspace');
			if AVERAGE > 1,
				TIM = mcgetavg(path_file,'TIM',corner(1),edges(1),AVERAGE)';	
			else
				TIM = mcgetvar(path_file,'TIM',corner(1),edges(1))';	
			end
			% do some massage
			missing_ens=find(TIM < -99999);
			TIM(missing_ens)=[];
			nens=length(TIM);
		end
		if exist('vel')~=1,
			disp('Loading velocity data into MATLAB workspace')
			% build a single holder for all the beam data
			k=1;	% select velocity
			%eval([ DNAMES(k,:), '=zeros(NBEAMS*nbins,edges(1));' ]);
			eval([ DNAMES(k,:), '=zeros(NBEAMS*nbins,nens);' ]);
			if AVERAGE > 1, eval([ DNAMES(k,:), '_sts=zeros(NBEAMS*nbins,nens);' ]); end
			j=1;
			for i=1:nbins:NBEAMS*nbins,
				if AVERAGE > 1,
					eval(['[',DNAMES(k,:),'(i:i+nbins-1,:),',DNAMES(k,:),'_name,',...
						DNAMES(k,:),'_units, ',DNAMES(k,:),'_sts]=mcgetavg(path_file,''',...
						[DNAMES(k,:) int2str(j)],''',corner,edges,AVERAGE);']);
				else
					eval(['[',DNAMES(k,:),'(i:i+nbins-1,:),',DNAMES(k,:),'_name,',...
						DNAMES(k,:),'_units]=mcgetvar(path_file,''',...
						[DNAMES(k,:) int2str(j)],''',corner,edges);']);
				end
				j=j+1;
			end
			% remove bad data for matlab convenience
			eval([DNAMES(k,:),'(:,missing_ens)=[];']);
			eval(['flags=find(',DNAMES(k,:),'==-999999);']);
			eval([DNAMES(k,:),'(flags) = ones(size(flags)).*NaN;']);
			clear flags j k
		end	
		if exist('Hdg')~=1 | exist('Ptch')~=1 | exist('Roll')~=1 | exist('sv')~=1,
			disp('Loading variable leader data (compass, pitch, roll) into MATLAB workspace')
			% get orientation data
			if AVERAGE > 1,
				[Hdg, Hdg_name, Hdg_units, Hdg_sts] = mcgetavg(path_file,'Hdg',corner(1),edges(1), AVERAGE);
				[Ptch, Ptch_name, Ptch_units, Hdg_sts] = mcgetavg(path_file,'Ptch',corner(1),edges(1), AVERAGE);
				[Roll, Roll_name, Roll_units, Hdg_sts] = mcgetavg(path_file,'Roll',corner(1),edges(1), AVERAGE);
			else
				[Hdg, Hdg_name, Hdg_units] = mcgetvar(path_file,'Hdg',corner(1),edges(1));
				[Ptch, Ptch_name, Ptch_units] = mcgetvar(path_file,'Ptch',corner(1),edges(1));
				[Roll, Roll_name, Roll_units] = mcgetvar(path_file,'Roll',corner(1),edges(1));
			end
			% do some massage
			Hdg(missing_ens)=[];
			Ptch(missing_ens)=[];
			Roll(missing_ens)=[];
			% get the speed of sound data
			if AVERAGE > 1,
				[sv, sv_name, sv_units, sv_sts] = mcgetavg(path_file,'sv',corner(1),edges(1), AVERAGE);
			else
				[sv, sv_name, sv_units] = mcgetvar(path_file,'sv',corner(1),edges(1));
			end
			% do some massage
			sv(missing_ens)=[];
		end
		disp('Converting velocity to earth coordinates.')
		% Use data(nbins*n+1:(nbins*n)+nbins,:) to recover a single ensemble
		% of data (nbins, nbeams) where n=nens-1.  
		earth=zeros(size(vel));
		for n=1:nens,
			junk=bm2earth(reshape(vel(:,n),nbins,NBEAMS),...
				Hdg(n), heading_bias, Ptch(n), Roll(n),...
				sv(n), 1500, beams_up, xfreq, convex,...
				sensor_config, beam_angle);
			earth(:,n)=junk(:);
		end
		vsave=vel;
		vel=earth;
		disp('The vel variable now contains data in earth cordinates')
		disp('Beam 1 = east')
		disp('Beam 2 = north')
		disp('Beam 3 = vertical')
		disp('Beam 4 = error')
		disp(sprintf('Heading Bias used was %4.2f',heading_bias))
		transform = 3;	% set to earth
		clear earth
	else
		disp('Data is in ship or instrument coordinates.')
	end
end
%
% ------ recover beam data
%
if cop == 9,
	if exist('vsave')~=1,
		disp('Old beam data not available');
	else
		vel=vsave;
		disp('The vel variable now contains data in beam cordinates')
		transform = 0;
	end
end
%
% ----- save workspace to *****.mat
%
if cop==10,
	% construct a name
	matfile=lower(path_file);
	n=max(findstr(matfile,'.'));
	matfile(n+1:n+3)='mat';
	k=1;
	if exist(matfile)==2,
		k=menu(['Overwrite existing ',matfile],'yes','no');
	end
	if k==1,
		eval(['save ' matfile])
		disp(['Workspace saved as ',matfile])
	end
	clear k n
end
%
% ------ save the workspace
%
if cop==11,
	disp(['Saving workspace to ',pwd,'\matlab.mat'])
	save
 end

%
%	------- save the latest plot to disk
%
if cop==12,
	t=get(0,'children');
	if isempty(t) | length(t)<2,
		disp('No existing plots');
	else
		% get a name
		set(menuptr,'visible','off');
		buf = input('Enter a file name for this plot: ','s');
		eval(['print -f',int2str(max(t)),' -dps ',pth,buf]);
		set(menuptr,'visible','on');
	end		
end
%
% ------ toggle a diary
%
if cop==13,
	if exist('dfile')~=1,
		disp('Creating a new diary file');
		% construct a name
		dfile=lower(path_file);
		k=max(findstr(dfile,'.'));
		dfile(k-3:k-1)='log';
		dfile(k+1:k+3)='txt';
		eval(['diary ' dfile])
	else
		set(0,'diaryfile',dfile);
		diary
	end
	disp(['diary file = ',get(0,'diaryfile')])
	disp(['diary is now ',get(0,'diary')])
	t=clock;
	disp(sprintf('%s %i:%i:%i',date,t(4),t(5),floor(t(6))));
end

%
% -------- purge large matrices
%
if cop==14,
	clear sv Hdg Ptch Roll Tx TIM xmitc xmitv dac VDD1 VDD3 VDC
	clear sv_name Hdg_name Ptch_name Roll_name Tx_name TIM_name...
		xmitc_name xmitv_name dac_name VDD1_name VDD3_name VDC_name
	clear sv_units Hdg_units Ptch_units Roll_units Tx_units TIM_units...
		xmitc_units xmitv_units dac_units VDD1_units VDD3_units VDC_units
	for i=1:4, 
		eval(['clear ',DNAMES(i,:),' ',DNAMES(i,:),'_units ',...
			DNAMES(i,:),'_name',DNAMES(i,:),'_mean',DNAMES(i,:),'_std']); 
	end
	clear vsave
end









