function epname = mavs_bst(infile, outfile)
% function epname = mavs_bst(infile, outfile)
% calibrates MAVS variables and pressure from burst 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
%		Reynolds Stress and pressure standard deviation.
% Based on midas_bst as written by Fran Hotchkiss on 7 Mar 2000, 
% 		modified to add P_4027 on 13 Jun 2000.
% Fran Hotchkiss 20 Feb 2001.

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
chosen{:}%;

% Choose time parameters
if running(batch)
	str = get(batch);
   eval(['BURST.average_interval = ' str ';']);      
	str = get(batch);
   eval(['BURST.first_average_minutes = ' str ';']);      
	str = get(batch);
   eval(['BURST.first_average_seconds = ' str ';']);      
else
	BURST.average_interval = nc.average_interval(1);
	BURST.first_average_minutes = 0;
	BURST.first_average_seconds = 0;
	B= BURST;
	BURST = guido(B,'Related files and output time base.');

	if isempty(BURST)
	   disp ('File and Time Base Menu must be completed. Execution ends.')
	   return
   end
end
   
%Create time base for output file.
firstime =  nc{'time'}.base_date(:);
base = datenum(firstime(1),firstime(2),firstime(3),...
   firstime(4),firstime(5),firstime(6));
firstime(5) = BURST.first_average_minutes;
firstime(6) = BURST.first_average_seconds;
ntime = datenum(firstime(1),firstime(2),firstime(3),...
   firstime(4),firstime(5),firstime(6));
btime = nc{'time'}(:);
% Create a vector of indicating which average time step
% each of the bursts should be averaged into.
speravg = BURST.average_interval;
avgnum = round((btime+(base-ntime)*86400)./speravg);
% Find decimal day numbers for time data.
btime = base + btime./86400;
% Cycle through averaging intervals and average time.
num_changes=diff(avgnum);
num_changes(length(num_changes)+1)=99;
lastbst=find(num_changes);
n=length(lastbst);
first = 1;
tavg = NaN * lastbst;
navgd = NaN * lastbst;
for iavg = 1:n;
   last = lastbst(iavg);
   tavg(iavg) = mean(btime(first:last));
   navgd(iavg) = 1 + last - first;
	first = last +1;
end
nvarout = nvarout+1;
epname{nvarout} = 'NAV_105';
inname{nvarout} = 'time';
vardesc = [ vardesc ':NAV'];
datname{nvarout} = 'navgd';
% Call function to make EPIC time.
%tavg = tavg';
alltime = EP_Time(tavg);
time = alltime(:,1);
time2 = alltime(:,2);
stime = ep_datenum(alltime(1,:));
[n m] = size(alltime);
ltime = ep_datenum(alltime(n,:));
   
   
% Find and calibrate pressure.
for i = 1:length(chosen)
   if strncmp('pres',chosen{i},4);
      disp('calibrating pressure');
      ivar = nc{chosen{i}};
		% Input data and calibration constants for pressure.
		rawdat = ivar(:);
		clk = ivar.clock_hertz(1);
		cnt = ivar.period_count(1);
		prcode = ivar.EPIC_code(1);
		% Calculate pressure period.
      prper = (rawdat .* 1000000.) ./ (2. * clk .* cnt);
      if prcode == 4027
         % Temperature is not needed for this calibration.
 			pra = ivar.cal_A(1);
			prb = ivar.cal_B(1);
			prt0 = ivar.cal_T0(1);
         % Call a function to calibrate pressure.
         prmbar = p4027 (prper,pra,prb,prt0);
      else
	      % Load temperature data or obtain a constant value.
 			if running(batch)
	   	   str=get(batch);
   			eval(['UserValueT.use_file_temperature{2} = ' str ';']);      
		      str=get(batch);
   			eval(['UserValueT.temperature_file = ' str ';']);      
	   	   str=get(batch);
   			eval(['UserValueT.temperature_variable = ' str ';']);      
		      str=get(batch);
   			eval(['UserValueT.constant_temperature = ' str ';']);      
			else
   	   	UserValueT.use_file_temperature={'checkbox',1};
				UserValueT.temperature_file = [infile(1:5) 'tcp-a.nc'];
		      UserValueT.temperature_variable = 'T_20';
   		   UserValueT.constant_temperature = 10;
      		c = UserValueT;
	         UserValueT = guido(c,'Use file temperature, or constant value?');
   	   end
      	if UserValueT.use_file_temperature{2};
         	temc = netcdf(UserValueT.temperature_file);
            tvar = temc{UserValueT.temperature_variable};
            if isempty(tvar)
               disp ('input temperature variable is empty')
               ncclose
               return
            end
            
   	      tcel = timterp(tvar,ivar);
      	   close(temc)
         	calcomment = [calcomment 'Temperature variable '...
            	   UserValueT.temperature_variable ' from file '...
               	UserValueT.temperature_file ' used for pressure. '];
	      else
   	      tcel = UserValueT.constant_temperature;
      		calcomment = [calcomment 'Constant temperature of '...
         	      num2str(tcel) ' Celsius used for pressure. '];
	      end
			% Calibration procedure depends on EPIC code.
			if prcode == 4028;
   			prC1 = ivar.cal_C1(1);
	   		prC2 = ivar.cal_C2(1);
   			prC3 = ivar.cal_C3(1);
				prD1 = ivar.cal_D1(1);
				prD2 = ivar.cal_D2(1);
				prT1 = ivar.cal_T1(1);
				prT2 = ivar.cal_T2(1);
				prT3 = ivar.cal_T3(1);
				prT4 = ivar.cal_T4(1);
				prT5 = ivar.cal_T5(1);
				prCarray = [prC3 prC2 prC1];
				prDarray = [prD2 prD1];
				prTarray = [prT5 prT4 prT3 prT2 prT1];
				% Call a function to calibrate pressure.
				prmbar = p4028 (prper,tcel,prCarray,prDarray,prTarray);
			elseif prcode == 4026;
   			prah = ivar.cal_AH(1);
	   		prac = ivar.cal_AC(1);
   			prbh = ivar.cal_BH(1);
   			prbc = ivar.cal_BC(1);
   			prt0h = ivar.cal_T0H(1);
	   		prt0c = ivar.cal_T0C(1);
   			prtr = ivar.cal_TR(1);
				% Call a function to calibrate pressure.
   			prmbar = p4026 (prper,tcel,prah,prac,prbh,prbc,prt0h,prt0c,prtr);
         end
      end
      
 		% Cycle through averaging intervals,
		%		Making pressure statistics.
		n=length(lastbst);
		first = 1;
      pavg = NaN * lastbst;
      psdev = NaN * lastbst;
		for iavg = 1:n
         last = lastbst(iavg);
         disp (['averaging pressures ' num2str(first) ' thru ' num2str(last)]);
         pavg(iavg) = mean(prmbar(first:last));
         if last>first
            psdev(iavg) = std(prmbar(first:last));
         end
			first = last +1;
		end
      nvarout = nvarout+1;
      epname{nvarout} = 'P_4023';
      inname{nvarout} = name(ivar);
      vardesc = [ vardesc ':P'];
      datname{nvarout} = 'pavg';
		if prcode == 4028
	   	calatts.P_4023 = {'clock_hertz','period_count',...
       	    'cal_C1','cal_C2','cal_C3','cal_D1',...
        	    'cal_D2','cal_T1','cal_T2','cal_T3',...
        	    'cal_T3','cal_T4','cal_T5'};
		elseif prcode == 4026
	   	calatts.P_4023 = {'clock_hertz','period_count',...
       	    'cal_AH','cal_AC','cal_BH','cal_BC',...
              'cal_T0H','cal_T0C','cal_TR'};
		elseif prcode == 4027
	   	calatts.P_4023 = {'clock_hertz','period_count',...
       	    'cal_A','cal_B','cal_T0'};
      end
      nvarout = nvarout+1;
      epname{nvarout} = 'SDP_850';
      inname{nvarout} = name(ivar);
      vardesc = [ vardesc ':SDP'];
      datname{nvarout} = 'psdev';
      calatts.SDP_850 = calatts.P_4023;
   end
end

% Find and calibrate MAVS.
for i = 1:length(chosen)
   if strncmp('bass',chosen{i},4)
      disp ('calibrating MAVS');
    	bass1 = nc{'bass_1'}(:);
      bass2 = nc{'bass_2'}(:);
		mavsnum = nc{'bass_1'}.instrument_number(:);
		mavszeros = [0 0 0 0];
		ivar = nc{'bass_1'};
		valid_range = [-100000 100000];
		allgood = (min(bass1,[],3) >= valid_range(1)) & (max(bass1,[],3) <= valid_range(2));
      allgood = allgood & (min(bass2,[],3) >= valid_range(1)) & (max(bass2,[],3) <= valid_range(2));
		% Choose MAVS calibration routine zeros.
		MAVS.axes = {{'All','A omitted','B omitted','C omitted',...
      		         'D omitted'},1};
		MAVS.zeros={{'pre-deployment','post-deployment','provided'},2};
		MAVS.zero_A = 0.;
		MAVS.zero_B = 0.;
		MAVS.zero_C = 0.;
		MAVS.zero_D = 0.;
		if running(batch)
		   str = get(batch);
		   eval(['MAVS.axes{2} = ' str ';']);      
		   str = get(batch);
		   eval(['MAVS.zeros{2} = ' str ';']);      
		   str = get(batch);
		   eval(['MAVS.zero_A = ' str ';']);      
		   str = get(batch);
		   eval(['MAVS.zero_B = ' str ';']);      
		   str = get(batch);
		   eval(['MAVS.zero_C = ' str ';']);      
		   str = get(batch);
		   eval(['MAVS.zero_D = ' str ';']);      
		else
			B = MAVS;
		   MAVS = guido(B,'Select MAVS calibration parameters');
			if isempty(MAVS)
		      disp ('MAVS Parameter Menu must be completed. Execution ends.')
		      close (nc)
		 	  	return
		   end
		end
            
		switch MAVS.zeros{2}
		case 1
			mavszeros(1) = ivar.pre_zero_offset_A(1);
			mavszeros(2) = ivar.pre_zero_offset_B(1);
			mavszeros(3) = ivar.pre_zero_offset_C(1);
			mavszeros(4) = ivar.pre_zero_offset_D(1);
		   calcomment = [calcomment 'Pre zeroes used. ']%;
		case 2
			mavszeros(1) = ivar.post_zero_offset_A(1);
			mavszeros(2) = ivar.post_zero_offset_B(1);
			mavszeros(3) = ivar.post_zero_offset_C(1);
			mavszeros(4) = ivar.post_zero_offset_D(1);
		   calcomment = [calcomment 'Post zeroes used. ']%;
		otherwise
			mavszeros(1) = MAVS.zero_A;
			mavszeros(2) = MAVS.zero_B;
			mavszeros(3) = MAVS.zero_C;
			mavszeros(4) = MAVS.zero_D;
		   calcomment = [calcomment 'User-provided zeroes. ']%;
		end

		switch MAVS.axes{2}
      case 1
         disp ('Call mavs2earth.m to calibrate MAVS data.')
         place_holder = [1];
		   [Ve, Vn, Vup, H, P, R, mavszeros] = mavs2earth(bass1, bass2, place_holder,...
		   	place_holder, place_holder, mavsnum, mavszeros, 1);
		otherwise
		   disp ('MAVS calibration with less than 4 axes not implemented.')
		   close (nc)
		   return
		end

		% Put data in normal units.
		Ve = Ve/10;
		Vn = Vn/10;
		Vup = Vup/10;
      
 		% Cycle through averaging intervals
      disp ('Computing Reynolds Stresses.')
		bstcounter = [1:length(Ve)];
		bstcounter = bstcounter';
		n=length(lastbst);
		first = 1;
      uvariance = NaN * lastbst;
      vvariance = NaN * lastbst;
      wvariance = NaN * lastbst;
      uvcovar = NaN * lastbst;
      uwcovar = NaN * lastbst;
      vwcovar = NaN * lastbst;
		for iavg = 1:n
   		last = lastbst(iavg);
         if last>first
		      disp(['averaging velocity records ' num2str(first) ' through ' num2str(last)]); 
	   		this_avg = allgood & (bstcounter >= first) &  (bstcounter <= last);
  				u_period(iavg) = Zero_x(Ve(this_avg),speravg);
	   		v_period(iavg) = Zero_x(Vn(this_avg),speravg); 
      		Allvars = cov(Ve(this_avg),Vn(this_avg));
      		if (size(Allvars) == [2 2])
       		  uvariance(iavg) = Allvars(1,1);
       		  uvcovar(iavg) = Allvars(2,1);
       		  vvariance(iavg) = Allvars(2,2);
     			end
	  			Allvars = cov(Ve(this_avg),Vup(this_avg));
      		if (size(Allvars) == [2 2])
		   		uwcovar(iavg) = Allvars(2,1);
		   		wvariance(iavg) = Allvars(2,2);
      		end
	   		Allvars = cov(Vn(this_avg),Vup(this_avg));
      		if (size(Allvars) == [2 2])
		   		vwcovar(iavg) = Allvars(1,2);
      		end
         end
			first = last +1;
		end
      ebassname = {'UVAR_4050','UVCOV_4051','VVAR_4052',...
         'UWCOV_4053','VWCOV_4054','WVAR_4055'};
   	rbassname = {'uvariance','uvcovar','vvariance',...
         'uwcovar','vwcovar','wvariance'};
      for inc = 1:6;
         ivarout = inc + nvarout;
         inname{ivarout} = name(ivar);
         epname{ivarout} = ebassname{inc};
         datname{ivarout} = rbassname{inc};
      end
      vardesc = [ vardesc ':UVAR:UVCOV:VVAR:UWCOV:VWCOV:WVAR'];
      nvarout = nvarout+6;
   end
end

% 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.MAVS_pod_depth = ' str ';']);      
else
	epicat.experiment = 'New York Bight';
	epicat.project = 'WHFC';
	epicat.comment = 'MAVS MIDAS, New York Bight Site C';
	epicat.MAVS_pod_depth = 69.1;
	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 =['Calculated using mavs_bst.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);
try
   attstr = calatts.glob;
   for icalatt = 1:length(attstr);
      str = attstr{icalatt};
      eval(['outc.'str '= bassvar.'str '(1);'])
   end
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
end
close (eps)

for i = 1:nvarout
	ovar = outc{epname{i}};
   ivar = nc{inname{i}};
	try;ovar.sensor_type = ivar.instrument_type(:);end;
	ovar.sensor_depth = ncfloat(epicat.MAVS_pod_depth);
	try;ovar.serial_number = ivar.instrument_number(:);end;
	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

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.MAVS_pod_depth;
for i = 1:nvarout
   ovar = outc{epname{i}};
   eval (['ovar(1:n) = ' datname{i} ';' ]);
end
close (nc);
close (outc);