disp('rdson_07')   %('read_son_day_rhb_07')
fclose('all');
%Sonic model R3 used on AMMA06. Running at 10Hz
% S/N 000138

clear st1;clear st2;clear st3;clear jadson;clear sonraw;

plotit=1;
prtit=0;

hr='00';
jhold=0;
jax=1;
jaxx=1;
backchk=0;
st1=[];
st2=[];
st3=[];
jaxd=[];
%if 0 % DEW 5/15/2007
for jam=1:1   %23,						%cycle thru daily files
    clear sonraw

    e=[way_raw_data_sonic, 'SONIC_15n_20070820000000.dat'];
    %e=[way_raw_data_sonic, 'SONIC_6_20070829000000.dat'];
    lvl = e(26:28);  % save ht of sonic
    flist=fopen(e,'r');
    % U Mass sonic format
    %2007 09 08 00 00 00.024766    1.31   4.41  -0.20   23.75   0 10.9
    %2007 09 08 00 00 00.056767    1.28   4.47  -0.15   23.76   0 10.9
    % Old P0 format
    %Sonic
    %12823142400046,-917,-843,225,3375
    %12823142400156,-914,-837,216,3373

    if flist>0,	   				%if the file exists,
        disp(['Reading U Mass Sonic file ',e])
        ic = 0;
        for kk = 1:115200  % 32 * 60 * 60 test loop DEW
            ic = ic + 1;
            %while feof(flist)==0,
            yr(ic)=fscanf(flist,'%g',[1]);
            mm(ic)=fscanf(flist,'%g',[1]);
            dd(ic)=fscanf(flist,'%g',[1]);
            hh(ic)=fscanf(flist,'%g',[1]);
            mn(ic)=fscanf(flist,'%g',[1]);
            ss(ic)=fscanf(flist,'%g',[1]);
            w(ic)=fscanf(flist,'%g',[1]);
            u(ic)=fscanf(flist,'%g',[1]);
            v(ic)=fscanf(flist,'%g',[1]);
            tt(ic)=fscanf(flist,'%g',[1]);
            tmp1=fscanf(flist,'%s',[1]);
            tmp2=fscanf(flist,'%g',[1]);

            % despike data  DEW 5/6/2007 fm CWK
            sss = size(yr);
            %put in because file P0_12812 was bombing after it look like it
            %read the file once after 64004 pts and then tried to read
            %again but was empty 71857 pts in file DEW 5/9/2007
            if ss(1) == 0
                break
            end
        end;       %end while feof
        st1(1,:)=datenum(yr,mm(:,2),dd,hh,mn,ss)-datenum(yr-1,12,31);
        st1(2,:) = u;
        st1(3,:) = v;
        st1(4,:) = w;
        st1(5,:) = tt;

        ddd = fix(st1(1,1));
        size(st1)


        sonraw = st1;   % load sonic array into temp array
        sonfix  % run despiking routine
        clear st1
        st1 = sonraw_dep; %copy despiked data into working array

        jd=num2str(st1(1,1));

        if 0   %DEW 12182007
            if backchk==0,
                tzilch=jam;
                backchk=1;
                tzero=st1(1);
            end;   % end if
            [nr1,nl1]=size(st1);
            nll=nl1-1;
            if nr1==5,
                st2(:,jax:jax+nl1-1)=st1(:,1:nl1);
                jax=jax+nl1;
                size(st2)
            else
                x=fgets(flist);
            end;    %end if


            %      sonfac=1200 or 600; %sonic runs at 32Hz, 115200 values/hour...
            sonfac=3600*fsonic/60;  % calculate the no of pts in 1 min
            if ~isempty(st2);
                jk=min(60,ceil(jax/sonfac));
                jade = st1(1,:);
                for j=1:jk,
                    jl=sonfac*(j-1)+1;
                    jj=min(jax-1,sonfac*j);
                    st3(:,j+jhold)=median(st2(:,jl:jj)')';
                    jadson(j+jhold)=jade(jl);
                end;    %end for
                stx=st2;
                jhold=jhold+jk;
                st2=[]; st1=[];
                jaxd=[jaxd jax];
                jax=1;
            end;       %if isempty
        end   % skip DEW 12182007
    end;         %end if flist
end;            %end for
%clear u v w Tsonic U V dir
u=st1(2,:);   %/100;
v=st1(3,:);    %/100;
w=st1(4,:);    %/100;

if 0  % DEW 12182007     outputsonic==1
    Tsonic=st3(5,:)/100; %Temperature in degC
    %  else
    speedson=st3(5,:)/50;
    Tsonic=(speedson.*speedson+(u.*u + w.*w)/2. + v.*v)/403-273.15;
end;

if rotationsonic==1
    C=strcmp(sonicmodel,'R3');
    if C==1                                            %means sonicmodel=R3
        U=u*cos(30/180*pi)-v*sin(30/180*pi);		   %to North
        V=u*sin(30/180*pi)+v*cos(30/180*pi);  		   %to West
    else
        U=-u*cos(30/180*pi)-v*sin(30/180*pi);		   %to North
        V=-u*sin(30/180*pi)+v*cos(30/180*pi);  		   %to West
    end;
else
    U=u;
    V=v;
end;

S=sqrt(u.*u+v.*v+w.*w);
dir=atan2(-V,U)*180/pi+180;		%jeff coordinates, +boward, +portward, 180 means wind from
%dir=mod(dir+720,360);
np=length(st1);   %jhold;

if prtit==1;
    prt_jas_son_rhb_07;
end;

if plotit==1;
    %set some plot limits
    Tll=-25;
    hll=70;
    uul=10;
    vul=10;
    Tul=-0;
    hul=100;
    sul=12;
    Rlul=30;
    Rlll=-100;
    Rtul=300;
    Rtll=125;
    Rsul=1000;
    Rsll=0;

    %plot T profiles
    j_st=st1(1,1);
    j_en=st1(1,np);

    tit = cruise;
    lt = length(cruise);
    tit(lt+1) = '-';
    tit((lt+2):lt+5) = year;
    tit(lt+6) = '-';
    mms = num2str(mm(1));
    ls = length(mms);
    if ls == 2
        tit((lt+7):lt+8) = mms;
    else
        tit(lt+7) = '0';
        tit(lt+8) = mms;
    end

    tit(lt+9) = '-';
    dds = num2str(dd(1));
    ls = length(dds);
    if ls == 2
        tit((lt+10):lt+11) = dds;
    else
        tit(lt+10) = '0';
        tit(lt+11) = dds;
    end

    tit(lt+12) = '-';
    hhs = num2str(hh(1));
    ls = length(hhs);
    if ls == 2
        tit((lt+13):lt+14) = hhs;
    else
        tit(lt+13) = '0';
        tit(lt+14) = hhs;
    end

    tit(lt+15:lt+18) = '....';
    tit((lt+18):lt+20) = lvl;
    tit(lt+21) = 'm';
    figure;
    plot(st1(1,1:np),st1(5,1:np));
    xlabel(['Jday', year]);ylabel('Temperature (C)');
    title(tit);
    if pltit
        print_buffer = [way_images_sonic cruise '_' year '_SONIC_T_' mms '_' dds '_' hhs '_' lvl 'm.png'];
        print('-dpng ', print_buffer);
    end

    figure;
    plot(st1(1,1:np),U(1:np),st1(1,1:np),V(1:np),st1(1,1:np),w(1:np),st1(1,1:np),S(1:np));
    xlabel(['Jday', year]);ylabel('Sonic Winds (m/s)');title('Relative Wind Components:  U=blue,  V=green,  W=red,   Speed=cyan');
    if pltit
        print_buffer = [way_images_sonic cruise '_' year '_SONIC_UVWS_' mms '_' dds '_' hhs '_' lvl 'm.png'];
        print('-dpng ', print_buffer);
    end
    if 0  % DEW 12182007
        for i=1:np % Jeff Otten's mod to display direction -180 to +180.
            if dir(i) > 180
                dir(i) = dir(i) - 360;
            end
        end
    end  % DEW 12182007

    figure;
    plot(st1(1,1:np),dir(1:np))
    xlabel('Jday (2006)'); ...
        ylabel('Sonic Wind Direction (deg)');
    title(tit);
    if pltit
        print_buffer = [way_images_sonic cruise '_' year '_SONIC_DIR_' mms '_' dds '_' hhs '_' lvl 'm.png'];
        print('-dpng ', print_buffer);
    end
    %axis([j_st j_en -180 180]);
    figure;
    plot(st1(1,1:np),sqrt(U(1:np).^2+V(1:np).^2+w(1:np).^2));
    xlabel(['Jday', year]);
    ylabel('Wind Speed m/s');
    title(tit);%axis([j_st j_en 0 360]);
    if pltit
        print_buffer = [way_images_sonic cruise '_' year '_SONIC_SPD_' mms '_' dds '_' hhs '_' lvl 'm.png'];
        print('-dpng ', print_buffer);
    end

end;
if ~isempty(st1)

    Tsonic=st1(5,:); %Temperature in degC
    Tslow=Tsonic;

    if rotationsonic==1
        C=strcmp(sonicmodel,'R3');
        if C==1                                            %means sonicmodel=R3
            U=u*cos(30/180*pi)-v*sin(30/180*pi);		   %to North
            V=u*sin(30/180*pi)+v*cos(30/180*pi);  		   %to West
        else
            U=-u*cos(30/180*pi)-v*sin(30/180*pi);		   %to North
            V=-u*sin(30/180*pi)+v*cos(30/180*pi);  		   %to West
        end;
    else
        U=u;
        V=v;
    end;

    S=sqrt(u.*u+v.*v+w.*w);
    uu=st1(2,:);
    vv=st1(3,:);
    ww=st1(4,:);
    tt=st1(5,:);%-273.15;
    UU=uu;
    VV=vv;
    dir=atan2(-V,U)*180/pi+180;		%jeff coordinates, +boward, +portward, 180 means wind from
    dir=mod(dir+720,360);
    np=jhold;

    %clear w vul v uul u tzilch tzero sul st3 speedson prtit prt_it plotit nr1 npz np nl1 jl jk jj jhold jd jazx jaxx jax jam jadson jade j_st j_en j hul hr hll flist
    %clear f e dmp2 dmp dir backchk ans V U Tul Tsonic Tll S Rtul Rtll Rsul Rsll Rlul Rlll

    if 1   %hhh==0 | hhh==01 | hhh==12 || hhh==13 | hhh==15 | hhh==18 | hhh==22,
        disp(['Calculating spectra for ',hhs])
        SS=sqrt(UU.^2+VV.^2);
        [Su,F]=psd2(detrend(SS),length(SS),20.);
        [Sus,Fsu]=specsmoo(Su,20.);
        [Sw,F]=psd2(detrend(ww),length(ww),20.);
        [Sws,Fsw]=specsmoo(Sw,20.);

        figure;
        plot(sonraw(2,:),'k')
        hold
        plot(sonraw_dep(2,:),'r.')
        %plot(sonraw_dep(3,:),'g.')
        %plot(sonraw_dep(4,:),'r.')
        %plot(sonraw_dep(5,:),'k.')
        legend('Raw-U','Despiked-U')
        title(tit)
        if pltit
            print_buffer = [way_images_sonic cruise '_' year '_SONIC_raw_' mms '_' dds '_' hhs '_' lvl 'm.png'];
            print('-dpng ', print_buffer);
        end
        figure;
        loglog(Fsu(length(Fsu)-31:length(Fsu)),((Fsu(length(Fsu)-31:length(Fsu))).^(-2/3))./(Fsu(length(Fsu)-18).^(-2/3)).*Sus(length(Fsu)-18).*Fsu(length(Fsu)-18),'r');
        hold
        loglog(Fsu,Sus.*Fsu,'b');
        loglog(Fsw,Sws.*Fsw,'g');
        title(tit);
        legend('5/3','Spd','W' );
        xlabel('Frequency in Hz');
        ylabel('Spectrum * Frequency');
        if pltit
            print_buffer = [way_images_sonic cruise '_' year '_Sonic_Spectra_' mms '_' dds '_' hhs '_' lvl 'm.png'];
            print('-dpng ', print_buffer);
        end
        if 0  Print out spectra file for hr 12
            specfile=[way_spec_data, '\w_sp',jd,hhs,'.txt'];
            fspec=fopen(specfile,'w');
            specs=[Fsu
                Sus
                Fsw
                Sws];
            fprintf(fspec,'%11.5f %11.5f %11.5f %11.5f \r\n',specs);
            fclose(fspec);
            clear fspec specs specfile
        end;%end hhh=12
        clear SS UU VV Su F Sus Fsu Sw F Sws Fsw
    end;%end hhh=0 or 12 or 18
    clear ww Tsonic UU VV jade jadson tt uu vv ww SS
else
    U=[];V=[];w=[];dir=[];Tslow=[];S=[];
end;

%clear w vul v uul u tzilch tzero sul st3 speedson prtit prt_it plotit nr1 npz np nl1 jl jk jj jhold jd jazx jaxx jax jam jadson jade j_st j_en j hul hr hll flist
%clear f e dmp2 dmp dir backchk ans V U Tul Tsonic Tll S Rtul Rtll Rsul Rsll Rlul Rlll