% read_motion
% Xcruise 2009 :: 2009-09-07 :: LB
% Reads the ship motion data files for a day.

% SPdeS sped by preallocating arrays; textscan of 24*36000 rows is still
% slow.
disp(['read_motion: ', cruise]);
clearit=true;

% motion data is 10 Hz
nsampl=86400*10;
ncols=9;
jd=sprintf('%03i',ddd);
jd_mot=NaN+zeros(nsampl*1.01,1); % preallocate with tiny 1% buffer
mot1=NaN+zeros(nsampl*1.01,ncols);
%jhold=0; % never used
jax1=1;

for jam=00:23
    hr=sprintf('%02i',jam);
    dfl=[way_raw_data_flux, yyyy(3:4),jd,'\mot0' yyyy(3:4) jd hr '_raw.txt'];
    disp(['Reading motion file  ',dfl])
    flist=fopen(dfl,'r');         %disp(dfl)
    if flist>0
        clear temp z
        while ~feof(flist) % end-of-file not reached
            try
                temp =textscan(flist,'%2f%2f%3f %f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', 'Ship Motion');
                z=cell2mat(temp); 
            catch
                [lncols,cncols]=size(temp{1,ncols});
                if lncols==0 || cncols==0; continue; end % read next batch of lines
                for ll=1:ncols
                    if length(temp{1,ll})~=lncols
                        temp{1,ll}(lncols+1)=[];
                    end;
                end;
                if isempty(temp{1,1});
                    temp{1,ncols}=temp{1,ncols}';
                end;
                z=cell2mat(temp); 
            end; %end try/catch    
                                        
            [nl1,nr1]=size(z);
            mot1(jax1:jax1+nl1-1,:)=z; 
            jd1=ddd + (jam+(mot1(jax1:jax1+nl1-1,1)+(mot1(jax1:jax1+nl1-1,2)+mot1(jax1:jax1+nl1-1,3)/1000)/60)/60)/24;
            jd_mot(jax1:jax1+nl1-1,:) = jd1;
            jax1=jax1+nl1;
           
        end; %end ~feof(flist)
        fclose(flist);     
    else
        disp(sprintf('\b<--File not found.'))
    end
end

% trim arrays to their filled size
mot1=mot1(1:jax1-1,:);
jd_mot=jd_mot(1:jax1-1,:);

% % Time in flux data files is the number of 100-nano second intervals since
% % Jan 1, 1601 divided by 10000. It was chosen to maintain compatibility
% % with an HP-UNIX version of the code form long ago
% %  0             'dd-mmm-yyyy HH:MM:SS'  
% %01-Mar-2000 15:45:17
% dates = datestr(datenum(1601,1,1) + datenum(mot1(1,:)*10000*100e-9/86400));
% hh = str2num(dates(:,13:14));
% mm = str2num(dates(:,16:17));
% ss = str2num(dates(:,19:20));
% jd_mot= ddd + (hh+((mm+(ss/60))/60))/24;
% 
% % test to see if first point is from the previous day
% if str2num(dates(1,1:2)) ~= str2num(dates(2,1:2))
%     jd_mot(1) = jd_mot(1)-1;
% end    

%  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 May 2012.
% accx=(-z(7,:)-0.014)./3.9969293*9.80 ;         %forward is positive 
% accy=(-z(8,:)- 0.0000)./3.9604125.*9.80;          %to port is positive
% accz=(-z(9,:)+ 0.022)./ 3.58409.*9.80;        %up is positive
% ratex=(+z(5,:)+0.0172821)/0.025/180*pi;     %port up is positive phi
% ratey=(-z(4,:)-0.0082864)/0.025/180*pi;      %bow down is positive theta
% ratez=(-z(6,:)-0.035699)/0.025/180*pi;  %bow to port is positive psi (right-handed)
%offsets removed in Campbell program
accx=(-mot1(:,7))./3.9969293*9.80 ;         %forward is positive 
accy=(-mot1(:,8))./3.9604125.*9.80;          %to port is positive
accz=(-mot1(:,9))./ 3.58409.*9.80;        %up is positive
ratex=(+mot1(:,5))/0.025/180*pi;     %port up is positive phi
ratey=(-mot1(:,4))/0.025/180*pi;      %bow down is positive theta
ratez=(-mot1(:,6))/0.025/180*pi;  %bow to port is positive psi (right-handed)

% Note: 0.025 is 1 / (100deg/s / 2.5Vdc)  =  1 / 40 = 0.025


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

% trim off bad data at either end
for mm=length(accz):-1:1;
   if accz(mm)>50,
      accz(mm)=0;
   else
      break;
   end;
end;
for mm=1:1:length(accz);
   if accz(mm)>50,
      accz(mm)=0;
   else
      break;
   end;
end;

if clearit
    clear headr dfl G Filtn Cax Cax1 Cay Cay1 Caz Caz1 Crx Crx1 Cry Cry1 %hr jd
    clear Crz Crz1 FF ans ii %tmotion
end;

% zillions of points makes the display slow...
if plotit &&  ~isempty(mot1(~isnan(mot1)))
    figure;
    subplot(3,2,1);plot(jd_mot,accx,'.','markersize',2);ylabel('accx (m.s^-^2)','FontSize',8);%xlabel(['Yearday ', yyyy]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    subplot(3,2,2);plot(jd_mot,accy,'.','markersize',2);ylabel('accy (m.s^-^2)','FontSize',8);%xlabel(['Yearday ', yyyy]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    subplot(3,2,3);plot(jd_mot,accz,'.','markersize',2);ylabel('accz (m.s^-^2)','FontSize',8);%xlabel(['Yearday ', yyyy]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    subplot(3,2,4);plot(jd_mot,ratex,'.','markersize',2);ylabel('ratex (rad/s)','FontSize',8);%xlabel(['Yearday ', yyyy]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');ylim([-.2 .2])
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    subplot(3,2,5);plot(jd_mot,ratey,'.','markersize',2);ylabel('ratey (rad/s)','FontSize',8);xlabel('Hour (UTC)','FontSize',8);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');ylim([-.2 .2])
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    subplot(3,2,6);plot(jd_mot,ratez,'.','markersize',2);ylabel('ratez (rad/s)','FontSize',8);xlabel('Hour (UTC)','FontSize',8);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');ylim([-.1 .1])
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    annotation(gcf,'textbox',[0.35 0.96 0.8 0.02462],'String',sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD motion pack',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','FitBoxToText','off','LineStyle','none','Interpreter','none'); 
    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
    if prtit; print(graphdevice,[way_images_flux cruise '_PSD_motionpack_accs_rates_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end
end

if clearit
    clear jd_mot accx accy accz ratex ratey ratez mot1 temp st3 z
end;