% 			ADCP_QA.M
% This routine examines the ADCP velocity data and does first a range
% check on the data and then does a gap filling 
% CNF 12/14/98
% CNF 1/11/99 modified so that a red line is plotted along the top of 
%	plots to indicate when bottom tracking is off and those 
%	ensembles that will later be dropped by rfmt4quoddy.m


% Minimum and maximum depths of good data
mxgdz = nanmax(nanmax(zc.*(ABS_U./ABS_U)));
mngdz = nanmin(nanmin(zc.*(ABS_U./ABS_U)));
[NBINS NPRFS] = size(zc);


% Plot the positions of the ADCP ensembles

figure(1)
clf
pos1 = [ 0.3000    0.15    0.6500    0.7500];
set(1, 'units', 'normal', 'position', pos1)

plot(nav(1,:),nav(2,:),'.')
alat = nanmean(nav(2,:)); f = cos(pi*alat/180);
set(gca,'PlotBoxAspectRatio',[1 f*1 1])
grid
xlabel('Longitude')
ylabel('Latitude')
title([cruise,': ',starttime,' to ',endtime])

% Perform a range check of the data
maxv = 150;
i = find(sqrt(ABS_U.^2 + ABS_V.^2) > maxv);
ABS_U(i) = nan*ones(size(i)); ABS_V(i) = nan*ones(size(i));

% Find bins in which there is no velocity which are embedded
% within otherwise good data and interpolate across the gap

for i = 1:NPRFS
  j = find(~isnan(ABS_U(:,i)));
  if length(j) > 1
    if (max(j)-j(1) + 1 > length(j))
      jj = diff(j);
      jjj = find(jj > 1);
      for k = 1:length(jjj)
        for kk = 1:jj(jjj(k))
          ABS_U(jjj(k)+kk,i) = ABS_U(jjj(k),i) + ...
	  (kk/jj(jjj))*(ABS_U(jjj(k)+jj(jjj),i) - ABS_U(jjj(k),i));
          ABS_V(jjj(k)+kk,i) = ABS_V(jjj(k),i) + ...
	  (kk/jj(jjj))*(ABS_V(jjj(k)+jj(jjj),i) - ABS_V(jjj(k),i));
        end
      end
      %break
    end
  end
end

clear i j k jj jjj kk

% Plot the velocity vectors
pos2 = [ 0.2000    0.10    0.7500    0.7500];
figure(2)
clf
set(2, 'units', 'normal', 'position', pos2)
axis([-75 -65 28 42]);
plot(nav(1,:),nav(2,:),'m-');
hold on
%hv=vecplot2(nav(1,:),nav(2,:),ABS_U(1,:),ABS_V(1,:),0.5,'m/s',-74.,30.)
hv=quiver(nav(1,:),nav(2,:),ABS_U(1,:),ABS_V(1,:),2,'.')
set(hv,'color','blue');
hold off
grid
%title(verbatim(['Observations: RMS = ',num2str(rms(ABS_U(1,:)+sqrt(-1.0)*ABS_V(1,:))),' m/s']))
title(['Observations: ADCP velocity, first bin'])
% Plot the velocities and backscatter intensity

pos2 = [ 0.2000    0.10    0.7500    0.7500];
colors = 'jet'; %choose a colorbar

% Determine time when bottom tracking is off so that those periods
% can be indicated on the plots by red +'s
ibtoff = find(isnan(bt_depth));

figure(3)
clf
set(3, 'units', 'normal', 'position', pos2)

crange = [-0.75 0.75];   % velocity range for color plot
arange = [-110 -80];     % range corrected backscatter intensity range

h=axes('units', 'normal', 'position', [0.0641 0.06 0.83 0.26]);
pcolor(ones(NBINS,1)*time,zc,Sv)
set(gca,'ydir','reverse','ylim',[0 50*(1+ floor(mxgdz/50))]);
caxis(arange)
eval(['colormap(', colors, ')']);
shading flat
set(h,'fontsize',8);
set(get(h,'title'),'fontsize',8);
set(get(h,'xlabel'),'fontsize',8);
set(get(h,'ylabel'),'fontsize',8);set(1, 'units', 'normal', 'position', pos1)
hold on
hoff=plot(time(ibtoff),zeros(size(time(ibtoff))),'+r'); set(hoff,'markersize',4)
plot(time,bt_depth)
hold off

xlabel([num2str(year) ' days']);
ylabel(['depth'])
title(['Range Corrected Backscatter Intensity (dB)']);
%set(1, 'units', 'normal', 'position', pos1)
x0 = time(1) - 0.05*(time(NPRFS) - time(1));
y0 = 1.15*(50*(1+ floor(mxgdz/50)));
h = text(x0,y0, ['BNL ' date]);
set(h,'fontsize',8);

%colorbar
h=axes('units', 'normal', 'position', [0.935 0.06 0.02 0.26]);
dcolor = arange(2)-arange(1);
pcolor([1:2], [arange(1):dcolor/31:arange(2)],[1:32;1:32]')
set(h,'fontsize',8);
set(get(h,'ylabel'),'fontsize',8);
shading flat
eval(['colormap(', colors, ')']);
ylabel(['dB'])


% Plot meridional velocity component

% Debugging section that estimates ship speed to check for 
% correlation with absolute water velocity
%dnav = diff(nav')';
%dt = diff(time)*24;
%T_spd = time(1:length(time)-1) + 0.5*dt;
%f = cos(mean(nav(2,:))*pi/180);
%i = find(T_spd < 324);
%us = 1.852*60*sqrt(f*f*dnav(1,:).^2 + dnav(2,:).^2)./dt;

h=axes('units', 'normal', 'position', [0.0641 0.38 0.83 0.26]);

pcolor(ones(NBINS,1)*time,zc,ABS_V)
set(gca,'ydir','reverse','ylim',[0 50*(1+ floor(mxgdz/50))]);
caxis(crange)
eval(['colormap(', colors, ')']);
shading flat
set(h,'fontsize',8);
set(get(h,'title'),'fontsize',8);
set(get(h,'xlabel'),'fontsize',8);
set(get(h,'ylabel'),'fontsize',8);

ylabel(['depth'])
title(['ADCP velocity: meridional (V)']);
hold on
% plot(T_spd(i),-us(i),'.'); set(gca,'ylim',[-50 50*(1+ floor(mxgdz/50))])
hoff=plot(time(ibtoff),zeros(size(time(ibtoff))),'+r'); set(hoff,'markersize',4)
plot(time,bt_depth)
hold off

%colorbar
h=axes('units', 'normal', 'position', [0.935 0.38 0.02 0.26]);
dcolor = crange(2)-crange(1);
pcolor([1:2], [crange(1):dcolor/31:crange(2)],[1:32;1:32]')
set(h,'fontsize',8);
set(get(h,'ylabel'),'fontsize',8);
shading flat
eval(['colormap(', colors, ')']);
ylabel(['m/s'])


% Plot the zonal velocity component

h=axes('units', 'normal', 'position', [0.0641 0.70 0.83 0.26]);

pcolor(ones(NBINS,1)*time,zc,ABS_U)
set(gca,'ydir','reverse','ylim',[0 50*(1+ floor(mxgdz/50))]);
caxis(crange)
shading flat
eval(['colormap(', colors, ')']);
set(h,'FontSize',8);
set(get(h,'title'),'fontsize',8);
set(get(h,'xlabel'),'fontsize',8);
set(get(h,'ylabel'),'fontsize',8);
ylabel(['depth'])
title(['ADCP velocity: zonal (U)']);
hold on
hoff=plot(time(ibtoff),zeros(size(time(ibtoff))),'+r'); set(hoff,'markersize',4)
plot(time,bt_depth)
hold off

%colorbar
h=axes('units', 'normal', 'position', [0.935 0.70 0.02 0.26]);
dcolor = crange(2)-crange(1);
pcolor([1:2], [crange(1):dcolor/31:crange(2)],[1:32;1:32]')
set(h,'FontSize',8);
set(get(h,'ylabel'),'fontsize',8);
shading flat
eval(['colormap(', colors, ')']);
ylabel(['m/s'])

% VK. clear alat arange colors crange dcolor f h mngdz mxgdz pos1 pos2
% VK. addition 2005/02/08
clear alat arange colors crange dcolor f h mngdz mxgdz pos1 pos2

return
figure(4)
clf
set(4, 'units', 'normal', 'position', pos2)

%crange = [-0.75 0.75];   % velocity range for color plot
arange = [-100 -50];     % range corrected backscatter intensity range

hh1=axes('units', 'normal', 'position', [0.1 0.16 0.83 0.76]);
%pcolor(ones(NBINS,1)*time,zc,Sv)
contourf(ones(NBINS,1)*time,zc,Sv,20)
set(gca,'ydir','reverse','ylim',[0 50*(1+ floor(mxgdz/50))]);
caxis(arange)
%eval(['colormap(', colors, ')']);
shading flat
grid
set(hh1,'fontsize',8);
set(get(hh1,'title'),'fontsize',8);
set(get(hh1,'xlabel'),'fontsize',8);
set(get(hh1,'ylabel'),'fontsize',8);set(1, 'units', 'normal', 'position', pos1)
hold on
hoff=plot(time(ibtoff),zeros(size(time(ibtoff))),'+r'); set(hoff,'markersize',4)
plot(time,bt_depth)
hold off

xlabel([num2str(year) ' days']);
ylabel(['depth'])
title(['Range Corrected Backscatter Intensity (dB)']);
%set(1, 'units', 'normal', 'position', pos1)
x0 = time(1) - 0.05*(time(NPRFS) - time(1));
y0 = 1.15*(50*(1+ floor(mxgdz/50)));
h = text(x0,y0, ['BNL ' date]);
set(h,'fontsize',8);

hh2=colorbar
set(hh2,'position', [0.92 0.16 0.02 0.56]);
set(gca,'position', [0.1 0.16 0.83 0.76]);
box on
%dcolor = arange(2)-arange(1);
%pcolor([1:2], [arange(1):dcolor/31:arange(2)],[1:32;1:32]')
set(hh2,'fontsize',8);
set(get(hh2,'ylabel'),'fontsize',8);
%shading flat
%eval(['colormap(', colors, ')']);
ylabel(['dB'])


clear alat arange colors crange dcolor f h mngdz mxgdz pos1 pos2
