function epname = bass_avg(infile, outfile)
% function epname = bass_avg(infile, outfile)
% calibrates BASS variables from the average records of MIDAS data.
% epname is a cell array of the names of the variables created.
% infile is the name of a netcdf file of raw data, to which
%		calibration constants have been attached as attributes.
% outfile is the name of an EPIC compliant netcdf file containing
%		the variables listed in epname.
% Fran Hotchkiss 7 Mar 2000.
% Updated to make proper use of error count.  22 Jun 2000 FSH
nvarout = 0;
vardesc = '';
epname{1} = '';
calatts = '';
varcomment = '';
calcomment = '';

nc = netcdf(infile);
if isempty(nc)
   [infile,inpath] = uigetfile('*.cdf','Choose input file');
   nc = netcdf(infile);
   if isempty(nc)
      disp (['Cannot open input file ' infile]);
      disp (['Execution ends.']);
      return;
   end
end

% Which variables should be processed?
if running(batch)
   str = get(batch);
   eval(['nvarin = ' str ';']);      
   for i = 1:nvarin;
      str = get(batch);
      eval(['chosen{i} = ' str ';']);
   end
else
   
	invars = var(nc);
	for i = 1:length(invars)
	   eval (['varlist.',name(invars{i}),'={''checkbox'' 0};'])
	end
	vlist = varlist;
	varlist = guido(vlist,'Which variables should be processed?');
	nvarin = 0;
	for i = 1:length(invars)
	   eval (['tf = varlist.',name(invars{i}),'{2};'])
	   if tf
	      nvarin = nvarin+1;
	      chosen{nvarin}=name(invars{i});
	   end   
   end
end

% Find and calibrate tilt first.
for i = 1:length(chosen)
   if strncmp('tiltx',chosen{i},5)
      ivar = nc{chosen{i}};
		% Input data and calibration constants for x tilt.
		rtx = ivar(:);
		a2d = ivar.AtoD_factor(1);
		txf = ivar.voltage_divider_factor(1);
		txa = ivar.cal_A(1);
		txb = ivar.cal_B(1);
		txc = ivar.cal_C(1);
		txd = ivar.cal_D(1);
		% Calibrate x tilt.
		Cal = [txd txc txb txa];
		rtx = a2d * txf .* rtx;
		txdeg = polyval(Cal,rtx);
      nvarout = nvarout+1;
      epname{nvarout} = 'tiltx_4017';
      inname{nvarout} = name(ivar);
      vardesc = [ vardesc ':tiltx'];
      datname{nvarout} = 'txdeg';
      calatts.tiltx_4017 = {'AtoD_factor','voltage_divider_factor',...
            'cal_A','cal_B','cal_C','cal_D'};
   end
end

% Find and calibrate y tilt next.
for i = 1:length(chosen)
   if strncmp('tilty',chosen{i},5)
      ivar = nc{chosen{i}};
		% Input data and calibration constants for y tilt.
		rty = ivar(:);
		a2d = ivar.AtoD_factor(1);
		tyf = ivar.voltage_divider_factor(1);
		tya = ivar.cal_A(1);
		tyb = ivar.cal_B(1);
		tyc = ivar.cal_C(1);
		tyd = ivar.cal_D(1);
		% Calibrate y tilt.
		Cal = [tyd tyc tyb tya];
		rty = a2d * tyf .* rty;
		tydeg = polyval(Cal,rty);
      nvarout = nvarout+1;
      epname{nvarout} = 'tilty_4018';
      inname{nvarout} = name(ivar);
      vardesc = [ vardesc ':tilty'];
      datname{nvarout} = 'tydeg';
      calatts.tilty_4018 = {'AtoD_factor','voltage_divider_factor',...
            'cal_A','cal_B','cal_C','cal_D'};
   end
end

% Find and calibrate compass next.
for i = 1:length(chosen)
   if strncmp('comp',chosen{i},4)
      ivar = nc{chosen{i}};
		% Input data and calibration constants for compass.
		rc = ivar(:);
		a2d = ivar.AtoD_factor(1);
		ca0 = ivar.cal_A0(1);
		ca1 = ivar.cal_A1(1);
		ca2 = ivar.cal_A2(1);
		ca3 = ivar.cal_A3(1);
		mv = nc.magnetic_variation(1);
		% Calibrate compass.
		Cal = [ca3 ca2 ca1 ca0];
		rc = a2d .* rc;
		cdeg = polyval(Cal,rc);
		% Correct for magnetic variation
		cdeg = cdeg + mv;
		low_cdeg = find(cdeg < 0.);
		cdeg(low_cdeg) = cdeg(low_cdeg) + 360.;
		high_cdeg = find(cdeg > 360.);
		cdeg(high_cdeg) = cdeg(high_cdeg) - 360.;
      nvarout = nvarout+1;
      epname{nvarout} = 'comp_1406';
      inname{nvarout} = name(ivar);
      vardesc = [ vardesc ':comp'];
      datname{nvarout} = 'cdeg';
      calatts.comp_1406 = {'AtoD_factor',...
            'cal_A0','cal_A1','cal_A2','cal_A3'};
   	varcomment.comp_1406 = 'corrected for magnetic variation';
   end
end

% Find and calibrate BASS now.
for i = 1:length(chosen)
   if strncmp('bass',chosen{i},4)
      % Choose BASS calibration routine and
      %		allow use of constant compass and/or tilts.
      BASS.axes = {{'All','A omitted','B omitted','C omitted',...
               'D omitted'},1};
      BASS.zeroes={{'pre-deployment','post-deployment'},2};
      BASS.use_const_compass={'checkbox',0};
      BASS.compass_degrees = 0.;
      BASS.use_constant_x_tilt={'checkbox',0};
      BASS.x_tilt_degrees = 0.;
      BASS.use_constant_y_tilt={'checkbox',0};
      BASS.y_tilt_degrees = 0.;
		if running(batch)
      	str = get(batch);
   		eval(['BASS.axes{2} = ' str ';']);      
      	str = get(batch);
   		eval(['BASS.zeroes{2} = ' str ';']);      
      	str = get(batch);
   		eval(['BASS.use_const_compass{2} = ' str ';']);      
	      str = get(batch);
   		eval(['BASS.compass_degrees = ' str ';']);      
	      str = get(batch);
   		eval(['BASS.use_constant_x_tilt{2} = ' str ';']);      
	      str = get(batch);
   		eval(['BASS.x_tilt_degrees = ' str ';']);      
	      str = get(batch);
   		eval(['BASS.use_constant_y_tilt{2} = ' str ';']);      
	      str = get(batch);
   		eval(['BASS.y_tilt_degrees = ' str ';']);      
      else
	      B = BASS;
         BASS = guido(B,'Select BASS calibration parameters');
	      if isempty(BASS)
  		 		disp ('BASS Parameter Menu must be completed. Execution ends.')
 	  		return
			end
      end
      
      ivar = nc{chosen{i}};
		% Input data and calibration constants for BASS.
		a = ivar(:,1,1,1);
		b = ivar(:,1,2,1);
		c = ivar(:,1,3,1);
      d = ivar(:,1,4,1);
      errvarname = ['bass_error_' chosen{i}(6)]
      errvar = nc{errvarname};
		a_err = errvar(:,1,1,1);
		b_err = errvar(:,1,2,1);
		c_err = errvar(:,1,3,1);
      d_err = errvar(:,1,4,1);
		switch BASS.zeroes{2}
		case 1
			za = ivar.pre_zero_offset_A(1);
			zb = ivar.pre_zero_offset_B(1);
			zc = ivar.pre_zero_offset_C(1);
         zd = ivar.pre_zero_offset_D(1);
         calcomment = [calcomment 'Pre zeroes used. '];
		otherwise
			za = ivar.post_zero_offset_A(1);
			zb = ivar.post_zero_offset_B(1);
			zc = ivar.post_zero_offset_C(1);
      	zd = ivar.post_zero_offset_D(1);
         calcomment = [calcomment 'Post zeroes used. '];
      end
		nf = ivar.norm_factor(1);
		pma = ivar.pod_misalignment_angle(1);
		if isempty(pma) ;
		   pma = 0.;
		end
		% BASS data are sums instead of averages, so divide by number summed.
      nburst = nc{'burst_per_avg'}(:);
		a = a ./ (nburst-a_err);
		b = b ./ (nburst-b_err);
		c = c ./ (nburst-c_err);
		d = d ./ (nburst-d_err);
		% Normalize data using zero offsets.
		a = a-za;
		b = b-zb;
		c = c-zc;
		d = d-zd;
		% Call a function to compute velocity in pod coordinates.
		switch BASS.axes{2}
		case 2
			[upod,vpod,wpod,epod] = BASSuvw_a(a,b,c,d,nf);
         calcomment = [calcomment 'A axis omitted. '];
		case 3
			[upod,vpod,wpod,epod] = BASSuvw_b(a,b,c,d,nf);
         calcomment = [calcomment 'B axis omitted. '];
   	case 4
      	[upod,vpod,wpod,epod] = BASSuvw_c(a,b,c,d,nf);
         calcomment = [calcomment 'C axis omitted. '];
		case 5
         [upod,vpod,wpod,epod] = BASSuvw_d(a,b,c,d,nf);
         calcomment = [calcomment 'D axis omitted. '];
  		otherwise
      	[upod,vpod,wpod,epod] = BASSuvw(a,b,c,d,nf);
   	end
		% Call a function to rotate velocity to logger coordinates.
		[ulog,vlog] = uv_rotate(upod,vpod,pma);
		% Call a function to rotate velocity into earth coordinates.      BASS.use_compass_data={'checkbox',1};
      if BASS.use_const_compass{2};
         cdeg = BASS.compass_degrees;
         calcomment = [calcomment 'Constant value of '...
               num2str(cdeg) ' deg used for compass. '];
      end
      if BASS.use_constant_x_tilt{2};
         txdeg = BASS.x_tilt_degrees;
         calcomment = [calcomment 'Constant value of '...
               num2str(txdeg) ' deg used for x tilt. '];
      end
      if BASS.use_constant_y_tilt{2};
         tydeg = BASS.y_tilt_degrees;
         calcomment = [calcomment 'Constant value of '...
               num2str(tydeg) ' deg used for y tilt. '];
      end
		txrad = txdeg .* (pi/180.);
		tyrad = tydeg .* (pi/180.);
		crad = cdeg .* (pi/180.);
		[uearth,vearth,wearth] = pod2earth (ulog,vlog,wpod,txrad,tyrad,crad);
		% Now make speed and direction.
		[vdir,vspd] = cart2pol(vearth,uearth);
		ilow = find(vdir < 0);
		vdir(ilow) =vdir(ilow) + 2*pi;
		vdir = vdir .* 180 ./ pi;
      ebassname = {'u_1205','v_1206','CD_310','CS_300','w_1204','Werr_1201'};
      rbassname = {'uearth','vearth','vdir','vspd','wearth','epod'};
      for inc = 1:6;
         ivarout = inc + nvarout;
         inname{ivarout} = name(ivar);
         epname{ivarout} = ebassname{inc};
         datname{ivarout} = rbassname{inc};
      end
      vardesc = [ vardesc ':u:v:CD:CS:w:Werr'];
      nvarout = nvarout+6;
   end
end

% Make EPIC time from MIDAS raw data set time.
% Input data and base time for time.
ivar = nc{'time'};
rawdat = ivar(:);
yr = ivar.base_date(1);
mo = ivar.base_date(2);
da = ivar.base_date(3);
hr = ivar.base_date(4);
mi = ivar.base_date(5);
sc = ivar.base_date(6);
% Find decimal day numbers for base time and time data.
base = datenum(yr,mo,da,hr,mi,sc);
rtime = base + rawdat./86400;
% Call function to make EPIC time.
alltime = EP_Time(rtime);
time = alltime(:,1);
time2 = alltime(:,2);
stime = ep_datenum(alltime(1,:));
[n m] = size(alltime);
ltime = ep_datenum(alltime(n,:));

% Get mooring log parameters.
if running(batch)
	str = get(batch);
   eval(['epicat.experiment = ' str ';']);      
	str = get(batch);
   eval(['epicat.project = ' str ';']);      
	str = get(batch);
   eval(['epicat.comment = ' str ';']);      
	str = get(batch);
   eval(['epicat.BASS_pod_depth = ' str ';']);      
else
	epicat.experiment = 'ECOHAB';
	epicat.project = 'WHFC';
	epicat.comment = 'Boston BASS MIDAS';
	epicat.BASS_pod_depth = 2.9;
	e = epicat;
   epicat = guido(e,'Input mooring log parameters.');
end


% Open output cdf
outc = netcdf(outfile,'noclobber');
if isempty(outc)
   [outfile,outpath] = uiputfile('*.nc','Choose output file');
   outc = netcdf(outfile,'noclobber');
   if isempty(outc)
      disp (['Cannot open output file ' outfile]);
      disp (['Data will be parked in temporary.nc']);
      outc = netcdf('temporary.nc','clobber');
   end
end

% Copy global attributes from ep_standard.nc to output cdf.
eps = netcdf('ep_standard.nc');
if isempty(eps)
   [epsfile,epspath] = uigetfile('*.nc','Choose epic standard file');
   eps = netcdf(epsfile);
   if isempty(eps)
      disp (['Cannot open epic standard file ' epsfile]);
      disp (['Execution ends.']);
      ncclose;
      return;
   end
end
stdatts = att(eps);
for i = 1, length(stdatts);
   copy(stdatts(i),outc);
end

%% Modify Global attributes:
 
outc.DATA_ORIGIN = nc.data_source(:);
outc.EXPERIMENT = epicat.experiment;
outc.PROJECT = epicat.project;
outc.MOORING = nc.mooring_number(:);
outc.DELTA_T = ncchar(num2str(nc.average_interval(1)));
outc.DATA_CMNT = epicat.comment;
outc.WATER_DEPTH = nc.water_depth(1);
outc.DESCRIPT = nc.comment(:);
outc.FILL_FLAG = nclong(0);
outc.VAR_FILL = ncfloat(NaN);
history =['Calibrated using bass_avg.m. ' calcomment ':' nc.history(:)];  
outc.history = history;
outc.start_time = datestr(stime,0);
outc.stop_time = datestr(ltime,0);
outc.latitude = dms2deg(nc.lat_ddmmss(:));
outc.longitude = -dms2deg(nc.lon_ddmmss(:));
outc.water_depth = nc.water_depth(1);
outc.CREATION_DATE = ncchar(datestr(now,0));
outc.INST_TYPE = nc.instrument_type(:);
outc.instrument_number = nc.instrument_number(:);
outc.magnetic_variation = nc.magnetic_variation(1);
vardesc(1) = '';
outc.VAR_DESC = ncchar(vardesc);
i = strmatch('u',epname);
if ~isempty(i);
	ivar = nc{inname{i}};
	outc.BASS_zero_offset_A = za;
	outc.BASS_zero_offset_B = zb;
	outc.BASS_zero_offset_C = zc;
	outc.BASS_zero_offset_D = zd;
	outc.BASS_norm_factor = ivar.norm_factor(1);
   outc.pod_misalignment_angle = ivar.pod_misalignment_angle(1);
end

%% Dimensions:
outc('time') = 0;
outc('depth') = 1;
outc('lon') = 1;
outc('lat') = 1;
 
 %% Variables and attributes:
copy(eps{'time'},outc,0,1);
copy(eps{'time2'},outc,0,1);
copy(eps{'depth'},outc,0,1);
copy(eps{'lon'},outc,0,1);
copy(eps{'lat'},outc,0,1);

for i = 1: nvarout;
   ivar = eps{epname{i}};
   if ~isempty(ivar)
      copy(ivar,outc,0,1);
   else
      display (['Variable ' epname{i} ' not found in ep_standard.nc. Execution ends.'])
      ncclose 
      return
   end
   ivar = nc{inname{i}};
	ovar = outc{epname{i}};
	ovar.sensor_type = ivar.instrument_type(:);
	ovar.sensor_depth = ncfloat(epicat.BASS_pod_depth);
	ovar.serial_number = ivar.instrument_number(:);
	eval (['ovar.minimum = ncfloat(min(' datname{i} '));']);
   eval (['ovar.maximum = ncfloat(max(' datname{i} '));']);
   try
   	str = epname{i};
      eval (['attstr = calatts.' str ';']);
      for icalatt = 1:length(attstr);
         str = attstr{icalatt};
         eval(['ovar.'str '= ivar.'str '(1);'])
      end
   end
   try
   	str = epname{i};
      eval(['str=varcomment.'epname{i} ';'])
      ovar.comment = ncchar(str );
   end
end
close (eps)
endef(outc)

n=length(time);
outc{'time'}(1:n) = time;
outc{'time2'}(1:n) = time2;
outc{'lat'}(1) = nc{'lat'}(1);
outc{'lon'}(1) = nc{'lon'}(1);
outc{'depth'}(1) = epicat.BASS_pod_depth;
for i = 1:nvarout
	ovar = outc{epname{i}};
   eval (['ovar(1:n) = ' datname{i} ';' ]);
end
close (nc);
close (outc);