%disp('openmotion.m')
% 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;
%disp('Open Motion-Pak Campbell File')
%dfl=['/home/user1/bin/data/P5_',dfx];
%dfl=['/scratch2/day',dfx1,'/P5_',dfx];
%dfl=['/home/user1/bin/data_backup/day',dfx1,'/P5_',dfx];
%dfl=['/home/jeh/fastex/P5_',dfx,';1'];
%dfl=['D:\data\day' dfx1 '\P5_',dfx];      
dfl=['f:\data\cruises\stratus_2007\raw\day',dfx1,'\mot007', dfx, '_raw.txt'];

flist=fopen(dfl,'r');%disp(dfl)
headr=fgetl(flist);%disp(headr)
headr1=fgetl(flist);
zm=fscanf(flist,'%g %g,%g,%g,%g,%g,%g',[7,inf]);
%                1  2  3  4  5  6  7              
z=zm(:,1:length(zm)-1);
tm1=z(1,:);tmotion=tm1;%clear tm1;clear tm2;
tmotion=tmotion-tmotion(1);tmotion=0:1/10:(length(tmotion)-1)/10;tmotion=tmotion';
mms=floor(tm1/1e5);secs=(tm1-mms*1e5)/1e3;
jd_mot=ddd+hhh/24+mms/60/24+secs/3600/24;
%  X is toward bow
%  Y is toward port
%  Z is up

accx=+z(6,:)./3.747.*9.81;              %forward is positive
accy=+z(5,:)./3.749.*9.81;              %to port is positive
accz=-z(7,:)./3.81*9.78;              %up is positive
%ratex=+z(4,:)/0.049888/180*pi;          %port up is positive phia
%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(3,:)/0.024993/180*pi;          %port up is positive phi
ratey=+z(2,:)/0.025016/180*pi;          %bow down is positive theta
ratez=+z(4,:)/0.024985/180*pi;          %bow to starboard is positive psi
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')

%lowpass filter time series at 2 Hz for signal sampled at 10 Hz
wp=[.02/5 2/5];ws= [.01/5 3/5];rp=3;rs=15;[n,wmn]=cheb1ord(wp,ws,rp,rs);
[b,a]=cheby1(n,rp,wmn);
yx=filter(b,a,ratex-mean(ratex));%roll rate rad/sec
yy=filter(b,a,ratey-mean(ratey));%pitch rate rad/sec

wpa=[2/5];wsa= [3/5];rp=3;rs=7.5;[n,wmn]=cheb1ord(wpa,wsa,rp,rs);
acx=filter(b,a,accx);%roll rate rad/sec
acy=filter(b,a,accy);%pitch rate rad/sec



%**********  convert to fft
yyx=fft(yx);
yyy=fft(yy);
yyx=fft(ratex-mean(ratex));
yyy=fft(ratey-mean(ratey));

fny=5;
nn=length(yx);
kk=1:nn;
%************  set up integral
i_yyx=-sqrt(-1)*nn/2/pi/fny./(kk-1)'.*yyx;
i_yyx(1)=0;
i_yyy=-sqrt(-1)*nn/2/pi/fny./(kk-1)'.*yyy;
i_yyy(1)=0;
%*************  compute angles in deg by inverse fft

roll=ifft(i_yyx)*180/pi;
pitch=ifft(i_yyy)*180/pi;


figure;plot(tmotion(1000:2000),roll(1000:2000),tmotion(1000:2000),pitch(1000:2000));title(['JD',dfx1,' Hr',dfx2,' UTC  Pitch and Roll']);ylabel('Pitch/roll in deg');xlabel('time in seconds')



figit=1;
if figit
Cax=psd(accx(1:12000),8192);Cay=psd(accy(1:12000),8192);Caz=psd(accz(1:12000),8192);Crx=psd(ratex(1:12000),8192);Cry=psd(ratey(1:12000),8192);Crz=psd(ratez(1:12000),8192);
%Chh=psd(heading,2048);[Chh1,Fh]=specsmooth(Chh,1);
[Cax1,FF]=specsmooth(Cax,20);[Cay1,FF]=specsmooth(Cay,20);[Caz1,FF]=specsmooth(Caz,20);
[Crx1,FF]=specsmooth(Crx,20);[Cry1,FF]=specsmooth(Cry,20);[Crz1,FF]=specsmooth(Crz,20);
nf=length(FF);dff=diff(FF);fc=1/(2*pi*2*pi);fa=180/pi*180/pi;
figure;loglog(FF(2:nf),Cax1(2:nf).*dff,FF(2:nf),Cay1(2:nf).*dff,FF(2:nf),Caz1(2:nf).*dff);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Accelerations #1']);xlabel('Freq in Hz');ylabel('f*S(acc)')
figure;loglog(FF(2:nf),fa*Crx1(2:nf).*dff,FF(2:nf),fa*Cry1(2:nf).*dff,FF(2:nf),fa*Crz1(2:nf).*dff);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Rotation Rates #1']);xlabel('Freq in Hz');ylabel('f*S(degrate)')
figure;loglog(FF(2:nf),fc*Cax1(2:nf).*dff./FF(2:nf)./FF(2:nf),FF(2:nf),fc*Cay1(2:nf).*dff./FF(2:nf)./FF(2:nf),FF(2:nf),fc*Caz1(2:nf).*dff./FF(2:nf)./FF(2:nf));title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Velocities #1']);xlabel('Freq in Hz');ylabel('f*S(vel)')
figure;loglog(FF(2:nf),fa*fc*Crx1(2:nf).*dff./FF(2:nf)./FF(2:nf),FF(2:nf),fa*fc*Cry1(2:nf).*dff./FF(2:nf)./FF(2:nf),FF(2:nf),fa*fc*Crz1(2:nf).*dff./FF(2:nf)./FF(2:nf));title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Angles #1']);xlabel('Freq in Hz');ylabel('f*S(degangle)')
%Cax=psd(accx(12001:24000),8192);Cay=psd(accy(12001:24000),8192);Caz=psd(accz(12001:24000),8192);Crx=psd(ratex(12001:24000),8192);Cry=psd(ratey(12001:24000),8192);Crz=psd(ratez(12001:24000),8192);
%[Cax1,FF]=specsmooth(Cax,20);[Cay1,FF]=specsmooth(Cay,20);[Caz1,FF]=specsmooth(Caz,20);
%[Crx1,FF]=specsmooth(Crx,20);[Cry1,FF]=specsmooth(Cry,20);[Crz1,FF]=specsmooth(Crz,20);
%figure;loglog(FF,Cax1./FF./FF,FF,Cay1./FF./FF,FF,Caz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Accelerations #2']);xlabel('Freq in Hz');ylabel('S(acc)')
%figure;loglog(FF,Crx1./FF./FF,FF,Cry1./FF./FF,FF,Crz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Rotation Rates #2']);xlabel('Freq in Hz');ylabel('S(rate)')

sigu=sqrt(sum(fc*Cax1(2:nf).*dff./FF(2:nf)./FF(2:nf)));
sigv=sqrt(sum(fc*Cay1(2:nf).*dff./FF(2:nf)./FF(2:nf)));
sigw=sqrt(sum(fc*Caz1(2:nf).*dff./FF(2:nf)./FF(2:nf)));
sigrol=sqrt(sum(fa*fc*Crx1(2:nf).*dff./FF(2:nf)./FF(2:nf)));
sigpit=sqrt(sum(fa*fc*Cry1(2:nf).*dff./FF(2:nf)./FF(2:nf)));
sighed=sqrt(sum(fa*fc*Crz1(2:nf).*dff./FF(2:nf)./FF(2:nf)));

['sigu= ' num2str(sigu) '   sigv= ' num2str(sigv) '   sigw= ' num2str(sigw)]
['sigroll= ' num2str(sigrol) '   sigpitch= ' num2str(sigpit)  '   sighead= ' num2str(sighed)]

end;
%Cax=psd(accx(24001:36000),8192);Cay=psd(accy(24001:36000),8192);Caz=psd(accz(24001:36000),8192);Crx=psd(ratex(24001:36000),8192);Cry=psd(ratey(24001:36000),8192);Crz=psd(ratez(24001:36000),8192);
%[Cax1,FF]=specsmooth(Cax,20);[Cay1,FF]=specsmooth(Cay,20);[Caz1,FF]=specsmooth(Caz,20);
%[Crx1,FF]=specsmooth(Crx,20);[Cry1,FF]=specsmooth(Cry,20);[Crz1,FF]=specsmooth(Crz,20);
%figure;loglog(FF,Cax1./FF./FF,FF,Cay1./FF./FF,FF,Caz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Accelerations #3']);xlabel('Freq in Hz');ylabel('S(acc)')
%figure;loglog(FF,Crx1./FF./FF,FF,Cry1./FF./FF,FF,Crz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Rotation Rates #3']);xlabel('Freq in Hz');ylabel('S(rate)')

%Cax=psd(accx(36001:48000),8192);Cay=psd(accy(36001:48000),8192);Caz=psd(accz(36001:48000),8192);Crx=psd(ratex(36001:48000),8192);Cry=psd(ratey(36001:48000),8192);Crz=psd(ratez(36001:48000),8192);
%[Cax1,FF]=specsmooth(Cax,20);[Cay1,FF]=specsmooth(Cay,20);[Caz1,FF]=specsmooth(Caz,20);
%[Crx1,FF]=specsmooth(Crx,20);[Cry1,FF]=specsmooth(Cry,20);[Crz1,FF]=specsmooth(Crz,20);
%figure;loglog(FF,Cax1./FF./FF,FF,Cay1./FF./FF,FF,Caz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Accelerations #4']);xlabel('Freq in Hz');ylabel('S(acc)')
%figure;loglog(FF,Crx1./FF./FF,FF,Cry1./FF./FF,FF,Crz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Rotation Rates #4']);xlabel('Freq in Hz');ylabel('S(rate)')

%Cax=psd(accx(48001:60000),8192);Cay=psd(accy(48001:60000),8192);Caz=psd(accz(48001:60000),8192);Crx=psd(ratex(48001:60000),8192);Cry=psd(ratey(48001:60000),8192);Crz=psd(ratez(48001:60000),8192);
%[Cax1,FF]=specsmooth(Cax,20);[Cay1,FF]=specsmooth(Cay,20);[Caz1,FF]=specsmooth(Caz,20);
%[Crx1,FF]=specsmooth(Crx,20);[Cry1,FF]=specsmooth(Cry,20);[Crz1,FF]=specsmooth(Crz,20);
%figure;loglog(FF,Cax1./FF./FF,FF,Cay1./FF./FF,FF,Caz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Accelerations #5']);xlabel('Freq in Hz');ylabel('S(acc)')
%figure;loglog(FF,Crx1./FF./FF,FF,Cry1./FF./FF,FF,Crz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Rotation Rates #5']);xlabel('Freq in Hz');ylabel('S(rate)')

%Cax=psd(accx(60000:length(accx)),8192);Cay=psd(accy(60000:length(accx)),8192);Caz=psd(accz(60000:length(accx)),8192);Crx=psd(ratex(60000:length(accx)),8192);Cry=psd(ratey(60000:length(accx)),8192);Crz=psd(ratez(60000:length(accx)),8192);
%[Cax1,FF]=specsmooth(Cax,20);[Cay1,FF]=specsmooth(Cay,20);[Caz1,FF]=specsmooth(Caz,20);
%[Crx1,FF]=specsmooth(Crx,20);[Cry1,FF]=specsmooth(Cry,20);[Crz1,FF]=specsmooth(Crz,20);
%figure;loglog(FF,Cax1./FF./FF,FF,Cay1./FF./FF,FF,Caz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Accelerations #6']);xlabel('Freq in Hz');ylabel('S(acc)')
%figure;loglog(FF,Crx1./FF./FF,FF,Cry1./FF./FF,FF,Crz1./FF./FF,Fh,Chh1);title(['JD',dfx1,' Hr',dfx2,' UTC  Spectra of Rotation Rates #6']);xlabel('Freq in Hz');ylabel('S(rate)')

%clear flist dfx2 dfx1 dfx dfl G Filtn Cax Cax1 Cay Cay1 Caz Caz1 Crx Crx1 Cry Cry1 Crz Crz1 FF ans ii
