%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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];
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 id to work on: ');
month=input('Enter the month (number) to work on: ');
year=input('Enter the year to work on: ');
fin=input('Enter the data file to append: ','s');
% create filename for the current month
pathm1=['/data/span9/bahr/m',num2str(mid),'/tstring/',num2str(year),'/'];
patho=pathm1;
fnm1=['m',num2str(mid),'_tstring_',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.
% set up seperate arrays for each t string sensor so we can hang onto
% time changes....
if fid > 2,
	fclose(fid);
	eval(['load ',pathm1,fnm1,'.mat']);
	for in=1:10,
		eval(['t',num2str(in),'=t',num2str(in),'n;']);
		eval(['f',num2str(in),'=f',num2str(in),'n;']);
	end
	eval(['timet=timen;']);
else
	for in=1:10,
		eval(['t',num2str(in),'=[];']);
		eval(['f',num2str(in),'=[];']);
	end
	eval(['timet=[];']);
end
% we should now have the old matlab process tstring files loaded.
% The next step is to load the new files with data to be processed.
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?
% the format for tstring files will be time, id, depth, temp, cond
% for the 10 sensors.  We have individual times for each microcat.
% do a loop for the sensors and sort each loop of data.
% only do the following if we have nonzero m1
if (isempty(m1) == 0),
% sort oasis time first
	dt=diff(m1(:,1));
	[ineg]=find(dt < 0);
	if isempty(ineg)==0,
		[newtt index]=sort(m1(:,1));
		m1=[newtt m1(index,[2:size(m1,2)])];
	end
% At this point we have sorted the whole mess by the oasis times.
% Now since at each oasis time microcats transmit latest value
% we might assume that the for microcats are sorted except for 
% bad values.
% now check times for each microcat and against the oasis time.
	for muk=1:10,
		dot=abs(m1(:,1)-m1(:,(2+(muk-1)*5)));
%# allow a period of twenty minutes to be between microcat and oasis time.
		idot=find(dot > 0.014);
		m1(idot,(2+(muk-1)*5))=m1(idot,1);
	end
	for itc=1:10,
		clear dt; clear ineg;
		ict=2+(itc-1)*5;
		dt=diff(m1(:,ict));
%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
% just do this for the time being.  If things are the way we think
% they are this should nolonger be here.....
			[newt1 idex1]=sort(m1(:,ict));
			m1(:,[ict:ict+4])=[newt1 m1(idex1,(ict+1):(ict+4))];
		end
	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
% Now we have to loop for each microcat...
% No! We only do it to remove the oasis times....
	itc=1;
	tmin=min(m1(:,1));
	lt=length(timet);
	if lt > 0,
		tmax=max(timet);
		cmin=min(timet);
	else,
		cmin=370;
		tmax=0;
	end
	lm1=size(m1(:,1),1);
	im=zeros(lm1,1);
	lt=size(timet,1);
	if lt > 0,
	        for i=1:lm1,
			if m1(i,1) > tmax,
				i=lm1;
			end
			id=find(timet==m1(i,1));
			if isempty(id)==0,
				im(i)=1;
			end
		end
	end
	m1(logical(im),:)=[];
end
%# After we change the size of m1 we need to check if it is empty for
%# continuation.
if (isempty(m1)==0),
%# Now that we have removed all the surpuflus points we go about setting up
%# the new flag arrays.  We also do a bounds check on the data.
	tft=zeros(size(m1));	
% Now set the flags to 1 for out of bounds data.;
% currently only checking temperatures
	for ib=1:10,
%# First define out of bounds data and flag that data.
		ll=((m1(:,(5+(ib-1)*5)) > 22)|(m1(:,(5+(ib-1)*5)) < 5));
		tft(ll,(5+(ib-1)*5))=1;
%# Then define data which is missing (Temperature).
		ll=(m1(:,(5+(ib-1)*5)) == -99);
		tft(ll,(5+(ib-1)*5))=3;
%# Then define data which is missing (conductivity).
		ll=((m1(:,(6+(ib-1)*5)) > 6)|(m1(:,(6+(ib-1)*5)) < 3));
		tft(ll,(6+(ib-1)*5))=1;
		ll=(m1(:,(6+(ib-1)*5)) == -99);
		tft(ll,(6+(ib-1)*5))=3;
%# Second define out of bounds data.
	end
		
% create an array without points we already have.....
% We need to worry that for times we may have blanks in the original but we want
% to replace them with new data (Do we store *.mat files without blanks?).
% we can't do the next step for the t-string we must also worry if the
% above eval's will work.
% 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.
	i=0; j=0;
% all the t'x' array will be the same size so we only need one to check...
	lt1=size(t1,1);
	bmask=zeros(lt1,1);
	lnew=size(m1,1);
	k=1;
	if(cmin==370),
		timet=m1(:,1);
		bmask=ones(lnew,1);
		k=lnew+1;
		for mm=1:10,
			st=2+(mm-1)*5;
			sp=6+(mm-1)*5;
			eval(['t',num2str(mm),'=[m1(:,[st:sp])',...
			      ' -99.*ones(size(m1,1),1)];']);
			eval(['f',num2str(mm),'=[tft(:,[st+1:sp])',...
			      ' tft(:,st+3)];']);
		end
	else,
		if(m1(k,1) < cmin),
			while(m1(k,1) < cmin),
				k=k+1;
			end
			timet=[m1([1:k-1],1); timet];
			for mm=1:10,
				st=2+(mm-1)*5;
				sp=6+(mm-1)*5;
			      eval(['t',num2str(mm),'=[m1([1:k-1],[st:sp]) ',...
				    ' -99.*ones(length([1:k-1]),1);',...
				    ' t',num2str(mm),'];']);
			     eval(['f',num2str(mm),'=[tft([1:k-1],[st+1:sp])',...
				   ' tft([1:k-1],[st+3]);',...
				   ' f',num2str(mm),'];']);
				bmask=[ones(length([1:k-1]),1)
				       bmask];
			end
		end
	end			
	while(k <= lnew),
		lt1=size(t1,1);
		if lt1==0,
			timet=m1(:,1);
			k=lnew+1;
			for mm=1:10,
				st=2+(mm-1)*5;
				sp=6+(mm-1)*5;
				eval(['t',num2str(mm),'=[m1(:,[st:sp])',...
					' -99.*ones(size(m1,1),1)];']);
%# Add in flags.  Note: set salt flag to temp flag for new data at this
%# point.  Also, no flag for time.
				eval(['f',num2str(mm),'=[tft(:,[st+1:sp])',...
				      ' tft(:,st+3)];']);
			end
			bmask=[bmask
		               ones(lnew,1)];
		else,
			[index]=find(timet(:,1) < m1(k,1));
			idmz=max(index);
			if idmz==lt1,
				timet=[timet
                                       m1([k:lnew],1)];
				for mm=1:10,
					st=2+(mm-1)*5;
					sp=6+(mm-1)*5;
				eval(['t',num2str(mm),'=[t',num2str(mm),';',...
 				      ' m1([k:lnew],[st:sp]) ',... 
				      ' -99.*ones(length(k:lnew),1)];']);
				eval(['f',num2str(mm),'=[f',num2str(mm),';',...
				      ' tft([k:lnew],[st+1:sp]) ',...
				      ' tft([k:lnew],[st+3])];']); 
				end
				bmask=[bmask
			       	       ones(length(k:lnew),1)];
				k=lnew;
			else,
				timet=[timet([1:idmz])
				      m1(k,1)
				      timet([idmz+1:lt1])];
				for mm=1:10,
					st=2+(mm-1)*5;
					sp=6+(mm-1)*5;
			eval(['t',num2str(mm),'=[t',num2str(mm),...
			      '([1:idmz],:); m1(k,[st:sp]) -99;',... 
			      ' t',num2str(mm),'([idmz+1:lt1],:)];']);
			eval(['f',num2str(mm),'=[f',num2str(mm),...
			      '([1:idmz],:); tft(k,[st+1:sp]) tft(k,st+3);',...
			      ' f',num2str(mm),'([idmz+1:lt1],:)];']);
				end
				bmask=[bmask([1:idmz])
				       1
			               bmask([idmz+1:lt1])];
		 	end
		end
		k=k+1;
	end
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We have a much bigger problem...How to coordinate all of the
% data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we can apply our gradient filter to the data (?)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(m1)==0,
	for i=1:10,
		eval(['ck=t',num2str(i),';']);
		eval(['cf=f',num2str(i),';']);
		nwhat=ctdfilt2(5,ck(:,5),bmask,ck(:,1),cf(:,4));
	end
end
% get the indecies of the insertion points
%for kk=3:4,
%	nwhat=ctdfilt2(kk,ctd(:,kk),bmask,ctd(:,1));
%	ctd(:,kk)=nwhat';
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Need to add a calculation of salinity from the conductivity 
% here after we filter the conductivity.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==>ADD SALT CALC HERE<==
c35=4.2914;	% S/m
% Why compute salt here?  Why not compute it earlier and then put the
% values in?  Want to filter before we compute salt?
%##########################################################
%# Only do the following if there is something to plot...
%# So, we need to check the size of t1 all arrays are the
%# defined to be the same size even if they are empty or full
%# of missing data.
%##########################################################
if(isempty(t1)==0),
for mm=1:10,
	clear crat; clear ccf;
	eval(['crat=t',num2str(mm),'(:,5)./c35;']);
%# grab the flag for salts (same as temp at this point)
	eval(['ccf=f',num2str(mm),'(:,4);']);
%# set flags and watch for computing bounds
	ll=(ccf == 0); %# bad values are flagged.
	if isempty(crat(ll))==0,	
	eval(['t',num2str(mm),'(ll,6)=sw_salt(crat(ll),t',num2str(mm),...
              '(ll,4),t',num2str(mm),'(ll,3));']);
%# Now we want to check the values and change bound flags accordingly.
       eval(['ls=((t',num2str(mm),'(:,6) > 32)&(t',num2str(mm),'(:,6) < 35));']);
	eval(['f',num2str(mm),'(~ls,5)=1;']);
%# Check gradient salt flags
	clear kko; clear kkf; clear jjf;
	eval(['kko=t',num2str(mm),';']);
	eval(['jjf=f',num2str(mm),';']);
	kkf=ctdfilt2(6,kko(:,6),ones(size(kko,1),1),kko(:,1),jjf(:,5));
	eval(['f',num2str(mm),'(:,5)=kkf(:);']);
%# Flag ones where we have no data as missing not out of bounds.
	eval(['ls=((f',num2str(mm),'(:,3)==3)|(f',num2str(mm),'(:,4)==3));']);
	eval(['f',num2str(mm),'(ls,5)=3;']);
	else,
	  eval(['f',num2str(mm),'(:,5)=1;']);
	end %# end of if isempty(crat)
end	%# end of for loop
end 	%# end of if isempty(t1)
% 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.
if isempty(t1)==0,
	for mm=1:10,
		eval(['t',num2str(mm),'n=t',num2str(mm),';']);
		eval(['f',num2str(mm),'n=f',num2str(mm),';']);
	end
	timen=timet;
	eval(['save ',pathm1,fnm1,' timen t1n t2n t3n t4n t5n ',...
              't6n t7n t8n t9n t10n f1n f2n f3n f4n f5n f6n f7n f8n ',...
	      'f9n f10n;']);
end %# end of if isempty(t1)
