disp('read_son_neaqs_GOMECC.m')
fclose('all');
clear st1;clear st2;clear st3;clear jadson;

read_parameters_rhb_GOMECC;
disp('read_son_day_rhb_07')

jd=num2str(ddd);

jhold=0; 
jax=1;
jaxx=1;
st1=zeros(6,1);
st2=zeros(6,75000);

jam=hhh;
if jam<10,
   hr=['0',num2str(jam)];
else
   hr=num2str(jam);
end;

%e=['/home/jeh/data/P0_' jd hr];
dfl=[way_rawdata ,'\GOMECC\day',jd,'\P0_' jd hr];

disp(dfl)
flist=fopen(dfl,'r');

st2=[];st1=[];
if flist>0,	   				%if the file exists, 
   dmp=fgets(flist);		%read header lines and skip blanks
   dmp2=fgets(flist); 		%read header lines and skip blanks
   while feof(flist)==0,
      st2=fscanf(flist,'%g,%g,%g,%g,%g',[5,inf]);
      st1=[st1 st2];
      dmp3=fgets(flist); 
   end;
end;   

tsons=st1(1,:);
tsons=0:1/20:(length(tsons)-1)/20;
tsons=tsons';

tmotion=0:1/20:71999/20;
tmotion=tmotion';

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

        time=st1(1,:)/1000;          %time in "s"
        U=st1(2,:)/100;
        V=st1(3,:)/100;
        W=st1(4,:)/100;

        u=despike_son2(U);
        v=despike_son2(V);
        w=despike_son2(W);
        Tsonic=despike_son2(Tson);
       
        %Despike 1. Dedtrend, 2. use despike, 3. retrend
        %[P,S]=polyfit(tsons',u,1);
        %uu=despike_son2(u-polyval(P,tsons'));
        %u=uu+polyval(P,tsons');
        
        %[P,S]=polyfit(tsons',v,1);
        %vv=despike_son2(v-polyval(P,tsons'));
        %v=vv+polyval(P,tsons');
        
        %[P,S]=polyfit(tsons',w,1);
        %ww=despike_son2(w-polyval(P,tsons'));
        %w=ww+polyval(P,tsons');
        
        
        %Find Nan and replace with the closest neighbour
        
        ii=find(isnan(u)); 
        for jj=1:length(ii)
            if isnan(u(1));
               u(ii(jj))=nanmean(U); 
            else
               u(ii(jj))=u(ii(jj)-1);
            end;
        end;
        
               
        ii=find(isnan(v)); 
        for jj=1:length(ii);
            if isnan(v(1))                
            v(ii(jj))=nanmean(V);
            else
            v(ii(jj))=v(ii(jj)-1);
            end;
        end;
       
        
        ii=find(isnan(w)); 
        for jj=1:length(ii)
            if isnan(w(1));
               w(ii(jj))=nanmean(W); 
            else
               w(ii(jj))=w(ii(jj)-1);
            end;
        end;
                
        
        ii=find(isnan(Tsonic)); 
        for jj=1:length(ii)
            if isnan(Tsonic(1));
               Tsonic(ii(jj))=nanmean(Tsonic); 
            else
               Tsonic(ii(jj))=Tsonic(ii(jj)-1);
            end;
        end;
        
                
        W=w;
        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;
        
mU=mean(U);
mV=mean(V);
mW=mean(W);

for ii=1:length(U);
   if U(ii)>35,
      U(ii)=mU;
   end;
   if V(ii)>35,
      V(ii)=mV;
   end;
   if W(ii)>35,
      W(ii)=mW;
   end;
end;
   
clear vul v uul u tzilch tzero sul st3 speedson prtit prt_it plotit nr1 npz np nl1 jl jk jj jhold jazx jaxx jax jam jadson jade j_st j_en j hul hll flist %jd hr
clear f e dmp2 dmp backchk ans Tul Tll S Rtul Rtll Rsul Rsll Rlul Rlll st1 st2 u v
clear mU mV mW

su=(1:72000).*0;
sv=(1:72000).*0;
sw=(1:72000).*0;
st=(1:72000).*0;
su(1)=U(1);
sv(1)=V(1);
sw(1)=W(1);
st(1)=Tsonic(1);

jj=2;
for ii=2:72000;
   while tsons(jj)<tmotion(ii) & jj < length(U);
      jj=jj+1;
   end;
   su(ii)=(U(jj)-U(jj-1))/(tsons(jj)-tsons(jj-1))*(tmotion(ii)-tsons(jj-1))+U(jj-1);
   sv(ii)=(V(jj)-V(jj-1))/(tsons(jj)-tsons(jj-1))*(tmotion(ii)-tsons(jj-1))+V(jj-1);
   sw(ii)=(W(jj)-W(jj-1))/(tsons(jj)-tsons(jj-1))*(tmotion(ii)-tsons(jj-1))+W(jj-1);
   st(ii)=(Tsonic(jj)-Tsonic(jj-1))/(tsons(jj)-tsons(jj-1))*(tmotion(ii)-tsons(jj-1))+Tsonic(jj-1);
end;

U=su;
V=sv;
W=sw;
Tsonic=st;

U=interpft(U,36000);
V=interpft(V,36000);
W=interpft(W,36000);
Tsonic=interpft(Tsonic,36000);

if hhh==0 | hhh==12 | hhh==18
   speed=[U V W Tsonic];
   save([way_rawdata ,'\GOMECC\flux_data\sonic_',jd, hr,'.txt'], 'speed','-ascii');
end;
   
%Un=mean(U).*cos(Mhead)-mean(-V).*sin(Mhead) + Msog.*cos(Mcog);   % + From South to North
%Ue=mean(U).*sin(Mhead)+mean(-V).*cos(Mhead) + Msog.*sin(Mcog);   % + From West to East

clear su sv sw st ii jj tsons tmotion

%disp('pause');pause
