disp('read_sonic_hr')
fclose('all');
clear all;

jd=input('julian day=ddd','s');
hr=input('gmt hr=hh','s');
sonic=input('sonic PORT#(2,3,4,5,6)','s');
%tic;
nx=10*(3600-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\98' jd hr '.p0'sonic];%note data found on I,D,and E also

flist=fopen(e,'r');
stuff1=fgets(flist);
stuff2=fgets(flist);
stuff3=fgets(flist);
disp(stuff1);
disp(stuff2);
zu=9;

%tic;
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
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);


%toc;
%tic;
t=.1:.1:nx/10;
%np=nx;

um=mean(u);
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);
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;
nt=10*60*10;%plot first ten minutes time series
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(['Rotated sonic components: u=blu; v=grn; w=orng; T=cyan. ' stuff1(7:12) ' '  stuff2(1:14)])  
%toc;
%tic;
%select fft segment
nf=2^15;
us=ur(1:nf);
vs=vr(1:nf);
ws=wr(1:nf);
Ts=T(1:nf);
usm=mean(us);
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;

%smooth spectra
tic;
if nf==2^16;
   sus=smth_mat*su;
	svs=smth_mat*sv;
	sws=smth_mat*sw;
	sTs=smth_mat*sT;
	cwus=smth_mat*swu;
	cwvs=smth_mat*swv;
	cwTs=smth_mat*swT;
else
   
	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;
   	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;
	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);
end;
toc;
%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(['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(['Cospectra: <wu>=blu; <wv>=grn; <wT>=orng. ' stuff1(7:12) ' '  stuff2(1:14)]) 

%save c:/users/cwf/matlab/sessions/30302_06 ur_6 vr_6 wr_6 T_6 sus_6 svs_6 sws_6 sTs_6 cwus_6 cwvs_6 cwTs_6	
	
	
	
	
	
	
	
	
	
	
	