function epname = midas_bst(infile, outfile)
% function epname = midas_bst(infile, outfile)
% calibrates BASS 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, zero crossing period, pressure standard deviation.
% Fran Hotchkiss 7 Mar 2000.
% P_4027 added 13 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
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 BASS.
for i = 1:length(chosen)
   if strncmp('bass',chosen{i},4)
      disp ('calibrating BASS');
    	BASS.axes = {{'All','A omitted','B omitted','C omitted',...
               'D omitted'},1};
      BASS.zeroes={{'pre-deployment','post-deployment'},2};
		BASS.reference_file = [infile(1:5) 'v-a_d1.nc'];
      BASS.use_file_compass={'checkbox',1};
      BASS.compass_variable = 'comp_1406';
      BASS.constant_compass = 0.;
      BASS.use_file_x_tilt={'checkbox',1};
      BASS.x_tilt_variable = 'tiltx_4017';
      BASS.constant_x_tilt_deg = 0.;
      BASS.use_file_y_tilt={'checkbox',1};
      BASS.y_tilt_variable = 'tilty_4018';
      BASS.constant_y_tilt_deg = 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.reference_file = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.use_file_compass{2} = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.compass_variable = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.constant_compass = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.use_file_x_tilt{2} = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.x_tilt_variable = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.constant_x_tilt_deg = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.use_file_y_tilt{2} = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.y_tilt_variable = ' str ';']);      
	      str=get(batch);
   		eval(['BASS.constant_y_tilt_deg = ' 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);
		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. '];
      	calatts.glob = {'pre_zero_offset_A',...
               'pre_zero_offset_B','pre_zero_offset_C',...
               'pre_zero_offset_D'};
         bassvar = ivar;
		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. '];
      	calatts.glob = {'post_zero_offset_A',...
               'post_zero_offset_B','post_zero_offset_C',...
               'post_zero_offset_D'};
         bassvar = ivar;
      end
		nf = ivar.norm_factor(1);
		pma = ivar.pod_misalignment_angle(1);
		if isempty(pma) ;
		   pma = 0.;
		end
		valid_range = ivar.valid_range(:);
		allgood = (min(a,b) >= valid_range(1)) & (max(a,b) <= valid_range(2));
		allgood = allgood & (min(c,d) >= valid_range(1)) & (max(c,d) <= valid_range(2));
		% 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);
      if BASS.use_file_compass{2}
         temc = netcdf(BASS.reference_file);
         tvar = temc{BASS.compass_variable};
         cdeg = timterp(tvar,ivar);
         close(temc)
         calcomment = [calcomment 'Compass variable '...
               BASS.compass_variable ' from file '...
               BASS.reference_file ' used for BASS. '];
      else
         cdeg = BASS.constant_compass;
      	calcomment = [calcomment 'Constant compass of '...
               num2str(cdeg) ' degrees used for BASS. '];
      end
      if BASS.use_file_x_tilt{2};
         temc = netcdf(BASS.reference_file);
         tvar = temc{BASS.x_tilt_variable};
         txdeg = timterp(tvar,ivar);
         close(temc)
         calcomment = [calcomment 'X-tilt variable '...
               BASS.x_tilt_variable ' from file '...
               BASS.reference_file ' used for BASS. '];
      else
         txdeg = BASS.constant_x_tilt_deg;
      	calcomment = [calcomment 'Constant x-tilt of '...
               num2str(txdeg) ' degrees used for BASS. '];
      end
      if BASS.use_file_y_tilt{2}
         temc = netcdf(BASS.reference_file);
         tvar = temc{BASS.y_tilt_variable};
         tydeg = timterp(tvar,ivar);
         close(temc)
         calcomment = [calcomment 'Y-tilt variable '...
               BASS.y_tilt_variable ' from file '...
               BASS.reference_file ' used for BASS. '];
      else
         tydeg = BASS.constant_y_tilt_deg;
      	calcomment = [calcomment 'Constant y-tilt of '...
               num2str(tydeg) ' degrees used for BASS. '];
      end
		txrad = txdeg .* (pi/180.);
		tyrad = tydeg .* (pi/180.);
		crad = cdeg .* (pi/180.);
      [uearth,vearth,wearth] = pod2earth (ulog,vlog,wpod,txrad,tyrad,crad);
      
 		% Cycle through averaging intervals,
      %		Making zero crossing periods and Reynolds stresses.
		bstcounter = [1:length(a)];
		bstcounter = bstcounter';
		n=length(lastbst);
		first = 1;
      u_period = NaN * lastbst;
      v_period = NaN * lastbst;
      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(uearth(this_avg),speravg);
	   		v_period(iavg) = Zero_x(vearth(this_avg),speravg); 
      		Allvars = cov(uearth(this_avg),vearth(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(uearth(this_avg),wearth(this_avg));
      		if (size(Allvars) == [2 2])
		   		uwcovar(iavg) = Allvars(2,1);
		   		wvariance(iavg) = Allvars(2,2);
      		end
	   		Allvars = cov(vearth(this_avg),wearth(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',...
         'peru_4056','perv_4057'};
   	rbassname = {'uvariance','uvcovar','vvariance',...
         'uwcovar','vwcovar','wvariance',...
      	'u_period','v_period'};
      for inc = 1:8;
         ivarout = inc + nvarout;
         inname{ivarout} = name(ivar);
         epname{ivarout} = ebassname{inc};
         datname{ivarout} = rbassname{inc};
      end
      vardesc = [ vardesc ':UVAR:UVCOV:VVAR:UWCOV:VWCOV:WVAR:peru:perv'];
      nvarout = nvarout+8;
   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.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 =['Calculated using midas_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.BASS_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.BASS_pod_depth;
for i = 1:nvarout
   ovar = outc{epname{i}};
   eval (['ovar(1:n) = ' datname{i} ';' ]);
end
close (nc);
close (outc);