function [Ve, Vn, Vup, H, P, R, mavszeros] = mavs2earth(bass1, bass2, burst_per_avg,...
   bass_error1, bass_error2, mavsnum, mavszeros, isburst);
%function [Ve, Vn, Vup, H, P, R, mavszeros] = mavs2earth(bass1, bass2, burst_per_avg,...
%   bass_error1, bass_error2, mavszeros);
% where 	bass1 = [A, B, C, D] are the four MAVS axis currents
%			bass2 = [My, Mx, P, R] 
%				My = the magnetic cosine in the Y direction * 1000
%				Mx = the magnetic cosine in the X direction * 1000
%				P, R   are the pitch and roll in 100ths of degrees
%			Ve, Vn, Vup velocity vectors in earth coordinates, mm/s magnetic.
%			H is the measured heading, magnetic, in degrees, computed from Mx and My
%			P is the measured pitch, in hundredths of degrees
%			R is the measured roll, in hundredths of degrees
%			mavszeros are the applied zero offsets, if provided
%			isburst is a logical variable which is 1 if the data are 
%				individual burst samples instead of averages.
%
% All inputs are required except mavszeros.  User will be prompted if mavszeros
%			are not provided.

%	written by Marinna Martini for the U.S. Geological Survey
%			Woods Hole Field Center with help from Sand Williams at WHOI.
% 	mavs2earth version 0.0	5/18/00

% 	modified for burst data, 20-Feb-2001 Fran Hotchkiss

disp ('mavs2earth.m is running')

% compute average from accumulated sums
if isburst
   b1norm = bass1;
   b2norm = bass2;
else
   b1norm = normbass(bass1, burst_per_avg, bass_error1);
   b2norm = normbass(bass2, burst_per_avg, bass_error2);
end

if exist('mavszeros') ~= 1, 
	prompt = {'Axis A:','Axis B:','Axis C:','Axis D:'};
	def = {'0','0','0','0'};
	dlgTitle = 'Enter MAVS zero offsets';
	lineNo = 1;
	mavszeroscell = inputdlg(prompt,dlgTitle,lineNo,def);
	mavszeros = str2num(char(mavszeroscell)); 
end
A = b1norm(:,1)-mavszeros(1);
B = b1norm(:,2)-mavszeros(2);
C = b1norm(:,3)-mavszeros(3);
D = b1norm(:,4)-mavszeros(4);
My = b2norm(:,1);
Mx = b2norm(:,2);
P = b2norm(:,3);	% pitch is in 100ths of degrees and needs no more massage
R = b2norm(:,4);	% roll is in 100ths of degrees and needs no more massage
clear bass1 bass2

% Process MAVS
MAVSfactor = 0.046;	% MAVS count * MAVSfactor = current cm/s

% remove points where data transfer failed.
% experiments show that this occurs 0.44% of the time
% and that it occurs through all the data, so test one axis
% remember which they are to flag them again after the calculations
BadFlag = -999999;
idxbad = find(A==-999999);

% the geometry of the axes are that each is inclined 45 degrees to 
% the horizontal (the plane perpendicular to the sensor tube) and 
% they are arranged 90 degrees in azimuth clockwise looking down
% (away from the instrument case).

% convert to instrument coordinates
% units of result are in mm/s
G = 4193;	% 
% U = velocity in the X direction (EAST)
U = (A+D-B-C).*(1000/G);
% V = velocity in Y direction (NORTH)
V = (A+B-C-D).*(1000/G);
% W = velocity in the Z direction
W = (A+B+C+D).*(707/G);

% apply pitch and roll
% convert pitch and roll, given in 100ths of degrees, into radians
Prad = (P./100).*(pi/180);
Rrad = (R./100).*(pi/180);
Vp = V.*cos(Prad)-W.*sin(Prad);
Wp = W.*cos(Prad)+V.*sin(Prad);
Ur = U.*cos(Rrad)-Wp.*sin(Rrad);
Wpr = Wp.*cos(Rrad)+U.*sin(Rrad);

% apply the compass rotation
Ve = (My.*Ur-Mx.*Vp)./1000;
Vn = (My.*Vp+Mx.*Ur)./1000;
Vup = Wpr;

if 0,
% do it without pitch and roll
Ve = (My.*U-Mx.*V)./1000;
Vn = (My.*V+Mx.*U)./1000;
Vup = W;
end

% insert the bad data flags
Baddata = zeros(size(idxbad))+BadFlag;
Ve(idxbad) = Baddata;
Vn(idxbad) = Baddata;
Vup(idxbad) = Baddata;

% metadata
% heading
% Mx and My are the direction cosines
   H=atan(-Mx./My).*(180/pi);
   % Quad I and III:  H>0
   % Quad I OK as is, find Quad III Mx<0
   idxquad = find(H>0 & Mx<0);
   H(idxquad) = H(idxquad)+180;
   % Quad II and IV:  H<0
   % Quad IV Mx>0
   idxquad = find(H<0 & Mx>0);
   H(idxquad) = H(idxquad)+180;
   % Quad II
   idxquad = find(H<0 & Mx<0);
   H(idxquad) = 360+H(idxquad);
   
   disp ('mavs2earth.m completed')