fclose all;
path='C:\Program Files\MATLAB71\work';

read_parameters_rhb_GOMECC;
disp('read_son_day_rhb_07')

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

plotit=1;
prtit=1;

jhold=0;
jax=1;
jaxx=1;
backchk=0;
st1=[];
st2=[];
st3=[];
jaxd=[];
w=[];
W=[];


%for ddd=196;
    jd=num2str(ddd);
    if  ddd<100;
        jd=['0' num2str(ddd)];	
        end;
    %for hhh=0:23;
    %    clear sonraw
        hr=num2str(hhh);
        if hhh<10;
           hr=['0' num2str(hhh)];	
        end;


flist=fopen([wat_rawdata, '\GOMECC\day' ,jd, '\P0_',jd,hr]);
disp(['P0_', jd,hr]);
   
    if flist>0;                                     %if file exists
       dump=fgetl(flist);                           %read header lines and skip blanks
       dump2=fgetl(flist);                          %read header lines and skip blanks
       st1=fscanf(flist,'%g,%g,%g,%g,%g',[5,inf]);    
    end;
    
    ss = size(st1);
    if ss(1)==0;
       break;
    end;
    
    time=st1(1,:)/1000;
    u=st1(2,:)/100;
    v=st1(3,:)/100;
    w=st1(4,:)/100;
    
    
    if outputsonic==1
       Tsonic=st1(5,:)/100; %Temperature in degC
    else
        speedson=st1(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;
    
    W=[W;w'];    
    sonic=[time', u', v', w', Tsonic'];
    save ([path, '\GOMECC\Analysis\sonic_' , jd, hr, '.txt'], 'sonic', '-ascii')
   
    
    
    sonraw = st1;                    % load sonic array into temp array
    sonfix                           % run despiking routine
    clear st1
    st1 = sonraw_dep;                %copy despiked data into working array

    if backchk==0,
       tzilch=hhh;
       backchk=1;
       tzero=st1(1);
    end;   
         
    [nr1,nl1]=size(st1);
    nll=nl1-1;
            
    if nr1==5,
       st2(:,jax:jax+nl1-1)=st1(:,1:nl1);
       jax=jax+nl1;
       else
       x=fgets(flist);
    end;    %end if
        
                                     %sonfac=1200 or 600; 
        fsonic=20;          %20Hz    %sonic runs at 10Hz, 36000 values/hour...
        sonfac=3600*fsonic/60;
        if ~isempty(st2);
            jk=min(60,ceil(jax/sonfac));
            jade=str2num(jd)+((st2(1,:)-tzero)/1000/60/60+tzilch)/24;
            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;         %end if flist
end;            %end for

clear u v w Tsonic U V dir
u=st3(2,:)/100;
v=st3(3,:)/100;
w=st3(4,:)/100;

if 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;

if length(w)> 72000 | length(w)< 72000 
   time=interpft(time,72000);
   w=interpft(w,72000);
   v=interpft(V,72000);
   u=interpft(U,72000);
   Tsonic=interpft(Tsonic,72000);
end;
        time=(1:71880);
        u=u(1:71880);
        v=v(1:71880);
        w=w(1:71880); 
        Tsonic=Tsonic(1:71880);
        W=[W;w];

sonic=[time', u, v, w, Tsonic'];
save ([path, '\GOMECC\Analysis\sonic_' , jd, hr, '.txt'], 'sonic', '-ascii')


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=jhold;

if prtit==1;
   prt_jas_son_rhb_07_GOMECC;
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=jadson(1);
    j_en=jadson(np);
    figure;plot(jadson(1:np),Tsonic(1:np));xlabel(['Jday', year]);ylabel('Temperature (C)');title('Sonic Temperature');
    figure;plot(jadson(1:np),U(1:np),jadson(1:np),V(1:np),jadson(1:np),w(1:np),jadson(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');

    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
    figure;plot(jadson(1:np),dir(1:np),'b', ... 		%data points
        [ddd+4/24 ddd+4/24],[-180 180],'k', ...		%4 hour watch dividers
        [ddd+8/24 ddd+8/24],[-180 180],'k', ...
        [ddd+12/24 ddd+12/24],[-180 180],'k', ...
        [ddd+16/24 ddd+16/24],[-180 180],'k', ...
        [ddd+20/24 ddd+20/24],[-180 180],'k', ...
        [ddd ddd+1],[45 45],'g', ...			%45 degree good zone
        [ddd ddd+1],[-45 -45],'g', ...
        [ddd ddd+1],[60 60],'r', ...			%60 degree questionable zone
        [ddd ddd+1],[-60 -60],'r'); ...
        xlabel('Jday (2006)'); ...
        ylabel('Sonic Wind Direction (deg)'); ...
        title('Relative Wind Direction'); ...
        axis([j_st j_en -180 180]);
        year=2007;
    figure;plot(jadson(1:np),sqrt(U(1:np).^2+V(1:np).^2+w(1:np).^2));
    xlabel(['Jday', year]);ylabel('Wind Speed m/s');title('Sonic Relative Wind Speed');%axis([j_st j_en 0 360]);
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


time=0:0.1:86399.9;
W=interpft(W,864000);

%wmean=tenminmean_GOMECC(sonic_real(:,4)/100,10,10);
%W=[W;wmean];


figure;plot(time,W,'.'); %plot 10 min mean ozone data
ylabel('m / s ');title(['10min mean, Julian Day:',jd]);

s=find(W>-20 & W<20);
size(W)-size(s)

%time=sonic_data(1,:)/1000; 
%u=sonic_data(:,2)/100;
%v=sonic_data(:,3)/100;
%w=sonic_data(:,4)/100;

%time_real=0:0.1:3599.9;                  %shorten vector to 36000 values
%sonic_real=interpft(sonic_data,36000);

%wmean=tenminmean_GOMECC(sonic_real(:,4)/100,10,10);
%W=[W;wmean];

%sonic_day=[sonic_day;sonic_data];
%sonic_day_corr=[sonic_day_corr;sonic_real];

%   end;
%end;

%time=sonic_day(:,1)/1000; 
%time_uncorr=time-time(1);
%time=0:0.1:86399.9;
%w_raw=sonic_day(:,4)/100;
%w_corr=sonic_day_corr(:,4)/100;


%figure;plot(time,w_corr);xlabel('time / s');ylabel('m / s');title(['vertical wind raw, Julian Day:',jd]);
%print([path, '\GOMECC\Auswertung\w_raw_' jd  '.jpg'],'-djpeg');

%figure;plot(time,w_corr,'.');hold on;plot(time_uncorr,w_raw,'r.'); %red points are the uncorrected, original data points
%xlabel('time / s');ylabel('m / s');title(['w_raw and w_corr, Julian Day:',jd]);

%figure;plot(W,'.'); %plot 10 min mean ozone data
%ylabel('m / s ');title(['10min mean, Julian Day:',jd]);
%print([path, '\GOMECC\Auswertung\W_' jd  '.jpg'],'-djpeg');