disp('read_sonics_hr')
fclose('all');
clear all;
jd=input('julian day','s');
hr=input('gmt hr','s');
%sonic=input('sonic COM port #(2,3,4,5,6)','s');
%sn=str2num(sonic)-1;
hrs=str2num(hr)+1;
   
nx=10*(3600-1);
ucor=[+.4 -.1 -.05  +.00 +.3];
vcor=[-.2 .00 .00 +.00 -.3];

plotit=0;
plotit=input('Spectral plots=1, no plots=0');
tic;

for sn=1:5;%read thru each sonic
   sonic=num2str(sn+1);
	ur=zeros(1,36000);
	vr=zeros(1,36000);
	wr=zeros(1,36000);
	Tr=zeros(1,36000);
	stux=zeros(4,36000);
   
	e=['F:\users\data\97' jd hr '.p0'sonic];
   
   flist=fopen(e,'r');
   if flist>0
      
		stuff1=fgets(flist);%read 3 headers
		stuff2=fgets(flist);
		stuff3=fgets(flist);
		disp(stuff1);
		disp(stuff2);

		stux=fscanf(flist,'%*1s%g%*1s%g%*1s%g%*1s%g',[4,36000]);
		[a,np]=size(stux);
		uraw=stux(1,1:np);
		vraw=stux(2,1:np);
		wraw=stux(3,1:np);
		Traw=stux(4,1:np);

	   %despike
   	uz=zeros(np,2);
   	on=1:1:np;
		uz(1:np,1)=uraw';
		uz(1:np,2)=on';
		u=despik(uz);
		uz(1:np,1)=vraw';
		v=despik(uz);
		uz(1:np,1)=wraw';
		w=despik(uz);
		uz(1:np,1)=Traw';
		T=despik(uz);

		t=.1:.1:nx/10;
   	u=u+ucor(sn);
   	v=v+vcor(sn);
		um=mean(u);%means for entire time block (1 hr)
		vm=mean(v);
		wm=mean(w);
		Tm=mean(T);

	%rotate to streamwise flow
		s=sqrt(um^2+vm^2);
		thet=atan2(vm,um);
   	phi=atan2(wm,s);
   	diry(sn)=thet*180/pi;
   	tilt(sn)=phi*180/pi;
		rot(1,1)=cos(phi)*cos(thet);
		rot(1,2)=cos(phi)*sin(thet);
		rot(1,3)=sin(phi);
		rot(2,1)=-sin(thet);
		rot(2,2)=cos(thet);
		rot(2,3)=0;
		rot(3,1)=-sin(phi)*cos(thet);
		rot(3,2)=-sin(phi)*sin(thet);
		rot(3,3)=cos(phi);
		for j=1:np;
   		x=[u(j) v(j) w(j)];
   		xr=rot*x';
   		ur(j)=xr(1);
   		vr(j)=xr(2);
   		wr(j)=xr(3);
		end; %end for j
		nt=10*60*10;%plot first ten minutes time series
   	if plotit==1%plot if plots enabled
      	figure;
			plot(t(1:nt),ur(1:nt),t(1:nt),vr(1:nt),t(1:nt),wr(1:nt),t(1:nt),T(1:nt)-Tm+2);
			xlabel('time (sec)');
			ylabel('u,v,w,T');
     	 %	title([stuff1(1:12)  '  '  stuff2]);
   		title(['Rotated sonic components: u=blu; v=grn; w=orng; T=cyan. ' stuff1(7:12) ' '  stuff2(1:14)])  
   
   	end;%end plotit

	%select fft segment
	   nf=2^15;
   	while np<nf;
      	nf=nf/2;
   	end;   
		us=ur(1:nf);
		vs=vr(1:nf);
		ws=wr(1:nf);
		Ts=T(1:nf);
		usm=mean(us);%means for 2^15 point segment
		vsm=mean(vs);
		wsm=mean(ws);
		Tsm=mean(Ts);

	%do fft/power spectr/cospectra
		[su,F]=psd(us-usm,nf,10,hamming(nf),'none');
		[sv,F]=psd(vs-vsm,nf,10,hamming(nf),'none');
		[sw,F]=psd(ws-wsm,nf,10,hamming(nf),'none');
		[sT,F]=psd(Ts-Tsm,nf,10,hamming(nf),'none');
		[swu,F]=csd(ws-wsm,us-usm,nf,10,hamming(nf),'none');
		[swv,F]=csd(ws-wsm,vs-vsm,nf,10,hamming(nf),'none');
		[swT,F]=csd(ws-wsm,Ts-Tsm,nf,10,hamming(nf),'none');

	%normalize spectra properly (we think)
		su=su*2/10;sv=sv*2/10;sw=sw*2/10;sT=sT*2/10;swu=swu*2/10;swv=swv*2/10;swT=swT*2/10;

	
		n=nf/2;
		c=.1;
		jx=1;
		dx=1;
		ix=1;
		inx=1;
		while jx<n,
   		dx=(dx*exp(c));
   		d=floor(dx);
   		acu=0;acv=0;acw=0;acT=0;acwu=0;acwv=0;acwT=0;ac2=0;
   		k=0;
   		for jx=ix:ix+d-1;
   	   	if jx==n,
      	   	break;
      		end;
      		acu=acu+su(jx);
      		acv=acv+sv(jx);
      		acw=acw+sw(jx);
      		acT=acT+sT(jx);
      		acwu=acwu+swu(jx);
      		acwv=acwv+swv(jx);
      		acwT=acwT+swT(jx);
      		ac2=ac2+F(jx);
      		k=k+1;
   		end; %end for  jx
   		sus(inx)=acu/k;
   		svs(inx)=acv/k;
   		sws(inx)=acw/k;
   		sTs(inx)=acT/k;
   		cwus(inx)=acwu/k;
   		cwvs(inx)=acwv/k;
   		cwTs(inx)=acwT/k;
   		fs(inx)=ac2/k;
   		dfs(inx)=F(jx)-F(ix)+F(2);
   		ix=jx+1;
   		inx=inx+1; 
		end;%end while jx


		if plotit==1
   	%plot smoothed variance spectra
			figure;
			loglog(fs,sus.*fs,fs,svs.*fs,fs,sws.*fs,fs,sTs.*fs);
			xlabel('f (Hz)');
			ylabel('f*Sx(f)');
      	%title([stuff1(1:12)  '  '  stuff2]);
     	 title(['Variance spectra: u=blu; v=grn; w=orng; T=cyan. ' stuff1(7:12) ' '  stuff2(1:14)]) 

      %plot smoothed cospectra
			figure;
			semilogx(fs,cwus.*fs,fs,cwvs.*fs,fs,cwTs.*fs);
			xlabel('f (Hz)');
			ylabel('f*Cwx(f)');
      	%title([stuff1(1:12)  '  '  stuff2]);
      	title(['Cospectra: <wu>=blu; <wv>=grn; <wT>=orng. ' stuff1(7:12) ' '  stuff2(1:14)]) 

		end;%end while plot

		cvs(hrs,sn,1)=sum(cwus.*dfs);
		cvs(hrs,sn,2)=sum(cwvs.*dfs);
		cvs(hrs,sn,3)=sum(cwTs.*dfs);
		sigu(sn)=sqrt(sum(sus.*dfs));
		sigv(sn)=sqrt(sum(svs.*dfs));
		sigw(sn)=sqrt(sum(sws.*dfs));
		sigT(sn)=sqrt(sum(sTs.*dfs));
		sigs(hrs,sn,1)=sqrt(sum(sus.*dfs));
		sigs(hrs,sn,2)=sqrt(sum(svs.*dfs));
		sigs(hrs,sn,3)=sqrt(sum(sws.*dfs));
		sigs(hrs,sn,4)=sqrt(sum(sTs.*dfs));

		umen(sn)=usm;
		vmen(sn)=vsm;
		wmen(sn)=wsm;
		Tmen(sn)=Tsm;
		sp(sn)=sqrt(usm^2+vsm^2+wsm^2);
   	gfac=4*(2*pi/sp(sn))^.667;
   	jacx=inx-1;
   	jacy=jacx-4;
   
		cu(sn)=gfac*median(sus(jacy:jacx).*fs(jacy:jacx).^1.667);
		cv(sn)=gfac*median(svs(jacy:jacx).*fs(jacy:jacx).^1.667);
		cw(sn)=gfac*median(sws(jacy:jacx).*fs(jacy:jacx).^1.667);
		cT(sn)=gfac*median((sTs(jacy:jacx)-min(sTs)).*fs(jacy:jacx).^1.667);
		%cx(hrs,sn,1)=gfac*median(sus(65:69).*fs(65:69).^1.667);
		%cx(hrs,sn,2)=gfac*median(svs(65:69).*fs(65:69).^1.667);
		%cx(hrs,sn,3)=gfac*median(sws(65:69).*fs(65:69).^1.667);
		%cx(hrs,sn,4)=gfac*median((sTs(65:69)-min(sTs)).*fs(65:69).^1.667);
		ust(sn)=sqrt(abs(real(cvs(hrs,sn,1))));
		wt(sn)=real(cvs(hrs,sn,3));
		zsn(sn)=str2num(stuff1(7:10));
		L_mo(sn)=(Tsm+273.16)/9.83/.4/wt(sn)*ust(sn)^3;
	end;%end if flist loop
end; %end sn loop
[P,S]=polyfit(log(zsn),sp,1);
zpl=[1 20];
uplx=polyval(P,log(zpl));
figure;
semilogy(sp,zsn,'x',uplx,zpl);
axis([0 .7 1 20]);
xlabel('Wind Speed (m/s)');
ylabel('z (m)');
title([stuff2]);
axis([0 median(sp)+3 1 20]);
figure;
semilogy(ust,zsn,'x');
axis([0 .7 1 20]);
xlabel('Friction velocity, u* (m/s)');
ylabel('z (m)');
title([stuff2]);
figure;
semilogy(wt,zsn,'o');
axis([-.02 .02 1 20]);
xlabel('wt (K m/s)');
ylabel('z (m)');
title([stuff2]);
figure;
semilogy(sigu,zsn,'x',sigv,zsn,'.',sigw,zsn,'o');
axis([0 1.5 1 20]);
ylabel('z (m)');
xlabel('sigma_x (m/s)');
Title(['sig_u=x;  sig_v=dots; sig_w=o. 'stuff2(1:14)])
figure;
loglog(cu,zsn,'x',cv,zsn,'.',cw,zsn,'o');
axis([.01 1 1 20]);
ylabel('z (m)');
xlabel('C_x^2 (m/s)^2/m^{2/3}');
Title(['Cu^2=x;  Cv^2=dots; Cw^2=o. 'stuff2(1:14)])

toc;






	
	
	
	
	
	
	
	
	
	