%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
fpath=['/data/span9/bahr/perlprogs/'];
mond=[0,31,59,90,120,151,181,212,243,273,304,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
eval(['load ',fpath,'mooringid.dat']);
eval(['mid=mooringid;']);
% grab time
eval(['load ',fpath,'time.dat']);
eval(['gtime=time;']);
year=gtime(1);
month=gtime(2);
day=gtime(3);
hour=gtime(4);
minu=gtime(5);
sec=0;
yday=gtime(6);
yearo=year;
% create filename for the current month
pathm1=['/data/span9/bahr/m',num2str(mid),'/adcp/',num2str(year),'/'];
patho=pathm1;
fnm1=['m',num2str(mid),'_gps_',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.
%# Note change to add a flag array for the data.  This new array will
%# contains flags defining how good the data used is for both plotting
%# and distribution.  The flag array will contain 1 less element than
%# the main array.  We will not flag the time, but pressure, temperature
%# conductivity and salinity will be flaged.
if fid > 2,
	fclose(fid);
	eval(['load ',pathm1,fnm1,'.mat']);
	eval(['b=gpsn;']);
	eval(['fb=gpsf;']);
else
	b=[];
	fb=[];
end
% Check to see if we want to load the old file.  We assume that data
% from the oasis mooring will not be older than 7 days.
if day > 8,
	a=[];
	fa=[];
else,
	if month==1,
		yearo=year-1;
		patho=['/data/span9/bahr/m',num2str(mid),'/adcp/',num2str(yearo),'/'];
		fnm1o=['m',num2str(mid),'_gps_dec_',num2str(yearo)];
	else,
		yearo=year;
		patho=pathm1;
		fnm1o=['m',num2str(mid),'_gps_',monc(month-1,:),'_',num2str(year)];
	end
% check to see if old file is there.  If it is load it.
	fido=fopen([patho,fnm1o,'.mat'],'r+');	
	if fido > 2,
		fclose(fido);
		eval(['load ',patho,fnm1o,'.mat']);
		eval(['a=gpsn;']);
		eval(['fa=gpsf;']);
	else
		a=[];
		fa=[];
	end
end
% create one array of data.  We have to be careful doing this at the
% year change.  Since the year is not explicitely stated in the file.
% Do a trick to handle year changes.  If we have a year change yearo!=year
% then we subtract the max days from the previous year?
%# Note we use 366 as December 31 at 12pm is denoted 365.5 days.
maxd=366;
if yearo~=year,
	maxd=366;
	if rem(yearo,4)==0,
		maxd=maxd+1;
	end
	if isempty(a)==0,
		a(:,1)=a(:,1)-maxd;
	end
end
gps=[a
     b];
%# set up the array flag for the flag information on the main array.
flag=[fa
      fb];
% we should now have the old matlab process gps sys files loaded.
% The next step is to load the new files with data to be processed.
ff1=['m',num2str(mid),'_gps.dat'];
ff2=['m',num2str(mid),'_gpso.dat'];
if (month==1)&(day >=2),
	fod1=fopen([pathm1,ff1],'r');
	fod2=fopen([pathm1,ff2],'r');
else,
	fod1=fopen([pathm1,ff1],'r');
	fod2=fopen([patho,ff2],'r');
end
if fod1 > 2,
	fclose(fod1);
	eval(['load ',pathm1,'m',num2str(mid),'_gps.dat']);
	eval(['m1=m',num2str(mid),'_gps;']);
else,
	m1=[];
end
if fod2 > 2,
	fclose(fod2);
	if(month==1)&(day >=2),
		eval(['load ',pathm1,'m',num2str(mid),'_gpso.dat']);
	else,
		eval(['load ',patho,'m',num2str(mid),'_gpso.dat']);
	end
	eval(['m1o=m',num2str(mid),'_gpso;']);
else,
	m1o=[];
end
% we only have to worry about day ==1 on this part.  But that doesn't
% cover the posibility that data from before may be in one of the old
% files?
if ((day==1)&(month==1)&(size(m1o,1)>0)),
	m1o(:,1)=m1o(:,1)-maxd;
end
% we have a problem here at the year change.  We need to worry about.
% since the files only look at year day and not have the year in it.
m1=[m1o
    m1];

% we now have the data.  Now we need to check times and sort times?
% format of gps files is time, press, temp, cond, salt
ydnow=mond(month)+day+hour/24+minu/(24*60);
if (rem(year,4)==0)&(month > 2),
	ydnow=ydnow+1;
end
ydnowl=ydnow-1;
if size(m1,1) > 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
	dt=diff(m1(:,1));
	idi=find(dt==0);
	m1(idi,:)=[];
% 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(gps,1) > 0,
		cmin=min(gps(:,1));
		tmax=max(gps(:,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(gps,1) > 0,
		for i=1:lm1,
			if m1(i,1) > tmax,
				i=lm1;
			end
			id=find(gps(:,1)==m1(i,1));
			if isempty(id)==0,
				im(i)=1;
			end
		end
	end
	m1(logical(im),:)=[];
end
if(isempty(m1) ==0),
% Remove points we already have in the metsys files.
%# Define our flag array for the new data.
	cfc=zeros(size(m1));
	cfc(:,1)=[];
%# Check bounds on latitude
	if mid==1,
		llat=((m1(:,2) > 36.72)&(m1(:,2) < 36.76));
		llon=((m1(:,3) > -122.035)&(m1(:,3) < -122.0));
	end
	if mid==2,
		llat=((m1(:,2) > 36.68)&(m1(:,2) < 36.71));
		llon=((m1(:,3) > -122.42)&(m1(:,3) < -122.35));
	end
%# Check bounds on pressure first
%	lp=(m1(:,2) > 0)&(m1(:,2) < 20);
%	lkp=(m1(:,2) > 0)&(m1(:,2) < 2);
%	kjunk=m1(~lkp,2);
%	if(isempty(kjunk)== 0),
%		lkj=kjunk(:)==-99;
%		kjunk=kjunk(~lkj);
%	   if (isempty(kjunk) == 0),
%		ffff=fopen('gps_message.txt','w+');
%		fprintf(ffff,'Mooring %d gps depth exceeded 1 m\n',mid);
%		fprintf(ffff,'A value of %f was recieved\n',kjunk(1,1));
%		fprintf(ffff,'Please check the values\n');
%		fprintf(ffff,'Fred\n');
%		fclose(ffff);
%		unix('mail bahr@oc.nps.navy.mil < gps_message.txt');
%	   end
%	end 
	cfc(~llat,1)=1;
	cfc(~llon,2)=1;
%# Check bounds on temperature next
%# Note: if Temperature is bad we set the salinity value to bad
%# also, as it is computed from conductivity using the temperature.
%	lt=(m1(:,3) > 5)&(m1(:,3) < 30);
%	cfc(~lt,[2 4])=1;
%# Check the bounds on the conductivity
%# If conductivity is bad set it and salinity to bad.
%	lc=(m1(:,4) > 2)&(m1(:,4) < 6);
%	cfc(~lc,[3 4])=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.
%######################################################################
%# when we add points we will need to expand the flag array accordingly
%# here.  This will be the first step.  We might also change how bmask
%# is used.  The mask, the array, and the flag array will be necessary
%# for any filter to work properly.  The filter will have to know if
%# a point is new and if the points around it can be used in the filter
%# or if it will have to expand to look for other points.
	i=0; j=0;
	lgps=size(gps,1);
	bmask=zeros(lgps,1);
	lnew=size(m1,1);
	k=1;
	if(cmin==370),
		gps=m1;
		flag=cfc;
		bmask=ones(lnew,1);
		k=lnew+1;
	else,
	   if(m1(k,1) < cmin),
		while(m1(k,1) < cmin),
			k=k+1;
		end
		gps=[m1(1:k-1,:); gps];
		flag=[cfc(1:k-1,:); flag];
		bmask=[ones(length([1:k-1]),1)
		       bmask];
	  end
	end
	while(k <= lnew),
		lgps=size(gps,1);
		if lgps==0,
			gps=m1;
			flag=cfc;
			k=lnew+1;
			bmask=[bmask
		       	       ones(lnew,1)];
		else,
			[index]=find(gps(:,1) < m1(k,1));
			idmz=max(index);
			if idmz==lgps,
				gps=[gps
                     	     	     m1([k:lnew],:)];
				flag=[flag
				      cfc([k:lnew],:)];
				bmask=[bmask
			               ones(length(k:lnew),1)];
			 	k=lnew;
			else
				gps=[gps([1:idmz],:)
		     	     	     m1(k,:)
		     	             gps([idmz+1:lgps],:)];
				flag=[flag([1:idmz],:)
				      cfc(k,:)
 				      flag([idmz+1:lgps],:)];
				bmask=[bmask([1:idmz])
				       1
			               bmask([idmz+1:lgps])];
			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.
%# We should also have the flag array the same size with 0's
%# for the new values.  We assume all data is good at first.
%# we will then remove it as we progress.
%# Now to compute bounds checks on data.
%#
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we can apply our gradient filter to the data (?)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get the indecies of the insertion points
%# Here we have a major change.  We need to have the filter return
%# the flag value and leave the variable alone.
%	for kk=3,
%		nwhat=sigfilt(gps(:,kk),gps(:,1),flag(:,kk-1),1,2.25,9,6,0.005);
%		if kk==4,
%			nwhat=sigfilt(gps(:,kk),gps(:,1),nwhat(:),2,2.25,9);
%		end
%		nwhat=gpsfiltf(kk,gps(:,kk),bmask,gps(:,1),flag(:,kk-1));
%		flag(:,kk-1)=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
%	crat=gps(:,4)./c35;
%	gps(:,5)=sw_salt(crat,gps(:,3),gps(:,2));
% Add a bounds check here to capture anything that escaped earlier.
%	ll=(gps(:,5) > 32)&(gps(:,5) < 35);
%	flag(ll,4)==1;
%	nwhat=sigfilt(gps(:,5),gps(:,1),flag(:,4),1,2.25,9,6,0.005);
%	nwhat=sigfilt(gps(:,5),gps(:,1),nwhat(:),2,2.25,9,6,0.005);
%	flag(:,4)=nwhat(:);
% 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 day > 8,
		gpsn=gps;
		gpsf=flag;
		eval(['save ',pathm1,fnm1,' gpsn gpsf']);
	else
		if month==1,
			ll=gps(:,1) < 0;
			gpsn=gps(ll,:);
			gpsn(:,1)=gpsn(:,1)+maxd;
			gpsf=flag(ll,:);
			eval(['save ',patho,fnm1o,' gpsn gpsf']);
			clear gpsn;
			clear gpsf;
			gpsn=gps(~ll,:);
			gpsf=flag(~ll,:);
			eval(['save ',pathm1,fnm1,' gpsn gpsf']);
		else
			dmax=mond(month)+1;
			if (month > 2)&(rem(year,4)==0),
				dmax=dmax+1;
			end
			ll=gps(:,1) < dmax;
			gpsn=gps(ll,:);
			gpsf=flag(ll,:);
			eval(['save ',patho,fnm1o,' gpsn gpsf']);
			clear gpsn;
			clear gpsf;
			gpsn=gps(~ll,:);
			gpsf=flag(~ll,:);
			eval(['save ',pathm1,fnm1,' gpsn gpsf']);
		end
	end
end % end of im m1 has data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now process the last days worth of data for display
% by various web pages.
% 1) compute 24 hour means and deviations along with
% current values.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(size(gps(:,1),1) > 0),
	ydp=ydnow-7;
	if ydp < 0,
		ydp=ydp-1;
	end
	ll=(gps(:,1) > ydnowl)&(gps(:,1) <= ydnow);
	lp=(gps(:,1) > ydp)&(gps(:,1) <= ydnow);
	gpsd=gps(ll,:);
	flagd=flag(ll,:);
	gpsp=gps(lp,:);
	flagp=flag(lp,:);
	lm=size(gpsd,1);
	lw=size(gpsd,2);
% now to write out a simple table of numbers for a perl
% script to load and create an htlm table from
	bad=-99;
% we need to check for the output file if it exists we remove it
	gpsplot(gpsp(:,1),gpsp(:,[2 3]),year,mid,flagp(:,[1 2]));
	dtplot(gpsp(:,1),year,mid,4);
	ffn=[fpath,'m',num2str(mid),'_nowgps.dat'];
	if(isempty(gpsd)==0),
		fid=fopen(ffn,'w');
		zt=ydtostr(year,gpsd(lm,1));
		zts=datestr(zt,0);
		fprintf(fid,'%s\n',zts(1,:));
		for jj=[3 5],
			agps=gpsd(:,jj);
			ll=flagd(:,jj-1)==0;
			agps=agps(ll);
		if isempty(agps)==0,
			fprintf(fid,'%8.4f %d %8.4f ',gpsd(lm,jj),flagd(lm,jj-1),mean(agps));
			fprintf(fid,'%8.4f %8.4f ',std(agps),min(agps));
			fprintf(fid,'%8.4f \n',max(agps));
		else
			fprintf(fid,'%8.4f %d %8.4f ',bad,flagd(lm,jj-1),bad);
			fprintf(fid,'%8.4f %8.4f ',bad,bad);
			fprintf(fid,'%8.4f \n',bad);
		end
	end		
	end
end
% We could add plots of the data here and save them to standard images.
% matlab would be good for all but stick and contour plots.
%
quit
