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

jd=sprintf('%03i',ddd);

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

st1=ones(8,36100)*NaN;
jam=hhh;
hr=sprintf('%02i',jam);

%dfl=[way_raw_data_flux, year(3:4),jd,'\son0' year(3:4) jd hr '_raw.txt'];
dfl=[way_raw_data_flux, year(3:4),jd,'/son0' year(3:4) jd hr '_raw.txt'];% Mac file system

disp(dfl)
flist=fopen(dfl,'r');
st1=[];temp=[];
if flist>0,	   				%if the file exists,
    while feof(flist)==0,  
        
        try
%             temp=textscan(flist,'%2f%2f%3f %*3c %f %f %f %*c %f %f %*3c %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', 'Sonic');
% temp=textscan(flist,'%2f%2f%3f %*3c %f %f %f %*c %f %f %*3c %*[^\n]','delimiter', 'T,M ', 'headerlines', 1,  'emptyvalue', NaN);
temp=textscan(flist,'%2f%2f%3f %*c %f %f %f %*c %f %f %*3c %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN);
% temp=textscan(flist,'%2f%2f%3f %*c %f %f %f %*c %f %f %*3c %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', {'T','TT'});
            st1=[st1 cell2mat(temp)'];
        catch
         for ll=1:8
             [l8,c8]=size(temp{1,8});
             if length(temp{1,ll})~=l8
                 temp{1,ll}(l8+1)=[];
             end;
         end;
        if ~isempty(temp{1,1});
             st1=[st1 cell2mat(temp)'];
        end;
        end;     
        
    end;  
end;   

% st1(:,any(isnan(st1),1))=[]; %remove lines containing NaNs (bad lines, too long...)
if isempty(st1)
    st1=NaN+zeros(8,36000);
end 

tref=0:1/10:35999/10;
tsons=st1(1,:)*60+st1(2,:)+st1(3,:)/1000; %real raw time in s

if length(st1(1,isnan(st1(1,:))))==36000
    tsons=tref;
end  

[b, m, n] = unique(tsons);

tsons=tsons(m);
u=st1(4,m);  %m/s
v=st1(5,m);  %m/s
W=st1(6,m);  %m/s
Tsonic=st1(7,m); %Temperature in degC

% try
    u=despike_dynamo2011(u);  %m/s
    v=despike_dynamo2011(v);  %m/s
    W=despike_dynamo2011(W);  %m/s
% catch
%     u=u;v=v;W=W;
% end

% [Su,Fson_u]=psd2(detrend(u),36000,10,hamming(36000));[Su,Fson_u]=specsmoo(Su,10);
% [Sv,Fson_v]=psd2(detrend(v),36000,10,hamming(36000));[Sv,Fson_v]=specsmoo(Sv,10);
% [Sw,Fson_w]=psd2(detrend(W),36000,10,hamming(36000));[Sw,Fson_w]=specsmoo(Sw,10);

% [Su,Fson_u]=psd2(detrend(u),36000,10,hamming(36000));[Su,Fson_u,df_u]=specsmooludo(Su,10);
% [Sv,Fson_v]=psd2(detrend(v),36000,10,hamming(36000));[Sv,Fson_v,df_v]=specsmooludo(Sv,10);
% [Sw,Fson_w]=psd2(detrend(W),36000,10,hamming(36000));[Sw,Fson_w,df_w]=specsmooludo(Sw,10);
% figure;loglog(Fson_w,Fson_w.*Sw,'--')

% % % **** Sonic configuration for VOCALS 2008****************************
% % sonicmodel='VOCALS'; %'R3'        %Sonic model
% % fsonic=10;                %frequency (Hz)
% % rotationsonic=false;      %rotate 30 deg = true
% % outputsonic=1; 		%Sonic output. 0=speed of sound, 1=Temperature (degC)
% % %  ********************************************************************
% Tsonic=st1(7,:); %Temperature in degC
% 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
		  
if rotationsonic
    switch sonicmodel           
        case {'R3','WindMasterPro'}     %omnidirectional  (N symbol pointing backward!)
            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
        case 'R3A'     %asymmetric      (N symbol pointing forward!)
            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
        case 'R2'      %omnidirectional  (N symbol pointing backward!)
            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
        case 'R2A'     %asymmetric      (N symbol pointing forward!)
            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
    switch sonicmodel
        case 'WindMasterPro'         %(N symbol pointing backward!)
            U=-u; %to North
            V=-v; %to West 
        case 'R3'                    %omnidirectional  (N symbol pointing backward!)
            U=-u; %to North
            V=-v; %to West  
        case 'R3A'                    %asymmetric      (N symbol pointing forward!)
            U=u; %to North
            V=v; %to West      
    end      
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+360,360);

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
clear f e dmp2 dmp backchk ans Tul Tll S Rtul Rtll Rsul Rsll Rlul Rlll st1 st2 u v
clear mU mV mW

% interpolation=1;
% if interpolation
% if length(U)< 36000
%         disp('<<<<<<<<<<<<<<<<<<  achtung >>>>>>>>>>>>>>>>>>>>>>>>')

% HiWIngs 2013. - Interpolate no matter what length of data set. Fluctuating time series seems to create distortion issues and reduce effectiveness of motion correction.     

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

    su=interp1(tsons,U,tref);su=su';ii=find(~isnan(su));msu=mean(su(ii));ii=find(isnan(su));su(ii)=msu;   %interpolation on the time base tmotion  
    sv=interp1(tsons,V,tref);sv=sv';ii=find(~isnan(sv));msv=mean(sv(ii));ii=find(isnan(sv));sv(ii)=msv;      %and replace the NaN by the mean.
    sw=interp1(tsons,W,tref);sw=sw';ii=find(~isnan(sw));msw=mean(sw(ii));ii=find(isnan(sw));sw(ii)=msw;
    st=interp1(tsons,Tsonic,tref);st=st';ii=find(~isnan(st));mst=mean(st(ii));ii=find(isnan(st));st(ii)=mst;
%     su=interp1(tsons,U,tref,'spline');su=su';ii=find(~isnan(su));msu=mean(su(ii));ii=find(isnan(su));su(ii)=msu;   %interpolation on the time base tmotion  
%     sv=interp1(tsons,V,tref,'spline');sv=sv';ii=find(~isnan(sv));msv=mean(sv(ii));ii=find(isnan(sv));sv(ii)=msv;      %and replace the NaN by the mean.
%     sw=interp1(tsons,W,tref,'spline');sw=sw';ii=find(~isnan(sw));msw=mean(sw(ii));ii=find(isnan(sw));sw(ii)=msw;
%     st=interp1(tsons,Tsonic,tref,'spline');st=st';ii=find(~isnan(st));mst=mean(st(ii));ii=find(isnan(st));st(ii)=mst;
% 
    % [Ssu,Fson_su]=psd2(detrend(su),36000,10,hamming(36000));[Ssu,Fson_su]=specsmoo(Ssu,10);
    % [Ssv,Fson_sv]=psd2(detrend(sv),36000,10,hamming(36000));[Ssv,Fson_sv]=specsmoo(Ssv,10);
    % [Ssw,Fson_sw]=psd2(detrend(sw),36000,10,hamming(36000));[Ssw,Fson_sw]=specsmoo(Ssw,10);
    % figure;loglog(Fson_w,Fson_w.*Sw,'--',Fson_sw,Fson_sw.*Ssw,'r--')

    % [Ssu,Fson_su]=psd2(detrend(su),36000,10,hamming(36000));[Ssu,Fson_su,df_su]=specsmooludo(Ssu,10);
    % [Ssv,Fson_sv]=psd2(detrend(sv),36000,10,hamming(36000));[Ssv,Fson_sv,df_sv]=specsmooludo(Ssv,10);
    % [Ssw,Fson_sw]=psd2(detrend(sw),36000,10,hamming(36000));[Ssw,Fson_sw,df_sw]=specsmooludo(Ssw,10);
    % [sum(Sw.*df_w) sum(Ssw.*df_sw) var(detrend(W))]

    % % %%%%%%%%%  file hhh=oo always starts late by 6-7s. Fix the interpolation issues at beginning of hhh 00 %%%%%%%%%%%%%%%%%%%%%%
    % % if hhh==00
    % %     t1=floor(tsons(1)*10);
    % %     su(1:t1+5)=su(t1+10);
    % %     sv(1:t1+5)=sv(t1+10);
    % %     sw(1:t1+5)=sw(t1+10);
    % %     st(1:t1+5)=st(t1+10);
    % % end;    
    % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    %--> does not seem to be present for VOCALS08

    U=su';
    V=sv';
    W=sw';
    Tsonic=st';
% elseif length(U)< 36000
% % if length(U)< 36000
%     U=[U (1:36000-length(U)).*0+U(length(U))];
%     V=[V (1:36000-length(V)).*0+V(length(V))];
%     W=[W (1:36000-length(W)).*0+W(length(W))];
%     Tsonic=[Tsonic (1:36000-length(Tsonic)).*0+U(length(Tsonic))];
%     tsons=[tsons (1:36000-length(tsons)).*0+tsons(length(tsons))];
% elseif length(U)> 36000    
% % if length(U)> 36000    
%     U=U(1:36000);
%     V=V(1:36000);
%     W=W(1:36000);
%     Tsonic=Tsonic(1:36000);
%     tsons=tsons(1:36000);
% % end
% end
% 
tmotion=0:1/10:35999/10;

