disp('read_motion_HiWinGS2013.m')

hr=sprintf('%02i',hhh);
jd=sprintf('%03i',ddd);
dfx=[jd hr];
temp=[];%z=[];
% z=NaN+zeros(9,36000);
z=[];
% dfl=[drive,'\day',jd,'\mot0' year jd hr '_raw.txt'];
% dfl=[way_raw_data_flux, yyyy(3:4),jd,'\mot0' yyyy(3:4) jd hr '_raw.txt'];
dfl=[way_raw_data_flux, yyyy(3:4),jd,'/mot0' yyyy(3:4) jd hr '_raw.txt']; % Mac file system
disp(dfl)
flist=fopen(dfl,'r');         %disp(dfl)
if flist>0
   while feof(flist)==0,

        try    
        temp =textscan(flist,'%2f%2f%3f %f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', 'Ship Motion');
        z=[z cell2mat(temp)'];
%         z(:,1:1+length(temp{1,1})-1)=cell2mat(temp)';
%        z(4:9,any(isnan(z),1))=NaN;      %put NaN into any rows containing NaNs (to remove bad lines, gap...)
        catch 
%          for ll=2:9
%              if length(temp{1,ll})~=length(temp{1,1})
%                  temp{1,ll}(length(temp{1,1}))=NaN;
%              end;
%          end;
         for ll=1:9
             if length(temp{1,ll})~=length(temp{1,9})
                 temp{1,ll}(length(temp{1,9})+1)=[];
             end;
         end;
         if ~isempty(temp{1,1});      
%            z(:,1:1+length(temp{1,1})-1)=cell2mat(temp)';
           z=[z cell2mat(temp)'];
         end;
        end; 
   end;
else
    disp('file not found')
end;

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

% % ttest=z(1,:)*60+z(2,:)+z(3,:)/1000; %real raw time in s
% tmotion=0:1/10:(length(z)-1)/10;  
% tmotion=tmotion';
tmotion=z(1,:)*60+z(2,:)+z(3,:)/1000; %real raw time in s
tref=0:1/10:35999/10;

if length(z(1,isnan(z(1,:))))==36000
    tmotion=tref';
end  
%  As mounted:  motion-pak
%  Xm is toward port
%  Ym is toward bow
%  Zm is down

[b, m, n] = unique(tmotion);
tmotion=tmotion(m);
mot1=z';
%  We desire the coordinate system such that
%  +X is toward bow
%  +Y is toward port
%  +Z is up
%  +rotation about the x-axis is port up (phi = roll) 
%  +rotation about the y-axis is bow down (theta = pitch)
%  +rotation about the z-axis is bow to port

%RHB motion pack tested in lab in September 2010.
accx=(+mot1(m,8)+0)./3.7756.*9.80;          %forward is positive 
accy=(+mot1(m,7)+0)./3.8215.*9.80;          %to port is positive
accz=(-mot1(m,9)+0)./3.7005.*9.80;              %up is positive
% %LB test when comparing to UConn and Wband Xbow --> does not seem to make
% %much difference on spectra on case study doy316. gets the pitch/roll
% %closer to Wband_Xbow too, but are they tilted the same way?
% accx=(+mot1(m,8)+0.09)./ 3.7756.*9.80;          %forward is positive 
% accy=(+mot1(m,7)-0.01)./3.8215.*9.80; 
% accz=(-mot1(m,9)+0.01)./3.692.*9.80;              %up is positive

ratex=(+mot1(m,5)+0)/0.025/180*pi;         %port up is positive phi
ratey=(+mot1(m,4)+0)/0.025/180*pi;          %bow down is positive theta
ratez=(-mot1(m,6)+0)/0.025/180*pi;      %bow to port is positive psi (right-handed)
% ratez=(-mot1(m,6)+0)/0.025023/180*pi;      %bow to port is positive psi (right-handed)


% try   
%     accx=despike_stratus07(accx);
%     accy=despike_stratus07(accy);
%     accz=despike_stratus07(accz);
%     ratex=despike_stratus07(ratex);
%     ratey=despike_stratus07(ratey);
%     ratez=despike_stratus07(ratez);
% catch
%     accx=accx;
%     accy=accy;
%     accz=accz;
%     ratex=ratex;
%     ratey=ratey;
%     ratez=ratez;
% end;

% figure;
% subplot(3,2,1);plot(accx);title('accx');
% subplot(3,2,2);plot(accy);title('accy');
% subplot(3,2,3);plot(accz);title('accz'); 
% subplot(3,2,4);plot(ratex);title('ratex');
% subplot(3,2,5);plot(ratey);title('ratey');
% subplot(3,2,6);plot(ratez);title('ratez');

G=sqrt(mean(accx)^2+mean(accy)^2+mean(accz)^2);

accx=accx';
accy=accy';
accz=accz';
ratex=ratex';
ratey=ratey';
ratez=ratez';
tmotion=tmotion';

% tref=0:1/10:35999.9/10;%tmotion=tmotion';
% saccx2=interp1(ttest,accx,tref);


% fclose(flist);

% interpolation=1;
% if interpolation
% figure;plot(diff(tmotion),'.')
%only interpolate if too big od data gap, otherwise wireless modem create
%delays (which compensate at the end to 0.1s)
% if length(z)< 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.     
    
    % Lenmotion=36000;
    % %interpolation at 10Hz
%     saccx=interp1(tmotion,accx,tref,'spline');ii=find(~isnan(saccx));msaccx=mean(saccx(ii));ii=find(isnan(saccx));saccx(ii)=msaccx;   
%     saccy=interp1(tmotion,accy,tref,'spline');ii=find(~isnan(saccy));msaccy=mean(saccy(ii));ii=find(isnan(saccy));saccy(ii)=msaccy;
%     saccz=interp1(tmotion,accz,tref,'spline');ii=find(~isnan(saccz));msaccz=mean(saccz(ii));ii=find(isnan(saccz));saccz(ii)=msaccz;
%     sratex=interp1(tmotion,ratex,tref,'spline');ii=find(~isnan(sratex));msratex=mean(sratex(ii));ii=find(isnan(sratex));sratex(ii)=msratex;   
%     sratey=interp1(tmotion,ratey,tref,'spline');ii=find(~isnan(sratey));msratey=mean(sratey(ii));ii=find(isnan(sratey));sratey(ii)=msratey;
%     sratez=interp1(tmotion,ratez,tref,'spline');ii=find(~isnan(sratez));msratez=mean(sratez(ii));ii=find(isnan(sratez));sratez(ii)=msratez;
%     removed spline here as it creates some glitches when data gaps are too big
    saccx=interp1(tmotion,accx,tref);ii=find(~isnan(saccx));msaccx=mean(saccx(ii));ii=find(isnan(saccx));saccx(ii)=msaccx;   
    saccy=interp1(tmotion,accy,tref);ii=find(~isnan(saccy));msaccy=mean(saccy(ii));ii=find(isnan(saccy));saccy(ii)=msaccy;
    saccz=interp1(tmotion,accz,tref);ii=find(~isnan(saccz));msaccz=mean(saccz(ii));ii=find(isnan(saccz));saccz(ii)=msaccz;
    sratex=interp1(tmotion,ratex,tref);ii=find(~isnan(sratex));msratex=mean(sratex(ii));ii=find(isnan(sratex));sratex(ii)=msratex;   
    sratey=interp1(tmotion,ratey,tref);ii=find(~isnan(sratey));msratey=mean(sratey(ii));ii=find(isnan(sratey));sratey(ii)=msratey;
    sratez=interp1(tmotion,ratez,tref);ii=find(~isnan(sratez));msratez=mean(sratez(ii));ii=find(isnan(sratez));sratez(ii)=msratez;
    
    accx=saccx;
    accy=saccy;
    accz=saccz;
    ratex=sratex;
    ratey=sratey;
    ratez=sratez;
% elseif length(z)< 36000
% % if length(z)< 36000
%     accx=[accx (1:36000-length(z)).*0+accx(length(z))];
%     accy=[accy (1:36000-length(z)).*0+accy(length(z))];
%     accz=[accz (1:36000-length(z)).*0+accz(length(z))];
%     ratex=[ratex (1:36000-length(z)).*0+ratex(length(z))];
%     ratey=[ratey (1:36000-length(z)).*0+ratey(length(z))];
%     ratez=[ratez (1:36000-length(z)).*0+ratez(length(z))];
% elseif length(z)> 36000
% % if length(z)> 36000
%     accx=accx(1:36000);
%     accy=accy(1:36000);
%     accz=accz(1:36000);
%     ratex=ratex(1:36000);
%     ratey=ratey(1:36000);
%     ratez=ratez(1:36000);
% end;

slopeaccx=polyfit(tref,accx,1);
% end
% % %%%%%%%%%  file hhh=oo always starts late by 6-7s. Fix the interpolation issues at beginning of hhh 00 %%%%%%%%%%%%%%%%%%%%%%
% % if hhh==00
% %     t1=floor(tmotion(1)*10);
% %     saccx(1:t1+5)=saccx(t1+10);
% %     saccy(1:t1+5)=saccy(t1+10);
% %     saccz(1:t1+5)=saccz(t1+10);
% %     sratex(1:t1+5)=sratex(t1+10);
% %     sratey(1:t1+5)=sratey(t1+10);
% %     sratez(1:t1+5)=sratez(t1+10);
% % end;    
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%--> does not seem to be present for VOCALS08


% slopeaccx=polyfit(tref,accx,1);
% % f = polyval(slopeaccx,tref); 
% % figure;plot(tref,accx,'.',tref,f,'r')

% try   
%     accx=despike_vocals08(accx);
%     accy=despike_vocals08(accy);
%     accz=despike_vocals08(accz);
%     ratex=despike_vocals08(ratex);
%     ratey=despike_vocals08(ratey);
%     ratez=despike_vocals08(ratez);
% catch
%     accx=accx;
%     accy=accy;
%     accz=accz;
%     ratex=ratex;
%     ratey=ratey;
%     ratez=ratez;
% end;