%function [ua,va,julian_day,pg,depth,echo_inten]=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

%eval(['load ' file]);

        if exist('pg')==0;
          disp('WARNING: percent good criteria not possible')
        else;
          disp('Note: percent good criteria  possible')
        end;

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


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
disp(['coord_system = ',coord_system]);

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
disp(['avg_method = ',avg_method]);

% 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.);

%%%%%%%%%%%%%%%%%%%%
%%%% JManning modification (so not to confuse w/other "julian" subroutine calls)
%[year,julian_day]=julian(century_days,day_seconds);
%year=input('Enter a 2-digit year: ');
year=97;
day_seconds2 = day_seconds/(100.*86400.);
julian_day=century_days-year*365-floor(year/4)+day_seconds2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% correct time to be UTC from EDT (add 4 hours)
% Note: the following three corrections should be made for al9805

disp('Warning: 4 hours NOT added to get a GMT time.')
disp('         We are assuming adcp time is  already gmt');
%julian_day=julian_day+(4./24.)*ones(size(julian_day));

disp('WARNING.... WARNING....NOT.. addding 1 hour fix to al9805 time')
%julian_day=julian_day+1/24;

disp('WARNING.... WARNING....NOT..[A subtracting 1 day from al9805 time')
%julian_day=julian_day-1;


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(error<-32000);
error(posbad) = NaN*posbad;
error=error/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);
wa=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);
                error(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');

%%%%%%% JManning modification (since I didn't have "union")
%badbt=union(badub,badvb);
% find the elements of bad vbars
all=[badub;badvb];

if all~=[];
  badbt=distinct(all);
  u_bottom(badbt) = NaN*badbt;
  v_bottom(badbt) = NaN*badbt;
  w_bottom(badbt) = NaN*badbt;
end;

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

%  reject data based on percent good and bottom tracking

%%%% JPM modification

%    note some cruises, like al9801, may not have a "pg" variable
%        if exist('pg')==0;
%          disp('WARNING: percent good criteria not possible')
%        else;
%          disp('Note: percent good criteria  possible')
%        end;
	posrej=find(pg(:,i)<PG_THRESH);
		
	if((abs(w_bottom(i))<WBOT_THRESH) & (abs(error_bottom(i))<ERRBOT_THRESH))
  		ua(1:jmin,i) = u(1:jmin,i) - u_bottom(i)*ones(jmin,1);
		ua(posrej,i) = NaN*posrej;		
	end
	if((abs(w_bottom(i))<WBOT_THRESH) & (abs(error_bottom(i))<ERRBOT_THRESH)) 
		va(1:jmin,i) = v(1:jmin,i) - v_bottom(i)*ones(jmin,1);
		va(posrej,i) = NaN*posrej;
        end;
%%%% JPM addition of "wa" vertical velocity???
	if((abs(w_bottom(i))<WBOT_THRESH) & (abs(error_bottom(i))<ERRBOT_THRESH)) 
		wa(1:jmin,i) = w(1:jmin,i);
		wa(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(error)>ERR_THRESH);
ua(posbad) = NaN*posbad;
va(posbad) = NaN*posbad;
wa(posbad) = NaN*posbad;
error(posbad) = NaN*posbad;

