%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M-file to process buoy m"x" metsys data.
% Simple bounds checks have already been
% applied to the data which have been output
% from the perl script.  This step is
% to concattenate the data and to create a
% much large files and do more sophisticated (sp?)
% processing of the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define an array
mond=[0,31,59,90,120,151,181,212,243,273,305,334];
monn=[31,28,31,30,31,30,31,31,30,31,30,31];
% minimum standard deviations for AT, RH, Ap, W/s
zsig=[0.1 0.05 0.25 1];
monc=['jan',
      'feb',
      'mar',
      'apr',
      'may',
      'jun',
      'jul',
      'aug',
      'sep',
      'oct',
      'nov',
      'dec'];
dtt=10.0/60.0/24.0;
% open file with mooring number which we will be working on
mid=input('Enter the mooring number to work on: ');
month=input('Enter the month of data to work on: ');
year=input('Enter the year of data to work on: ');
fin=input('Enter the data filename to append: ','s');
% create filename for the current month
pathm1=['/data/span9/bahr/m',num2str(mid),'/metsys/',num2str(year),'/'];
patho=pathm1;
fnm1=['m',num2str(mid),'_met_',monc(month,:),'_',num2str(year)];
% see if the file exists
fid=fopen([pathm1,fnm1,'.mat'],'r+');
% if the file exists then fid will be greater than 2.  If the file
% doesn't exist we want to create it but we also want to load
% the previous months file.  We also have to careful around the year
% change.
%# Add change for a flag array for each variable.
if fid > 2,
	fclose(fid);
	eval(['load ',pathm1,fnm1,'.mat']);
	eval(['met=metn;']);
	eval(['flag=metf;']);
else
	met=[];
	flag=[];
end
ff1=[fin(1,:)];
ff2=[fin(1,[1:length(ff1)-4])];
fod1=fopen([pathm1,ff1],'r');
if fod1 > 2,
	fclose(fod1);
	eval(['load ',pathm1,ff1]);
	eval(['m1=',ff2,';']);
else,
	m1=[];
end
% we now have the data.  Now we need to check times and sort times?
% format of met files loaded is time, airt, rh, press, wind sp, wind dir
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only go through the trouble if we have some data in m1.  At this point
% we can compute it's size.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(m1) ==0,
	dt=diff(m1(:,1));
%check for neg values to see if we have times out of order.
	[ineg]=find(dt < 0);
% we want to do all the time check stuff as vetor/matrix checks and
% try to avoid costly loops and redoing stuff that we don't need to
% do more than once.
	if isempty(ineg)==0,
% we need to do a sort of the file
		[newt1 idex1]=sort(m1(:,1));
		m1=[newt1 m1(idex1,2:size(m1,2))];
	end
% again below we have to worry about the year change.  This could
% be a recipe for disaster.  But try to get something that runs now
% and patch up the year bug later.
% Data is now sorted.  Now we need to find which values are in the
% cleaned data file already and remove get send them to the nowhere.
% we only want to process what we have to.  If I can think up a vectorized
% method for doing this, then this should be done before the sort to make
% the sort faster (fewer points sort).
% compute current year day time
% find the minimum time in the currently captured files
	tmin=min(m1(:,1));
	if size(met,1) > 0,
		tmax=max(met(:,1));
		cmin=min(met(:,1));
	else,
		cmin=370;
		tmax=0;
	end
% Note: Only m1 is beyond tmax of file.
% Now a hard check to do.  Look for equal times and remove them
% from the m1 file because we have already dealt with them.
	lm1=size(m1(:,1),1);
	im=zeros(lm1,1);
	if size(met,1) > 0,
		for i=1:lm1,
			if m1(i,1) > tmax,
				i=lm1;
			end
			id=find(met(:,1)==m1(i,1));
			if isempty(id)==0,
				im(i)=1;
			end
		end
	end
% Remove points we already have in the metsys files.
	m1(logical(im),:)=[];
%####################################################
% Now for the new data we do a simple bulk bounds check
% and set the flag array accordingly.
%####################################################
	mfm=zeros(size(m1));
%# Bounds for air temperature.
	lt=(m1(:,2) > 0)&(m1(:,2) < 40);
	mfm(~lt,2)=1;
%# Bounds for Relative Humidity.
	lh=(m1(:,3) > 20)&(m1(:,3) <= 100);
	mfm(~lh,3)=1;
%# Bound for Air Pressure.
	lp=(m1(:,4) > 950)&(m1(:,4) < 1040);
	mfm(~lp,4)=1;
%# Wind Speed Check.  If the wind speed is bad then we
%# must also assume that the wind direction is bad.
	lw=(m1(:,5) >=0)&(m1(:,5) < 100);
	mfm(~lw,[5 6])=1;
% Now we need to stuff our values in.  We also want to keep a "binary"
% mask of where the values are put in for computing our filters later.
% compare against current day -1 and see how much of old file
% we need to retain before stuffing in our values in a slow loop.
% NOTE: ABOVE needs to be modified for the transition that occurs at
% the new year.
% now we have paired down the file somewhat.  The next step is to
% actually loop and insert the data.  We set up some indicies first
% before actually doing the loop.
%# Now we have to worry about keeping track of the flags here.
	i=0; j=0;
	lmet=size(met,1);
	bmask=zeros(lmet,1);
	lnew=size(m1,1);
	k=1;
	if(cmin==370),
		met=m1;
		flag=mfm(:,[2:8]);
		bmask=ones(lnew,1);
		k=lnew+1;
	else,
		if(m1(k,1) < cmin),
			while(m1(k,1) < cmin),
				k=k+1;
			end
			met=[m1([1:k-1],:); met];
			flag=[mfm([1:k-1],[2:8]); flag];
			bmask=[ones(length([1:k-1]),1)
			       bmask];
		end
	end
	while(k <= lnew),
		lmet=size(met,1);
		if lmet==0,
			met=m1;
			flag=mfm(:,[2:8]);
			k=lnew+1;
			bmask=[bmask
		       	       ones(lnew,1)];
		else,
			[index]=find(met(:,1) < m1(k,1));
			idmz=max(index);
			if idmz==lmet,
				met=[met
                     	             m1([k:lnew],:)];
				flag=[flag
				     mfm([k:lnew],[2:8])];
				bmask=[bmask
			               ones(length(k:lnew),1)];
					k=lnew;
			else
				met=[met([1:idmz],:)
		     	     	     m1(k,:)
		     	     	     met([idmz+1:lmet],:)];
				flag=[flag([1:idmz],:)
				      mfm(k,[2:8]);
				      flag([idmz+1:lmet],:)];
				bmask=[bmask([1:idmz])
				       1
			               bmask([idmz+1:lmet])];
			end
		end
		k=k+1;
	end
% If the above works, bmask should be a file with 1's where
% we have inserted data and should be the same length as
% the array met with the inserted data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%# Next we apply any of the filters we want to the data and
%# set the flags to 2's to indicate the data was removed by
%# the filter and not because if failed a bounds check.
% Now we can apply our gradient filter to the data (?)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if size(met,2)==9,
		knk=[2:4 9];
	else,
		knk=[2:4];
	end
	for kk=knk,
		if kk <=4,
			ysig=zsig(kk-1);
		else,
			ysig=0;
		end
		if ((mid==3)&(kk==2)),
			cc=1.5;
		else,
			cc=2.25;
		end
		nwhat=sigfilt(met(:,kk),met(:,1),flag(:,kk-1),1,cc,9,6,ysig);
		if kk==2,
			nwhat=sigfilt(met(:,kk),met(:,1),nwhat(:),2,cc,9,6,ysig);
		end
%		nwhat=metfiltf(kk,met(:,kk),bmask,met(:,1),flag(:,kk-1));
		flag(:,kk-1)=nwhat(:);
	end
%# Find wind speeds flaged as bad and flag the corresponding wind 
%# directions as bad.
	ln=flag(:,4) > 2;
	flag(ln,5)=2;
% Need to save the appropriate data here!
% we need to be aware of year and month breaks
% in our data if we want to save it correctly.
	metn=met;
	metf=flag;
	eval(['save ',pathm1,fnm1,' metn metf']);
end % end of if m1 has data.
