%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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];
% minimum standard deviations? 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;
% add arrays of nominal lat and lon of moorings for output 
ylat=[36.745
      36.666
      36.567
      36.150];
ylon=[-122.022
      -122.400
      -123.972
      -122.433];
% open file with mooring number which we will be working on
eval(['load ',fpath,'mooringid.dat']);
eval(['mid=mooringid;']);
% load the gmt time that getdata.perl has output.
% time file format: year,month,day of m, hour, min, year day
eval(['load ',fpath,'time.dat']);
eval(['gtime=time;']);
% grab 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),'/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(['b=metn;']);
	eval(['fb=metf;']);
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),'/metsys/',num2str(yearo),'/'];
		fnm1o=['m',num2str(mid),'_met_dec_',num2str(yearo)];
	else,
		yearo=year;
		patho=pathm1;
		fnm1o=['m',num2str(mid),'_met_',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=metn;']);
		eval(['fa=metf;']);
	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?
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
met=[a
     b];
flag=[fa
      fb];
% we should now have the old matlab process met sys files loaded.
% The next step is to load the new files with data to be processed.
% add error checking should either of the file below not exist...
ff1=['m',num2str(mid),'_metsys.dat'];
ff2=['m',num2str(mid),'_metsyso.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),'_metsys.dat']);
	eval(['m1=m',num2str(mid),'_metsys;']);
else,
	m1=[];
end
if fod2 > 2,
	fclose(fod2);
	if(month==1)&(day >= 2),
	eval(['load ',pathm1,'m',num2str(mid),'_metsyso.dat']);
	else,
	eval(['load ',patho,'m',num2str(mid),'_metsyso.dat']);
	end
	eval(['m1o=m',num2str(mid),'_metsyso;']);
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? Only want to do this for january first.
if ((day==1)&(month==1)&(size(m1o,1)> 0)),
	m1o(:,1)=m1o(:,1)-maxd;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We need to worry about the arrays being empty because no data have
% been read or there is no data coming from the mooring.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
% Now remove any duplicate times
	dt=diff(m1(:,1));
	idi=find(dt==0);
	m1(idi,:)=[];
% Also remove any times within 5 seconds of the others (we had this fail before)
	dt=diff(m1(:,1)).*24*60;
	idi=find(dt < 0.08);
	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(met,1) > 0,
		tmax=max(met(:,1));
	else,
		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),:)=[];
end
if isempty(m1)==0,
%####################################################
% 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;
%# Speed of Sound in Air check if we have that data...
	lwid=size(m1,2);
	if lwid == 9,
		lw=(m1(:,9) >= 200)&(m1(:,9) < 400);
		mfm(~lw,9)=1;
	end
% 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.
%#************************************************************
%# Assuming we don't change the size of the metfile in the middle
%# this should work...
%#***************************************************************
	i=0; j=0;
	lmet=size(met,1);
	lwid2=size(met,2);
	bmask=zeros(lmet,1);
	lnew=size(m1,1);
	k=1;
	while(k <= lnew),
		lmet=size(met,1);
		if lmet==0,
			met=m1;
			flag=mfm(:,[2:lwid]);
			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:lwid])];
				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:lwid]);
				      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 lwid==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.
	if day > 8,
		metn=met;
		metf=flag;
		eval(['save ',pathm1,fnm1,' metn metf']);
	else
		if month==1,
			ll=met(:,1) < 0;
			metn=met(ll,:);
			metn(:,1)=metn(:,1)+maxd;
			metf=flag(ll,:);
			eval(['save ',patho,fnm1o,' metn metf']);
			clear metn;
			clear metf;
			metn=met(~ll,:);
			metf=flag(~ll,:);
			eval(['save ',pathm1,fnm1,' metn metf']);
		else
			dmax=mond(month)+1;
			if (month > 2)&(rem(year,4)==0),
				dmax=dmax+1;
			end
			ll=met(:,1) < dmax;
			metn=met(ll,:);
			metf=flag(ll,:);
			eval(['save ',patho,fnm1o,' metn metf']);
			clear metn; clear metf;
			metn=met(~ll,:);
			metf=flag(~ll,:);
			eval(['save ',pathm1,fnm1,' metn metf']);

		end
	end
end % end of if 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(met,1) > 0),
	ydp=ydnow-7;
	if ydp < 0,
		ydp=ydp-1;
	end
	ll=(met(:,1) > ydnowl)&(met(:,1) <= ydnow);
	lp=(met(:,1) > ydp)&(met(:,1) <= ydnow);
	metd=met(ll,:);
	metp=met(lp,:);
%# Add truncate flag values for plotting and stat computation
%# We don't want to include the values we have deemed bad.
	flagd=flag(ll,:);
	flagp=flag(lp,:);
	lm=size(metd,1);
	lw=size(metd,2);
	lp=size(metp,1);
%# We have expanded the array metd to include the buoy heading and 
%# the relative wind direction so we need to subtract 2 to get the
%# width where the "true" wind direction is.
	lw=lw-2;
% compute u and v components of wind for plotting.
	wwd=metp(:,6);
	wwd=wwd+180;
	lwd=wwd(:) > 360;
	wwd(lwd)=wwd(lwd)-360;
	up=metp(:,5).*sin(wwd(:).*(pi/180));
	vp=metp(:,5).*cos(wwd(:).*(pi/180));
	nu=[];
	nv=[];
	nt=[];
	tinc=10.0/60.0/24.0;
	if(lp > 0),
	wid=fopen([fpath,'wind',num2str(mid),'.dat'],'w');
% output data for inclusion in CODAR plots....
% one value will be an average of the last hour of data.
% the other value will be an 24 hour average of the winds.
% 1) Get data for past 24 hours
	wh=hour-0.5;
	wd=yday;
	wh=ydnow-(wd+wh/24);
	lih=((met(:,1) > wh)&(met(:,1) <= ydnow));
	metwh=met(lih,:);
	flagh=flag(lih,:);
	lfg=flagh(:,4)==0;
	metwh=metwh(lfg,:);
	wiz=fopen([fpath,'wind24.dat'],'a');
	fprintf(wiz,'%% yyyy mm dd hh mid lat lon u v\n');
% 2) Now compute the hour interval for computation
	if isempty(metwh)==0,
		hb=[hour:24];
		ha=[1:hour-1];
		db=ones(length(hb),1).*wd;
		da=ones(length(ha),1).*yday+1;
		hh=[hb(:);ha(:)];
		dd=[db(:);da(:)];
		hm=0.5;
		lih=length(hh);
% 3) Now to compute the actual values
		for iu=1:lih,
			hlow=dd(iu)+(hh(iu)-hm)/24;
			hhig=dd(iu)+(hh(iu)+hm)/24;
			loh=((metwh(:,1) > hlow)&(metwh(:,1) <= hhig));
			wwa=metwh(loh,5);
			wwb=metwh(loh,6);
			if isempty(wwb)==0,
				wwb=wwb+180;
				lll=wwb > 360;
				wwb(lll)=wwb(lll)-360;
				uua=wwa.*sin(wwb.*(pi/180));
				vva=wwa.*cos(wwb.*(pi/180));
				mud=mean(uua);
				mvd=mean(vva);
			else,
				mud=-99.0;
				mvd=-99.0;
			end
			if dd(iu) < 0,
				mon=month-1;
				if mon < 0,
					you=year-1;
					mon=12;
					dou=31;
				else,
					dou=monn(mon);
					you=year;
				end
			else,
				mon=month;
				you=year;
				dou=dd(iu)-mond(mon);
			end
			fprintf(wiz,'%4.4d %2.2d %2.2d ',you,mon,dou);
			fprintf(wiz,'%2.2d %d ',hh(iu),mid);
			fprintf(wiz,'%8.3f %8.3f ',ylat(mid),ylon(mid));
			fprintf(wiz,'%6.2f %6.2f \n',mud,mvd);
		end % end of for loop
	end % end of isempty
% now to create get just unflaged data for 24 hours
%	lid=flagd(:,4)==0;
%	widd=metd(lid,5);
%	widh=metd(lid,6);
%	wihd=metd;
%	wihf=flagd;
% now to get just the last hour of data...
%	ydh=ydnow-(1/24);
%	lh=(wihd(:,1)> ydh)&(wihd(:,1) <= ydnow);
%	wihd=wihd(lh,:);
%	wihf=wihf(lh,:);
%	lh=wihf(:,4)==0;
%	wihh=wihd(lh,6);
%	wihd=wihd(lh,5);
% now to compute components
%	if isempty(widh)==0,
%		widh=widh+180;
%		ldh=widh > 360;
%		widh(ldh)=widh(ldh)-360;
%		ud=widd.*sin(widh.*(pi/180));
%		vd=widd.*cos(widh.*(pi/180));
%		mud=mean(ud);
%		mvd=mean(vd);
%	else,
%		mud=-99;
%		mvd=-99;
%	end
%	if isempty(wihd)==0,
%		wihh=wihh+180;
%		lhh=wihh > 360;
%		wihh(lhh)=wihh(lhh)-360;
%		uh=wihd.*sin(wihh.*(pi/180));
%		vh=wihd.*cos(wihh.*(pi/180));
%		muh=mean(uh);
%		mvh=mean(vh);
%	else,
%		muh=-99;
%		mvh=-99;
%	end	
%	wiz=fopen([fpath,'wind24.dat'],'a');
%	wic=fopen([fpath,'windhr.dat'],'a');
%	if wic > 2,
%	 	fprintf(wic,'%8.3f %8.3f %6.2f %6.2f\n',ylat(mid),ylon(mid),muh,mvh);
%		fclose(wic);
%	end
%	if wiz > 2,
%		fprintf(wiz,'%8.3f %8.3f %6.2f %6.2f\n',ylat(mid),ylon(mid),mud,mvd);
%		fclose(wiz);
%	end
% find gaps and fill them with -99's for plot plus
		metpp=metp;
		lppp=metpp(:,1) < 0;
		metpp(lppp,1)=metpp(lppp,1)+1;
	for ig=1:size(metpp(:,1),1)-1,
		ddt=metpp(ig+1,1)-metpp(ig,1);
		if (ddt > 0.008),
			gap=[metpp(ig,1):tinc:metpp(ig+1,1)];
			gap(length(gap))=[];
			gap(1)=[];
			bad=ones(length(gap),1).*NaN;
			nt=[nt
			    metpp(ig,1)
			    gap'];
			nv=[nv
			    vp(ig)
			    bad];
			nu=[nu
			    up(ig)
			    bad];
		else,
			nt=[nt
			    metpp(ig,1)];
			nv=[nv
                            vp(ig)];
			nu=[nu
			    up(ig)];
		end
	end
	nt=[nt
	    metpp(size(metp,1),1)];
	nu=[nu
	    up(size(metp,1))];
	nv=[nv
      	    vp(size(metp,1))];	
	bad=-99;
	lwu=isnan(nu)==1;
	wu(lwu)=-99;
	wu(~lwu)=nu(~lwu);
	lwu=isnan(nv)==1;
	wv(lwu)=-99;
	wv(~lwu)=nv(~lwu);
	wi=[nt(:) wu(:) wv(:)];
	wo=[fpath,'wind',num2str(mid),'.dat'];
	writedata(wi,wo);
	ll=nt(:) > ydnow-2;
	ww=wi(ll,:);
	wwf=[fpath,'wind1',num2str(mid),'.dat'];
	writedata(ww,wwf);
	windplot(nt,nu,nv,year);
	fclose(wid);
	end
% now to write out a simple table of numbers for a perl
% script to load and create an htlm table from
bad=-99;
%# Do we want regular plots of buoy heading and relative wind direction?
for i=2:8,
	metplot(metp(:,1),metp(:,i),year,mid,i,flagp(:,i-1));
end
dtplot(metp(:,1),year,mid,1);
ffn=[fpath,'m',num2str(mid),'_nowmet.dat'];
if isempty(metd)==0,
zt=ydtostr(year,metd(lm,1));
zts=datestr(zt,0);
fid=fopen(ffn,'w');
fprintf(fid,'%s\n',zts(1,:));
if lw >= 7,
	lw=lw-1;
end
for jj=2:lw-2,
	amet=metd(:,jj);
	ll=flagd(:,jj-1)==0;
	amet=amet(ll);
	if isempty(amet)==0,
		fprintf(fid,'%7.2f %d %7.2f ',metd(lm,jj),flagd(lm,jj-1),mean(amet));
		fprintf(fid,'%7.2f %7.2f ',std(amet),min(amet));
		fprintf(fid,'%7.2f \n',max(amet));
	else
		fprintf(fid,'%7.2f %d %7.2f ',bad,flagd(lm,jj-1),bad);
		fprintf(fid,'%7.2f %7.2f ',bad,bad);
		fprintf(fid,'%7.2f \n',bad);
	end		
end
	amet=metd(:,5);
	amed=metd(:,6);
	ll=flagd(:,4)==0;
	amed=amed(ll);
	amet=amet(ll);
	if isempty(amet)==0,
		y=vectavg(amet,amed);
		fprintf(fid,'%7.2f %d %7.2f ',metd(lm,5),flagd(lm,5-1),y(1));
		fprintf(fid,'%7.2f %7.2f ',std(amet),min(amet));
		fprintf(fid,'%7.2f \n',max(amet));
		fprintf(fid,'%7.2f %d %7.2f ',metd(lm,6),flagd(lm,6-1),y(2));
		fprintf(fid,'%7.2f %7.2f ',bad,min(amed));
		fprintf(fid,'%7.2f \n',max(amed));
	end
	for gi=7:8,
		amed=metd(:,gi);
%		ll=flagd(:,gi-1)==0;
		ll=(isnan(amed)==0)&(amed~=-99);
		amed=amed(ll);
		amet=ones(size(amed,1),1);
		if isempty(amed)==0,
			y=vectavg(amet,amed);
			fprintf(fid,'%7.2f %d %7.2f ',metd(lm,gi),flagd(lm,gi-1),y(2));
			fprintf(fid,'%7.2f %7.2f ',bad,min(amed));
			fprintf(fid,'%7.2f \n',max(amed));
		end
	end
end
end
%
 quit
