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')
way_rawdata='C:\Data';

jhold=0; 
jax=1;
jaxx=1;


for ddd=192;
    jd=num2str(ddd);
for hhh=1:23
    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;
           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;

        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);
        
       
        %Despike 1. Dedtrend, 2. use despike, 3. retrend
        [P,S]=polyfit(tsons',U,1);
        uu=despike_son2(U-polyval(P,tsons'));
        ud=uu+polyval(P,tsons');
        
        [P,S]=polyfit(tsons',V,1);
        vv=despike_son2(V-polyval(P,tsons'));
        vd=vv+polyval(P,tsons');
        
        [P,S]=polyfit(tsons',W,1);
        ww=despike_son2(w-polyval(P,tsons'));
        wd=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(ud)); 
        for jj=1:length(ii)
            if isnan(ud(1));
               ud(ii(jj))=nanmean(ud); 
            else
               ud(ii(jj))=ud(ii(jj)-1);
            end;
        end;
        
        ii=find(isnan(vd)); 
        for jj=1:length(ii);
            if isnan(vd(1))                
               vd(ii(jj))=nanmean(vd);
            else
               vd(ii(jj))=vd(ii(jj)-1);
            end;
        end;
        
        ii=find(isnan(wd)); 
        for jj=1:length(ii)
            if isnan(wd(1));
               wd(ii(jj))=nanmean(wd); 
            else
               wd(ii(jj))=wd(ii(jj)-1);
            end;
        end;
        
        Wd=wd
        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;
        
        
        if rotationsonic==1
           C=strcmp(sonicmodel,'R3');
           if C==1                                              %means sonicmodel=R3
              Ud=ud*cos(30/180*pi)-vd*sin(30/180*pi);              %to North
              Vd=ud*sin(30/180*pi)+vd*cos(30/180*pi);              %to West
           else
              Ud=-ud*cos(30/180*pi)-vd*sin(30/180*pi);             %to North
              Vd=-ud*sin(30/180*pi)+vd*cos(30/180*pi);             %to West
           end;
        else
            Ud=ud;
            Vd=vd;
        end;
        
mU=mean(U);
mV=mean(V);
mW=mean(W);

mUd=mean(Ud);
mVd=mean(Vd);
mWd=mean(Wd);


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;
   
for ii=1:length(Ud);
   if Ud(ii)>35,
      Ud(ii)=mU;
   end;
   if Vd(ii)>35,
      Vd(ii)=mV;
   end;
   if Wd(ii)>35,
      Wd(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);

%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