disp('ICEALOT_allcruise_mb');
%process mailbox .raw files
clear all;

s1=['m:\*.raw'];
s='m:\';
%s1='F:\AMMA_2007\RHB\radiometer\Mailbox\Raw\*.raw';
%s='F:\AMMA_2007\RHB\radiometer\Mailbox\Raw\';
dr=dir(s1);
[n m]=size(dr);
mx=n;
j=1;
nj=1;
kx=1;					  % counts number of 10-min blocks of data
min=0:6;				  % 10-min parts of an hour
jdmx=min/6/24;			  % 10-min fraction of hour
%mx=2*24;
for i=1:mx-1 %don't try to read the current active file
    k=1;
    clear dath;
    flnm=dr(i).name

    %    if flnm(10:12)=='raw';
    fid=fopen([s flnm],'r');

    while feof(fid)==0 ;
        %                   1   2   3   4   5  6    7  8

        stux1=fscanf(fid,'%02d%02d%02d%02d%02d%02d %d %f %f %f',[10,1]);
        stux2=fscanf(fid,'%g %g %g %g %g %g %g %g %g',[9,1]);
        stux3=fscanf(fid,'%g %g %g %g %g %g %g %g %g',[9,1]);
        if (isempty(stux1))
            break
        else
            datx(j,1)=datenum(2000+stux1(1),stux1(2),stux1(3),stux1(4),stux1(5),stux1(6))-datenum(2000+stux1(1)-1,12,31);
            if (length(stux2)==9) &(length(stux3)==9)
                datx(j,2)=stux2(9);%vapor
                datx(j,3)=stux3(9);%liquid
                datx(j,4)=stux2(6);%tau 20
                datx(j,5)=stux2(8);%tb 20
                datx(j,6)=stux3(6);%tau 30
                datx(j,7)=stux3(8);%tb 30
                datx(j,8)=stux1(7);%mode
                datx(j,9)=stux1(8);%angle
                datx(j,10)=stux2(7);%cf20
                datx(j,11)=stux3(7);%cf30
                datx(j,12)=stux2(1);%sky20
                datx(j,13)=stux2(2);%nsky20
                datx(j,14)=stux2(3);%bb20
                datx(j,15)=stux2(4);%nbb20
                datx(j,16)=stux2(5);%tbb20
                datx(j,17)=stux3(1);%sky30
                datx(j,18)=stux3(2);%nsky30
                datx(j,19)=stux3(3);%bb30
                datx(j,20)=stux3(4);%nbb30
                datx(j,21)=stux3(5);%tbb30
                datx(j,22)=stux1(9);%pitch
                datx(j,23)=stux1(10);%roll
            else
                datx(j,1:23)=NaN;
            end;

            dath(k,1)=datx(j,1);
            dath(k,2)=datx(j,2);
            dath(k,3)=datx(j,3);

            j=j+1;
            k=k+1;
        end;%isempty
    end;%while feof
    fclose(fid);

    if k>1
        [a b]=size(dath);
        ii=find(dath(:,2)<-1 | dath(:,2)>8);
        [nz mz]=size(ii)
        dath(ii,2)=NaN;
        dath(ii,3)=NaN;
        jdd=dath(:,1);
        jd0=floor(jdd(1));  % first Julian Day of data series
        jj=find(isfinite(dath(:,2)));

        if isempty(jj)
            dathr(nj,1)=datx(j-a,1);
            dathr(nj,2)=NaN;
            dathr(nj,3)=NaN;
            dathr(nj,4)=0;
            dathr(nj,5)=NaN;
            dathr(nj,6)=NaN;
        else
            kk=1;			 % if hourly data exists start 10-min calc.
            pp=jdd-floor(jdd); % fractional part of day

            while kk<7
                kj=find(pp>(pp(1)+jdmx(kk)) & pp<(pp(1)+jdmx(kk+1)));  % select 10-min section
                if isempty(kj)
                    datmn(kx,1)=dath(1,1)+jdmx(kk);
                    datmn(kx,2)=NaN;
                    datmn(kx,3)=NaN;
                    datmn(kx,4)=NaN;
                    datmn(kx,5)=NaN;
                    datmn(kx,6)=NaN;
                else
                    datmn(kx,1)=dath(kj(1),1);
                    kkj=find(~isnan(dath(kj,2)));%skip tip cals
                    if isempty(kkj)
                        datmn(kx,2)=NaN;
                        datmn(kx,3)=NaN;
                    else
                        datmn(kx,2)=median(dath(kj(kkj),2));
                        datmn(kx,3)=median(dath(kj(kkj),3));
                    end;
                    y1=dath(kj,2);
                    y=sort(y1);
                    [ny my]=size(y);
                    j1=max(floor(ny*.16),1);j2=max(floor(ny*.84),1);
                    datmn(kx,4)=ny;
                    datmn(kx,5)=(y(j2)-y(j1))/2;
                    y2=dath(kj,3);
                    y=sort(y2);
                    datmn(kx,6)=(y(j2)-y(j1))/2;
                end; %end if isempty kj
                %jdmin(kx)=jd0+jdx(jj)+jdmx(kk); % integrate timestep
                kk=kk+1;
                kx=kx+1;
            end  %while kk
            dathr(nj,1)=dath(1,1);
            dathr(nj,2)=median(dath(jj,2));
            dathr(nj,3)=median(dath(jj,3));
            y1=dath(jj,2);
            y=sort(y1);
            [ny my]=size(y);
            j1=max(floor(ny*.16),1);j2=max(floor(ny*.84),1);
            dathr(nj,4)=ny;
            dathr(nj,5)=(y(j2)-y(j1))/2;
            y2=dath(jj,3);
            y=sort(y2);
            dathr(nj,6)=(y(j2)-y(j1))/2;
        end;
        nj=nj+1;
    end;
    % end;%if filename

end;%for i
%NaN out the tip cals
ii=find(datx(:,2)<-1 | datx(:,2)>8);
datx(ii,2)=NaN;
datx(ii,3)=NaN;
%
daty=sortrows(datx,1);
dathry=sortrows(dathr,1);
datmny=sortrows(datmn,1);
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

    ss='c:\data\icealot\data\mailbox';

    %ss=s;
    e=[ss 'ICEALOT_08_micro_sb.txt'];
    flist1=fopen([ss 'ICEALOT_08_micro_s.txt'],'w');
    flist2=fopen([ss 'ICEALOT_08_micro_h.txt'],'w');
    flist3=fopen([ss 'ICEALOT_08_mn.txt'],'w');
else
    flist1=1;flist2=1;
end;

%fprintf(flist1,'%12.6f %10.5f %10.5f \r',daty');
fprintf(flist1,'%12.6f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f \r',daty');
%				jd	     vap   liq    tau20  tb30   tau30  tb30   mode   angle  cf20   cf30   sky20  nsky20 bb20   nbb20  tbb20  sky30  nsky30 bb30   nbb30  tbb30 pitch    roll
%                1        2      3      4      5      6      7      8      9     10     11     12     13     14     15     16     17     18     19    20      21     22      23
fprintf(flist2,'%12.6f %10.5f %10.5f %5d %10.5f %10.5f \r',dathry');
%						jd	     vap  liq   #obs  sigv   sigl

fprintf(flist3,'%12.6f %10.5f %10.5f %5d %10.5f %10.5f \r',datmny');
%						jd	     vap  liq   #obs  sigv   sigl


fclose('all');


%[jdav,y,yvar,nn]=timave_rawmicro(daty(:,1),10,daty);%average to 10-min
%y(:,1)=y(:,1)-5/24/60;%change mean time to beginning of record

%ee=[ss 'micro_10m1.txt'];
%save(ee,'y','-ascii','-tabs');


%save e datx -ascii;
ld = length(dathry);
ejd = max(dathry(:,1));
bjd = dathry(1,1);
% Plot last day
%figure;plot(dathry(:,1),dathry(:,2),dathry(:,1),dathry(:,5));axis([ejd-1 ejd 0 8]);
%xlabel('DOY (2007)');ylabel('Integrated atmospheric water vapor (cm)');
%figure;plot(dathry(:,1),dathry(:,3)*1e4,dathry(:,1),dathry(:,6)*1e4);axis([ejd-1 ejd -10 400]);
%xlabel('DOY (2007)');ylabel('Integrated cloud liquid content (g/m^2)');
%Plot multiple days
%figure;plot(dathry(:,1),dathry(:,2),dathry(:,1),dathry(:,5));axis([bjd ejd 0 8]);
way_images_mb = 'c:\data\icealot\data\mailbox\'
figure(1);
plot(dathry(:,1),dathry(:,2),'k');axis([bjd ejd 0 9]);
xlabel('Year Day (2008)');ylabel('Integrated atmospheric water vapor (cm)');
title('ICEALOT 2008 MWR')
print_buffer = [way_images_mb 'ICEALOT_MWR_all_vap.png'];
print('-dpng ', print_buffer);
%figure;plot(dathry(:,1),dathry(:,3)*1e4,dathry(:,1),dathry(:,6)*1e4);axis([bjd ejd -10 400]);

figure(2);
plot(dathry(:,1),dathry(:,3)*1e4,'k');axis([bjd ejd -10 800]);
xlabel('Year Day (2008)');ylabel('Integrated cloud liquid content (g/m^2)');
title('ICEALOT_2008 MWR')
print_buffer = [way_images_mb 'ICEALOT_MWR_all_liq.png'];
print('-dpng ', print_buffer);

save f:/AMMA_2007/RHB/radiometer/mailbox/raw/MBpwall

zilch=0;
if zilch
    %******************************   LWP zero bias calculations
    %datmny(:,3) is LWP
    %datmny(:,6) is stand dev of LWP in interval


    hg=0:23;
    hg=.0004+.04167*hg;
    dx=[];for i=1:24;ii=find(abs(datmny(:,1)-floor(datmny(:,1))-hg(i))<.002);dx=[dx ii'];end;%find indices with tip cals
    dy=datmny;
    dy(dx,:)=[];%time series without tip cals
    dz=sortrows(dx',1);%put in time order
    kk=[];
    liq0=zeros(length(datmny(:,1)),1);%set up radiometer bias correction
    for i=1:length(dx)-1
        kk=find(datmny(:,1)>datmny(dz(i),1) & datmny(:,1)<datmny(dz(i+1),1));%find 10-min values in 1 hr interval

        if ~isempty(kk)
            jj=find(1e4*datmny(kk,6)<25);%check for low std of LWP retrieval
            if ~isempty(jj);liq0([kk(1)-1 kk'],1)=min(datmny(kk(jj),3));end;%pick min value for zero bias
        end;
    end;

    ik=find( ~isnan(liq0) & ~isnan(datmny(:,3)));%valid microwave values
    ij=find((datmny(:,6))*1e4>15 & ~isnan(liq0) & ~isnan(datmny(:,3)));%presence of clouds
    disp('Mean LWP when cloud present and mean cloud fraction')
    [mean((datmny(ij,3)-liq0(ij))*1e4) length(ij)/length(ik)]

    figure;plot(datmny(:,1),datmny(:,3)*1e4,datmny(:,1),(datmny(:,3)-liq0)*1e4,'-x',dy(:,1),dy(:,6)*1e4);axis([276 293 -100 500]);
    xlabel('Julian Day (2007)');ylabel('LWP (g/m^2)');
end;

if 0
% Need to run amma07_read_mw21_composite and then pwv_sonde_07 to calculate sonde PWV
load f:/AMMA_2007/RHB/balloon/raw/AMMA_PWV_sonde

figure
plot(jd,pw,'r*')
hold
%plot(jd,pw,'r*')
plot(dathry(:,1),dathry(:,2),'k');axis([bjd ejd 0 9]);
xlabel('Julian Day (UTC)');
ylabel('PWV (cm)');
title('AMMA 2007 MWR and Sonde')
legend('Sonde','MWR')
print_buffer = ['F:\AMMA_2007\RHB\radiometer\Mailbox\Raw_Images\ICEALOT_PWV_sonde_MWR.jpg'];
print('-djpeg90 ', print_buffer);

end