disp('nau_motion')
% Inputs campbell data accelerations and angular rates 
% Computes one-minute averages of these variables and saves these to a file.
%  Note:  hhh=hour

fclose('all'); 
if hhh<10,
   dfx2=['0',num2str(hhh)];
else,
   dfx2=num2str(hhh);
end;

if ddd<10,
   dfx1=['00',num2str(ddd)];
elseif ddd<100,
   dfx1=['0',num2str(ddd)];
else,
   dfx1=num2str(ddd);
end;

dfx=[dfx1 dfx2];

Filtn=0;
dfl=['D:\hp_data\day',num2str(ddd),'\P5_',dfx];

flist=fopen(dfl,'r');%disp(dfl)
headr=fgetl(flist);%disp(headr)
headr1=fgetl(flist);
z=fscanf(flist,'%g%*c%g%*c%g%*c%g%*c%g%*c%g%*c%g%*c%g',[8,inf]);
%                1     2   3    4     5   6     7   8                  
tm1=z(1,:);tm2=z(2,:);tm2=tm2./10000;tmotion=tm1+tm2;%clear tm1;clear tm2;
tmotion=tmotion-tmotion(1);tmotion=0:1/20:(length(tmotion)-1)/20;tmotion=tmotion';

%  X is toward bow
%  Y is toward port
%  Z is up

accx=+z(7,:)./3.747.*9.81;              %forward is positive
accy=+z(6,:)./3.749.*9.81;              %to port is positive
accz=-z(8,:)./3.754.*9.81;              %up is positive
%ratex=+z(4,:)/0.049888/180*pi;          %port up is positive phi
%ratey=+z(3,:)/0.050343/180*pi;          %bow down is positive theta
%ratez=+z(5,:)/0.050084/180*pi;          %bow to starboard is positive psi
ratex=+z(4,:)/0.024993/180*pi;          %port up is positive phi
ratey=+z(3,:)/0.025016/180*pi;          %bow down is positive theta
ratez=+z(5,:)/0.024985/180*pi;          %bow to starboard is positive psi
clear z
G=sqrt(mean(accx)^2+mean(accy)^2+mean(accz)^2);
accx=accx';accy=accy';accz=accz';ratex=ratex';ratey=ratey';ratez=ratez';
if accz(1) > 100, accz(1)=0; end;
for ii=2:length(accz);
   if accz(ii) > 100, 
      accz(ii) = accz(ii-1);
   end;
end;
%clear z headr;

figure;plot(tmotion(1000:2000),accx(1000:2000),tmotion(1000:2000),accy(1000:2000));title(['JD',dfx1,' Hr',dfx2,' UTC  Horizontal Accelerations']);ylabel('acceleration in m/s/s');xlabel('time in seconds')
figure;plot(tmotion(1000:2000),accz(1000:2000));title(['JD',dfx1,' Hr',dfx2,' UTC  Acceleration Vertical']);ylabel('acceleration in m/s/s');xlabel('time in seconds')
figure;plot(tmotion(1000:2000),ratex(1000:2000),tmotion(1000:2000),ratey(1000:2000),tmotion(1000:2000),ratez(1000:2000));title(['JD',dfx1,' Hr',dfx2,' UTC  Rotation Rates']);ylabel('rate of rotation in rad/s');xlabel('time in seconds')
pause
clear tm1 tm2 tmotion
disp('computing spectra')
[Cax,F]=psd2(detrend(accx),length(accx));[Cax1,FF]=specsmooth(Cax,20);clear Cax F
[Cay,F]=psd2(detrend(accy),length(accy));[Cay1,FF]=specsmooth(Cay,20);clear Cay F
[Caz,F]=psd2(detrend(accz),length(accz));[Caz1,FF]=specsmooth(Caz,20);clear Caz F
[Crx,F]=psd2(detrend(ratex),length(ratex));[Crx1,FF]=specsmooth(Crx,20);clear Crx F
[Cry,F]=psd2(detrend(ratey),length(ratey));[Cry1,FF]=specsmooth(Cry,20);clear Cry F
[Crz,F]=psd2(detrend(ratez),length(ratez));[Crz1,FF]=specsmooth(Crz,20);clear Crz F

clear F Cax Cay Caz Crx Cry Crz z tm1 tm2 ii headr1 headr flist

figure;loglog(FF,Cax1.*FF,FF,Cay1.*FF,FF,Caz1.*FF);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Accelerations ']);xlabel('Freq in Hz');ylabel('S(acc)*F')
figure;loglog(FF,Crx1.*FF,FF,Cry1.*FF,FF,Crz1.*FF);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Rotation Rates ']);xlabel('Freq in Hz');ylabel('S(rate)*F')

clear FF Crx1 Cry1 Crz1 Cax1 Cay1 Caz1 dfx1 dfx2 dfx dfl ans G Filtn tmotion

