disp('read_sonics_day_mod4')
fclose('all');
clear all;
sav_spec=1;%save spectra =1
for jdy=107:124, %loop through all days on the CD
jd=num2str(jdy);
%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;
blk=1;
hr='00';
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 jam=1:24,%cycle thru 24 hourly stats files
   if jam<11,
      hr=['0'num2str(jam-1)];
   else
      hr=num2str(jam-1);
   end;
   hrs=str2num(hr)+1;

	for sn=1:5;%read thru each sonic
   	sonic=num2str(sn+1);
      zsn(hrs,sn)=-sn;
      sp(hrs,sn)=999.;
      dir(hrs,sn)=999.;
      Tmen(hrs,sn)=999.;
      wmen(hrs,sn)=9.99;
      ust(hrs,sn)=9.99;
      wt(hrs,sn)=999.;
      L_mo(hrs,sn)=999.;
      sigu(hrs,sn)=9.99;
      sigv(hrs,sn)=9.99;
      sigw(hrs,sn)=9.99;
      sigT(hrs,sn)=9.99;
      cu(hrs,sn)=9.99;
      cv(hrs,sn)=9.99;
      cw(hrs,sn)=9.99;
      cT(hrs,sn)=9.99;
      fl(hrs,sn)=0; %good(0)/bad(1) data flag
      npts(hrs,sn)=0;
      np=0;
      ur=zeros(1,36000);
      vr=zeros(1,36000);
      wr=zeros(1,36000);
      Tr=zeros(1,36000);
      T=zeros(1,36000);
      stux=zeros(4,36000);
   
jdy
if jdy>=107 & jdy<=109 
  e=['D:\DasWin32\Data\ascii\98' jd hr '.p0'sonic]
elseif jdy>= 110 & jdy<=112
  e=['E:\DasWin32\Data\ascii\98' jd hr '.p0'sonic]
elseif jdy>=113 & jdy<=115
  e=['F:\Users\pguest\98' jd hr '.p0'sonic]
elseif jdy>=116 & jdy<= 118
  e=['G:\Users\pguest\98' jd hr '.p0'sonic]
elseif jdy >=119 & jdy <=121
  e=['H:\Users\pguest\98' jd hr '.p0'sonic]
elseif jdy >= 122 & jdy <=124
  e=['i:\DasWin32\Data\ascii\98' jd hr '.p0'sonic]
end;
	   flist=fopen(e,'r');
      if flist>0,   %if the file exists
			
                   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)
         nfw=2^13;
         if np<2*nfw, %if there is a break early in the hour, rescan data after break
                   stux=zeros(4,36000);
                   stuff3=fgets(flist);%read 3 headers
	     	
                   stuff4=fgets(flist);
	      
                   stuff5=fgets(flist);
                   stuff6=fgets(flist);
	           disp(stuff3);
		   disp(stuff4);
                   disp(stuff5);
                   disp(stuff6);
                   stux=fscanf(flist,'%*1s%g%*1s%g%*1s%g%*1s%g',[4,36000]);
                   [a,np2]=size(stux)
                   np=max([np np2]);
         end;
         if np<2*nfw, %if there is a second break early in the hour, rescan data after break
                   stux=zeros(4,36000);
                   stuff3=fgets(flist);%read 3 headers
	     	
                   stuff4=fgets(flist);
	      
                   stuff5=fgets(flist);
                   stuff6=fgets(flist);
	           disp(stuff3);
		   disp(stuff4);
                   disp(stuff5);
                   disp(stuff6);
                   stux=fscanf(flist,'%*1s%g%*1s%g%*1s%g%*1s%g',[4,36000]);
                   [a,np3]=size(stux)
                   np=max([np np3]);
         end;
        if np>=2*nfw,%if there is at least 27 min of continuous data, process fluxes
		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=zeros(np,2);
	    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);
            dirx=atan2(-vm,um)*180/pi;
            if (dirx<0) 
            	dirx=dirx+360;
            end
	%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-use only the first nf points in time series
	%nf=2^15;
        novrlp=nfw/2;
        nfinc=nfw-novrlp;
        mlt=(np-nfw)/nfinc;
        mlt=mlt-mod(mlt,1);
        nf=mlt*nfinc+nfw
%        if np<2^15
%          nf=2*nfw;
%        else
%          nf=2^15;
%        end;
        %nf=np;

		us=ur(1:nf);
		vs=vr(1:nf);
		ws=wr(1:nf);
		Ts=T(1:nf);
		usm=mean(us);%means for nf-point segment
		vsm=mean(vs);
		wsm=mean(ws);
		Tsm=mean(Ts);

	%do fft/power spectr/cospectra
		[su,F]=psd(us,nfw,10,hamming(nfw),novrlp,'linear');
		[sv,F]=psd(vs,nfw,10,hamming(nfw),novrlp,'linear');
		[sw,F]=psd(ws,nfw,10,hamming(nfw),novrlp,'linear');
		[sT,F]=psd(Ts,nfw,10,hamming(nfw),novrlp,'linear');
		[swu,F]=csd(ws,us,nfw,10,hamming(nfw),novrlp,'linear');
		[swv,F]=csd(ws,vs,nfw,10,hamming(nfw),novrlp,'linear');
		[swT,F]=csd(ws,Ts,nfw,10,hamming(nfw),novrlp,'linear');

	%normalize spectra properly (we think)
		%fct=norm(hamming(nfw))^2/sum(hamming(nfw))^2;
      	  	fct=1;
		su=fct*su*2/10;sv=fct*sv*2/10;sw=fct*sw*2/10;sT=fct*sT*2/10;
       		swu=fct*swu*2/10;swv=fct*swv*2/10;swT=fct*swT*2/10;

		n=nfw/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 sav_spec==1; % save spectra to file
                  clear spec
                  g=['c:\Sheba\data\Summer_98\son_spec' jd '.txt'];
                  gidt=fopen(g,'a');
                  jdate=str2num(jd)+(jam-1)/24;                    
                  datcol(1,1)=jdate;
                  datcol(2,1)=jdy;
                  datcol(3,1)=jam-1;
                  datcol(4,1)=str2num(stuff1(7:10));
                  datcol(5,1)=sn;
                  datcol(6,1)=nf;
                  datcol(7,1)=inx-1;
                  fprintf(gidt,'%9.4f %4i %3i %5.2f %3i %7i %4i\r',datcol);                  
                  for spct=1:inx-1
                    spec(1,spct)=fs(spct);
                    spec(2,spct)=sus(spct)*fs(spct);
                    spec(3,spct)=svs(spct)*fs(spct);
                    spec(4,spct)=sws(spct)*fs(spct);
                    spec(5,spct)=sTs(spct)*fs(spct);
                    spec(6,spct)=cwus(spct)*fs(spct);
                    spec(7,spct)=cwvs(spct)*fs(spct);
                    spec(8,spct)=cwTs(spct)*fs(spct);
                  end;%for spct
                %vectorized print, 8 columns, spct rows
                  fprintf(gidt,'%7.4f %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e\r',spec);
                  clear spec datcol
                  fclose(gidt)
                end; %end sav_spec           

		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 if plot

		cwux(hrs,sn)=sum(cwus(1:60).*dfs(1:60));
		cwvx(hrs,sn)=sum(cwvs(1:60).*dfs(1:60));
		cwTx(hrs,sn)=sum(cwTs(1:60).*dfs(1:60));
		sigu(hrs,sn)=sqrt(sum(sus.*dfs));
		sigv(hrs,sn)=sqrt(sum(svs.*dfs));
		sigw(hrs,sn)=sqrt(sum(sws.*dfs));
        	sigT(hrs,sn)=sqrt(sum(sTs.*dfs));
         	ssm=sqrt(usm*usm+vsm*vsm);
			
 		if sigu(hrs,sn)/ssm > 0.4 | sigv(hrs,sn)/ssm > 0.4 
            		fl(hrs,sn)=1
         	end ; %end test for high std dev.
         	umen(hrs,sn)=usm;
		vmen(hrs,sn)=vsm;
		wmen(hrs,sn)=wsm;
		Tmen(hrs,sn)=Tsm;
         	dir(hrs,sn)=dirx;
		sp(hrs,sn)=sqrt(usm^2+vsm^2+wsm^2);
		gfac=4*(2*pi/sp(hrs,sn))^.667;
		cu(hrs,sn)=gfac*median(sus(51:55).*fs(51:55).^1.667);
		cv(hrs,sn)=gfac*median(svs(51:55).*fs(51:55).^1.667);
		cw(hrs,sn)=gfac*median(sws(51:55).*fs(51:55).^1.667);
		cT(hrs,sn)=gfac*median((sTs(51:55)-min(sTs)).*fs(51:55).^1.667);
			
		ust(hrs,sn)=-sign(cwux(hrs,sn))*sqrt(abs(real(cwux(hrs,sn))));
		rho=100*1013/((273.15+Tsm)*287);
        	wt(hrs,sn)=1004*rho*real(cwTx(hrs,sn));
		wtx(hrs,sn)=real(cwTx(hrs,sn));
               zsn(hrs,sn)=str2num(stuff1(7:10));
               if wtx(hrs,sn)==0,
                  L_mo(hrs,sn)=-999
               else
                  L_mo(hrs,sn)=-(Tsm+273.15)/9.83/.4/wtx(hrs,sn)*ust(hrs,sn)^3;
               end;
	 else   %conditional on np (array size)
           zsn(hrs,sn)=str2num(stuff1(7:10));
           fl(hrs,sn)=1;
         end;%end if np     
       else %flist conditional
         fl(hrs,sn)=1;
       end;%end if flist loop
       npts(hrs,sn)=np
     end; %end sn loop
     if plotit==1,
   	[P,S]=polyfit(log(zsn(hrs,1:sn)),sp(hrs,1:sn),1);
		zpl=[1 20];	
		uplx=polyval(P,log(zpl));
		figure;
		semilogy(sp(hrs,1:sn),zsn(hrs,1:sn),'x',uplx,zpl);
		axis([0 .7 1 20]);
		xlabel('Wind Speed (m/s)');
		ylabel('z (m)');
		title([stuff2]);
		axis([0 median(sp(hrs,1:sn))+3 1 20]);
		figure;
		semilogy(ust(hrs,1:sn),zsn(hrs,1:sn),'x');
		axis([0 .7 1 20]);
		xlabel('Friction velocity, u* (m/s)');	
		ylabel('z (m)');
		title([stuff2]);
		figure;
		semilogy(wt(hrs,1:sn),zsn(hrs,1:sn),'o');
		axis([-.02 .02 1 20]);
		xlabel('wt (K m/s)');
		ylabel('z (m)');
		title([stuff2]);
		figure;
		semilogy(sigu(hrs,1:sn),zsn(hrs,1:sn),'x',sigv(hrs,1:sn),zsn(hrs,1:sn),'.',sigw(hrs,1:sn),zsn(hrs,1:sn),'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)])
   end;%end plotit
   jz(jam)=str2num(jd)+(jam-1)/24;
end; %end jam loop
%The following takes data computed by the read_sonics_hr progam and creates an ascii
%file appropriate for reading into a spreadsheet (e.g., lotus, quatropro, etc)

prt_it=1;  %set to 0 to print to screen, 1 to print to a file (i.e., make ascii file)
if prt_it==1
%	f='D:\matlabstf\cwf\sheba_proc\proc_son3.txt'
%      f='E:\Users\opersson\proc_son_sum.txt'
      f='c:\Sheba\data\Summer_98\proc_son_sp_e.txt';
	fidt=fopen(f,'a');
else
   fidt=1;
end; 
%take sonic data and stick into matrix
for j=1:24; %loop over hours
for sn=1:5; %print out each sonic
jaz(1,sn)=jz(j);%Julian Day
jaz(2,sn)=j-1;%hour
jaz(3,sn)=zsn(j,sn);%instrument heights
jaz(4,sn)=min(99.9,sp(j,sn));%streamwise winds
jaz(5,sn)=min(999,dir(j,sn));%wind direction
jaz(6,sn)=sign(Tmen(j,sn))*min(99.9,abs(Tmen(j,sn)));%mean temperature
jaz(7,sn)=sign(wmen(j,sn))*min(9.99,abs(wmen(j,sn)));%mean vertical velocity
jaz(8,sn)=sign(ust(j,sn))*min(9.99,abs(ust(j,sn)));%friction velocity
jaz(9,sn)=sign(wt(j,sn))*min(999,abs(wt(j,sn)));%heat flux
jaz(10,sn)=sign(L_mo(j,sn))*min(999,abs(L_mo(j,sn)));%Monin-Obuhkov length
jaz(11,sn)=min(99.9,sigu(j,sn));%u variance
jaz(12,sn)=min(99.9,sigv(j,sn));%v variance
jaz(13,sn)=min(99.9,sigw(j,sn));%w variance
jaz(14,sn)=min(99.9,sigT(j,sn));%T variance
jaz(15,sn)=min(99.9,cu(j,sn));%u structure function
jaz(16,sn)=min(99.9,cv(j,sn));%v structure function
jaz(17,sn)=min(99.9,cw(j,sn));%w structure function
jaz(18,sn)=min(99.9,cT(j,sn));%T structure function
jaz(19,sn)=fl(j,sn);%good/bad data flag
jaz(20,sn)=npts(j,sn);
end; %end instrument loop through sonics
%vectorized print, 19 columns, 5 rows
fprintf(fidt,'%9.4f %4i %7.2f %6.1f %5.0f %7.2f %9.4f %7.3f %7.1f %8.1f %7.2f %7.2f %7.2f %7.2f %7.3f %7.3f %7.3f %7.3f %3.0f %7i\r',jaz);
clear jaz
end; %end hour loop
fclose('all')
clear jz zsn sp dir Tmen wmen ust wt L_mo sigu sigv sigw sigT cu cv cw cT fl npts 
end; %end jdy loop

toc;
