From - Tue May 26 10:03:07 1998
Received: from midori.gso.uri.edu (midori.gso.uri.edu [198.116.10.79]) by whsun1.wh.whoi.edu (8.6.12/8.6.12) with ESMTP id JAA01884 for <jmanning@whsun1.wh.whoi.edu>; Tue, 26 May 1998 09:27:51 -0400
Received: (from john@localhost) by midori.gso.uri.edu (950413.SGI.8.6.12/950213.SGI.AUTOCF) id GAA24658 for jmanning@whsun1.wh.whoi.edu; Tue, 26 May 1998 06:32:15 -0700
Date: Tue, 26 May 1998 06:32:15 -0700
From: john@midori.gso.uri.edu (John Ryan)
Message-Id: <199805261332.GAA24658@midori.gso.uri.edu>
Apparently-To: jmanning@whsun1.wh.whoi.edu
X-Mozilla-Status: 8001

function [julian_day,mean_latitude,mean_longitude,ua,va,trans_temp,depth,w_bottom,error_bottom]=toabs(file)
%  script to convert output from bb2mat to absolute velocities using
%  bottom-track velocities.  This script reads in the output from bb2mat
% and scales all variables accordingly.
%  convert configuration data to meaningful expressions

%% Changed error variable name (output from bb2mat) to herror to avoid
%% confusion by matlab with the function error 
%% This also has to be changed in the matlab output files
%% JR 5/98

eval(['load ' file]);

% PG_THRESH=70;  % set percent good threshold at 70
% ERR_THRESH=0.1; % set error velocity threshold at 0.1 m/s
% ERRBOT_THRESH=0.02; % set bottom error velocity threshold at 0.02 m/s
% WBOT_THRESH=0.15; % set bottom w velocity threshold at 0.15 m/s
load bbadcp_thresholds

if (coord_system == 0) 
	coord_system='beam coordinates';
elseif(coord_system == 1)
	coord_system='earth coordinates';
elseif(coord_system == 2)
	coord_system='ship coordinates';
elseif(coord_system == 3)
	coord_system='instrument coordinates';
else
	coord_system='unknown';
end

if (orientation == 0)
	orientation='up';
elseif(orientation == 1)
	orientation='down';
else
	orientation='unknown';
end

if (beam_pattern == 0)
	beam_pattern='convex';
elseif(beam_pattern == 1)
	beam_pattern='concave';
else
	beam_pattern='unknown';
end

if (avg_method == 0)
	avg_method='temporal';
elseif(avg_method == 1)
	avg_method='spatial';
else
	avg_method='unknown';
end

% sound speed flag

if (ss_flag == 0)
	ss_flag='t/s';
elseif(ss_flag == 1)
	ss_flag='no';
else
	ss_flag='unknown';
end



% time between pings, , avgeraging interval converted to seconds (or meters
% if spatial averaging is used)

delta_t = delta_t/100.;
avg_interval = avg_interval/100.;

% bin length (WS command), blank after transit (WF command),  and
% adcp depth, converted to meters

bin_length = bin_length/100;
blank_length = blank_length/100;
adcp_depth = adcp_depth/100;

% convert magnetic variation, compass offset, and transducer misalignment
%  to degrees

magnetic_offset = magnetic_offset/1000.;
compass_offset = compass_offset/1000.;
trans_offset = trans_offset/1000.;

% intensity conversion factor, in dB/count, and absorption, in dB/m

intensity_scale = intensity_scale/1000;
adsorption=absorption/1000;

% salinity, in psu

salinity=salinity/1000.;

% sound speed, in m/s

sound_speed=sound_speed/1000.;

% convert bin 1 distance and transmit pulse length to meters

bin1_distance=bin1_distance/100.;
xmit_pulse_len=xmit_pulse_len/100.;

% convert distance measures to m

%east_ship_displ=east_ship_displ/(100.);
%north_ship_displ=north_ship_displ/(100.);


dist_travelled=dist_travelled/(100.);
x_displacements=x_displacements/(100.);
y_displacements=y_displacements/(100.);
path_length=path_length/100.;

% convert lat, lon to decimal degrees

mean_latitude=mean_latitude/(1000.*3600.);
mean_longitude=mean_longitude/(1000.*3600.);

% redo displacement calculation, RD says they are suspect

scale=1.1112e5;  % meters per degree of latitude
fact=cos((pi/180)*(mean(mean_latitude)));  % factor to scale longitude
east_ship_displ=fact*scale*(mean_longitude-...
	mean_longitude(1)*ones(size(mean_longitude)))+x_displacements(1);
north_ship_displ=scale*(mean_latitude-...
	mean_latitude(1)*ones(size(mean_latitude)))+y_displacements(1);

% convert day_sconds, elapsed time to days

elapsed_time = elapsed_time/(100.*86400.);
[year,julian_day]=julian(century_days,day_seconds);
%day_seconds = day_seconds/(100.*86400.);

% correct time to be UTC from EDT (add 4 hours)
% julian_day=julian_day+(4./24.)*ones(size(julian_day));

% convert pitch, roll, heading and transducer tempertature to degrees

pitch=pitch/100.;
roll=roll/100.;
heading=heading/100.;
trans_temp=trans_temp/100.;

% convert bottom track velocities to  m/s after converting 
% badflags to NaN
badpos=find(u_bottom<-32000);
u_bottom(badpos) = NaN*badpos;
u_bottom=u_bottom/1000.;
v_bottom(badpos) = NaN*badpos;
v_bottom=v_bottom/1000.;
w_bottom(badpos) = NaN*badpos;
w_bottom=w_bottom/1000.;
east_ship_vel=east_ship_vel/1000.;
north_ship_vel=north_ship_vel/1000.;

% convert displacements, ranges, to meters

range_1(badpos) = NaN*badpos;
range_1=range_1/100.;
range_2(badpos) = NaN*badpos;
range_2=range_2/100.;
range_3(badpos) = NaN*badpos;
range_3=range_3/100.;
range_4(badpos) = NaN*badpos;
range_4=range_4/100.;

% treat error_bottom separately, seems to allow 3 beam solutions

badpos=find(error_bottom<-32000);
error_bottom(badpos) = NaN*badpos;
error_bottom=error_bottom/1000.;

%find size of arrays

[m,n] = size(u);

% set up depth array

for i=1:m
	depth(i) = -(adcp_depth+bin1_distance+bin_length*(i-1));
end;

% first replace rd badflags with nans, convert to m/s
posbad=find(u<-32000);
u(posbad) = NaN*posbad;
u=u/1000.;
v(posbad) = NaN*posbad;
v=v/1000.;
w(posbad) = NaN*posbad;
w=w/1000.;

% do error velocity separately in case 3 beam solutions were allowed
posbad=find(herror<-32000);
herror(posbad) = NaN*posbad;
herror=herror/1000.;

% loop over u and v and  subtract off bottom velocity, badflag
% data below depth range to create absolute velocity estimates

ua=NaN*ones(m,n);
va=NaN*ones(m,n);

% despike w_bottom

%[badbt,goodub]=w_despike(julian_day',w_bottom');
%badbt=union(badub,badvb,badwb);

%u_bottom(badbt) = NaN*badbt;
%v_bottom(badbt) = NaN*badbt;
%w_bottom(badbt) = NaN*badbt;
% try a different approach - will median filter the smooth velocity records
% makes sense that changes in speed should occur smoothly

%u_bottom=medfilt1(u_bottom,3);
%v_bottom=medfilt1(v_bottom,3);

for i=1:n

	range=0.25*(range_1(i)+range_2(i)+range_3(i)+range_4(i));

	bottom_depth(i)=adcp_depth+range;

%   decrease range to account for bottom interference

	range = cos(beam_angle/180*pi)*range;
	if(range <= 0)  % presume this means we are too deep to bottom track
		jmin=m;
	else
		jmin=round((range-bin1_distance)/bin_length);
		if(jmin <= 0) jmin =1; end;
	end
	if(jmin >= m | isnan(jmin)) 
		jmin=m;  % also assuming we are deep
	else   % badflag existing arrays below jmin depth
		u(jmin+1:m,i)=NaN*ones(m-jmin,1);
		v(jmin+1:m,i)=NaN*ones(m-jmin,1);
                w(jmin+1:m,i)=NaN*ones(m-jmin,1);
                herror(jmin+1:m,i)=NaN*ones(m-jmin,1);
	end;
	ubar(i)=mean(u(1:jmin,i)) - u_bottom(i);
	vbar(i)=mean(v(1:jmin,i)) - v_bottom(i);
end;
% despike ubar and vbar 

[badub,goodub]=harv_de_spike(julian_day',ubar');
[badvb,goodvb]=harv_de_spike(julian_day',vbar');

badbt=union(badub,badvb);

u_bottom(badbt) = NaN*badbt;
v_bottom(badbt) = NaN*badbt;
w_bottom(badbt) = NaN*badbt;

for i=1:n
% fill in absolute velocities where valid data exist

%  reject data based on percent good and bottom tracking

	posrej=find(pg(:,i)<PG_THRESH);
		
	if((abs(w_bottom(i))<WBOT_THRESH) & (abs(error_bottom(i))<ERRBOT_THRESH))
  		ua(:,i) = u(:,i) - u_bottom(i)*ones(m,1);
		ua(posrej,i) = NaN*posrej;		
	end
	if((abs(w_bottom(i))<WBOT_THRESH) & (abs(error_bottom(i))<ERRBOT_THRESH)) 
		va(:,i) = v(:,i) - v_bottom(i)*ones(m,1);
		va(posrej,i) = NaN*posrej;
 	end;
end;

% convert echo intensity to dB

% e1=e1*intensity_scale;
% e2=e2*intensity_scale;
% e3=e3*intensity_scale;
% e4=e4*intensity_scale;

% echo_inten=ones(m,n);
% echo_inten = 0.25*(e1+e2+e3+e4);

% do these for error checking on horizontal velocities, currently
% hardwired to 10 cm/s

posbad=find(abs(herror)>ERR_THRESH);
ua(posbad) = NaN*posbad;
va(posbad) = NaN*posbad;
herror(posbad) = NaN*posbad;


