%======================================================================
% L O A D S A D C P . M
% doc: Sun Jun 27 23:42:04 2004
% dlm: Fri Mar 5 15:49:11 2010
% (c) 2004 ladcp@
% uE-Info: 133 14 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================
% CHANGES BY ANT:
% Jun 30, 2004: - added Figure 9 title
% Jul 9, 2004: - made lack of SADCP data in time window into real warning
% Jan 7, 2009: - tightened use of exist()
% Mar 4, 2010: - changed default of p.sadcp_dtok from 5min to zero
function [di,p]=loadsadcp(f,di,p)
% function [di,p]=loadsadcp(f,di,p)
%
%
%
% need to make array di.svel=[depth,U,V,Vel_error];
% will use array di.sadcp_lon
% di.sadcp_lat
%
%
% find SADCP profiles within time of station +/- slack (in fractional days)
p = setdefv(p,'sadcp_dtok',0);
if existf(f,'sadcp')==1
if exist(f.sadcp,'file')
disp(['LOADSADCP: load SADCP data file ',f.sadcp])
load(f.sadcp);
% By Now you should have:
% tim_sadcp(t) - Julian Days
% lat_sadcp(t) - Degrees N
% lon_sadcp(t) - Degrees E
% u_sadcp(z,t) - m/s
% v_sadcp(z,t) - m/s
% z_sadcp(z,1) - Meter (Positive Depth)
ii=find(tim_sadcp>(julian(p.time_start)-p.sadcp_dtok) & tim_sadcp<(julian(p.time_end)+p.sadcp_dtok));
if length(ii)<1
warn = ' no SADCP data found in time window';
p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
disp(warn)
return;
end
% check if position matches
if abs(p.lon)+abs(p.lat)~=0 & (abs(mean(lon_sadcp(ii))-p.lon)>0.1/cos(p.lat*pi/360) |...
abs(mean(lat_sadcp(ii))-p.lat)>0.1)
disp(' position of SADCP data is more than 0.1 degree away from LADCP')
disp(' no SADCP data used')
figure(9)
clf
orient tall
subplot(211)
text(0,0.5,' SADCP data is more than 0.1 degree away from LADCP',...
'color','r','fontsize',16,'fontweight','bold')
axis off
subplot(212)
plot(lon_sadcp(ii),lat_sadcp(ii),'rp')
hold on
if existf(di,'slon')
plot(di.slon,di.slat,'.g')
end
if abs(p.lon)+abs(p.lat)~=0
plot(p.lon,p.lat,'pr')
plot(p.poss(3)+p.poss(4)/60,p.poss(1)+p.poss(2)/60,'bp','markersize',15)
plot(p.pose(3)+p.pose(4)/60,p.pose(1)+p.pose(2)/60,'kp','markersize',15)
end
title('ship nav (g.) start (bp) end (kp) SADCP (rp)')
xlabel('longitude')
ylabel('latitude')
streamer([p.name,' Figure 9']);
pause(0.001)
return;
end
% interpolate SADCP navigation to LADCP time series
di.sadcp_lon=interp1(tim_sadcp,lon_sadcp(:),di.time_jul);
di.sadcp_lat=interp1(tim_sadcp,lat_sadcp(:),di.time_jul);
% set position from SADCP nav
if abs(p.lon)+abs(p.lat)==0
slat=di.sadcp_lat(1);
slon=di.sadcp_lon(1);
elat=di.sadcp_lat(end);
elon=di.sadcp_lon(end);
p.poss=[fix(slat), (slat-fix(slat))*60, fix(slon), (slon-fix(slon))*60];
p.pose=[fix(elat), (slat-fix(elat))*60, fix(elon), (slon-fix(elon))*60];
end
% if no other ship navigation exists, use SADCP navigation
if existf(di,'slon')==0
di.slon=di.sadcp_lon;
di.slat=di.sadcp_lat;
else
if sum(isfinite(di.slon+di.slat))==0
di.slon=di.sadcp_lon;
di.slat=di.sadcp_lat;
end
end
if length(find(isfinite(u_sadcp(:,ii))))<1,
disp(' no finite SADCP data found ')
return;
end
disp([' found ',int2str(length(ii)),' SADCP profiles '])
u=squeeze(u_sadcp(:,ii));
v=squeeze(v_sadcp(:,ii));
% compute velocity standard deviation
if numel(u) == length(u)
v_err = u*0+0.1;
else
v_err=(nstd(u')+nstd(v'))';
u=meannan(u')';
v=meannan(v')';
end
ij=find(v_err==0);
v_err(ij)=0.1;
nvel=sum(isfinite(u+v)')'+v_err*0;
v_err=v_err*max(nvel)./nvel;
z=z_sadcp;
% make output array
izok=find(isfinite(u+v) & z<(p.maxdepth));
if isfinite(p.zbottom),
izok=izok(find(z(izok)<(p.zbottom-30)));
end
di.svel=[z(izok),u(izok),v(izok),v_err(izok)];
else
disp([' can not find SADCP data file:',f.sadcp])
return
end
else
disp(' no SADCP data file given')
return
end
% plot some results
figure(9)
clf
subplot(211)
plot(squeeze(u_sadcp(:,ii)),-z_sadcp,'-r')
hold on
plot(squeeze(v_sadcp(:,ii)),-z_sadcp,'--g')
ax=axis;
axis(ax)
grid
streamer(p.name);
ylabel('depth [m]')
xlabel('velocity [m/s]')
title('SADCP U(-r) and V(--g) profiles')
subplot(212)
plot(di.slon,di.slat,'.g')
hold on
plot(lon_sadcp(ii),lat_sadcp(ii),'rp')
plot(p.poss(3)+p.poss(4)/60,p.poss(1)+p.poss(2)/60,'bp','markersize',15)
plot(p.pose(3)+p.pose(4)/60,p.pose(1)+p.pose(2)/60,'kp','markersize',15)
title('ship nav (g.) start (bp) end (kp) SADCP (rp)')
xlabel('longitude')
ylabel('latitude')
pause(0.001)
% ------------------------------------------------------------
function y = nstd(x,flag,dim)
%NSTD Standard deviation, ignoring NaN.
% Same as STD, but NaN's are ignored.
%
% Copyright (c) 1997 by Toby Driscoll.
% Adapted from STD.M, written by The MathWorks, Inc.
% added backward compatibility G.Krahmann, LODYC Paris
if nargin<2, flag = 0; end
if nargin<3,
dim = min(find(size(x)~=1));
if isempty(dim), dim = 1; end
end
% Avoid divide by zero.
if size(x,dim)==1, y = zeros(size(x)); return, end
tile = ones(1,max(ndims(x),dim));
tile(dim) = size(x,dim);
xc = x - repmat(meannan(x),tile); % Remove mean
mask = isnan(xc);
xc(mask) = 0;
s = sum(~mask,dim);
s(s==0) = NaN;
if flag,
y = sqrt(sum(conj(xc).*xc,dim)./s);
else
z = (s==1);
s(z) = Inf;
y = sqrt(sum(conj(xc).*xc,dim)./(s-1));
end