disp('nau_ophir')
% 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 OPHIR file
dfl=['D:\hp_data\day',num2str(ddd),'\P4_',dfx];       

flist=fopen(dfl,'r');
headr=fgetl(flist)
headr1=fgetl(flist);
%disp('     ###############     Ophir Evaluation     ###############     ')
z=fscanf(flist,'%g,%g,%g%g%g',[5,inf]);
if ~isempty(z);
   z=z(:,1:length(z)-1);
   t1=z(1,:)-z(1,1);t2=z(2,:)/10000-z(2,1)/10000;
   Oclear=z(4,:);Oabsorb=z(3,:);
%   ophx=input('Ophir number?  [89111 = 1;   89113 = 2;   andreas = 3]   ');
   ophx=3;
   tophir=t1+t2;clear t1 t2;
   if ophx==1,
      a0=1.039;
      a1=-.001364;
      ki=1/0.05215;
      bi=1/0.7288;
      disp('ophir 89111');%installed 126 1212
   elseif ophx==2,
      a0=.7643;
      a1=-.001664;
      ki=1/0.05161;
      bi=1/0.74;
      disp('ophir 89113');%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
   end;
   a=a0*(1+a1*45);%45=Tcan+15
   x=(302.15/273.15)^1.54*(ki)^bi;%3012.15=Ta+273.15
   
   mm=1:12000;

   for pp=1:6;
      ii=mm+12000*(pp-1);
      if pp==6;
         ii=length(Oclear)-12000+1:length(Oclear);
      end;
      Oclear_10=detrend(Oclear(ii));Oabsorb_10=detrend(Oabsorb(ii));
      med_Oabsorb=median(Oabsorb_10);med_Oclear=median(Oclear_10);
      sort_Oabsorb=sort(Oabsorb_10);sort_Oclear=sort(Oclear_10);
      jk=floor(0.5*length(sort_Oabsorb));
      jkp=floor(.667*length(sort_Oabsorb));
		jkm=floor(.333*length(sort_Oabsorb));
      
      sig_Oabsorb=0.5*(sort_Oabsorb(jkp)-sort_Oabsorb(jkm));
      sig_Oclear=0.5*(sort_Oclear(jkp)-sort_Oclear(jkm));
      if abs(Oclear_10(1)-med_Oclear) > 7*sig_Oclear,
         Oclear_10(1)=med_Oclear;
      end;
      if abs(Oabsorb_10(1)-med_Oabsorb) > 7*sig_Oabsorb,
         Oabsorb_10(1)=med_Oabsorb;
      end;
      for kk=2:length(Oclear_10);
         if abs(Oclear_10(kk)-med_Oclear) > 7*sig_Oclear,
            Oclear_10(kk)=Oclear_10(kk-1);
         end;   
         if abs(Oabsorb_10(kk)-med_Oabsorb) > 7*sig_Oabsorb,
            Oabsorb_10(kk)=Oabsorb_10(kk-1);
         end;
      end;

      rc=(Oabsorb_10+median(Oabsorb(ii)))./(Oclear_10+median(Oclear(ii)));
      abshum=x*(-log(rc/a)).^bi;

      Ophir_stat_clear=std(Oclear_10);
      Ophir_stat_absorb=std(Oabsorb_10);
      abshum_mean=mean(abshum);

      disp(['Standard Deviation ',num2str(pp),'/6    Clear = ',num2str(Ophir_stat_clear),'    Absorb = ',num2str(Ophir_stat_absorb),'   Mean abs hum = ',num2str(abshum_mean),'  [g/m/m/m]'])
      disp(['     Median Raw Counts    Clear = ',num2str(median(Oclear(ii))),'     Absorb = ',num2str(median(Oabsorb(ii)))])
	end;
   rc=Oabsorb./Oclear;
   abshum=x*(-log(rc/a)).^bi;
else,
   abshum=(1:72000).*0;
end;

specsizeo=length(abshum);
abm=mean(abshum);
Ch2O=psd(abshum-abm,specsizeo);[Ch2O1,FFo]=specsmooth(Ch2O,20);

figure;plot(tophir,abshum);title(['JD',dfx1,' Hr',dfx2,' UTC  Ophir abshum']);xlabel('time in sec');ylabel('abshum');
figure;loglog(FFo,Ch2O1.*FFo);title(['JD',dfx1,'   Hr',dfx2,'  Ophir spectrum']);xlabel('freq in Hz');ylabel('F*S_a_b_s_h_u_m');

clear Ch Ch1 FF z x tophir specsizeo sort_Oclear sort_Oabsorb sig_Oclear sig_Oabsorb pp ophx mm
clear med_Oclear med_Oabsorb kk ki jkp jkm jk ii headr1 headr flist Ch2O Ch2O1 DD FFo Oabsorb_10
clear Oclear_10 Ophir_stat_absorb Ophir_stat_clear a a0 a1 abm abshum_mean ans bi dfl dfx dfx1 dfx2 rc