%disp('opensons')
% Inputs data from ASCII data files
%  Note:  hhh=hour

fclose('all'); 
if hhh<10,
   dfx2=['0',num2str(hhh)];
else,
   dfx2=num2str(hhh);
end;

if ddd<10,
   dfx1=['00',num2str(ddd)];
elseif ddd<100,
   dfx1=['0',num2str(ddd)];
else,
   dfx1=num2str(ddd);
end;

dfx=[dfx1 dfx2];

%  Open Sonic file		
%dfl=['/home/user1/bin/data/P0','_',dfx];%disp(dfl)
%dfl=['/scratch2/day',dfx1,'/P0','_',dfx];%disp(dfl)
%dfl=['/home/user1/bin/data_backup/day',dfx1,'/P0','_',dfx];%disp(dfl)
%dfl=['/home/jeh/fastex/P0','_',dfx,';1'];%disp(dfl)
%dfl=['/home/jeh/fastex/P0_',dfx,';1'];
dfl=['D:\hp_data\P0_',dfx];       

clear dfx
flist=fopen(dfl,'r');
headr=fgetl(flist);%disp(headr)
%headr1=fgetl(flist);
%disp('Open Sonic File')

z=fscanf(flist,'%g%*c%g%*c%g%*c%g%*c%g%*c%g%*c',[6,inf]);
%                1    2    3    4    5    6   
z=z(:,1:length(z)-1);%disp(['length of sonic file is ',num2str(length(z))])
t1=z(1,:)-z(1,1);t2=z(2,:)/10000-z(2,1)/10000;
tsons=t1+t2;clear t1 t2;
tsons=tsons-tsons(1);tsons=0:1/20.8325:(length(tsons))/20.8325;tsons=tsons';

u=z(3,:)/100;v=z(4,:)/100;w=z(5,:)/100;Tsonic=z(6,:)/50;%clear z;

Tsonic=(Tsonic.*Tsonic+(u.*u + w.*w)/2. + v.*v)/403-273.15;

sx=-u*cos(30/180*pi)-v*sin(30/180*pi);     %to North
sy=-u*sin(30/180*pi)+v*cos(30/180*pi);     %to West

%Un=mean(sx).*cos(Mhead)-mean(-sy).*sin(Mhead)-Msog.*cos(Mcog);   %to north
%Ue=mean(sx).*sin(Mhead)+mean(-sy).*cos(Mhead)-Msog.*sin(Mcog);   %to east
%Twindir=atan2(-Ue,-Un).*180/pi;if Twindir<0,Twindir=Twindir+360;end;

%ssx=(1:length(tt)).*0;ssx(1)=u(1);
%ssy=(1:length(tt)).*0;ssy(1)=v(1);
%ssz=(1:length(tt)).*0;ssz(1)=w(1);
%sst=(1:length(tt)).*0;sst(1)=Tsonic(1);
%jj=2;
%for ii=2:length(tt);
%   while tb(jj) < tt(ii) & jj < length(u);
%      jj=jj+1;
%   end;
%   ssx(ii)=(u(jj)+u(jj-1))/2;
%   ssy(ii)=(v(jj)+v(jj-1))/2;
%   ssz(ii)=(w(jj)+w(jj-1))/2;
%   sst(ii)=(Tsonic(jj)+Tsonic(jj-1))/2;
%end;
%clear u v w  ii jj
u=sx;v=sy;%Tsonic=sst;
clear sx sy sz ;

if length(u) > 65535,
   specsize=65536;
elseif length(u) > 32767,
   specsize=32768;   
elseif length(u) > 16383,
   specsize=16384;   
elseif length(u) > 8191,
   specsize=8192;   
elseif length(u) > 4095,
   specsize=4096;   
elseif length(u) > 2047,
   specsize=2048;   
elseif length(u) > 1023,
   specsize=1024;   
elseif length(u) > 511,
   specsize=512;   
elseif length(u) > 255,
   specsize=256;   
elseif length(u) > 127,
   specsize=128;   
elseif length(u) > 63,
   specsize=64;   
elseif length(u) > 31,
   specsize=32;   
elseif length(u) > 15,
   specsize=16;   
end;
DD=20.832;
tm=mean(Tsonic);
[Cts,F]=psd(Tsonic-tm,specsize,20.83,hamming(specsize),'none');[Cts1,FF]=specsmooth(Cts,DD);clear Cts

%su=sx;sv=sy;sw=w;clear u v w;
%tsons=tb;clear tb;
mm=1:18000;
jv=length(Tsonic);
figure;plot(tsons(1:length(Tsonic)),Tsonic);title(['JD',dfx1,' Hr',dfx2,' UTC  Sonic Temperature']);ylabel('deg');xlabel('time in sec');
figure;loglog(FF,Cts1.*FF);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectrum Sonic Temperature']);ylabel('F*S(x,y,z)');xlabel('Freq in Hz');clear Cts1
jx=min(jv,length(u));
um=mean(u);vm=mean(v);wm=mean(w);
[Cus,F]=psd(u-um,specsize,20.83,hamming(specsize),'none');[Cus1,FF]=specsmooth(Cus,DD);clear Cus;
[Cvs,F]=psd(v-vm,specsize,20.83,hamming(specsize),'none');[Cvs1,FF]=specsmooth(Cvs,DD);clear Cvs;
[Cws,F]=psd(w-wm,specsize,20.83,hamming(specsize),'none');[Cws1,FF]=specsmooth(Cws,DD);clear Cws;

figure;plot(tsons(1:jx),u(1:jx));title(['JD',dfx1,' Hr',dfx2,' UTC  U wind component']);ylabel('m/s');xlabel('time in sec');
figure;loglog(FF,Cus1.*FF,FF,Cvs1.*FF,FF,Cws1.*FF);title(['JD',dfx1,' Hr',dfx2,' UTC  Uncorrected U,V, W Spectra']);ylabel('F*S_{xx} (m/s)^2/Hz');xlabel('Freq (Hz)');
%for pp=1:4;
%   ii=mm+18000*(pp-1);
%   if pp==4;
%      ii=length(su)-18000+1:length(su);
%   end;
%   u=su(ii);
%   v=sv(ii);
%   w=sw(ii);
%   tb=tsons(ii);
%   Ts=Tsonic(ii);
%   rotcorr;
%   figure;plot(tb,u,tb,v,tb,w);title(['JD',dfx1,' Hr',dfx2,' UTC  Raw Sonic Velocities']);ylabel('m/s');xlabel('time in sec');
%   figure;plot(tb,Ustream,tb,Vstream,tb,Wstream);title(['JD',dfx1,' Hr',dfx2,' UTC  Rotated Sonic Velocities']);ylabel('m/s');xlabel('time in sec');

%   Cu=psd(Ustream,specsize);[Cu1,FF]=specsmooth(Cu);clear Cu
%   Cv=psd(Vstream,specsize);[Cv1,FF]=specsmooth(Cv);clear Cv
%   Cw=psd(Wstream,specsize);[Cw1,FF]=specsmooth(Cw);clear Cw
%   figure;loglog(FF,Cu1.*FF,FF,Cv1.*FF,FF,Cw1.*FF);title(['JD',dfx1,' Hr',dfx2,' #',num2str(pp),' UTC  Uncorrected Sonic Velocities']);ylabel('F*S(x,y,z)');xlabel('Freq in Hz');clear Cu1 Cv1 Cw1
%   disp(['size of FF =',num2str(size(FF))])
%   clear alpha beta Sb Ub Vb Wb FF
%   uu=cov(Ustream);disp(['uu = ',num2str(uu),'   streamwise wind   variance']),clear uu
%   vv=cov(Vstream);disp(['vv = ',num2str(vv),'   transverse wind   variance']),clear vv
%   ww=cov(Wstream);disp(['ww = ',num2str(ww),'   vertical velocity variance']),clear ww
%   tts=cov(Ts);disp(['tts = ',num2str(tts),'   sonic temperature variance']),clear tts

%   disp('')
%   uv=cov(Ustream,Vstream);uv=uv(1,2);disp(['uv = ',num2str(uv),'   U-V   covariance']),clear uv
%   uw=cov(Ustream,Wstream);uw=uw(1,2);disp(['uw = ',num2str(uw),'   U-W   covariance']),clear uw
   %ut=cov(Ustream,Tsonic); ut=ut(1,2);disp(['ut = ',num2str(ut),'   U-T   covariance']),clear ut
%   vw=cov(Vstream,Wstream);vw=vw(1,2);disp(['vw = ',num2str(vw),'   V-W   covariance']),clear vw
   %vt=cov(Vstream,Tsonic); vt=vt(1,2);disp(['vt = ',num2str(vt),'   V-T   covariance']),clear vt
%   wt=cov(Wstream,Ts); wt=wt(1,2);disp(['wt = ',num2str(wt),'   W-T   covariance']),clear wt
   %Cuw=csd(Ustream,Wstream,specsize);[Cuw1,FF]=specsmooth(Cuw);clear Cuw;qx=Cuw1.*FF;figure;semilogx(FF(5:length(FF)),real(qx(5:length(FF))),FF(5:length(FF)),imag(qx(5:length(FF))));title(['JD',dfx1,' Hr',dfx2,' UTC  Co- & Qu-Spectra U-W']);ylabel('F*Cuw');xlabel('Freq in Hz');clear Cuw1
   %Cvw=csd(Vstream,Wstream,specsize);[Cvw1,FF]=specsmooth(Cvw);clear Cvw;qx=Cvw1.*FF;figure;semilogx(FF(5:length(FF)),real(qx(5:length(FF))),FF(5:length(FF)),imag(qx(5:length(FF))));title(['JD',dfx1,' Hr',dfx2,' UTC  Co- & Qu-Spectra V-W']);ylabel('F*Cvw');xlabel('Freq in Hz');clear Cvw1
   %Cwt=csd(Wstream,Ts,specsize);[Cwt1,FF]=specsmooth(Cwt);clear Cwt; qx=Cwt1.*FF;figure;semilogx(FF(5:length(FF)),real(qx(5:length(FF))),FF(5:length(FF)),imag(qx(5:length(FF))));title(['JD',dfx1,' Hr',dfx2,' UTC  Co- & Qu-Spectra W-Tsonic']);ylabel('F*Cwt');xlabel('Freq in Hz');clear Cwt1

%   clear x vtemp tcan lv ktemp headr flist a ts 
%end;

%clear specsize u v w tb pp mm ii ans Wstream Ts Vstream Ustream dfl