disp('eval_ophir_sjohnson_05')
% 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;
clear t1 t2;
dfx=[dfx1 dfx2];

%   Open OPHIR file
dfl=['d:\data\day',dfx1,'\P4_',dfx];    
Tmean=avgTvais;
qmean=avgQvais;
tref=Tmean; %15;
Rgas=287.1;
rhoa=1010*100/(Rgas*(Tmean+273.16)*(1+0.61*qmean/1000));%air density
flist=fopen(dfl,'r');
headr=fgetl(flist);
headr1=fgetl(flist);
%disp('     ###############     Ophir Evaluation     ###############     ')
z=fscanf(flist,'%g,%g%g%g',[4,inf]);
if ~isempty(z);
   z=z(:,1:length(z)-1);
   t1=z(1,:)-z(1,1);
   Oclear=z(3,:);Oabsorb=z(2,:);
   %ophx=input('Ophir number?  [89111 = 1;   89113 = 2;   andreas = 3; 98157 = 4]   ');
   ophx=4;  %Used on Rico 
   tophir=t1/1000;
   if ophx==1,
      a0=.7403;
      a1=-.00172;
      ki=1/0.0399;
      bi=1/0.7848;
      disp('ophir 98160');%installed start epic2001 230
   elseif ophx==2,
      a0=.7813;
      a1=-.00118;
      ki=1/0.0423;
      bi=1/0.7648;
      disp('ophir 98161');%installed 119, taken down 126 1212
   elseif ophx==3,
      a0=0.8095;
      a1=-.00051;
      ki=1/0.047;
      bi=1/0.7543;
%      disp('ophir andreas');%installed for nauru........unit # 93149
   elseif ophx==4,
      a0=0.7447;
      a1=-.00150;
      ki=1/0.0416;
      bi=1/0.7710;
      disp('ophir 98157');%installed for weller 03
   end;%end ophx
   a=a0*(1+a1*(tref+15));%45=Tcan+15
   x=((tref+273.15)/273.15)^1.54*(ki)^bi;%3012.15=Ta+273.15
   if Oclear(1)<=0,
      Oclear(1)=10;
   end;
   if Oabsorb(1)<=0,
      Oabsorb(1)=10;
   end;
   for ii=2:length(Oclear);
      if Oclear(ii)<=0,
         Oclear(ii)=Oclear(ii-1);
      end
      if Oabsorb(ii)<=0,
         Oabsorb(ii)=Oabsorb(ii-1);
      end;
   end;%end for ii
   rc=Oabsorb./Oclear;
   abshum=x*(-log(rc/a)).^bi;
   abshum_mean=mean(abshum);

if isreal(abshum)==0,
   abshum=real(abshum);
end;

moh=median(abshum);
ii=find(abshum<0);abshum(ii)=moh;

fudg_oph=1.05;

spechum=abshum/rhoa*fudg_oph;
[Oph,F]=psd2(detrend(spechum),length(abshum),20);
[Ophs,Fso]=specsmoo(Oph,20);

licorit=1;
if licorit
    read_co2_rhb_04
    spechum_l=Licor_H2O_a;
    [Lic,F]=psd2(detrend(spechum_l),length(spechum_l),20);
    [Lics,Fsl]=specsmoo(Lic,20);
end;

if hhh==12,
   specfile=['d:\data\spec_files\q_sp',num2str(ddd),num2str(hhh),'.txt'];
   fspec=fopen(specfile,'w');
   if licorit
    specs=[Fso Ophs Fsl Lics];
    fprintf(fspec,'%11.5f %11.5f %11.5f %11.5f \r\n',specs);
   else
    specs=[Fso Ophs];
    fprintf(fspec,'%11.5f %11.5f \r\n',specs);
   end;

   fclose(fspec);
   clear fspec specs specfile
end;
if plot_f & hhh==12
    if licorit
        figure;loglog(Fsl,Lics.*Fsl,Fso,Ophs.*Fso,Fso(length(Fso)-31:length(Fso)),((Fso(length(Fso)-31:length(Fso))).^(-2/3))./(Fso(length(Fso)-18).^(-2/3)).*Ophs(length(Fso)-18).*Fso(length(Fso)-18));
    else
        figure;loglog(Fso,Ophs.*Fso,Fso(length(Fso)-31:length(Fso)),((Fso(length(Fso)-31:length(Fso))).^(-2/3))./(Fso(length(Fso)-18).^(-2/3)).*Ophs(length(Fso)-18).*Fso(length(Fso)-18));
    end;

    title(['Hygrometer Spectrum:   Day ',num2str(ddd),'   Hour ',num2str(hhh)]);
    xlabel('Frequency in Hz');
    ylabel('Spectrum * Frequency (g/kg)^2)');

    if licorit;if length(Fso)==length(Fsl);figure;semilogx(Fso,sqrt(Ophs./Lics),[1e-4 10],[1 1]);end;end;
end;

end;           %if ~isempty z