%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M-file to process buoy m"x" ADCP 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.  At this point we will
% rotate the data for magnetic declination and
% attempt to remove mooring motion (lateral) from
% the computed velocities.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;
ylat=[36.717 36.783
      36.630 36.717
      36.517 36.617
      36.100 36.200];
ylon=[-122.05 -122.00
      -122.45 -122.37
      -123.00 -122.95
      -122.47 -122.40];
% 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=time(1);
month=time(2);
day=time(3);
hour=time(4);
minu=time(5);
sec=0;
yday=time(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),'_adcp_uv_',monc(month,:),'_',num2str(year)];
fnm2=['m',num2str(mid),'_adcp_sc_',monc(month,:),'_',num2str(year)];
% see if the file exists
fid1=fopen([pathm1,fnm1,'.mat'],'r+');
fid2=fopen([pathm1,fnm2,'.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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The adcp_uv file will contain the following arrays:
% un-->U array vn-->V array en-->Error velocity array
% cmn-> Correlation magnitude array.
% gn--> AGC array gpsn-->GPS data array (time,lat,lon,stat)
% tn--> Time vector for ADCP data (NOT GPS data).
% fu--> Flag array for U velocities
% fv--> Flag array for V velocities
% fe--> Flag array for error velocity
% fg--> Flag array for AGC
% fgps--> Flag array for GPS info.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if fid1 > 2,
	fclose(fid1);
	eval(['load ',pathm1,fnm1,'.mat']);
% load velocity arrays
	eval(['bu=un;']);
	eval(['bv=vn;']);
	eval(['be=en;']);
% load correlation magnitude arrays (multidimensional)
	eval(['bc=cmn;']);
% load echo intensity (agc) array (multidimensional)
	eval(['bg=gn;']);
% load percent 4 beam good solution array
	eval(['bg4=pg4n;']);
% load percent 3 beam good solution array
	eval(['bg3=pg3n;']);
% load time vector
	eval(['bt=tn;']);
% load gps data array
	eval(['bgps=gpsn;']);
% load velocity flags
	eval(['bfu=fu;']);
	eval(['bfv=fv;']);
% load gps data flags
	eval(['bfgps=fgps;']);
else,
	bu=[];
	bv=[];
	be=[];
	bc=[];
	bg=[];
	bg4=[];
	bg3=[];
	bt=[];
	bgps=[];
	bfu=[];
	bfv=[];
	bfgps=[];
end
% Now load the scalar information
% add flag to this for transducer temperature
if fid2 > 2,
	fclose(fid2);
	eval(['load ',pathm1,fnm2,'.mat']);
	eval(['c=adcpscn;']);
	eval(['cf=adcpfsc;']);
else,
	c=[];
	cf=[];
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),
	au=[];
	av=[];
	ae=[];
	ac=[];
	ag=[];
	ag4=[];
	ag3=[];
	at=[];
	agps=[];
	afu=[];
	afv=[];
	afgps=[];
	d=[];
	df=[];
else,
	if (month==1),
		yearo=year-1;
		patho=['/data/span9/bahr/m',num2str(mid),'/adcp/',num2str(yearo),'/'];
		fnm1o=['m',num2str(mid),'_adcp_uv_dec_',num2str(yearo)];
		fnm2o=['m',num2str(mid),'_adcp_sc_dec_',num2str(yearo)];
	else,
		yearo=year;
		patho=pathm1;
		fnm1o=['m',num2str(mid),'_adcp_uv_',monc(month-1,:),'_',num2str(year)];
		fnm2o=['m',num2str(mid),'_adcp_sc_',monc(month-1,:),'_',num2str(year)];
	end
% check to see if old file is there.  If it is load it.
	fid1o=fopen([patho,fnm1o,'.mat'],'r+');
	fid2o=fopen([patho,fnm2o,'.mat'],'r+');	
	if (fid1o > 2),
		fclose(fid1o);
		eval(['load ',patho,fnm1o,'.mat']);
		eval(['au=un;']);
		eval(['av=vn;']);
		eval(['ae=en;']);
		eval(['ac=cmn;']);
		eval(['ag=gn;']);
		eval(['ag4=pg4n;']);
		eval(['ag3=pg3n;']);
		eval(['at=tn;']);
		eval(['agps=gpsn;']);
		eval(['afu=fu;']);
		eval(['afv=fv;']);
		eval(['afgps=fgps;']);
	else,
		au=[];
		av=[];
		ae=[];
		ac=[];
		ag=[];
		ag4=[];
		ag3=[];
		at=[];
		agps=[];
		afu=[];
		afv=[];
		afgps=[];
	end
	if (fid2o > 2),
		fclose(fid2o);
		eval(['load ',patho,fnm2o,'.mat']);
		eval(['d=adcpscn;']);
		eval(['df=adcpfsc;']);
	else,
		d=[];
		df=[];
	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(at)==0,
		at(:,1)=at(:,1)-maxd;
	         d(:,1)=d(:,1)-maxd;
	end
end
% Set up the vector arrays.  Note we need to check sizes to see if
% we need to augment an array in shape or size.
% we also need to worry about depths.
% If we pad velocity and other arrays we also need to fill in the 
% flag arrays.  We fill in the flag arrays with the 3 value (NO DATA).
wa=size(au,2); wb=size(bu,2);
if wa~=wb,
	if wa < wb,
		pad=ones(size(au,1),length([1:(wb-wa)])).*NaN;
		pad3=ones(size(au,1),length([1:(wb-wa)]),4).*NaN;
		padf=ones(size(au,1),length([1:(wb-wa)])).*3;
		au=[au pad];
		av=[av pad];
		ae=[ae pad];
% padding the mulitdimensional arrays will be fun...
		ac=[ac pad3];
		ag=[ag pad3];
% percent good array can be padded easily
		ag4=[ag4 pad];
		ag3=[ag3 pad];
		afu=[afu padf];
		afv=[afv padf];
	else,
		pad=ones(size(bu,1),length([1:(wa-wb)])).*NaN;
		pad3=ones(size(bu,1),length([1:(wa-wb)]),4).*NaN;
		padf=ones(size(bu,1),length([1:(wa-wb)])).*3;
		bu=[bu pad];
		bv=[bv pad];
		be=[be pad];
% multi-dim arrays
		bc=[bc pad3];
		bg=[bg pad3];
% percent good arrays
		bg4=[bg4 pad];
		bg3=[bg3 pad];		
		bfu=[bfu padf];
		bfv=[bfv padf];
	end
end
% Now we can merge the arrays. 
adt=[at
     bt];
u=[au
   bu];
uf=[afu
    bfu];
v=[av
   bv];
vf=[afv
    bfv];
e=[ae
   be];
cm=[];
% I hope this concattenation of multi-dim arrays works
% Don't want to concattenate empty multidimensional arrays.
if isempty(ac)==0,
	cm=ac;
end
if isempty(bc)==0,
	cm=[cm
    	    bc];
end
g=[];
if isempty(ag)==0,
	g=ag;
end
if isempty(bg)==0,
	g=[g
	   bg];
end
pg4=[ag4
     bg4];
pg3=[ag3
     bg3];
gps=[agps
     bgps];
gpsf=[afgps
      bfgps];
% set up the corresponding adcp scalar data arrays.
adcpsc=[d
        c];
fsc=[df
     cf];
% we should now have the old matlab process adcp files loaded.
% The next step is to load the new files with data to be processed.
% up to this point we have not decide how the arrays are stored in
% the ADCP files.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% So we can handle changes in number of velocities put out by the
% ADCP's we will create arrays which grow in time and are n columns
% wide (where n is the greater bin number) so we will have:
% Depth vector       U array       V array     E array      AGC array by 4
% D1--->D(n)      U1--->U(n)      V1--->V(n)   E1--->E(n)   A1--->A(n)
%                 |               |            |            |  
%                 |               |            |            |
%                 v               v            v            v
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now also load GPS data for mooring motion if the files exist......
% What do we do if there is no GPS data for a mooring....
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ff1=['m',num2str(mid),'_adcp_uv.dat'];
fs1=['m',num2str(mid),'_adcp_sc.dat'];
ff2=['m',num2str(mid),'_adcpo_uv.dat'];
fs2=['m',num2str(mid),'_adcpo_sc.dat'];
fg1=['m',num2str(mid),'_gps.dat'];
fg2=['m',num2str(mid),'_gpso.dat'];
fod1=fopen([pathm1,ff1],'r');
fsd1=fopen([pathm1,fs1],'r');
fog1=fopen([pathm1,fg1],'r');
if(month==1)&(day >=2),
	fod2=fopen([pathm1,ff2],'r');
	fsd2=fopen([pathm1,fs2],'r');
	fog2=fopen([pathm1,fg2],'r');
else,
	fod2=fopen([patho,ff2],'r');
	fsd2=fopen([patho,fs2],'r');
	fog2=fopen([patho,fg2],'r');
end
if fod1 > 2,
	fclose(fod1);
	eval(['load ',pathm1,'m',num2str(mid),'_adcp_uv.dat']);
	eval(['m1uv=m',num2str(mid),'_adcp_uv;']);
else,
	m1uv=[];
end
if fsd1 > 2,
	fclose(fsd1);
	eval(['load ',pathm1,'m',num2str(mid),'_adcp_sc.dat']);
	eval(['m1sc=m',num2str(mid),'_adcp_sc;']);
else,
	m1sc=[];
end
if fod2 > 2,
	fclose(fod2);
	if(month==1)&(day>=2),
	eval(['load ',pathm1,'m',num2str(mid),'_adcpo_uv.dat']);
	else,
	eval(['load ',patho,'m',num2str(mid),'_adcpo_uv.dat']);
	end
	eval(['m1ouv=m',num2str(mid),'_adcpo_uv;']);
else,
	m1ouv=[];
end
if fsd2 > 2,
	fclose(fsd2);
	if(month==1)&(day >=2),
	eval(['load ',pathm1,'m',num2str(mid),'_adcpo_sc.dat']);
	else,
	eval(['load ',patho,'m',num2str(mid),'_adcpo_sc.dat']);
	end
	eval(['m1osc=m',num2str(mid),'_adcpo_sc;']);
else,
	m1osc=[];
end
if fog1 > 2,
	fclose(fog1);
	eval(['load ',pathm1,'m',num2str(mid),'_gps.dat']);
	eval(['mgps=m',num2str(mid),'_gps;']);
else,
	mgps=[];
end
if fog2 > 2,
	fclose(fog2);
	if(month==1)&(day >=2),
	eval(['load ',pathm1,'m',num2str(mid),'_gpso.dat']);
	else,
	eval(['load ',patho,'m',num2str(mid),'_gpso.dat']);
	end
	eval(['mgpso=m',num2str(mid),'_gpso;']);
else,
	mgpso=[];
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)&(isempty(m1ouv)==0)),
	m1ouv(:,1)=m1ouv(:,1)-maxd;
	m1osc(:,1)=m1osc(:,1)-maxd;
end
% can't just add as aline above because we may not have a gps signal but
% have adcp data.
if ((day==1)&(month==1)&(isempty(mgpso)==0)),
	mgpso(:,1)=mgpso(:,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=[m1ouv
    m1uv];
% m2 contains the ADCP scalar data corresponding to the vector data in m1.
m2=[m1osc
    m1sc];
% mg will contain the gps information
mg=[mgpso
    mgps];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next we want to get the GPS data sorted and checked out
% At this point we want to get GPS data cleaned up and ready 
% to compute velocities from.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (isempty(mg) == 0),
	dtg=diff(mg(:,1));
	[inegg]=find(dtg < 0);
	if isempty(inegg)==0,
		[newtg indg]=sort(mg(:,1));
		mg=[newtg mg(indg,2:size(mg,2))];
	end
% remove duplicate times
	dt=diff(mg(:,1));
	idi=find(dt==0);
	mg(idi,:)=[];
	tmin=min(mg(:,1));
	if isempty(gps)==0,
		tmax=max(gps(:,1));
	else
		tmax=0;
	end
	lgp=size(mg,1);
	ig=zeros(lgp,1);
	if isempty(gps)==0,
		for i=1:lgp,
			if mg(i,1) > tmax,
				i=lgp;
			end
			igz=find(gps(:,1)==mg(i,1));
			if isempty(igz)==0,
				ig(i)=1;
			end
		end
	end
	mg(logical(ig),:)=[]; % remove identical times.
end
mymax=ylat(mid,2);
mymin=ylat(mid,1);
mxmax=ylon(mid,2);
mxmin=ylon(mid,1);
if(isempty(mg)==0),
% Now we eliminate data based upon some criteria from the above data.
% Add arrays earlier defining the max/min lat/lon for the mooring positions
% Then flag the data.
% Define a flag array for new gps data.
	fmg=zeros(size(mg));
% Check which values are out of bounds and flag them.
	ll=((mg(:,2) < mymin)|(mg(:,2) > mymax));
	fmg(ll,2)=1;
	ll=((mg(:,3) < mxmin)|(mg(:,3) > mxmax));
	fmg(ll,3)=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==>fill in here determination if data is good.<==
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now fill in the gps array with the new data....
	i=0; j=0;
	lgp=size(gps,1);
	gmask=zeros(lgp,1);
	lnew=size(mg,1);
	k=1;
	while(k <= lnew),
		lgp=size(gps,1);
% matlab matrix is empty so all data is new and fill it in
		if lgp==0,
			gps=mg;
			gpsf=fmg(:,[2 3]);
			k=lnew+1;
			gmask=[ones(lnew,1)];
		else,
% matlab matix is not empty and fill it in correct spots
			[index]=find(gps(:,1) < mg(k,1));
			idmz=max(index);
% check if we are at the end of the array if so fill to the
% end
			if idmz==lgp,
				gps=[gps
				     mg([k:lnew],:)];
				gpsf=[gpsf
				      fmg([k:lnew],[2:3])];
				gmask=[gmask
				       ones(length([k:lnew]),1)];
% we just have one point in the middle to put in
			else,
				gps=[gps([1:idmz],:)
		                     mg(k,:)
                                     gps([idmz+1:lgp],:)];
				gpsf=[gpsf([1:idmz],:)
				      fmg(k,[2:3])
				      gpsf([idmz+1:lgp],:)];
				gmask=[gmask([1:idmz])
 				       1
                                       gmask([idmz+1:lgp])];
			end
		end
		k=k+1;
	end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%>>>>****Ended here until figured out what to do with 
% GPS data (keep it in a file with ADCP or seperate
%%%%%>>>>
% compute current year day time
ydnow=mond(month)+day+hour/24+minu/(24*60);
if (rem(year,4)==0)&(month > 2),
	ydnow=ydnow+1;
end
ydnowl=ydnow-1;
%% Sort out the ADCP data now.
% we now have the data.  Now we need to check times and sort times?
% since the scalar data comes from the same as the vector data we
% should process them at the same time here.  Both should have the
% same time stamps.
% ADCP array formats are given above.  We will have a Time vector and
% Depth vector with arrays of U,V, Error Vel, AGC
% Note m1 and m2 should be the same length.  They are both computed from the
% same files.  If for some wierd reason they are not the same we need to worry.
if ((isempty(m1) == 0)&(isempty(m2) == 0)),
	dt1=diff(m1(:,1));
	dt2=diff(m2(:,1));
%check for neg values to see if we have times out of order.
	[ineg1]=find(dt1 < 0);
	[ineg2]=find(dt2 < 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(ineg1)==0,
% we need to do a sort of the file
		[newt1 idex1]=sort(m1(:,1));
		m1=[newt1 m1(idex1,2:size(m1,2))];
	end
% remove duplicate times (both arrays should be the same size
	dt=diff(m1(:,1));
	idi=find(dt==0);
	m1(idi,:)=[];
	if isempty(ineg2)==0,
		[newt2 idex2]=sort(m2(:,1));
		m2=[newt2 m2(idex2,2:size(m2,2))];
	end
	dt=diff(m2(:,1));
	idi=find(dt==0);
	m2(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).
% find the minimum time in the currently captured files
	tmin1=min(m1(:,1));
	tmin2=min(m2(:,1));
	if isempty(adt)==0,
		tmax1=max(adt);
	else,
		tmax1=0;
	end
	if isempty(adcpsc)==0,
		tmax2=max(adcpsc(:,1));
	else,
		tmax2=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);
	lm2=size(m2(:,1),1);
	im1=zeros(lm1,1);
	im2=zeros(lm2,1);
	if isempty(adt)==0,
		for i=1:lm1,
			if m1(i,1) > tmax1,
				i=lm1;
			end
			id=find(adt(:,1)==m1(i,1));
			if isempty(id)==0,
				im1(i)=1;
			end
		end
	end
	if isempty(adcpsc)==0,
		for i=1:lm2,
			if m2(i,1) > tmax2,
				i=lm2;
			end
			id=find(adcpsc(:,1)==m2(i,1));
			if isempty(id)==0,
				im2(i)=1;
			end
		end
	end
	m1(logical(im1),:)=[]; % Remove points we already have in the files.
	m2(logical(im2),:)=[];
end
if ((isempty(m1) == 0)&(isempty(m2) == 0)),
%
% 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.
% 
	i=0; j=0;
% adt=acoustic doppler time (vector)
	ladt=size(adt,1);
	ladcpsc=size(adcpsc,1);
% even though we have seperate arrays for U,V,and such we only need one
% mask for point insertion since we insert for all arrays at the same point.
% we do however need a seperate mask for masking the insertion of good 
% adcp x-ducer temps.
% (we should only need one mask, theoretically the sc file and the uv file
% should be the same sizes, use 2 for now and change later?)
	bmask1=zeros(ladt,1);
	bmask2=zeros(ladcpsc,1);
	lnew1=size(m1,1);
	lnew2=size(m2,1);
	k=1;
% Format of m1 data file is:
% time,u(1)-->u(n),9999,v(1)-->v(n),9999,e(1)-->e(n),9999,cm1(1)-->cm1(n),9999,
% cm2(1)->cm2(n),9999,cm3(1)->cm3(n),9999,...9999,agc1(1)->agc1(n),9999,
% agc2(1)->agc2(n),9999,...,9999,pg3(1)->pg3(n),9999,pg4(1)->pg4(n).
% Format of m2 data file is:
% time,head,stdhead,pitch,stdpitch,roll,stdroll,xtemp,nbin,nping,
% bin_length,blank_length,depth,ver,rev,sspd,tpmin,tpsec,tphun,
% ex,ez,wa
% We only want to flag heading, stdhead, pitch,stdpitch, roll, stdroll, and xtemp.
%
% before we start the while loop we need to know which array will be
% bigger n or the old arrays.
	nbin=m2(1,9);
	narr=size(u,2);
% This is the usual case arrays are the same size and we don't
% have to pad either to get them to fit.
	if nbin==narr,
		lpad=0;
	end
% This case old array had more bins so we need to pad new data to
% fill in old array size.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: At this point we need to recompute the depth vector
% if nbin < narr we don't need a new vector, BUT it could be
% possible that the deployment has changed the depth at which
% the ADCP is at.  Which will change all the bins.
% if nbin > narr: we need to add to the depth vector.  Also
% the depth of the ADCP could change and cause all types of
% problems.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: Proper depths can be gotten for all bins by computing 
% using the values in the scalar ADCP data file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if nbin < narr,
		lpad=narr-nbin;
	end
% This case old arrays are too small and we need to expand to fill
% in the data.
	if nbin > narr,
		lpad=nbin-narr;
		lx=size(u,1);
		pad=ones(lx,lpad).*NaN;
		pad3=ones(lx,lpad,4).*NaN;
		u=[u pad];
		v=[v pad];
		e=[e pad];
		g=[g pad3];
		cm=[cm pad3];
		pg3=[pg3 pad];
		pg4=[pg4 pad];
		uf=[uf pad];
		vf=[vf pad];
		lpad=0;
	end
%*************************************************************
%************************************************************
% Warning!! We can't just stuff velocities in.  The stored
% velocities will be corrected for mooring motion AND rotated
% to degrees true (take into account magnetic declination)
%************************************************************
%************************************************************
% Unfortuneatly, we can't quite use multidimensional arrays like
% the other arrays.  If we append to 1 "page" the others get a
% row of zeros which we don't want.  So we pull out the pages
% temporarily as we append data and then recreate our multi-dim
% array.
	if(isempty(cm)==0),
		c1=squeeze(cm(:,:,1));
		c2=squeeze(cm(:,:,2));
		c3=squeeze(cm(:,:,3));
		c4=squeeze(cm(:,:,4));
		g1=squeeze(g(:,:,1));
		g2=squeeze(g(:,:,2));
		g3=squeeze(g(:,:,3));
		g4=squeeze(g(:,:,4));
	else,
		c1=[]; c2=[]; c3=[]; c4=[];
		g1=[]; g2=[]; g3=[]; g4=[];
	end
	clear cm; clear g;
	while(k <= lnew1),
		ladt=size(adt,1);
		if ladt==0,
% for this case all arrays are empty and we just fill in
			adt=m1(:,1);
			k=lnew1+1;
			bmask1=[bmask1
		       	       ones(lnew1,1)];
			u=m1(:,[2:(nbin+1)]).*100;
			v=m1(:,[(nbin+3):(2*nbin+2)]).*100;
			e=m1(:,[(2*nbin+4):(3*nbin+3)]).*100;
			c1=m1(:,[(3*nbin+5):(4*nbin+4)]);
			c2=m1(:,[(4*nbin+6):(5*nbin+5)]);
			c3=m1(:,[(5*nbin+7):(6*nbin+6)]);
			c4=m1(:,[(6*nbin+8):(7*nbin+7)]);
			g1=m1(:,[(7*nbin+9):(8*nbin+8)]);
			g2=m1(:,[(8*nbin+10):(9*nbin+9)]);
			g3=m1(:,[(9*nbin+11):(10*nbin+10)]);
			g4=m1(:,[(10*nbin+12):(11*nbin+11)]);
			pg3=m1(:,[(11*nbin+13):(12*nbin+12)]);
			pg4=m1(:,[(12*nbin+14):(13*nbin+13)]);
			uf=zeros(lnew1,nbin);
			vf=zeros(lnew1,nbin);
		else,
% here the arrays are not empty and we have to be carefull.
			[index]=find(adt(:,1) < m1(k,1));
			idmz=max(index);
			if idmz==ladt,
				adt=[adt
                     	     	     m1([k:lnew1],1)];
				bmask1=[bmask1
			               ones(length(k:lnew1),1)];
				pad=ones(length(k:lnew1),lpad).*NaN;
				uf=[uf
				    zeros(length(k:lnew1),nbin) pad];
				vf=[vf
				    zeros(length(k:lnew1),nbin) pad];
				u=[u
                                   m1([k:lnew1],[2:(nbin+1)]).*100 pad];
				v=[v
				   m1([k:lnew1],[(nbin+3):(2*nbin+2)]).*100 pad];
				e=[e
                                   m1([k:lnew1],[(2*nbin+4):(3*nbin+3)]).*100 pad];
				c1=[c1
				   m1([k:lnew1],[(3*nbin+5):(4*nbin+4)]) pad];
				c2=[c2
				   m1([k:lnew1],[(4*nbin+6):(5*nbin+5)]) pad];
				c3=[c3
				   m1([k:lnew1],[(5*nbin+7):(6*nbin+6)]) pad];
				c4=[c4
				   m1([k:lnew1],[(6*nbin+8):(7*nbin+7)]) pad];
				g1=[g1
				   m1([k:lnew1],[(7*nbin+9):(8*nbin+8)]) pad];
				g2=[g2
				   m1([k:lnew1],[(8*nbin+10):(9*nbin+9)]) pad];
				g3=[g3
				   m1([k:lnew1],[(9*nbin+11):(10*nbin+10)]) pad];
				g4=[g4
				   m1([k:lnew1],[(10*nbin+12):(11*nbin+11)]) pad];
				pg3=[pg3
				   m1([k:lnew1],[(11*nbin+13):(12*nbin+12)]) pad];
				pg4=[pg4
				   m1([k:lnew1],[(12*nbin+14):(13*nbin+13)]) pad];
			 	k=lnew1;
			else
				adt=[adt([1:idmz],:)
		     	     	     m1(k,1)
		     	             adt([idmz+1:ladt],:)];
				bmask1=[bmask1([1:idmz])
				       1
			               bmask1([idmz+1:ladt])];
				pad=ones(1,lpad).*NaN;
				uf=[uf([1:idmz],:)
				    zeros(1,nbin) pad
				    uf([idmz+1:ladt],:)];
				vf=[vf([1:idmz],:)
				    zeros(1,nbin) pad
				    vf([idmz+1:ladt],:)];
				u=[u([1:idmz],:)
				   m1(k,[2:(nbin+1)]).*100 pad
			           u([idmz+1:ladt],:)];
				v=[v([1:idmz],:)
				   m1(k,[(nbin+3):(2*nbin+2)]).*100 pad
			           v([idmz+1:ladt],:)];
				e=[e([1:idmz],:)
				   m1(k,[(2*nbin+4):(3*nbin+3)]).*100 pad
				   e([idmz+1:ladt],:)];
				c1=[c1([1:idmz],:)
				   m1(k,[(3*nbin+5):(4*nbin+4)]) pad
				   c1([idmz+1:ladt],:)];
				c2=[c2([1:idmz],:)
				   m1(k,[(4*nbin+6):(5*nbin+5)]) pad
				   c2([idmz+1:ladt],:)];
				c3=[c3([1:idmz],:)
				   m1(k,[(5*nbin+7):(6*nbin+6)]) pad
				   c3([idmz+1:ladt],:)];
				c4=[c4([1:idmz],:)
				   m1(k,[(6*nbin+8):(7*nbin+7)]) pad
				   c4([idmz+1:ladt],:)];
				g1=[g1([1:idmz],:)
				   m1(k,[(7*nbin+9):(8*nbin+8)]) pad
				   g1([idmz+1:ladt],:)];
				g2=[g2([1:idmz],:)
				   m1(k,[(8*nbin+10):(9*nbin+9)]) pad
				   g2([idmz+1:ladt],:)];
				g3=[g3([1:idmz],:)
				   m1(k,[(9*nbin+11):(10*nbin+10)]) pad
				   g3([idmz+1:ladt],:)];
				g4=[g4([1:idmz],:)
				   m1(k,[(10*nbin+12):(11*nbin+11)]) pad
				   g4([idmz+1:ladt],:)];
				pg3=[pg3([1:idmz],:)
				   m1(k,[(11*nbin+13):(12*nbin+12)]) pad
				   pg3([idmz+1:ladt],:)];
				pg4=[pg4([1:idmz],:)
				   m1(k,[(12*nbin+14):(13*nbin+13)]) pad
				   pg4([idmz+1:ladt],:)];
			end
		end
		k=k+1;
	end
	cm(:,:,1)=c1;
	cm(:,:,2)=c2;
	cm(:,:,3)=c3;
	cm(:,:,4)=c4;
	clear c1; clear c2; clear c3; clear c4;
	g(:,:,1)=g1;
	g(:,:,2)=g2;
	g(:,:,3)=g3;
	g(:,:,4)=g4;
	clear g1; clear g2; clear g3; clear g4;
	k=1;
	while(k <= lnew2),
		ladcpsc=size(adcpsc,1);
		if ladcpsc==0,
			adcpsc=m2;
			k=lnew2+1;
			bmask2=[bmask2
		       	       ones(lnew2,1)];
			fsc=zeros(lnew2,size(adcpsc,2)-1);
		else,
			[index]=find(adcpsc(:,1) < m2(k,1));
			idmz=max(index);
			if idmz==ladcpsc,
				adcpsc=[adcpsc
                     	     	        m2([k:lnew2],:)];
				bmask2=[bmask2
			               ones(length(k:lnew2),1)];
				fsc=[fsc
				     zeros(length(k:lnew2),size(m2,2)-1)];
			 	k=lnew2;
			else
				adcpsc=[adcpsc([1:idmz],:)
		     	     	     	m2(k,:)
		     	             	adcpsc([idmz+1:ladcpsc],:)];
				bmask2=[bmask2([1:idmz])
				       1
			               bmask2([idmz+1:ladcpsc])];
				fsc=[fsc([1:idmz],:)
				     zeros(1,size(m2,2)-1)
				     fsc([idmz+1:ladcpsc],:)];
			end
		end
		k=k+1;
	end
% Do a simple test to put marked bad velocities (by the adcp software)
% to have flags of 3 (no data).
%	if mid==2,
% for a narrow band we need to fake things abit, rely on pg3 pg4
%
%		cmnb=65.*ones(size(cm));
%		uf=uvflag_m2(adt,cm,g,pg3,pg4,nbin);
%		uf=uvflag(adt,cmnb,g,pg3,pg4,nbin);
%	else,
		uf=uvflag(adt,cm,g,pg3,pg4,nbin);
%	end
	vf=uf;
	for yi=1:nbin,
		yu=u(:,yi)<=-3276.8;
		yv=v(:,yi)<=-3276.8;
		uf(yu,yi)=3;
		vf(yv,yi)=3;
	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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we can apply our gradient filter to the data (?)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For now, no filter on adcp velocities........
%====> put filter here <====%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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,
		tn=adt;
		un=u;
		vn=v;
		en=e;
		cmn=cm;
		gn=g;
		pg3n=pg3;
		pg4n=pg4;
		fu=uf;
		fv=vf;
		gpsn=gps;
		fgps=gpsf;
		adcpscn=adcpsc;
		adcpfsc=fsc;
		eval(['save ',pathm1,fnm1,' tn un vn en cmn gn pg3n pg4n fu fv gpsn fgps']);
		eval(['save ',pathm1,fnm2,' adcpscn adcpfsc']);
	else
		if month==1,
			ll=adt(:,1) < 0;
			tn=adt(ll)+maxd;
			un=u(ll,:);
			vn=v(ll,:);
			en=e(ll,:);
			cmn=cm(ll,:,:);
			gn=g(ll,:,:);
			pg3n=pg3(ll,:);
			pg4n=pg4(ll,:);
			fu=uf(ll,:);
			fv=vf(ll,:);
% GPS data has it's own time column.  So we need to do it seperate
			lg=gps(:,1) < 0;
			gpsn=gps(lg,:);
			gpsn(:,1)=gpsn(:,1)+maxd;
%			ilg=find(gps(:,1) < 0);
%			gps(ilg,1)=gps(ilg,1)+maxd;
%			gpsn=gps(lg,:);
			fgps=gpsf(lg,:);
			adcpscn=adcpsc(ll,:);
			adcpscn(:,1)=adcpscn(:,1)+maxd;
			adcpfsc=fsc(ll,:);
			eval(['save ',patho,fnm1o,' tn un vn en cmn gn pg3n pg4n fu fv gpsn fgps']);
			eval(['save ',patho,fnm2o,' adcpscn adcpfsc']);
			clear tn; clear un; clear vn; clear en;
			clear cmn; clear pg3n; clear pg4n; clear fu; clear fv;
			clear fgps; clear adcpfsc;
			clear gn; clear gpsn; clear adcpscn;
			tn=adt(~ll);
			un=u(~ll,:);
			vn=v(~ll,:);
			en=e(~ll,:);
			cmn=cm(~ll,:,:);
			gn=g(~ll,:,:);
			pg3n=pg3(~ll,:);
			pg4n=pg4(~ll,:);
			fu=uf(~ll,:);
			fv=vf(~ll,:);
			gpsn=gps(~lg,:);
			fgps=gpsf(~lg,:);
			adcpscn=adcpsc(~ll,:);
			adcpfsc=fsc(~ll,:);
			eval(['save ',pathm1,fnm1,' tn un vn en cmn gn pg3n pg4n fu fv gpsn fgps']);
			eval(['save ',pathm1,fnm2,' adcpscn adcpfsc']);
		else
			dmax=mond(month)+1;
			if (month > 2)&(rem(year,4)==0),
				dmax=dmax+1;
			end
			ll=adt(:,1) < dmax;
			tn=adt(ll);
			un=u(ll,:);
			vn=v(ll,:);
			en=e(ll,:);
			cmn=cm(ll,:,:);
			gn=g(ll,:,:);
			pg3n=pg3(ll,:);
			pg4n=pg4(ll,:);
			fu=uf(ll,:);
			fv=vf(ll,:);
% GPS has it's own time stamp.
			lg=gps(:,1) < dmax;
			gpsn=gps(lg,:);
			fgps=gpsf(lg,:);
			adcpscn=adcpsc(ll,:);
			adcpfsc=fsc(ll,:);
%			adcpscn(:,1)=adcpscn(:,1)+maxd;
			eval(['save ',patho,fnm1o,' tn un vn en cmn gn pg3n pg4n fu fv gpsn fgps']);
			eval(['save ',patho,fnm2o,' adcpscn adcpfsc']);
			clear tn; clear un; clear vn; clear en;
			clear cmn; clear fu; clear fv; clear pg3n; clear pg4n;
			clear gn; clear gpsn; clear adcpscn;
			clear fgps; clear adcpfsc;
			tn=adt(~ll);
			un=u(~ll,:);
			vn=v(~ll,:);
			en=e(~ll,:);
			cmn=cm(~ll,:,:);
			gn=g(~ll,:,:);
			pg3n=pg3(~ll,:);
			pg4n=pg4(~ll,:);
			fu=uf(~ll,:);
			fv=vf(~ll,:);
			gpsn=gps(~lg,:);
			fgps=gpsf(~lg,:);
			adcpscn=adcpsc(~ll,:);
			adcpfsc=fsc(~ll,:);
			eval(['save ',pathm1,fnm1,' tn un vn en cmn gn pg3n pg4n fu fv gpsn fgps']);
			eval(['save ',pathm1,fnm2,' adcpscn adcpfsc']);

		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 isempty(adt) == 0,
	ydp=ydnow-7;
	if ydp < 0,
		ydp=ydp-1;
	end
	ll=(adt(:,1) > ydnowl)&(adt(:,1) <= ydnow);
	ls=(adcpsc(:,1) > ydnowl)&(adcpsc(:,1) <= ydnow);
	lp=(adt(:,1) > ydp)&(adt(:,1) <= ydnow);
	lq=(adcpsc(:,1) > ydp)&(adcpsc(:,1) <= ydnow);
%	max(adt)
% set up arrays of 24 hour data
	adtd=adt(ll,:);
	ud=u(ll,:);
	vd=v(ll,:);
	ufd=uf(ll,:);
	vfd=vf(ll,:);
% set up arrays of 7 days of data
	adtp=adt(lp,:);
	up=u(lp,:);
	vp=v(lp,:);
	ufp=uf(lp,:);
	vfp=vf(lp,:);
	adcpscd=adcpsc(ls,:);
	adcpfscd=fsc(ls,:);
	adcpscp=adcpsc(lq,:);
	adcpfscp=fsc(lq,:);
	lm1=size(adtd,1);
	lm2=size(adcpscd,1);
	lw1=size(adtd,2);
	lw2=size(adcpscd,2);
% The algorithm to compute transducer depths is
% depth=bin#*bin+blank+depth of Transducer+(transmit pulse width- bin width)/2
% where we have assumed that transmit pulse width=bin width
	for lkl=1:lm2,
     depd(lkl,:)=[1:adcpscd(lkl,9)].*adcpscd(lkl,11)+adcpscd(lkl,12)+adcpscd(lkl,13);
	end
	lmm2=size(adcpscp,1);
	for qs=1:lmm2,
     depp(qs,:)=[1:adcpscp(qs,9)].*adcpscp(qs,11)+adcpscp(qs,12)+adcpscp(qs,13);
	end
% now to write out a simple table of numbers for a perl
% script to load and create an htlm table from
	bad=-99;
% plot ADCP trancducer Temp
%	xdtplot(adcpscp(:,1),adcpscp(:,8),year,mid);
% plot UV profile
% we need to get the most current profile?
	hu=size(ud,1);
	if hu >= 1,
		uvprofplot(ud(hu,:),vd(hu,:),depd(hu,:),mid,ufd(hu,:),vfd(hu,:));
	end
	hv=size(up,1);
	if hv >=1,
		uvplot(adtp,up,vp,mid,ufp,vfp,year);
	end
% create uv profile information for plots
%%%%%%%%new code%%%%%%%%
	jjj=adtp < 0;
	adtt=adtp;
	adtt(jjj)=adtt(jjj)+1;
%%%%%%%%end new code %%%%%%%%
	if isempty(up)==0,
		kku=ufp > 0;
		kkv=vfp > 0;
		ucp=up;
		vcp=vp;
		ucp(kku)= -3200;
		vcp(kkv)= -3200;
	contouruv(adtt,depp(1,:),ucp,vcp,year,mid);
	[lu wu]=size(up);
% redo this bit below slightly.
% Below is data to writeout the uv velocity for contouring.  Ignore for
% now.
	ftemp=[];
	for ou=1:wu,
		ftemp=[ftemp
			adtt(:) ones(lu,1).*depp(:,ou) ucp(:,ou) vcp(:,ou)];
	end
	ov=(ftemp(:,3) < -3000)|(ftemp(:,4) < -3000);
	ftemp(ov,[3 4])=1e+35;
	fff=[fpath,'uv.dat'];
	writeuvdat(ftemp,fff);
	end
% plot U,V feathers for a set depth
%		uvfeather(adtp,up,vp,year,mid);
	aend=size(adtp,1);
	if aend~=0,
		ffn=[fpath,'m',num2str(mid),'_nowadcpsc.dat'];
		fid=fopen(ffn,'w');
		ffm=[fpath,'m',num2str(mid),'_nowadcpuv.dat'];
		fid1=fopen(ffm,'w');
		zt=ydtostr(year,adtp(aend));
		zts=datestr(zt,0);
		fprintf(fid,'%s\n',zts(1,:));
		fprintf(fid1,'%s\n',zts(1,:));
		jj=8;
		aadcpscd=adcpscd(:,jj);
		ll=adcpfscd(:,jj-1)==0;
		aadcpscd=aadcpscd(ll);
		if isempty(aadcpscd)==0,
		  fprintf(fid,'%7.2f %7.2f ',adcpscd(lm1,jj),mean(aadcpscd));
		  fprintf(fid,'%7.2f %7.2f ',std(aadcpscd),min(aadcpscd));
		  fprintf(fid,'%7.2f \n',max(aadcpscd));
		else
		  fprintf(fid,'%7.2f %7.2f ',bad,bad);
		  fprintf(fid,'%7.2f %7.2f ',bad,bad);
		  fprintf(fid,'%7.2f \n',bad);
		end
		if isempty(adcpscd)==0,
		  nbin=adcpscd(1,9);
		else,
		  nbin=0;
		end		
		for jj=1:nbin,
		  uu=ud(:,jj);
		  vv=vd(:,jj);
		  ll=ufd(:,jj)==0;
		  vl=vfd(:,jj)==0;
		  uu=uu(ll);
		  vv=vv(vl);
		  if isempty(uu)==0,
		    fprintf(fid1,'%6.1f %9.3f %7.3f ',depd(lm1,jj),ud(hu,jj),mean(uu));
			fprintf(fid1,'%7.3f %7.3f ',std(uu),min(uu));
			fprintf(fid1,'%7.3f ',max(uu));
		  else,
			fprintf(fid1,'%6.1f %9.3f %7.3f ',bad,bad,bad);
			fprintf(fid1,'%7.3f %7.3f ',bad,bad);
			fprintf(fid1,'%7.3f ',bad);
		  end
		  if isempty(vv)==0,
			fprintf(fid1,'%9.3f %7.3f ',vd(hu,jj),mean(vv));
			fprintf(fid1,'%7.3f %7.3f ',std(vv),min(vv));
			fprintf(fid1,'%7.3f\n',max(vv));
		  else,
			fprintf(fid1,'%9.3f %7.3f ',bad,bad);
			fprintf(fid1,'%7.3f %7.3f ',bad,bad);
			fprintf(fid1,'%7.3f\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
