disp('fastexrate.m')
%clear;
%openmotion;opengyro;openship;opensons;openophir
%

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];

sfl=['c:\hare\fastex\P_',dfx,'.dat'];

tgyro(length(tgyro))=tmotion(length(tmotion));
qq=1;shead=(1:length(tmotion)).*0;chead=(1:length(tmotion)).*0;heading=[heading heading(length(heading))];

for pp=1:length(tmotion);
   if tmotion(pp)==tgyro(qq),
      shead(pp)=sin(heading(qq)./180.*pi);
      chead(pp)=cos(heading(qq)./180.*pi);
   elseif tmotion(pp) > tgyro(qq),
      shead(pp)=(sin(heading(qq+1)./180.*pi)-sin(heading(qq)./180.*pi))/(tgyro(qq+1)-tgyro(qq))*(tmotion(pp)-tgyro(qq))+sin(heading(qq)./180.*pi);
      chead(pp)=(cos(heading(qq+1)./180.*pi)-cos(heading(qq)./180.*pi))/(tgyro(qq+1)-tgyro(qq))*(tmotion(pp)-tgyro(qq))+cos(heading(qq)./180.*pi);
      if tmotion(pp) > tgyro(qq+1), qq=qq+1;end;
   end;
end;

clear heading;
heading=atan2(shead,chead).*180./pi;
clear tgyro qq pp;
for ppp=1:length(heading);
   if heading(ppp)<0,
      heading(ppp)=heading(ppp)+360;
   end;
end;
heading=heading';

G=sqrt(mean(accx)^2+mean(accy)^2+mean(accz)^2);
ratex=detrend(ratex);ratey=detrend(ratey);ratez=detrend(ratez);

ssx=(1:length(accx)).*0;ssx(1)=su(1);
ssy=(1:length(accx)).*0;ssy(1)=sv(1);
ssz=(1:length(accx)).*0;ssz(1)=sw(1);
sst=(1:length(accx)).*0;sst(1)=Tsonic(1);
jj=2;

for ii=2:length(accx);
   while tsons(jj)<tmotion(ii) & jj < length(su);
      jj=jj+1;
   end;
   ssx(ii)=(su(jj)-su(jj-1))/(tsons(jj)-tsons(jj-1))*(tmotion(ii)-tsons(jj-1))+su(jj-1);
   ssy(ii)=(sv(jj)-sv(jj-1))/(tsons(jj)-tsons(jj-1))*(tmotion(ii)-tsons(jj-1))+sv(jj-1);
   ssz(ii)=(sw(jj)-sw(jj-1))/(tsons(jj)-tsons(jj-1))*(tmotion(ii)-tsons(jj-1))+sw(jj-1);
   sst(ii)=(Tsonic(jj)-Tsonic(jj-1))/(tsons(jj)-tsons(jj-1))*(tmotion(ii)-tsons(jj-1))+Tsonic(jj-1);
end;

clear su sv sw headr tsons flist tmotion Tsonic;
fclose('all');

maccx=mean(accx);
maccy=mean(accy);
maccz=mean(accz);

%disp('Complementary filtering')

Fc=1/30;                 %60 seconds
%Fc=1/300;                  %10 seconds

tau=1/(Fc*2*pi);bc=Fc/10;[bf af]=butter(1,bc);dt=1/20;
om_c=tan(dt/2/tau);tau_c=dt/2/om_c;
clear Fc;

%load c:\hare\fxyz.prn;
load /home/user1/fxyz.prn;

fx=fxyz(1);
fy=fxyz(2);
fz=fxyz(3);
px=fxyz(4);
py=fxyz(5);
pz=fxyz(6);

if fx==0 | isempty(fx),
   fx=maccy/G;
end;
if fy==0 | isempty(fy),
   fy=maccx/G;
end;
if fz==0 | isempty(fz),
   fz=mean(heading);
end;

thetarate=ratey*tau_c+accx/G;                              %pitchrate
theta=filter(bf,af,thetarate,maccx/G);
%theta=filter(bf,af,thetarate,fy);

phirate=ratex*tau_c+(accy/G)./cos(theta);
phi=filter(bf,af,phirate,maccy/G);
%phi=filter(bf,af,phirate,fx);

heading=heading/180*pi;
psirate=ratez*tau_c+heading;                              %no turns thru zero allowed
psi=filter(bf,af,psirate,mean(heading));
%psi=filter(bf,af,psirate,fz);

fx=filtic(bf,af,[phi(length(phi)) phi(length(phi)-1)],[phirate(length(phirate)) phirate(length(phirate)-1)]);
fy=filtic(bf,af,[theta(length(theta)) theta(length(theta)-1)],[thetarate(length(thetarate)) thetarate(length(thetarate)-1)]);
fz=filtic(bf,af,[psi(length(psi)) psi(length(psi)-1)],[psirate(length(psirate)) psirate(length(psirate)-1)]);

%Compute Omega X R for sonics

Rz=1.25;Rx=0;Ry=0;
OXRx=ratey*Rz-ratez*Ry;
OXRy=-ratex*Rz+ratez*Rx;
OXRz=+ratex*Ry-ratey*Rx;
OM=[OXRx OXRy OXRz];

%TOM=0.*OM';

clear ratex ratey ratez maccx maccy maccz thetarate phirate psirate OXRx;
clear OXRy OXRz Rx Ry Rz;

Plat=[accx'
      accy'
      accz'];

%clear accz accy;

Vel=[ssy              %Matrix defined as [V U W]
     ssx               %
     ssz];             %

Platacc=Plat;
Velrot=Vel;
Plat=Plat';
Vel=Vel';

%disp('Rotating coordinate system');

cos_theta=cos(theta);sin_theta=sin(theta);
cos_phi=cos(phi);sin_phi=sin(phi);
cos_psi=cos(psi);sin_psi=sin(psi);

clear theta;clear psi;clear phi;

%T11=cos_psi.*cos_theta;                           %JEH Masters
%T12=sin_psi.*cos_phi-cos_psi.*sin_theta.*sin_phi;
%T13=-sin_psi.*sin_phi-cos_psi.*sin_theta.*cos_phi;
%T21=-sin_psi.*cos_theta;
%T22=cos_psi.*cos_phi+sin_psi.*sin_theta.*sin_phi;
%T23=-cos_psi.*sin_phi+sin_psi.*sin_theta.*cos_phi;
%T31=sin_theta;
%T32=cos_theta.*sin_phi;
%T33=cos_theta.*cos_phi;

T11=cos_psi.*cos_theta;                           %JBE
T12=sin_psi.*cos_phi+cos_psi.*sin_theta.*sin_phi;
T13=-sin_psi.*sin_phi+cos_psi.*sin_theta.*cos_phi;
T21=-sin_psi.*cos_theta;
T22=cos_psi.*cos_phi-sin_psi.*sin_theta.*sin_phi;
T23=-cos_psi.*sin_phi-sin_psi.*sin_theta.*cos_phi;
T31=-sin_theta;
T32=cos_theta.*sin_phi;
T33=cos_theta.*cos_phi;

clear cos_theta;clear sin_theta;clear cos_phi;
clear sin_phi;clear cos_psi;clear sin_psi;
accx=accx';accy=accy';accz=accz';

%ssx=ssx';ssy=ssy';ssz=ssz';

T11=T11';T12=T12';T13=T13';
T21=T21';T22=T22';T23=T23';
T31=T31';T32=T32';T33=T33';

%for zz=1:length(accx);
%   T=[T11(zz) T12(zz) T13(zz)
%      T21(zz) T22(zz) T23(zz)
%      T31(zz) T32(zz) T33(zz)];
%   Platacc(:,zz)=T*Plat(zz,:)';     %Rotated accelerations
%    Velrot(:,zz)=T*Vel(zz,:)';      %Rotated sonic velocities
%    TOM(:,zz)=T*OM(zz,:)';          %Rotated  Om X R
   
Platacc(1,:) = T11.*accx + T12.*accy + T13.*accz;
Platacc(2,:) = T21.*accx + T22.*accy + T23.*accz;
Platacc(3,:) = T31.*accx + T32.*accy + T33.*accz;

Velrot(1,:) = T11.*ssx + T12.*ssy + T13.*ssz;
Velrot(2,:) = T21.*ssx + T22.*ssy + T23.*ssz;
Velrot(3,:) = T31.*ssx + T32.*ssy + T33.*ssz;

%end;

clear Plat Vel T OM ssx ssy ssz accy accz;

maccx=mean(Platacc(1,:));
maccy=mean(Platacc(2,:));
maccz=mean(Platacc(3,:));

Platacc(1,:)=Platacc(1,:)-maccx;
Platacc(2,:)=Platacc(2,:)-maccy;
Platacc(3,:)=Platacc(3,:)-maccz;

clear maccx maccy maccz;

Platvel=(1:length(accx)).*0;
Platvel=[Platvel
         Platvel
         Platvel];

Acclow=Platacc.*0;

%Fc=1/30;                 %60 seconds
Fc=1/120;                 %10 seconds
tau=1/(Fc*2*pi);bc=Fc/10;[bf af]=butter(1,bc);dt=1/20;
om_c=tan(dt/2/tau);tau_c=dt/2/om_c;
clear Fc;

Acclow(1,:)=filtfilt(bf,af,Platacc(1,:));
Acclow(2,:)=filtfilt(bf,af,Platacc(2,:));
Acclow(3,:)=filtfilt(bf,af,Platacc(3,:));

Platform=Platacc.*0;

for zz=1:length(accx);
   for kk=1:3;
      Platform(kk,zz)=Platacc(kk,zz)-Acclow(kk,zz);      %Highpassed Rotated Accelerations
   end;
end;

clear Platacc Acclow kk;

%disp('Integrating')

Fc=1/30;           %seconds
%Fc=1/150;                 %10 seconds
tau=1/(Fc*2*pi);bc=Fc/10;[bf af]=butter(1,bc);dt=1/20;
om_c=tan(dt/2/tau);tau_c=dt/2/om_c;
clear Fc;

Platrate(1,:)=Platform(1,:)*tau_c;
Platrate(2,:)=Platform(2,:)*tau_c;
Platrate(3,:)=Platform(3,:)*tau_c;

if isempty(px),
   px=0;
end;
if isempty(py),
   py=0;
end;
if isempty(pz),
   pz=0;
end;

[Platvel(1,:),Zf]=filter(bf,af,Platrate(1,:),px);
[Platvel(2,:),Zf]=filter(bf,af,Platrate(2,:),py);
[Platvel(3,:),Zf]=filter(bf,af,Platrate(3,:),pz);

px=filtic(bf,af,[Platvel(1,length(Platvel)) Platvel(1,length(Platvel)-1)],[Platrate(1,length(Platrate)) Platrate(1,length(Platrate)-1)]);
py=filtic(bf,af,[Platvel(2,length(Platvel)) Platvel(2,length(Platvel)-1)],[Platrate(2,length(Platrate)) Platrate(2,length(Platrate)-1)]);
pz=filtic(bf,af,[Platvel(3,length(Platvel)) Platvel(3,length(Platvel)-1)],[Platrate(3,length(Platrate)) Platrate(3,length(Platrate)-1)]);

fxyz=[fx fy fz px py pz];
save c:\hare\fxyz.prn fxyz -ascii

clear dt ii jj om_c zz tau_c bc G ans T11 T12 T13 T21 T22 T23 T31 T32;
clear T33 Bf0 Bf1 Af0 tau lv accx;
clear Platrate Platform;

Result=Velrot;

Result(1,:)=Velrot(1,:)-Platvel(1,:);
Result(2,:)=Velrot(2,:)-Platvel(2,:);
Result(3,1:length(Result(3,:))-5)=Velrot(3,1:length(Result(3,:))-5)+Platvel(3,6:length(Result(3,:)));

Result(3,length(Result(1,:))-4:length(Result(1,:)))=Velrot(3,length(Result(1,:))-4:length(Result(1,:)))+Platvel(3,length(Result(1,:))-4:length(Result(1,:)));

figure;plot(1:length(Result),Result(1,:),1:length(Velrot),Velrot(1,:),1:length(Platvel),Platvel(1,:));title('axis 1')
figure;plot(1:length(Result),Result(2,:),1:length(Velrot),Velrot(2,:),1:length(Platvel),Platvel(2,:));title('axis 2')
figure;plot(1:length(Result),Result(3,:),1:length(Velrot),Velrot(3,:),1:length(Platvel),Platvel(3,:));title('axis vert')

[Sv1,FF]=psd(detrend(Velrot(1,:)),65536,20);
[Sv2,FF]=psd(detrend(Velrot(2,:)),65536,20);
[Sv3,FF]=psd(detrend(Velrot(3,:)),65536,20);
[Sp1,FF]=psd(detrend(Platvel(1,:)),65536,20);
[Sp2,FF]=psd(detrend(Platvel(2,:)),65536,20);
[Sp3,FF]=psd(detrend(Platvel(3,:)),65536,20);
[Sr1,FF]=psd(detrend(Result(1,:)),65536,20);
[Sr2,FF]=psd(detrend(Result(2,:)),65536,20);
[Sr3,FF]=psd(detrend(Result(3,:)),65536,20);

figure;loglog(FF,Sv1.*FF,FF,Sp1.*FF,FF,Sr1.*FF)
figure;loglog(FF,Sv2.*FF,FF,Sp2.*FF,FF,Sr2.*FF)
figure;loglog(FF,Sv3.*FF,FF,Sp3.*FF,FF,Sr3.*FF)

pause
%[A,Fa]=psd(Result(1,:),65536,20);[sA,sFa]=specsmoo(A,20);
%[B,Fb]=psd(Result(2,:),65536,20);[sB,sFb]=specsmoo(B,20);
%[C,Fc]=psd(Result(3,:),65536,20);[sC,sFc]=specsmoo(C,20);
%figure;loglog(sFa,sA.*sFa,'r-',sFb,sB.*sFb,'g-',sFc,sC.*sFc,'b-');title('Result')
%[A,Fa]=psd(Platvel(1,:),65536,20);[sA,sFa]=specsmoo(A,20);
%[B,Fb]=psd(Platvel(2,:),65536,20);[sB,sFb]=specsmoo(B,20);
%[C,Fc]=psd(Platvel(3,:),65536,20);[sC,sFc]=specsmoo(C,20);
%figure;loglog(sFa,sA.*sFa,'r-',sFb,sB.*sFb,'g-',sFc,sC.*sFc,'b-');title('Platvel')
%[A,Fa]=psd(Velrot(1,:),65536,20);[sA,sFa]=specsmoo(A,20);
%[B,Fb]=psd(Velrot(2,:),65536,20);[sB,sFb]=specsmoo(B,20);
%[C,Fc]=psd(Velrot(3,:),65536,20);[sC,sFc]=specsmoo(C,20);
%figure;loglog(sFa,sA.*sFa,'r-',sFb,sB.*sFb,'g-',sFc,sC.*sFc,'b-');title('Velrot')

%rotcorr;           %converts [u,v,w] to [Ustream, Vstream, Wstream]
%clear Uu Vv Ww ans tm u v w;
%clear Umean Vmean Wmean ss ut uu uv uw vt vv vw wt ww;
%clear yy alpha beta Ub Vb Wb Sb FF Svu1 Spu1 Sru1 Svv1 Spv1 Srv1 Svw1 Spw1 Srw1
%b4a='Before';
%sonicplots;
%

u=Result(1,:);
v=Result(2,:);
w=Result(3,:);

rotcorr;           %converts [u,v,w] to [Ustream, Vstream, Wstream]
clear Uu Vv Ww ans tm u v w;
clear Umean Vmean Wmean ss ut uu uv uw vt vv vw wt ww;
clear yy alpha beta Ub Vb Wb Sb FF Svu1 Spu1 Sru1 Svv1 Spv1 Srv1 Svw1 Spw1 Srw1

firstq=fix(length(Wstream)/6*1);
seconq=fix(length(Wstream)/6*2);
thirdq=fix(length(Wstream)/6*3);
fourtq=fix(length(Wstream)/6*4);
fifthq=fix(length(Wstream)/6*5);
sixthq=fix(length(Wstream)/6*6);
Us=0;
zi=600;P=1013;
zu=19.2;zt=19.2;zq=17.5;

if length(ship_Tsea)<60,
   ship_Tsea=[ship_Tsea (1:60-length(ship_Tsea)).*0+ship_Tsea(length(ship_Tsea))];
end;   
if length(Tsea_snake)<60,
   Tsea_snake=[Tsea_snake (1:60-length(Tsea_snake)).*0+Tsea_snake(length(Tsea_snake))];
end;   
if length(ship_Tair)<60,
   ship_Tair=[ship_Tair (1:60-length(ship_Tair)).*0+ship_Tair(length(ship_Tair))];
end;   
if length(ship_hum)<60,
   ship_hum=[ship_hum (1:60-length(ship_hum)).*0+ship_hum(length(ship_hum))];
end;   
if length(rain)<60,
   rain=[rain (1:60-length(rain)).*0+rain(length(rain))];
end;   

wu_1=cov(Ustream(1:firstq),Wstream(1:firstq));
wv_1=cov(Vstream(1:firstq),Wstream(1:firstq));
wt_1=cov(Wstream(1:firstq),sst(1:firstq));
Ubar_1=sqrt(mean(Ustream(1:firstq)).^2+mean(Vstream(1:firstq)).^2+mean(Wstream(1:firstq)).^2);
Tbar_1=mean(ship_Tair(1:40));
rho_1=P*100./(273.15+Tbar_1)./287;
wa_1=cov(Wstream(1:firstq),abshum(1:firstq)./rho_1);
wh_1=cov(Wstream(1:firstq),H2O(1:firstq)./rho_1);
wc_1=cov(Wstream(1:firstq),CO2(1:firstq));
%Qbar_1=0.622*mean(ship_hum(1:40))/100*6.112*exp((17.67*Tbar_1)/(Tbar_1+243.5))/1013.25*1000;
Qbar_1=0.622*mean(RHv(1:10))/100/1.2*6.112*exp((17.67*Tbar_1)/(Tbar_1+243.5))/1013.25*1000;
Tsea_1=mean(Tsea_snake(1:10));
Esat=6.112.*exp(17.67.*Tsea_1./(Tsea_1+243.5))*0.98;
Qsea_1=Esat*0.622/(P-0.378*Esat)*1000;
Ss_1=std(shead(1:firstq));Sc_1=std(chead(1:firstq));
Sp_1=std(Platvel(1:firstq));
rain_1=mean(rain(1:10));
Rl_1=(mean(LongIR1(1:10))+mean(LongIR2(1:10)))/2;
Rs_1=(mean(Solar1(1:10))+mean(Solar2(1:10)))/2;
x_1=[Ubar_1 Us Tsea_1 Tbar_1 Qsea_1 Qbar_1 Rs_1 Rl_1 rain_1 zi P zu zt zq];
y_1=coul_lkb(x_1);

wu_2=cov(Ustream(firstq+1:seconq),Wstream(firstq+1:seconq));
wv_2=cov(Vstream(firstq+1:seconq),Wstream(firstq+1:seconq));
wt_2=cov(Wstream(firstq+1:seconq),sst(firstq+1:seconq));
Ubar_2=sqrt(mean(Ustream(firstq+1:seconq)).^2+mean(Vstream(firstq+1:seconq)).^2+mean(Wstream(firstq+1:seconq)).^2);
Tbar_2=mean(ship_Tair(41:80));
rho_2=P*100./(273.15+Tbar_2)./287;
wa_2=cov(Wstream(firstq+1:seconq),abshum(firstq+1:seconq)./rho_2);
wh_2=cov(Wstream(firstq+1:seconq),H2O(firstq+1:seconq)./rho_2);
wc_2=cov(Wstream(firstq+1:seconq),CO2(firstq+1:seconq));
%Qbar_2=0.622*mean(ship_hum(41:80))/100*6.112*exp((17.67*Tbar_2)/(Tbar_2+243.5))/1013.25*1000;
Qbar_2=0.622*mean(RHv(11:20))/100/1.2*6.112*exp((17.67*Tbar_2)/(Tbar_2+243.5))/1013.25*1000;
Tsea_2=mean(Tsea_snake(11:20));
Esat=6.112.*exp(17.67.*Tsea_2./(Tsea_2+243.5))*0.98;
Qsea_2=Esat*0.622/(P-0.378*Esat)*1000;
Ss_2=std(shead(firstq+1:seconq));Sc_2=std(chead(firstq+1:seconq));
Sp_2=std(Platvel(firstq+1:seconq));
rain_2=mean(rain(11:20));
Rl_2=(mean(LongIR1(11:20))+mean(LongIR2(11:20)))/2;
Rs_2=(mean(Solar1(11:20))+mean(Solar2(11:20)))/2;
x_2=[Ubar_2 Us Tsea_2 Tbar_2 Qsea_2 Qbar_2 Rs_2 Rl_2 rain_2 zi P zu zt zq];
y_2=coul_lkb(x_2);

wu_3=cov(Ustream(seconq+1:thirdq),Wstream(seconq+1:thirdq));
wv_3=cov(Vstream(seconq+1:thirdq),Wstream(seconq+1:thirdq));
wt_3=cov(Wstream(seconq+1:thirdq),sst(seconq+1:thirdq));
Ubar_3=sqrt(mean(Ustream(seconq+1:thirdq)).^2+mean(Vstream(seconq+1:thirdq)).^2+mean(Wstream(seconq+1:thirdq)).^2);
Tbar_3=mean(ship_Tair(81:120));
rho_3=P*100./(273.15+Tbar_3)./287;
wa_3=cov(Wstream(seconq+1:thirdq),abshum(seconq+1:thirdq)./rho_3);
wh_3=cov(Wstream(seconq+1:thirdq),H2O(seconq+1:thirdq)./rho_3);
wc_3=cov(Wstream(seconq+1:thirdq),CO2(seconq+1:thirdq));
%Qbar_3=0.622*mean(ship_hum(81:120))/100*6.112*exp((17.67*Tbar_3)/(Tbar_3+243.5))/1013.25*1000;
Qbar_3=0.622*mean(RHv(21:30))/100/1.2*6.112*exp((17.67*Tbar_3)/(Tbar_3+243.5))/1013.25*1000;
Tsea_3=mean(Tsea_snake(21:30));
Esat=6.112.*exp(17.67.*Tsea_3./(Tsea_3+243.5))*0.98;
Qsea_3=Esat*0.622/(P-0.378*Esat)*1000;
Ss_3=std(shead(seconq+1:thirdq));Sc_3=std(chead(seconq+1:thirdq));
Sp_3=std(Platvel(seconq+1:thirdq));
rain_3=mean(rain(21:30));
Rl_3=(mean(LongIR1(21:30))+mean(LongIR2(21:30)))/2;
Rs_3=(mean(Solar1(21:30))+mean(Solar2(21:30)))/2;
x_3=[Ubar_3 Us Tsea_3 Tbar_3 Qsea_3 Qbar_3 Rs_3 Rl_3 rain_3 zi P zu zt zq];
y_3=coul_lkb(x_3);

wu_4=cov(Ustream(thirdq+1:fourtq),Wstream(thirdq+1:fourtq));
wv_4=cov(Vstream(thirdq+1:fourtq),Wstream(thirdq+1:fourtq));
wt_4=cov(Wstream(thirdq+1:fourtq),sst(thirdq+1:fourtq));
Ubar_4=sqrt(mean(Ustream(thirdq+1:fourtq)).^2+mean(Vstream(thirdq+1:fourtq)).^2+mean(Wstream(thirdq+1:fourtq)).^2);
Tbar_4=mean(ship_Tair(121:160));
rho_4=P*100./(273.15+Tbar_4)./287;
wa_4=cov(Wstream(thirdq+1:fourtq),abshum(thirdq+1:fourtq)./rho_4);
wh_4=cov(Wstream(thirdq+1:fourtq),H2O(thirdq+1:fourtq)./rho_4);
wc_4=cov(Wstream(thirdq+1:fourtq),CO2(thirdq+1:fourtq));
%Qbar_4=0.622*mean(ship_hum(121:160))/100*6.112*exp((17.67*Tbar_4)/(Tbar_4+243.5))/1013.25*1000;
Qbar_4=0.622*mean(RHv(31:40))/100/1.2*6.112*exp((17.67*Tbar_4)/(Tbar_4+243.5))/1013.25*1000;
Tsea_4=mean(Tsea_snake(31:40));
Esat=6.112.*exp(17.67.*Tsea_4./(Tsea_4+243.5))*0.98;
Qsea_4=Esat*0.622/(P-0.378*Esat)*1000;
Ss_4=std(shead(thirdq+1:fourtq));Sc_4=std(chead(thirdq+1:fourtq));
Sp_4=std(Platvel(thirdq+1:fourtq));
rain_4=mean(rain(31:40));
Rl_4=(mean(LongIR1(31:40))+mean(LongIR2(31:40)))/2;
Rs_4=(mean(Solar1(31:40))+mean(Solar2(31:40)))/2;
x_4=[Ubar_4 Us Tsea_4 Tbar_4 Qsea_4 Qbar_4 Rs_4 Rl_4 rain_4 zi P zu zt zq];
y_4=coul_lkb(x_4);

wu_5=cov(Ustream(fourtq+1:fifthq),Wstream(fourtq+1:fifthq));
wv_5=cov(Vstream(fourtq+1:fifthq),Wstream(fourtq+1:fifthq));
wt_5=cov(Wstream(fourtq+1:fifthq),sst(fourtq+1:fifthq));
Ubar_5=sqrt(mean(Ustream(fourtq+1:fifthq)).^2+mean(Vstream(fourtq+1:fifthq)).^2+mean(Wstream(fourtq+1:fifthq)).^2);
Tbar_5=mean(ship_Tair(161:200));
rho_5=P*100./(273.15+Tbar_5)./287;
wa_5=cov(Wstream(fourtq+1:fifthq),abshum(fourtq+1:fifthq)./rho_5);
wh_5=cov(Wstream(fourtq+1:fifthq),H2O(fourtq+1:fifthq)./rho_5);
wc_5=cov(Wstream(fourtq+1:fifthq),CO2(fourtq+1:fifthq));
%Qbar_5=0.622*mean(ship_hum(161:200))/100*6.112*exp((17.67*Tbar_5)/(Tbar_5+243.5))/1013.25*1000;
Qbar_5=0.622*mean(RHv(41:50))/100/1.2*6.112*exp((17.67*Tbar_5)/(Tbar_5+243.5))/1013.25*1000;
Tsea_5=mean(Tsea_snake(41:50));
Esat=6.112.*exp(17.67.*Tsea_5./(Tsea_5+243.5))*0.98;
Qsea_5=Esat*0.622/(P-0.378*Esat)*1000;
Ss_5=std(shead(fourtq+1:fifthq));Sc_5=std(chead(fourtq+1:fifthq));
Sp_5=std(Platvel(fourtq+1:fifthq));
rain_5=mean(rain(41:50));
Rl_5=(mean(LongIR1(41:50))+mean(LongIR2(41:50)))/2;
Rs_5=(mean(Solar1(41:50))+mean(Solar2(41:50)))/2;
x_5=[Ubar_5 Us Tsea_5 Tbar_5 Qsea_5 Qbar_5 Rs_5 Rl_5 rain_5 zi P zu zt zq];
y_5=coul_lkb(x_5);

wu_6=cov(Ustream(fifthq+1:sixthq),Wstream(fifthq+1:sixthq));
wv_6=cov(Vstream(fifthq+1:sixthq),Wstream(fifthq+1:sixthq));
wt_6=cov(Wstream(fifthq+1:sixthq),sst(fifthq+1:sixthq));
Ubar_6=sqrt(mean(Ustream(fifthq+1:sixthq)).^2+mean(Vstream(fifthq+1:sixthq)).^2+mean(Wstream(fifthq+1:sixthq)).^2);
Tbar_6=mean(ship_Tair(201:length(ship_Tair)));
rho_6=P*100./(273.15+Tbar_6)./287;
wa_6=cov(Wstream(fifthq+1:sixthq),abshum(fifthq+1:sixthq)./rho_6);
wh_6=cov(Wstream(fifthq+1:sixthq),H2O(fifthq+1:sixthq)./rho_6);
wc_6=cov(Wstream(fifthq+1:sixthq),CO2(fifthq+1:sixthq));
%Qbar_6=0.622*mean(ship_hum(201:length(RHv)))/100*6.112*exp((17.67*Tbar_6)/(Tbar_6+243.5))/1013.25*1000;
Qbar_6=0.622*mean(RHv(51:length(RHv)))/100/1.2*6.112*exp((17.67*Tbar_6)/(Tbar_6+243.5))/1013.25*1000;
Tsea_6=mean(Tsea_snake(51:length(Tsea_snake)));
Esat=6.112.*exp(17.67.*Tsea_6./(Tsea_6+243.5))*0.98;
Qsea_6=Esat*0.622/(P-0.378*Esat)*1000;
Ss_6=std(shead(fifthq+1:sixthq));Sc_6=std(chead(fifthq+1:sixthq));
Sp_6=std(Platvel(fifthq+1:sixthq));
rain_6=mean(rain(51:60));
Rl_6=(mean(LongIR1(51:60))+mean(LongIR2(51:60)))/2;
Rs_6=(mean(Solar1(51:60))+mean(Solar2(51:60)))/2;
x_6=[Ubar_6 Us Tsea_6 Tbar_6 Qsea_6 Qbar_6 Rs_6 Rl_6 rain_6 zi P zu zt zq];
y_6=coul_lkb(x_6);

clear H2O LongIR1 LongIR2 P Platvel LongIR1 LongIR2 RHv Result Solar1 Solar2
clear zi zq zt zu shead rain ppp px py pz heading fz
clear fx fy fxyz dfx1 dfx2 chead bf af Zf
clear Velrot Tsea_snake Tairv Esat ship_Tair ship_Tsea
clear ship_cog ship_gyro RHv ship_sog gpslat_sign gpslon_sign gpshour
clear gpsmin gpssec gpshead gpstime

%y_y=[y_1 Ubar_1 Tbar_1 Tsea_1 Qbar_1 wu_1(1,2) wv_1(1,2) wt_1(1,2) wa_1(1,2) wh_1(1,2) rho_1 stdC_1 stdA_1 medC_1 medA_1 Ss_1 Sc_1 Sp_1 rain_1
%     y_2 Ubar_2 Tbar_2 Tsea_2 Qbar_2 wu_2(1,2) wv_2(1,2) wt_2(1,2) wa_2(1,2) wh_2(1,2) rho_2 stdC_2 stdA_2 medC_2 medA_2 Ss_2 Sc_2 Sp_2 rain_2
%     y_3 Ubar_3 Tbar_3 Tsea_3 Qbar_3 wu_3(1,2) wv_3(1,2) wt_3(1,2) wa_3(1,2) wh_3(1,2) rho_3 stdC_3 stdA_3 medC_3 medA_3 Ss_3 Sc_3 Sp_3 rain_3
%     y_4 Ubar_4 Tbar_4 Tsea_4 Qbar_4 wu_4(1,2) wv_4(1,2) wt_4(1,2) wa_4(1,2) wh_4(1,2) rho_4 stdC_4 stdA_4 medC_4 medA_4 Ss_4 Sc_4 Sp_4 rain_4
%     y_5 Ubar_5 Tbar_5 Tsea_5 Qbar_5 wu_5(1,2) wv_5(1,2) wt_5(1,2) wa_5(1,2) wh_5(1,2) rho_5 stdC_5 stdA_5 medC_5 medA_5 Ss_5 Sc_5 Sp_5 rain_5
%     y_6 Ubar_6 Tbar_6 Tsea_6 Qbar_6 wu_6(1,2) wv_6(1,2) wt_6(1,2) wa_6(1,2) wh_6(1,2) rho_6 stdC_6 stdA_6 medC_6 medA_6 Ss_6 Sc_6 Sp_6 rain_6];

%clear y_1 y_2 y_3 y_4 y_5 y_6 Ubar_1 Ubar_2 Ubar_3 Ubar_4 Ubar_5 Ubar_6 Tbar_1 Tbar_2 Tbar_3 Tbar_4 Tbar_5 Tbar_6
%clear Tsea_1 Tsea_2 Tsea_3 Tsea_4 Tsea_5 Tsea_6 Qbar_1 Qbar_2 Qbar_3 Qbar_4 Qbar_5 Qbar_6
%clear rho_1 rho_2 rho_3 rho_4 rho_5 rho_6 zi zq zt zu
%   y=[hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF hl_webb Cd Ch Ce Cdn_10 Chn_10 Cen_10];
%       ********** output data  ***************
%	01 hsb=			sensible heat flux (w/m^2)
%	02 hlb=			latent heat flux (w/m^2)
%	03 tau=			stress (nt/m^2)
%	04 zo=			velocity roughness length (m)
%	05 zot			temperature roughness length (m)
%	06 zoq=			moisture roughness length (m)
%	07 L=			Monin_Obukhov stability length
%	08 usr=			turbulent friction velocity (m/s), including gustiness
%	09 tsr			temperature scaling parameter (K)
%	10 qsr			humidity scaling parameter (g/g)
%	11 dter=		cool skin temperature depression (K)
%	12 dqer=		cool skin humidity depression (g/g)
%	13 tkt=			cool skin thickness (m)
%	14 RF=			rain heat flux(w/m^2)
%	15 hl_webb=		webb addition to latent flux (w/m^2)
%	16 Cd=			velocity drag coefficient at zu, referenced to u
%	17 Ch=			heat transfer coefficient at zt
%	18 Ce=			moisture transfer coefficient at zq
%	19 Cdn_10=		10-m velocity drag coeeficient, including gustiness
%	20 Chn_10=		10-m heat transfer coeeficient, including gustiness
%	21 Cen_10=		10-m humidity transfer coeeficient, including gustiness

%ss=['save ',sfl,' y_y -ascii -double'];disp(ss)
%eval(ss);
