disp('read_ophir_hr')
fclose('all');
jd=input('julian day','s');
hr=input('gmt hr','s');

tic;
nx=10*(7200-1);

ar=zeros(1,72000);
cr=zeros(1,72000);
er=zeros(1,72000);
stux=zeros(3,72000);
%	e=['D:\DasWin32\data\ascii\97' jd hr '.p08'];
%	e=['E:\users\data\97' jd hr '.p08'];
	e=['K:\DasWin32\data\97' jd hr '.p08'];

flist=fopen(e,'r');
stuff1=fgets(flist);
ophir=fscanf(flist,'%*s%*s%*s%*s%g%g%g%g%g%g%g%g%g%g%g%g%*s%g%g%g%g%g%g%g%g',[20,1]);%stuff2=fgets(flist);
mo=ophir(12);
no=ophir(11);
ko=ophir(7);
bo_inv=1/ophir(6);
mo=ophir(12);
Ao=ophir(5);
A1=ophir(8);
rhq=ophir(14);
Rcor=ophir(17);
Tair=ophir(18);
Tcan=ophir(19);
qxq=rhq*((Tair+273.16)/273.16)^mo/(-log(Rcor/Ao))^bo_inv;
afac=Ao*(1+A1*Tcan);



disp(stuff1);
%disp(stuff2);
stuff4=fgets(flist);
stuff5=fgets(flist);
zu=9;

tic;
stux=fscanf(flist,'%*c%g%g%g',[3,72000]);
[a,npa]=size(stux);
np=floor(npa/2);
np2=2*np;
jam=stux(1:2,1:np2)';
jamx=reshape(jam,[2,np2]);
ah=jamx(1,1:np2);
bh=jamx(1,1:np2);
ch=(ah+bh)/2;
jamy=reshape(ch,[np,2]);
araw=jamy(1:np,1)';
craw=jamy(1:np,2)';

%despike
on=1:1:np;
uz=zeros(np,2);

uz(1:np,1)=araw';
uz(1:np,2)=on';
n_a=despik(uz);
uz(1:np,1)=craw';
n_c=despik(uz);

rcx=n_a./n_c;
rcx=rcx/afac;
qcx=-log(rcx);
qcxi=power(qcx,bo_inv);
rhow=qcxi*qxq;




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


rhom=mean(rhow);


nt=10*60*30;%plot first ten minutes time series
figure;
plot(t(1:nt),rhow(1:nt));
xlabel('time (sec)');
ylabel('Q (g/m^3)');

toc;
tic;
%select fft segment
nf=2^15;
rhos=rhow(1:nf);
rhosm=mean(rhos);

%do fft/power spectr/cospectra
[sq,F]=psd(rhos-rhosm,nf,10,hamming(nf),'none');
%normalize spectra properly (we think)
sq=sq*2/10;
%smooth spectra
n=nf/2;
c=.1;
jx=1;
dx=1;
ix=1;
inx=1;
while jx<n,
   dx=(dx*exp(c));
   d=floor(dx);
   acq=0;ac2=0;k=0;
   for jx=ix:ix+d-1;
      if jx==n,
         break;
      end;
      acq=acq+sq(jx);
      ac2=ac2+F(jx);
      k=k+1;
   end;
   sqs(inx)=acq/k;
   fs(inx)=ac2/k;
   dfs(inx)=F(jx)-F(ix)+F(2);
   ix=jx+1;
   inx=inx+1; 
end;
sqs(inx)=acq/k;
fs(inx)=ac2/k;  
dfs(inx)=F(jx)-F(ix)+F(2);
toc;
%plot smoothed variance spectra
figure;
loglog(fs,sqs.*fs);
xlabel('f (Hz)');
ylabel('f*Sq(f) (g/m^3)^2');

%
%save c:/users/cwf/matlab/sessions/30302_06 rhow sqs_6 	
	
	
	
	
	
	
	
	
	
	
	