disp('fastex07.m')
%  09/11/98   Multiply FASTEX abshum by 20% according to analysis of data set

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 Mcog Msog Mhead Ue Un
heading=atan2(shead,chead).*180./pi;
heading = unwrap(heading);
heading=heading';

clear tgyro qq pp;

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 tsons headr tmotion Tsonic
clear tsons headr tmotion Tsonic

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

%disp('Complementary filtering')

Fc=1/60;                 %60 seconds
bc=Fc/10;
[bf af]=butter(1,bc);

tau=1/(Fc*2*pi);
dt=1/20;
om_c=tan(dt/2/tau);tau_c=dt/2/om_c;

wp=1/120/(20/2);
ws=1/60/(20/2);   
[n,wn]=buttord(wp,ws,1,20);

[bhigh,ahigh]=butter(n,wn,'high');

clear Fc wn wp ws n

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)/180*pi;
end;

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

%thetaslow = atan2(-accx,G) - filtfilt(bhigh,ahigh,atan2(-accx,G));
%theta = thetaslow + filtfilt(bhigh,ahigh,(cumsum(ratey) -.5*ratey -.5*ratey(1))/20);

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

%phislow = atan2(accy,G) - filtfilt(bhigh,ahigh,atan2(accy,G));
%phi = phislow + filtfilt(bhigh,ahigh,(cumsum(ratex) - .5*ratex - .5*ratex(1))/20);

heading=heading/180*pi;                     %convert to radians
psirate=ratez*tau_c-heading;                              %+ to - 09/15/98
%psi=filter(bf,af,psirate,mean(heading));
psi=filter(bf,af,psirate,fz);

%psislow = (-heading - filtfilt(bhigh,ahigh,-heading));
%psi = psislow + filtfilt(bhigh,ahigh,(cumsum(ratez) -.5*ratez - .5*ratez(1))/20);

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];
OM=OM';

%TOM=0.*OM';

clear maccx maccy maccz ahigh bhigh
clear thetarate phirate psirate OXRx OXRy OXRz

%clear OXRy OXRz Rx Ry Rz;

Plat=[accx'
      accy'
      accz'];

Vel=[ssx              %
     ssy              %
     ssz];            %

Platacc=Plat;
Velrot=Vel;

%*****************************************************************
%            Start Edson correction
%*****************************************************************

deg_rate=[ratex' 
   ratey' 
   ratez'];

clear ratex ratey ratez

R=[Rx
   Ry
   Rz];

[bhi,ahi,bf2,af2,tau2,spinup2]=filtcoef(20,1/60);

[euler,dr]=angles(ahi,bhi,20,Plat,deg_rate,heading');
[uvwplat,xyzplat]=accels(bhi,ahi,20,Plat,euler);clear xyzplat

[uvw]=sonic(Vel,dr,euler,uvwplat,R);

Jcovtest=cov(uvw(3,:),uvwplat(3,:));
Jcovtest=Jcovtest(1,2);
lagflag=0;

for lag=1:8;
   Jcovtestlag=cov(uvw(3,1:length(uvw)-lag),uvwplat(3,lag+1:length(uvw)));
   %figure;plot([uvw(3,1:length(uvw)-lag)' -uvwplat(3,lag+1:length(uvw))']);axis([20000 22000 -5 5]);
   Jcovtestlag=Jcovtestlag(1,2);
   if Jcovtestlag < 0 & Jcovtest < 0,
      if abs(Jcovtestlag) > abs(Jcovtest),
         Jcovtest=Jcovtestlag;
         lagflag=lag;
      end;
   elseif Jcovtestlag > 0 & Jcovtest > 0,
      if abs(Jcovtestlag) < abs(Jcovtest),
         Jcovtest=Jcovtestlag;
         lagflag=lag;
      end;
   elseif Jcovtestlag < 0 & Jcovtest > 0,
         Jcovtest=Jcovtestlag;
         lagflag=lag;
   end;
end;
for lag=-1:-1:-8;
   Jcovtestlag=cov(uvw(3,-lag+1:length(uvw)),uvwplat(3,1:length(uvw)+lag));
   %figure;plot([uvw(3,-lag+1:length(uvw))' -uvwplat(3,1:length(uvw)+lag)']);axis([20000 22000 -5 5]);
   Jcovtestlag=Jcovtestlag(1,2);
   if Jcovtestlag < 0 & Jcovtest < 0,
      if abs(Jcovtestlag) > abs(Jcovtest),
         Jcovtest=Jcovtestlag;
         lagflag=lag;
      end;
   elseif Jcovtestlag > 0 & Jcovtest > 0,
      if abs(Jcovtestlag) < abs(Jcovtest),
         Jcovtest=Jcovtestlag;
         lagflag=lag;
      end;
   elseif Jcovtestlag < 0 & Jcovtest > 0,
         Jcovtest=Jcovtestlag;
         lagflag=lag;
   end;
end;

uvw(1,:)=uvw(1,:) + uvwplat(1,:) ;
uvw(2,:)=uvw(2,:) + uvwplat(2,:) ;

if lagflag > 0,
   uvw(3,1:length(uvw)-lagflag)=uvw(3,1:length(uvw)-lagflag) + uvwplat(3,lagflag+1:length(uvw));
   uvw(3,length(uvw)-lagflag+1:length(uvw))= uvw(3,length(uvw)-lagflag+1:length(uvw)) + uvwplat(3,1:lagflag);
else,
   uvw(3,-lagflag+1:length(uvw)) = uvw(3,-lagflag+1:length(uvw)) + uvwplat(3,1:length(uvw)+lagflag);
   uvw(3,1:-lagflag) = uvw(3,1:-lagflag) + uvwplat(3,length(uvw)+lagflag+1:length(uvw));
end;   

%*****************************************************************
%            End Edson correction
%*****************************************************************

%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 psi phi R Rx Ry Rz Vel af2 ahi bf2 bhi deg_rate dr euler

%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 sin_theta cos_phi sin_phi cos_psi sin_psi;

accx=accx';accy=accy';accz=accz';
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
%end;

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+OM(1,:)) + T12.*(ssy+OM(2,:)) + T13.*(ssz+OM(3,:));
Velrot(2,:) = T21.*(ssx+OM(1,:)) + T22.*(ssy+OM(2,:)) + T23.*(ssz+OM(3,:));
Velrot(3,:) = T31.*(ssx+OM(1,:)) + T32.*(ssy+OM(2,:)) + T33.*(ssz+OM(3,:));

clear Plat 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/60;                 %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/60;           %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;

Jcovtest=cov(Velrot(3,:),Platvel(3,:));
Jcovtest=Jcovtest(1,2);
lagflag=0;

for lag=1:8;
   Jcovtestlag=cov(Velrot(3,1:length(Velrot)-lag),Platvel(3,lag+1:length(Velrot)));
   %figure;plot([Velrot(3,1:length(Velrot)-lag)' - Platvel(3,lag+1:length(Velrot))']);axis([20000 22000 -5 5]);
   Jcovtestlag=Jcovtestlag(1,2);
   if Jcovtestlag < 0 & Jcovtest < 0,
      if abs(Jcovtestlag) > abs(Jcovtest),
         Jcovtest=Jcovtestlag;
         lagflag=lag;
      end;
   elseif Jcovtestlag > 0 & Jcovtest > 0,
      if abs(Jcovtestlag) < abs(Jcovtest),
         Jcovtest=Jcovtestlag;
         lagflag=lag;
      end;
   elseif Jcovtestlag < 0 & Jcovtest > 0,
         Jcovtest=Jcovtestlag;
         lagflag=lag;
   end;
end;
for lag=-1:-1:-8;
   Jcovtestlag=cov(Velrot(3,-lag+1:length(Velrot)),Platvel(3,1:length(Velrot)+lag));
   %figure;plot([Velrot(3,-lag+1:length(Velrot))' -Platvel(3,1:length(Velrot)+lag)']);axis([20000 22000 -5 5]);
   Jcovtestlag=Jcovtestlag(1,2);
   if Jcovtestlag < 0 & Jcovtest < 0,
      if abs(Jcovtestlag) > abs(Jcovtest),
         Jcovtest=Jcovtestlag;
         lagflag=lag;
      end;
   elseif Jcovtestlag > 0 & Jcovtest > 0,
      if abs(Jcovtestlag) < abs(Jcovtest),
         Jcovtest=Jcovtestlag;
         lagflag=lag;
      end;
   elseif Jcovtestlag < 0 & Jcovtest > 0,
         Jcovtest=Jcovtestlag;
         lagflag=lag;
   end;
end;

clear Jcovtest Jcovtestlag lag

Result(1,:)=Velrot(1,:)+Platvel(1,:);
Result(2,:)=Velrot(2,:)+Platvel(2,:);

if lagflag > 0,
   Result(3,1:length(Velrot)-lagflag)=Velrot(3,1:length(Velrot)-lagflag) + Platvel(3,lagflag+1:length(Velrot));
   Result(3,length(Velrot)-lagflag+1:length(Velrot))= Velrot(3,length(Velrot)-lagflag+1:length(Velrot)) + Platvel(3,1:lagflag);
else,
   Result(3,-lagflag+1:length(Velrot)) = Velrot(3,-lagflag+1:length(Velrot)) + Platvel(3,1:length(Velrot)+lagflag);
   Result(3,1:-lagflag) = Velrot(3,1:-lagflag) + Platvel(3,length(Velrot)+lagflag+1:length(Velrot));
end;

clear lagflag

%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,:)));

[Ustream,Vstream,Wstream,alpha,beta]=rotcorr(Result(1,:),Result(2,:),Result(3,:));           %converts [u,v,w] to [Ustream, Vstream, Wstream] for one hour

clear Uu Vv Ww ans tm
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

if length(ship_Tsea)<240,
   ship_Tsea=[ship_Tsea (1:240-length(ship_Tsea)).*0+ship_Tsea(length(ship_Tsea))];
end;   
if length(ship_wind)<240,
   ship_wind=[ship_wind (1:240-length(ship_wind)).*0+ship_wind(length(ship_wind))];
end;   
if length(ship_dir)<240,
   ship_dir=[ship_dir (1:240-length(ship_dir)).*0+ship_dir(length(ship_dir))];
end;   
if length(ship_Tair)<240,
   ship_Tair=[ship_Tair (1:240-length(ship_Tair)).*0+ship_Tair(length(ship_Tair))];
end;   
if length(ship_hum)<240,
   ship_hum=[ship_hum (1:240-length(ship_hum)).*0+ship_hum(length(ship_hum))];
end;   
if length(ship_gyro)<240,
   ship_gyro=[ship_gyro (1:240-length(ship_gyro)).*0+ship_gyro(length(ship_gyro))];
end;   
if length(rain)<60,
   rain=[rain (1:60-length(rain)).*0+rain(length(rain))];
end;
if length(RHv)<60,
   RHv=[RHv (1:60-length(RHv)).*0+RHv(length(RHv))];
end;
if length(Tairv)<60,
   Tairv=[Tairv (1:60-length(Tairv)).*0+Tairv(length(Tairv))];
end;
if length(LongIR1)<60,
   LongIR1=[LongIR1 (1:60-length(LongIR1)).*0+LongIR1(length(LongIR1))];
end;
if length(LongIR2)<60,
   LongIR2=[LongIR2 (1:60-length(LongIR2)).*0+LongIR2(length(LongIR2))];
end;
if length(Solar1)<60,
   Solar1=[Solar1 (1:60-length(Solar1)).*0+Solar1(length(Solar1))];
end;
if length(Solar2)<60,
   Solar2=[Solar2 (1:60-length(Solar2)).*0+Solar2(length(Solar2))];
end;
if length(Tsea_snake)<60,
   Tsea_snake=[Tsea_snake (1:60-length(Tsea_snake)).*0+Tsea_snake(length(Tsea_snake))];
end;   
if length(abshum)<length(Wstream),
   abshum=[abshum (1:length(Wstream)-length(abshum)).*0+abshum(length(abshum))];
end;   
if length(H2O)<length(Wstream),
   H2O=[H2O (1:length(Wstream)-length(H2O)).*0+H2O(length(H2O))];
end;   
if length(Quw)<length(Wstream),
   Quw=[Quw (1:length(Wstream)-length(Quw)).*0+Quw(length(Quw))];
end;   
if length(CO2)<length(Wstream),
   CO2=[CO2 (1:length(Wstream)-length(CO2)).*0+CO2(length(CO2))];
end;   
if length(Oabsorb)<length(Wstream),
   Oabsorb=[Oabsorb (1:length(Wstream)-length(Oabsorb)).*0+Oabsorb(length(Oabsorb))];
end;   
if length(Oclear)<length(Wstream),
   Oclear=[Oclear (1:length(Wstream)-length(Oclear)).*0+Oclear(length(Oclear))];
end;   
if length(Tdry)<length(Wstream),
   Tdry=[Tdry (1:length(Wstream)-length(Tdry)).*0+Tdry(length(Tdry))];
end;   

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

firstr=fix(length(ship_sog)/6*1);
seconr=fix(length(ship_sog)/6*2);
thirdr=fix(length(ship_sog)/6*3);
fourtr=fix(length(ship_sog)/6*4);
fifthr=fix(length(ship_sog)/6*5);
sixthr=fix(length(ship_sog)/6*6);

firsts=fix(length(LongIR1)/6*1);
secons=fix(length(LongIR1)/6*2);
thirds=fix(length(LongIR1)/6*3);
fourts=fix(length(LongIR1)/6*4);
fifths=fix(length(LongIR1)/6*5);
sixths=fix(length(LongIR1)/6*6);

firstt=fix(length(gpshead)/6*1);
secont=fix(length(gpshead)/6*2);
thirdt=fix(length(gpshead)/6*3);
fourtt=fix(length(gpshead)/6*4);
fiftht=fix(length(gpshead)/6*5);
sixtht=fix(length(gpshead)/6*6);

firstu=fix(length(su)/6*1);
seconu=fix(length(su)/6*2);
thirdu=fix(length(su)/6*3);
fourtu=fix(length(su)/6*4);
fifthu=fix(length(su)/6*5);
sixthu=fix(length(su)/6*6);

Us=0;
zi=600;P=1013;
zu=19.2;zt=19.2;zq=17.5;


%**************************************!!!!!!!!!!!!!!!!!!!!!!!!
%abshum=abshum.*1.2;        %Gain Correction Only
%**************************************!!!!!!!!!!!!!!!!!!!!!!!!


%*********************************************************************************
%       Compute products
%*********************************************************************************


%*************************Minute 0-10********************************************

mean_sog_1=mean(ship_sog(1:firstr));
if mean_sog_1 < 1.5, mean_sog_1=0;end;
std_sog_1=std(ship_sog(1:firstr));
min_sog_1=min(ship_sog(1:firstr));
max_sog_1=max(ship_sog(1:firstr));

mean_cog_1=atan2(mean(sin(ship_cog(1:firstr)/180*pi)),mean(cos(ship_cog(1:firstr)/180*pi)))*180/pi;
std_cog_1=std(ship_cog(1:firstr));
min_cog_1=min(ship_cog(1:firstr));
max_cog_1=max(ship_cog(1:firstr));
if mean_cog_1 < 0, mean_cog_1=mean_cog_1+360;end;
if min_cog_1 < 0, min_cog_1=min_cog_1+360;end;
if max_cog_1 < 0, max_cog_1=max_cog_1+360;end;

Result(1,1:firstq)=Result(1,1:firstq)+mean_sog_1.*cos(mean_cog_1/180*pi)+mean(Platvel(1,1:firstq));
Result(2,1:firstq)=Result(2,1:firstq)+mean_sog_1.*sin(mean_cog_1/180*pi)+mean(Platvel(2,1:firstq));
Result(3,1:firstq)=Result(3,1:firstq)+mean(Platvel(3,1:firstq));

uvw(1,1:firstq)=uvw(1,1:firstq)+mean_sog_1.*cos(mean_cog_1/180*pi)+mean(uvwplat(1,1:firstq));
uvw(2,1:firstq)=uvw(2,1:firstq)+mean_sog_1.*sin(mean_cog_1/180*pi)+mean(uvwplat(2,1:firstq));
uvw(3,1:firstq)=uvw(3,1:firstq)+mean(uvwplat(3,1:firstq));

[Ustream_1,Vstream_1,Wstream_1,alpha_1,beta_1]=rotcorr(Result(1,1:firstq),Result(2,1:firstq),Result(3,1:firstq));           %converts [u,v,w] to [Ustream, Vstream, Wstream] for ten minutes
[Uj_1,Vj_1,Wj_1,alphaj_1,betaj_1]=rotcorr(uvw(1,1:firstq),uvw(2,1:firstq),uvw(3,1:firstq));

Ubarj_1=mean(Uj_1);
Tbar_1=mean(ship_Tair(1:firstr));      %Best
Tbarv_1=mean(Tairv(1:firsts));
rho_1=P*100./(273.15+Tbar_1)./287;

wu_1=cov(Ustream_1,Wstream_1);
wv_1=cov(Vstream_1,Wstream_1);
wt_1=cov(Wstream_1,sst(1:firstq));
wtd_1=cov(Wstream_1,Tdry(1:firstq));
Ustr_1=sqrt(mean(su(1:firstu)).^2 + mean(sv(1:firstu)).^2 + mean(sw(1:firstu)).^2);

wuj_1=cov(Uj_1,Wj_1);wvj_1=cov(Vj_1,Wj_1);
wtj_1=cov(Wj_1,sst(1:firstq));
waj_1=cov(Wj_1,abshum(1:firstq)./rho_1);   % m/s * g/kg

mean_head_1=atan2(mean(sin(heading(1:firstq))),mean(cos(heading(1:firstq))))*180/pi;
std_head_1=std(heading(1:firstq))*180/pi;
min_head_1=min(heading(1:firstq))*180/pi;
max_head_1=max(heading(1:firstq))*180/pi;
if mean_head_1<0, mean_head_1=mean_head_1+360;end;
if min_head_1<0, min_head_1=min_head_1+360;end;
if max_head_1<0, max_head_1=max_head_1+360;end;

mean_wind_1=mean(ship_wind(1:firstr));
std_wind_1=std(ship_wind(1:firstr));
max_wind_1=max(ship_wind(1:firstr));
min_wind_1=min(ship_wind(1:firstr));

mean_dir_1=mean(ship_dir(1:firstr));
std_dir_1=std(ship_dir(1:firstr));
max_dir_1=max(ship_dir(1:firstr));
min_dir_1=min(ship_dir(1:firstr));

Reldir_1=atan2(mean(sv(1:firstu)),-mean(su(1:firstu)))*180/pi;
Un_1=mean(su(1:firstu)).*cos(mean_head_1/180*pi)-mean(-sv(1:firstu)).*sin(mean_head_1/180*pi)+mean_sog_1.*cos(mean_cog_1/180*pi);   %to north
Ue_1=mean(su(1:firstu)).*sin(mean_head_1/180*pi)+mean(-sv(1:firstu)).*cos(mean_head_1/180*pi)+mean_sog_1.*sin(mean_cog_1/180*pi);   %to east
Ubar_1=sqrt(Un_1.^2 + Ue_1.^2);

%disp(['Ubar_1 = ',num2str(Ubar_1)]);
%disp(['Ubarj_1 = ',num2str(Ubarj_1)]);

Twindir_1=atan2(-Ue_1,-Un_1).*180/pi;
if Reldir_1<0,Reldir_1=Reldir_1+360;end;
if Twindir_1<0,Twindir_1=Twindir_1+360;end;

wa_1=cov(Wstream_1,abshum(1:firstq)./rho_1);   % m/s * g/kg
wh_1=cov(Wstream_1,H2O(1:firstq)./rho_1);
wc_1=cov(Wstream_1,CO2(1:firstq));

RH_1=mean(ship_hum(1:firstr));
RHv_1=mean(RHv(1:firsts));        %Best
Qbar_1=0.622*RH_1/100*6.112*exp((17.67*Tbar_1)/(Tbar_1+243.5))/1013.25*1000;
Qbarv_1=0.622*RHv_1/100*6.112*exp((17.67*Tbar_1)/(Tbar_1+243.5))/1013.25*1000;

Tseasnake_1=mean(Tsea_snake(1:firsts));
Esatsnake_1=6.112.*exp(17.67.*Tseasnake_1./(Tseasnake_1+243.5))*0.98;
Qseasnake_1=Esatsnake_1*0.622/(P-0.378*Esatsnake_1)*1000;

Tsea_1=mean(ship_Tsea(1:firstr));
Esat_1=6.112.*exp(17.67.*Tsea_1./(Tsea_1+243.5))*0.98;
Qsea_1=Esat_1*0.622/(P-0.378*Esat_1)*1000;

Ssin_1=std(shead(1:firstq));
Scos_1=std(chead(1:firstq));
Sp1_1=std(Platvel(1,1:firstq));
Sp2_1=std(Platvel(2,1:firstq));
Sp3_1=std(Platvel(3,1:firstq));

rain_1=mean(rain(1:firsts));
Rl_1=(median(LongIR1(1:firsts))+median(LongIR2(1:firsts)))/2;
Rs_1=(median(Solar1(1:firsts))+median(Solar2(1:firsts)))/2;
std_Rl_1=std((LongIR1(1:firsts) + LongIR2(1:firsts))/2);
std_Rs_1=std((Solar1(1:firsts) + Solar2(1:firsts))/2);

mean_CO2_1=mean(CO2(1:firstq));
mean_H2Oa_1=mean(H2O(1:firstq)./rho_1);
mean_H2Op_1=mean(abshum(1:firstq)./rho_1);       %  g/kg
mean_H2Ow_1=mean(Quw(1:firstq));
med_CO2_1=median(CO2(1:firstq));
med_H2Oa_1=median(H2O(1:firstq)./rho_1);
med_H2Op_1=median(abshum(1:firstq)./rho_1);       % g/kg
med_H2Ow_1=median(Quw(1:firstq));
std_CO2_1=std(CO2(1:firstq));
std_H2Oa_1=std(H2O(1:firstq)./rho_1);
std_H2Op_1=std(abshum(1:firstq)./rho_1);       % g/kg
std_H2Ow_1=std(Quw(1:firstq));

std_absorb_1=std(Oabsorb(1:firstq));
std_clear_1=std(Oclear(1:firstq));
med_absorb_1=median(Oabsorb(1:firstq));
med_clear_1=median(Oclear(1:firstq));

mean_gpscog_1=atan2(mean(sin(gpshead(1:firstt)/180*pi)),mean(cos(gpshead(1:firstt)/180*pi)))*180/pi;
std_gpscog_1=std(gpshead(1:firstt));
min_gpscog_1=min(gpshead(1:firstt));
max_gpscog_1=max(gpshead(1:firstt));
if mean_gpscog_1<0, mean_gpscog_1=mean_gpscog_1+360;end;
if min_gpscog_1<0, min_gpscog_1=min_gpscog_1+360;end;
if max_gpscog_1<0, max_gpscog_1=max_gpscog_1+360;end;

mean_gpsspeed_1=mean(gpsspeed(1:firstt))/1.94;
min_gpsspeed_1=min(gpsspeed(1:firstt))/1.94;
max_gpsspeed_1=max(gpsspeed(1:firstt))/1.94;

Lat_1=gpslat(1);
Lon_1=gpslon(1);

mean_gyro_1=atan2(mean(sin(ship_gyro(1:firstr)/180*pi)),mean(cos(ship_gyro(1:firstr)/180*pi)))*180/pi;
std_gyro_1=std(ship_gyro(1:firstr));
min_gyro_1=min(ship_gyro(1:firstr));
max_gyro_1=max(ship_gyro(1:firstr));
if mean_gyro_1 < 0, mean_gyro_1=mean_gyro_1+360;end;
if min_gyro_1 < 0, min_gyro_1=min_gyro_1+360;end;
if max_gyro_1 < 0, max_gyro_1=max_gyro_1+360;end;

xship_1=[Ubar_1 Us Tsea_1 Tbar_1 Qsea_1 Qbar_1 Rs_1 Rl_1 rain_1 zi P zu zt zq];
xsnake_1=[Ubar_1 Us Tseasnake_1 Tbar_1 Qseasnake_1 Qbar_1 Rs_1 Rl_1 rain_1 zi P zu zt zq];
yship_1=cor26(xship_1);
ysnake_1=cor26(xsnake_1);

%*************************Minute 10-20*******************************************

mean_sog_2=mean(ship_sog(firstr+1:seconr));
if mean_sog_2 < 1.5, mean_sog_2=0;end;
std_sog_2=std(ship_sog(firstr+1:seconr));
min_sog_2=min(ship_sog(firstr+1:seconr));
max_sog_2=max(ship_sog(firstr+1:seconr));

mean_cog_2=atan2(mean(sin(ship_cog(firstr+1:seconr)/180*pi)),mean(cos(ship_cog(firstr+1:seconr)/180*pi)))*180/pi;
std_cog_2=std(ship_cog(firstr+1:seconr));
min_cog_2=min(ship_cog(firstr+1:seconr));
max_cog_2=max(ship_cog(firstr+1:seconr));
if mean_cog_2 < 0, mean_cog_2=mean_cog_2+360;end;
if min_cog_2 < 0, min_cog_2=min_cog_2+360;end;
if max_cog_2 < 0, max_cog_2=max_cog_2+360;end;

Result(1,firstq+1:seconq)=Result(1,firstq+1:seconq)+mean_sog_2.*cos(mean_cog_2/180*pi)+mean(Platvel(1,firstq+1:seconq));
Result(2,firstq+1:seconq)=Result(2,firstq+1:seconq)+mean_sog_2.*sin(mean_cog_2/180*pi)+mean(Platvel(2,firstq+1:seconq));
Result(3,firstq+1:seconq)=Result(3,firstq+1:seconq)+mean(Platvel(3,firstq+1:seconq));

uvw(1,firstq+1:seconq)=uvw(1,firstq+1:seconq)+mean_sog_2.*cos(mean_cog_2/180*pi)+mean(uvwplat(1,firstq+1:seconq));
uvw(2,firstq+1:seconq)=uvw(2,firstq+1:seconq)+mean_sog_2.*sin(mean_cog_2/180*pi)+mean(uvwplat(2,firstq+1:seconq));
uvw(3,firstq+1:seconq)=uvw(3,firstq+1:seconq)+mean(uvwplat(3,firstq+1:seconq));

[Ustream_2,Vstream_2,Wstream_2,alpha_2,beta_2]=rotcorr(Result(1,firstq+1:seconq),Result(2,firstq+1:seconq),Result(3,firstq+1:seconq));           %converts [u,v,w] to [Ustream, Vstream, Wstream] for ten minutes
[Uj_2,Vj_2,Wj_2,alphaj_2,betaj_2]=rotcorr(uvw(1,firstq+1:seconq),uvw(2,firstq+1:seconq),uvw(3,firstq+1:seconq));

Ubarj_2=mean(Uj_2);
Tbar_2=mean(ship_Tair(firstr+1:seconr));      %Best
Tbarv_2=mean(Tairv(firsts+1:secons));
rho_2=P*100./(273.15+Tbar_2)./287;

wu_2=cov(Ustream_2,Wstream_2);
wv_2=cov(Vstream_2,Wstream_2);
wt_2=cov(Wstream_2,sst(firstq+1:seconq));
wtd_2=cov(Wstream_2,Tdry(firstq+1:seconq));
Ustr_2=sqrt(mean(su(firstu+1:seconu)).^2 + mean(sv(firstu+1:seconu)).^2 + mean(sw(firstu+1:seconu)).^2);

wuj_2=cov(Uj_2,Wj_2);
wvj_2=cov(Vj_2,Wj_2);
wtj_2=cov(Wj_2,sst(firstq+1:seconq));
waj_2=cov(Wj_2,abshum(firstq+1:seconq)./rho_2);       % m/s * g/kg

mean_head_2=atan2(mean(sin(heading(firstq+1:seconq))),mean(cos(heading(firstq+1:seconq))))*180/pi;
std_head_2=std(heading(firstq+1:seconq))*180/pi;
min_head_2=min(heading(firstq+1:seconq))*180/pi;
max_head_2=max(heading(firstq+1:seconq))*180/pi;
if mean_head_2<0, mean_head_2=mean_head_2+360;end;
if min_head_2<0, min_head_2=min_head_2+360;end;
if max_head_2<0, max_head_2=max_head_2+360;end;

mean_wind_2=mean(ship_wind(firstr+1:seconr));
std_wind_2=std(ship_wind(firstr+1:seconr));
max_wind_2=max(ship_wind(firstr+1:seconr));
min_wind_2=min(ship_wind(firstr+1:seconr));

mean_dir_2=mean(ship_dir(firstr+1:seconr));
std_dir_2=std(ship_dir(firstr+1:seconr));
max_dir_2=max(ship_dir(firstr+1:seconr));
min_dir_2=min(ship_dir(firstr+1:seconr));

Reldir_2=atan2(mean(sv(firstu+1:seconu)),-mean(su(firstu+1:seconu)))*180/pi;
Un_2=mean(su(firstu+1:seconu)).*cos(mean_head_2/180*pi)-mean(-sv(firstu+1:seconu)).*sin(mean_head_2/180*pi)+mean_sog_2.*cos(mean_cog_2/180*pi);   %to north
Ue_2=mean(su(firstu+1:seconu)).*sin(mean_head_2/180*pi)+mean(-sv(firstu+1:seconu)).*cos(mean_head_2/180*pi)+mean_sog_2.*sin(mean_cog_2/180*pi);   %to east
Ubar_2=sqrt(Un_2.^2 + Ue_2.^2);

%disp(['Ubar_2 = ',num2str(Ubar_2)]);
%disp(['Ubarj_2 = ',num2str(Ubarj_2)]);

Twindir_2=atan2(-Ue_2,-Un_2).*180/pi;
if Reldir_2<0,Reldir_2=Reldir_2+360;end;
if Twindir_2<0,Twindir_2=Twindir_2+360;end;

wa_2=cov(Wstream_2,abshum(firstq+1:seconq)./rho_2);       % m/s * g/kg       
wh_2=cov(Wstream_2,H2O(firstq+1:seconq)./rho_2);
wc_2=cov(Wstream_2,CO2(firstq+1:seconq));

RH_2=mean(ship_hum(firstr+1:seconr));
RHv_2=mean(RHv(firsts+1:secons));        %Best
Qbar_2=0.622*RH_2/100*6.112*exp((17.67*Tbar_2)/(Tbar_2+243.5))/1013.25*1000;
Qbarv_2=0.622*RHv_2/100*6.112*exp((17.67*Tbar_2)/(Tbar_2+243.5))/1013.25*1000;

Tseasnake_2=mean(Tsea_snake(firsts+1:secons));
Esatsnake_2=6.112.*exp(17.67.*Tseasnake_2./(Tseasnake_2+243.5))*0.98;
Qseasnake_2=Esatsnake_2*0.622/(P-0.378*Esatsnake_2)*1000;

Tsea_2=mean(ship_Tsea(firstr+1:seconr));
Esat_2=6.112.*exp(17.67.*Tsea_2./(Tsea_2+243.5))*0.98;
Qsea_2=Esat_2*0.622/(P-0.378*Esat_2)*1000;

Ssin_2=std(shead(firstq+1:seconq));
Scos_2=std(chead(firstq+1:seconq));
Sp1_2=std(Platvel(1,firstq+1:seconq));
Sp2_2=std(Platvel(2,firstq+1:seconq));
Sp3_2=std(Platvel(3,firstq+1:seconq));

rain_2=mean(rain(firsts+1:secons));
Rl_2=(median(LongIR1(firsts+1:secons))+median(LongIR2(firsts+1:secons)))/2;
Rs_2=(median(Solar1(firsts+1:secons))+median(Solar2(firsts+1:secons)))/2;
std_Rl_2=std((LongIR1(firsts+1:secons) + LongIR2(firsts+1:secons))/2);
std_Rs_2=std((Solar1(firsts+1:secons) + Solar2(firsts+1:secons))/2);

mean_CO2_2=mean(CO2(firstq+1:seconq));
mean_H2Oa_2=mean(H2O(firstq+1:seconq)./rho_2);
mean_H2Op_2=mean(abshum(firstq+1:seconq)./rho_2);       % g/kg
mean_H2Ow_2=mean(Quw(firstq+1:seconq));
med_CO2_2=median(CO2(firstq+1:seconq));
med_H2Oa_2=median(H2O(firstq+1:seconq)./rho_2);
med_H2Op_2=median(abshum(firstq+1:seconq)./rho_2);       % g/kg
med_H2Ow_2=median(Quw(firstq+1:seconq));
std_CO2_2=std(CO2(firstq+1:seconq));
std_H2Oa_2=std(H2O(firstq+1:seconq)./rho_2);
std_H2Op_2=std(abshum(firstq+1:seconq)./rho_2);       % g/kg
std_H2Ow_2=std(Quw(firstq+1:seconq));

std_absorb_2=std(Oabsorb(firstq+1:seconq));
std_clear_2=std(Oclear(firstq+1:seconq));
med_absorb_2=median(Oabsorb(firstq+1:seconq));
med_clear_2=median(Oclear(firstq+1:seconq));

mean_gpscog_2=atan2(mean(sin(gpshead(firstt+1:secont)/180*pi)),mean(cos(gpshead(firstt+1:secont)/180*pi)))*180/pi;
std_gpscog_2=std(gpshead(firstt+1:secont));
min_gpscog_2=min(gpshead(firstt+1:secont));
max_gpscog_2=max(gpshead(firstt+1:secont));
if mean_gpscog_2<0, mean_gpscog_2=mean_gpscog_2+360;end;
if min_gpscog_2<0, min_gpscog_2=min_gpscog_2+360;end;
if max_gpscog_2<0, max_gpscog_2=max_gpscog_2+360;end;

mean_gpsspeed_2=mean(gpsspeed(firstt+1:secont))/1.94;
min_gpsspeed_2=min(gpsspeed(firstt+1:secont))/1.94;
max_gpsspeed_2=max(gpsspeed(firstt+1:secont))/1.94;

Lat_2=gpslat(firstt+1);
Lon_2=gpslon(firstt+1);

mean_gyro_2=atan2(mean(sin(ship_gyro(firstr+1:seconr)/180*pi)),mean(cos(ship_gyro(firstr+1:seconr)/180*pi)))*180/pi;
std_gyro_2=std(ship_gyro(firstr+1:seconr));
min_gyro_2=min(ship_gyro(firstr+1:seconr));
max_gyro_2=max(ship_gyro(firstr+1:seconr));
if mean_gyro_2 < 0, mean_gyro_2=mean_gyro_2+360;end;
if min_gyro_2 < 0, min_gyro_2=min_gyro_2+360;end;
if max_gyro_2 < 0, max_gyro_2=max_gyro_2+360;end;

xship_2=[Ubar_2 Us Tsea_2 Tbar_2 Qsea_2 Qbar_2 Rs_2 Rl_2 rain_2 zi P zu zt zq];
xsnake_2=[Ubar_2 Us Tseasnake_2 Tbar_2 Qseasnake_2 Qbar_2 Rs_2 Rl_2 rain_2 zi P zu zt zq];
yship_2=cor26(xship_2);
ysnake_2=cor26(xsnake_2);

%*************************Minute 20-30********************************************

mean_sog_3=mean(ship_sog(seconr+1:thirdr));
if mean_sog_3 < 1.5, mean_sog_3=0;end;
std_sog_3=std(ship_sog(seconr+1:thirdr));
min_sog_3=min(ship_sog(seconr+1:thirdr));
max_sog_3=max(ship_sog(seconr+1:thirdr));

mean_cog_3=atan2(mean(sin(ship_cog(seconr+1:thirdr)/180*pi)),mean(cos(ship_cog(seconr+1:thirdr)/180*pi)))*180/pi;
std_cog_3=std(ship_cog(seconr+1:thirdr));
min_cog_3=min(ship_cog(seconr+1:thirdr));
max_cog_3=max(ship_cog(seconr+1:thirdr));
if mean_cog_3 < 0, mean_cog_3=mean_cog_3+360;end;
if min_cog_3 < 0, min_cog_3=min_cog_3+360;end;
if max_cog_3 < 0, max_cog_3=max_cog_3+360;end;

Result(1,seconq+1:thirdq)=Result(1,seconq+1:thirdq)+mean_sog_3.*cos(mean_cog_3/180*pi)+mean(Platvel(1,seconq+1:thirdq));
Result(2,seconq+1:thirdq)=Result(2,seconq+1:thirdq)+mean_sog_3.*sin(mean_cog_3/180*pi)+mean(Platvel(2,seconq+1:thirdq));
Result(3,seconq+1:thirdq)=Result(3,seconq+1:thirdq)+mean(Platvel(3,seconq+1:thirdq));

uvw(1,seconq+1:thirdq)=uvw(1,seconq+1:thirdq)+mean_sog_3.*cos(mean_cog_3/180*pi)+mean(uvwplat(1,seconq+1:thirdq));
uvw(2,seconq+1:thirdq)=uvw(2,seconq+1:thirdq)+mean_sog_3.*sin(mean_cog_3/180*pi)+mean(uvwplat(2,seconq+1:thirdq));
uvw(3,seconq+1:thirdq)=uvw(3,seconq+1:thirdq)+mean(uvwplat(3,seconq+1:thirdq));

[Ustream_3,Vstream_3,Wstream_3,alpha_3,beta_3]=rotcorr(Result(1,seconq+1:thirdq),Result(2,seconq+1:thirdq),Result(3,seconq+1:thirdq));           %converts [u,v,w] to [Ustream, Vstream, Wstream] for ten minutes
[Uj_3,Vj_3,Wj_3,alphaj_3,betaj_3]=rotcorr(uvw(1,seconq+1:thirdq),uvw(2,seconq+1:thirdq),uvw(3,seconq+1:thirdq));

Ubarj_3=mean(Uj_3);
Tbar_3=mean(ship_Tair(seconr+1:thirdr));      %Best
Tbarv_3=mean(Tairv(secons+1:thirds));
rho_3=P*100./(273.15+Tbar_3)./287;

wu_3=cov(Ustream_3,Wstream_3);
wv_3=cov(Vstream_3,Wstream_3);
wt_3=cov(Wstream_3,sst(seconq+1:thirdq));
wtd_3=cov(Wstream_3,Tdry(seconq+1:thirdq));
Ustr_3=sqrt(mean(su(seconu+1:thirdu)).^2 + mean(sv(seconu+1:thirdu)).^2 + mean(sw(seconu+1:thirdu)).^2);

wuj_3=cov(Uj_3,Wj_3);
wvj_3=cov(Vj_3,Wj_3);
wtj_3=cov(Wj_3,sst(seconq+1:thirdq));
waj_3=cov(Wj_3,abshum(seconq+1:thirdq)./rho_3);       % m/s * g/kg

mean_head_3=atan2(mean(sin(heading(seconq+1:thirdq))),mean(cos(heading(seconq+1:thirdq))))*180/pi;
std_head_3=std(heading(seconq+1:thirdq))*180/pi;
min_head_3=min(heading(seconq+1:thirdq))*180/pi;
max_head_3=max(heading(seconq+1:thirdq))*180/pi;
if mean_head_3<0, mean_head_3=mean_head_3+360;end;
if min_head_3<0, min_head_3=min_head_3+360;end;
if max_head_3<0, max_head_3=max_head_3+360;end;

mean_wind_3=mean(ship_wind(seconr+1:thirdr));
std_wind_3=std(ship_wind(seconr+1:thirdr));
max_wind_3=max(ship_wind(seconr+1:thirdr));
min_wind_3=min(ship_wind(seconr+1:thirdr));

mean_dir_3=mean(ship_dir(seconr+1:thirdr));
std_dir_3=std(ship_dir(seconr+1:thirdr));
max_dir_3=max(ship_dir(seconr+1:thirdr));
min_dir_3=min(ship_dir(seconr+1:thirdr));

Reldir_3=atan2(mean(sv(seconu+1:thirdu)),-mean(su(seconu+1:thirdu)))*180/pi;
Un_3=mean(su(seconu+1:thirdu)).*cos(mean_head_3/180*pi)-mean(-sv(seconu+1:thirdu)).*sin(mean_head_3/180*pi)+mean_sog_3.*cos(mean_cog_3/180*pi);   %to north
Ue_3=mean(su(seconu+1:thirdu)).*sin(mean_head_3/180*pi)+mean(-sv(seconu+1:thirdu)).*cos(mean_head_3/180*pi)+mean_sog_3.*sin(mean_cog_3/180*pi);   %to east
Ubar_3=sqrt(Un_3.^2 + Ue_3.^2);

%disp(['Ubar_3 = ',num2str(Ubar_3)]);
%disp(['Ubarj_3 = ',num2str(Ubarj_3)]);

Twindir_3=atan2(-Ue_3,-Un_3).*180/pi;
if Reldir_3<0,Reldir_3=Reldir_3+360;end;
if Twindir_3<0,Twindir_3=Twindir_3+360;end;

wa_3=cov(Wstream_3,abshum(seconq+1:thirdq)./rho_3);       % m/s * g/kg
wh_3=cov(Wstream_3,H2O(seconq+1:thirdq)./rho_3);
wc_3=cov(Wstream_3,CO2(seconq+1:thirdq));

RH_3=mean(ship_hum(seconr+1:thirdr));
RHv_3=mean(RHv(secons+1:thirds));        %Best
Qbar_3=0.622*RH_3/100*6.112*exp((17.67*Tbar_3)/(Tbar_3+243.5))/1013.25*1000;
Qbarv_3=0.622*RHv_3/100*6.112*exp((17.67*Tbar_3)/(Tbar_3+243.5))/1013.25*1000;

Tseasnake_3=mean(Tsea_snake(secons+1:thirds));
Esatsnake_3=6.112.*exp(17.67.*Tseasnake_3./(Tseasnake_3+243.5))*0.98;
Qseasnake_3=Esatsnake_3*0.622/(P-0.378*Esatsnake_3)*1000;

Tsea_3=mean(ship_Tsea(seconr+1:thirdr));
Esat_3=6.112.*exp(17.67.*Tsea_3./(Tsea_3+243.5))*0.98;
Qsea_3=Esat_3*0.622/(P-0.378*Esat_3)*1000;

Ssin_3=std(shead(seconq+1:thirdq));
Scos_3=std(chead(seconq+1:thirdq));
Sp1_3=std(Platvel(1,seconq+1:thirdq));
Sp2_3=std(Platvel(2,seconq+1:thirdq));
Sp3_3=std(Platvel(3,seconq+1:thirdq));

rain_3=mean(rain(secons+1:thirds));
Rl_3=(median(LongIR1(secons+1:thirds))+median(LongIR2(secons+1:thirds)))/2;
Rs_3=(median(Solar1(secons+1:thirds))+median(Solar2(secons+1:thirds)))/2;
std_Rl_3=std((LongIR1(secons+1:thirds) + LongIR2(secons+1:thirds))/2);
std_Rs_3=std((Solar1(secons+1:thirds) + Solar2(secons+1:thirds))/2);

mean_CO2_3=mean(CO2(seconq+1:thirdq));
mean_H2Oa_3=mean(H2O(seconq+1:thirdq)./rho_3);
mean_H2Op_3=mean(abshum(seconq+1:thirdq)./rho_3);       % g/kg
mean_H2Ow_3=mean(Quw(seconq+1:thirdq));
med_CO2_3=median(CO2(seconq+1:thirdq));
med_H2Oa_3=median(H2O(seconq+1:thirdq)./rho_3);
med_H2Op_3=median(abshum(seconq+1:thirdq)./rho_3);       % g/kg
med_H2Ow_3=median(Quw(seconq+1:thirdq));
std_CO2_3=std(CO2(seconq+1:thirdq));
std_H2Oa_3=std(H2O(seconq+1:thirdq)./rho_3);
std_H2Op_3=std(abshum(seconq+1:thirdq)./rho_3);       % g/kg
std_H2Ow_3=std(Quw(seconq+1:thirdq));

std_absorb_3=std(Oabsorb(seconq+1:thirdq));
std_clear_3=std(Oclear(seconq+1:thirdq));
med_absorb_3=median(Oabsorb(seconq+1:thirdq));
med_clear_3=median(Oclear(seconq+1:thirdq));

mean_gpscog_3=atan2(mean(sin(gpshead(secont+1:thirdt)/180*pi)),mean(cos(gpshead(secont+1:thirdt)/180*pi)))*180/pi;
std_gpscog_3=std(gpshead(secont+1:thirdt));
min_gpscog_3=min(gpshead(secont+1:thirdt));
max_gpscog_3=max(gpshead(secont+1:thirdt));
if mean_gpscog_3<0, mean_gpscog_3=mean_gpscog_3+360;end;
if min_gpscog_3<0, min_gpscog_3=min_gpscog_3+360;end;
if max_gpscog_3<0, max_gpscog_3=max_gpscog_3+360;end;

mean_gpsspeed_3=mean(gpsspeed(secont+1:thirdt))/1.94;
min_gpsspeed_3=min(gpsspeed(secont+1:thirdt))/1.94;
max_gpsspeed_3=max(gpsspeed(secont+1:thirdt))/1.94;

Lat_3=gpslat(secont+1);
Lon_3=gpslon(secont+1);

mean_gyro_3=atan2(mean(sin(ship_gyro(seconr+1:thirdr)/180*pi)),mean(cos(ship_gyro(seconr+1:thirdr)/180*pi)))*180/pi;
std_gyro_3=std(ship_gyro(seconr+1:thirdr));
min_gyro_3=min(ship_gyro(seconr+1:thirdr));
max_gyro_3=max(ship_gyro(seconr+1:thirdr));
if mean_gyro_3 < 0, mean_gyro_3=mean_gyro_3+360;end;
if min_gyro_3 < 0, min_gyro_3=min_gyro_3+360;end;
if max_gyro_3 < 0, max_gyro_3=max_gyro_3+360;end;

xship_3=[Ubar_3 Us Tsea_3 Tbar_3 Qsea_3 Qbar_3 Rs_3 Rl_3 rain_3 zi P zu zt zq];
xsnake_3=[Ubar_3 Us Tseasnake_3 Tbar_3 Qseasnake_3 Qbar_3 Rs_3 Rl_3 rain_3 zi P zu zt zq];
yship_3=cor26(xship_3);
ysnake_3=cor26(xsnake_3);

%*************************Minute 30-40********************************************

mean_sog_4=mean(ship_sog(thirdr+1:fourtr));
if mean_sog_4 < 1.5, mean_sog_4=0;end;
std_sog_4=std(ship_sog(thirdr+1:fourtr));
min_sog_4=min(ship_sog(thirdr+1:fourtr));
max_sog_4=max(ship_sog(thirdr+1:fourtr));

mean_cog_4=atan2(mean(sin(ship_cog(thirdr+1:fourtr)/180*pi)),mean(cos(ship_cog(thirdr+1:fourtr)/180*pi)))*180/pi;
std_cog_4=std(ship_cog(thirdr+1:fourtr));
min_cog_4=min(ship_cog(thirdr+1:fourtr));
max_cog_4=max(ship_cog(thirdr+1:fourtr));
if mean_cog_4 < 0, mean_cog_4=mean_cog_4+360;end;
if min_cog_4 < 0, min_cog_4=min_cog_4+360;end;
if max_cog_4 < 0, max_cog_4=max_cog_4+360;end;

Result(1,thirdq+1:fourtq)=Result(1,thirdq+1:fourtq)+mean_sog_4.*cos(mean_cog_4/180*pi)+mean(Platvel(1,thirdq+1:fourtq));
Result(2,thirdq+1:fourtq)=Result(2,thirdq+1:fourtq)+mean_sog_4.*sin(mean_cog_4/180*pi)+mean(Platvel(2,thirdq+1:fourtq));
Result(3,thirdq+1:fourtq)=Result(3,thirdq+1:fourtq)+mean(Platvel(3,thirdq+1:fourtq));

uvw(1,thirdq+1:fourtq)=uvw(1,thirdq+1:fourtq)+mean_sog_4.*cos(mean_cog_4/180*pi)+mean(uvwplat(1,thirdq+1:fourtq));
uvw(2,thirdq+1:fourtq)=uvw(2,thirdq+1:fourtq)+mean_sog_4.*sin(mean_cog_4/180*pi)+mean(uvwplat(2,thirdq+1:fourtq));
uvw(3,thirdq+1:fourtq)=uvw(3,thirdq+1:fourtq)+mean(uvwplat(3,thirdq+1:fourtq));

[Ustream_4,Vstream_4,Wstream_4,alpha_4,beta_4]=rotcorr(Result(1,thirdq+1:fourtq),Result(2,thirdq+1:fourtq),Result(3,thirdq+1:fourtq));           %converts [u,v,w] to [Ustream, Vstream, Wstream] for ten minutes
[Uj_4,Vj_4,Wj_4,alphaj_4,betaj_4]=rotcorr(uvw(1,thirdq+1:fourtq),uvw(2,thirdq+1:fourtq),uvw(3,thirdq+1:fourtq));

Ubarj_4=mean(Uj_4);
Tbar_4=mean(ship_Tair(thirdr+1:fourtr));      %Best
Tbarv_4=mean(Tairv(thirds+1:fourts));
rho_4=P*100./(273.15+Tbar_4)./287;

wu_4=cov(Ustream_4,Wstream_4);
wv_4=cov(Vstream_4,Wstream_4);
wt_4=cov(Wstream_4,sst(thirdq+1:fourtq));
wtd_4=cov(Wstream_4,Tdry(thirdq+1:fourtq));
Ustr_4=sqrt(mean(su(thirdu+1:fourtu)).^2 + mean(sv(thirdu+1:fourtu)).^2 + mean(sw(thirdu+1:fourtu)).^2);

wuj_4=cov(Uj_4,Wj_4);
wvj_4=cov(Vj_4,Wj_4);
wtj_4=cov(Wj_4,sst(thirdq+1:fourtq));
waj_4=cov(Wj_4,abshum(thirdq+1:fourtq)./rho_4);       % m/s * g/kg

mean_head_4=atan2(mean(sin(heading(thirdq+1:fourtq))),mean(cos(heading(thirdq+1:fourtq))))*180/pi;
std_head_4=std(heading(thirdq+1:fourtq))*180/pi;
min_head_4=min(heading(thirdq+1:fourtq))*180/pi;
max_head_4=max(heading(thirdq+1:fourtq))*180/pi;
if mean_head_4<0, mean_head_4=mean_head_4+360;end;
if min_head_4<0, min_head_4=min_head_4+360;end;
if max_head_4<0, max_head_4=max_head_4+360;end;

mean_wind_4=mean(ship_wind(thirdr+1:fourtr));
std_wind_4=std(ship_wind(thirdr+1:fourtr));
max_wind_4=max(ship_wind(thirdr+1:fourtr));
min_wind_4=min(ship_wind(thirdr+1:fourtr));

mean_dir_4=mean(ship_dir(thirdr+1:fourtr));
std_dir_4=std(ship_dir(thirdr+1:fourtr));
max_dir_4=max(ship_dir(thirdr+1:fourtr));
min_dir_4=min(ship_dir(thirdr+1:fourtr));

Reldir_4=atan2(mean(sv(thirdu+1:fourtu)),-mean(su(thirdu+1:fourtu)))*180/pi;
Un_4=mean(su(thirdu+1:fourtu)).*cos(mean_head_4/180*pi)-mean(-sv(thirdu+1:fourtu)).*sin(mean_head_4/180*pi)+mean_sog_4.*cos(mean_cog_4/180*pi);   %to north
Ue_4=mean(su(thirdu+1:fourtu)).*sin(mean_head_4/180*pi)+mean(-sv(thirdu+1:fourtu)).*cos(mean_head_4/180*pi)+mean_sog_4.*sin(mean_cog_4/180*pi);   %to east
Ubar_4=sqrt(Un_4.^2 + Ue_4.^2);

%disp(['Ubar_4 = ',num2str(Ubar_4)]);
%disp(['Ubarj_4 = ',num2str(Ubarj_4)]);

Twindir_4=atan2(-Ue_4,-Un_4).*180/pi;
if Reldir_4<0,Reldir_4=Reldir_4+360;end;
if Twindir_4<0,Twindir_4=Twindir_4+360;end;

wa_4=cov(Wstream_4,abshum(thirdq+1:fourtq)./rho_4);       % m/s * g/kg
wh_4=cov(Wstream_4,H2O(thirdq+1:fourtq)./rho_4);
wc_4=cov(Wstream_4,CO2(thirdq+1:fourtq));

RH_4=mean(ship_hum(thirdr+1:fourtr));
RHv_4=mean(RHv(thirds+1:fourts));        %Best
Qbar_4=0.622*RH_4/100*6.112*exp((17.67*Tbar_4)/(Tbar_4+243.5))/1013.25*1000;
Qbarv_4=0.622*RHv_4/100*6.112*exp((17.67*Tbar_4)/(Tbar_4+243.5))/1013.25*1000;

Tseasnake_4=mean(Tsea_snake(thirds+1:fourts));
Esatsnake_4=6.112.*exp(17.67.*Tseasnake_4./(Tseasnake_4+243.5))*0.98;
Qseasnake_4=Esatsnake_4*0.622/(P-0.378*Esatsnake_4)*1000;

Tsea_4=mean(ship_Tsea(thirdr+1:fourtr));
Esat_4=6.112.*exp(17.67.*Tsea_4./(Tsea_4+243.5))*0.98;
Qsea_4=Esat_4*0.622/(P-0.378*Esat_4)*1000;

Ssin_4=std(shead(thirdq+1:fourtq));
Scos_4=std(chead(thirdq+1:fourtq));
Sp1_4=std(Platvel(1,thirdq+1:fourtq));
Sp2_4=std(Platvel(2,thirdq+1:fourtq));
Sp3_4=std(Platvel(3,thirdq+1:fourtq));

rain_4=mean(rain(thirds+1:fourts));
Rl_4=(median(LongIR1(thirds+1:fourts))+median(LongIR2(thirds+1:fourts)))/2;
Rs_4=(median(Solar1(thirds+1:fourts))+median(Solar2(thirds+1:fourts)))/2;
std_Rl_4=std((LongIR1(thirds+1:fourts) + LongIR2(thirds+1:fourts))/2);
std_Rs_4=std((Solar1(thirds+1:fourts) + Solar2(thirds+1:fourts))/2);

mean_CO2_4=mean(CO2(thirdq+1:fourtq));
mean_H2Oa_4=mean(H2O(thirdq+1:fourtq)./rho_4);
mean_H2Op_4=mean(abshum(thirdq+1:fourtq)./rho_4);       % g/kg
mean_H2Ow_4=mean(Quw(thirdq+1:fourtq));
med_CO2_4=median(CO2(thirdq+1:fourtq));
med_H2Oa_4=median(H2O(thirdq+1:fourtq)./rho_4);
med_H2Op_4=median(abshum(thirdq+1:fourtq)./rho_4);       % g/kg
med_H2Ow_4=median(Quw(thirdq+1:fourtq));
std_CO2_4=std(CO2(thirdq+1:fourtq));
std_H2Oa_4=std(H2O(thirdq+1:fourtq)./rho_4);
std_H2Op_4=std(abshum(thirdq+1:fourtq)./rho_4);       % g/kg
std_H2Ow_4=std(Quw(thirdq+1:fourtq));

std_absorb_4=std(Oabsorb(thirdq+1:fourtq));
std_clear_4=std(Oclear(thirdq+1:fourtq));
med_absorb_4=median(Oabsorb(thirdq+1:fourtq));
med_clear_4=median(Oclear(thirdq+1:fourtq));

mean_gpscog_4=atan2(mean(sin(gpshead(thirdt+1:fourtt)/180*pi)),mean(cos(gpshead(thirdt+1:fourtt)/180*pi)))*180/pi;
std_gpscog_4=std(gpshead(thirdt+1:fourtt));
min_gpscog_4=min(gpshead(thirdt+1:fourtt));
max_gpscog_4=max(gpshead(thirdt+1:fourtt));
if mean_gpscog_4<0, mean_gpscog_4=mean_gpscog_4+360;end;
if min_gpscog_4<0, min_gpscog_4=min_gpscog_4+360;end;
if max_gpscog_4<0, max_gpscog_4=max_gpscog_4+360;end;

mean_gpsspeed_4=mean(gpsspeed(thirdt+1:fourtt))/1.94;
min_gpsspeed_4=min(gpsspeed(thirdt+1:fourtt))/1.94;
max_gpsspeed_4=max(gpsspeed(thirdt+1:fourtt))/1.94;

Lat_4=gpslat(thirdt+1);
Lon_4=gpslon(thirdt+1);

mean_gyro_4=atan2(mean(sin(ship_gyro(thirdr+1:fourtr)/180*pi)),mean(cos(ship_gyro(thirdr+1:fourtr)/180*pi)))*180/pi;
std_gyro_4=std(ship_gyro(thirdr+1:fourtr));
min_gyro_4=min(ship_gyro(thirdr+1:fourtr));
max_gyro_4=max(ship_gyro(thirdr+1:fourtr));
if mean_gyro_4 < 0, mean_gyro_4=mean_gyro_4+360;end;
if min_gyro_4 < 0, min_gyro_4=min_gyro_4+360;end;
if max_gyro_4 < 0, max_gyro_4=max_gyro_4+360;end;

xship_4=[Ubar_4 Us Tsea_4 Tbar_4 Qsea_4 Qbar_4 Rs_4 Rl_4 rain_4 zi P zu zt zq];
xsnake_4=[Ubar_4 Us Tseasnake_4 Tbar_4 Qseasnake_4 Qbar_4 Rs_4 Rl_4 rain_4 zi P zu zt zq];
yship_4=cor26(xship_4);
ysnake_4=cor26(xsnake_4);

%*************************Minute 40-50********************************************

mean_sog_5=mean(ship_sog(fourtr+1:fifthr));
if mean_sog_5 < 1.5, mean_sog_5=0;end;
std_sog_5=std(ship_sog(fourtr+1:fifthr));
min_sog_5=min(ship_sog(fourtr+1:fifthr));
max_sog_5=max(ship_sog(fourtr+1:fifthr));

mean_cog_5=atan2(mean(sin(ship_cog(fourtr+1:fifthr)/180*pi)),mean(cos(ship_cog(fourtr+1:fifthr)/180*pi)))*180/pi;
std_cog_5=std(ship_cog(fourtr+1:fifthr));
min_cog_5=min(ship_cog(fourtr+1:fifthr));
max_cog_5=max(ship_cog(fourtr+1:fifthr));
if mean_cog_5 < 0, mean_cog_5=mean_cog_5+360;end;
if min_cog_5 < 0, min_cog_5=min_cog_5+360;end;
if max_cog_5 < 0, max_cog_5=max_cog_5+360;end;

Result(1,fourtq+1:fifthq)=Result(1,fourtq+1:fifthq)+mean_sog_5.*cos(mean_cog_5/180*pi)+mean(Platvel(1,fourtq+1:fifthq));
Result(2,fourtq+1:fifthq)=Result(2,fourtq+1:fifthq)+mean_sog_5.*sin(mean_cog_5/180*pi)+mean(Platvel(2,fourtq+1:fifthq));
Result(3,fourtq+1:fifthq)=Result(3,fourtq+1:fifthq)+mean(Platvel(3,fourtq+1:fifthq));

uvw(1,fourtq+1:fifthq)=uvw(1,fourtq+1:fifthq)+mean_sog_5.*cos(mean_cog_5/180*pi)+mean(uvwplat(1,fourtq+1:fifthq));
uvw(2,fourtq+1:fifthq)=uvw(2,fourtq+1:fifthq)+mean_sog_5.*sin(mean_cog_5/180*pi)+mean(uvwplat(2,fourtq+1:fifthq));
uvw(3,fourtq+1:fifthq)=uvw(3,fourtq+1:fifthq)+mean(uvwplat(3,fourtq+1:fifthq));

[Ustream_5,Vstream_5,Wstream_5,alpha_5,beta_5]=rotcorr(Result(1,fourtq+1:fifthq),Result(2,fourtq+1:fifthq),Result(3,fourtq+1:fifthq));           %converts [u,v,w] to [Ustream, Vstream, Wstream] for ten minutes
[Uj_5,Vj_5,Wj_5,alphaj_5,betaj_5]=rotcorr(uvw(1,fourtq+1:fifthq),uvw(2,fourtq+1:fifthq),uvw(3,fourtq+1:fifthq));

Ubarj_5=mean(Uj_5);
Tbar_5=mean(ship_Tair(fourtr+1:fifthr));      %Best
Tbarv_5=mean(Tairv(fourts+1:fifths));
rho_5=P*100./(273.15+Tbar_5)./287;

wu_5=cov(Ustream_5,Wstream_5);
wv_5=cov(Vstream_5,Wstream_5);
wt_5=cov(Wstream_5,sst(fourtq+1:fifthq));
wtd_5=cov(Wstream_5,Tdry(fourtq+1:fifthq));
Ustr_5=sqrt(mean(su(fourtu+1:fifthu)).^2 + mean(sv(fourtu+1:fifthu)).^2 + mean(sw(fourtu+1:fifthu)).^2);

wuj_5=cov(Uj_5,Wj_5);
wvj_5=cov(Vj_5,Wj_5);
wtj_5=cov(Wj_5,sst(fourtq+1:fifthq));
waj_5=cov(Wj_5,abshum(fourtq+1:fifthq)./rho_5);       % m/s * g/kg

mean_head_5=atan2(mean(sin(heading(fourtq+1:fifthq))),mean(cos(heading(fourtq+1:fifthq))))*180/pi;
std_head_5=std(heading(fourtq+1:fifthq))*180/pi;
min_head_5=min(heading(fourtq+1:fifthq))*180/pi;
max_head_5=max(heading(fourtq+1:fifthq))*180/pi;
if mean_head_5<0, mean_head_5=mean_head_5+360;end;
if min_head_5<0, min_head_5=min_head_5+360;end;
if max_head_5<0, max_head_5=max_head_5+360;end;

mean_wind_5=mean(ship_wind(fourtr+1:fifthr));
std_wind_5=std(ship_wind(fourtr+1:fifthr));
max_wind_5=max(ship_wind(fourtr+1:fifthr));
min_wind_5=min(ship_wind(fourtr+1:fifthr));

mean_dir_5=mean(ship_dir(fourtr+1:fifthr));
std_dir_5=std(ship_dir(fourtr+1:fifthr));
max_dir_5=max(ship_dir(fourtr+1:fifthr));
min_dir_5=min(ship_dir(fourtr+1:fifthr));

Reldir_5=atan2(mean(sv(fourtu+1:fifthu)),-mean(su(fourtu+1:fifthu)))*180/pi;
Un_5=mean(su(fourtu+1:fifthu)).*cos(mean_head_5/180*pi)-mean(-sv(fourtu+1:fifthu)).*sin(mean_head_5/180*pi)+mean_sog_5.*cos(mean_cog_5/180*pi);   %to north
Ue_5=mean(su(fourtu+1:fifthu)).*sin(mean_head_5/180*pi)+mean(-sv(fourtu+1:fifthu)).*cos(mean_head_5/180*pi)+mean_sog_5.*sin(mean_cog_5/180*pi);   %to east
Ubar_5=sqrt(Un_5.^2 + Ue_5.^2);

%disp(['Ubar_5 = ',num2str(Ubar_5)]);
%disp(['Ubarj_5 = ',num2str(Ubarj_5)]);

Twindir_5=atan2(-Ue_5,-Un_5).*180/pi;
if Reldir_5<0,Reldir_5=Reldir_5+360;end;
if Twindir_5<0,Twindir_5=Twindir_5+360;end;

wa_5=cov(Wstream_5,abshum(fourtq+1:fifthq)./rho_5);       % m/s * g/kg
wh_5=cov(Wstream_5,H2O(fourtq+1:fifthq)./rho_5);
wc_5=cov(Wstream_5,CO2(fourtq+1:fifthq));

RH_5=mean(ship_hum(fourtr+1:fifthr));
RHv_5=mean(RHv(fourts+1:fifths));        %Best
Qbar_5=0.622*RH_5/100*6.112*exp((17.67*Tbar_5)/(Tbar_5+243.5))/1013.25*1000;
Qbarv_5=0.622*RHv_5/100*6.112*exp((17.67*Tbar_5)/(Tbar_5+243.5))/1013.25*1000;

Tseasnake_5=mean(Tsea_snake(fourts+1:fifths));
Esatsnake_5=6.112.*exp(17.67.*Tseasnake_5./(Tseasnake_5+243.5))*0.98;
Qseasnake_5=Esatsnake_5*0.622/(P-0.378*Esatsnake_5)*1000;

Tsea_5=mean(ship_Tsea(fourtr+1:fifthr));
Esat_5=6.112.*exp(17.67.*Tsea_5./(Tsea_5+243.5))*0.98;
Qsea_5=Esat_5*0.622/(P-0.378*Esat_5)*1000;

Ssin_5=std(shead(fourtq+1:fifthq));
Scos_5=std(chead(fourtq+1:fifthq));
Sp1_5=std(Platvel(1,fourtq+1:fifthq));
Sp2_5=std(Platvel(2,fourtq+1:fifthq));
Sp3_5=std(Platvel(3,fourtq+1:fifthq));

rain_5=mean(rain(fourts+1:fifths));
Rl_5=(median(LongIR1(fourts+1:fifths))+median(LongIR2(fourts+1:fifths)))/2;
Rs_5=(median(Solar1(fourts+1:fifths))+median(Solar2(fourts+1:fifths)))/2;
std_Rl_5=std((LongIR1(fourts+1:fifths) + LongIR2(fourts+1:fifths))/2);
std_Rs_5=std((Solar1(fourts+1:fifths) + Solar2(fourts+1:fifths))/2);

mean_CO2_5=mean(CO2(fourtq+1:fifthq));
mean_H2Oa_5=mean(H2O(fourtq+1:fifthq)./rho_5);
mean_H2Op_5=mean(abshum(fourtq+1:fifthq)./rho_5);       % g/kg
mean_H2Ow_5=mean(Quw(fourtq+1:fifthq));
med_CO2_5=median(CO2(fourtq+1:fifthq));
med_H2Oa_5=median(H2O(fourtq+1:fifthq)./rho_5);
med_H2Op_5=median(abshum(fourtq+1:fifthq)./rho_5);       % g/kg
med_H2Ow_5=median(Quw(fourtq+1:fifthq));
std_CO2_5=std(CO2(fourtq+1:fifthq));
std_H2Oa_5=std(H2O(fourtq+1:fifthq)./rho_5);
std_H2Op_5=std(abshum(fourtq+1:fifthq)./rho_5);       % g/kg
std_H2Ow_5=std(Quw(fourtq+1:fifthq));

std_absorb_5=std(Oabsorb(fourtq+1:fifthq));
std_clear_5=std(Oclear(fourtq+1:fifthq));
med_absorb_5=median(Oabsorb(fourtq+1:fifthq));
med_clear_5=median(Oclear(fourtq+1:fifthq));

mean_gpscog_5=atan2(mean(sin(gpshead(fourtt+1:fiftht)/180*pi)),mean(cos(gpshead(fourtt+1:fiftht)/180*pi)))*180/pi;
std_gpscog_5=std(gpshead(fourtt+1:fiftht));
min_gpscog_5=min(gpshead(fourtt+1:fiftht));
max_gpscog_5=max(gpshead(fourtt+1:fiftht));
if mean_gpscog_5<0, mean_gpscog_5=mean_gpscog_5+360;end;
if min_gpscog_5<0, min_gpscog_5=min_gpscog_5+360;end;
if max_gpscog_5<0, max_gpscog_5=max_gpscog_5+360;end;

mean_gpsspeed_5=mean(gpsspeed(fourtt+1:fiftht))/1.94;
min_gpsspeed_5=min(gpsspeed(fourtt+1:fiftht))/1.94;
max_gpsspeed_5=max(gpsspeed(fourtt+1:fiftht))/1.94;

Lat_5=gpslat(fourtt+1);
Lon_5=gpslon(fourtt+1);

mean_gyro_5=atan2(mean(sin(ship_gyro(fourtr+1:fifthr)/180*pi)),mean(cos(ship_gyro(fourtr+1:fifthr)/180*pi)))*180/pi;
std_gyro_5=std(ship_gyro(fourtr+1:fifthr));
min_gyro_5=min(ship_gyro(fourtr+1:fifthr));
max_gyro_5=max(ship_gyro(fourtr+1:fifthr));
if mean_gyro_5 < 0, mean_gyro_5=mean_gyro_5+360;end;
if min_gyro_5 < 0, min_gyro_5=min_gyro_5+360;end;
if max_gyro_5 < 0, max_gyro_5=max_gyro_5+360;end;

xship_5=[Ubar_5 Us Tsea_5 Tbar_5 Qsea_5 Qbar_5 Rs_5 Rl_5 rain_5 zi P zu zt zq];
xsnake_5=[Ubar_5 Us Tseasnake_5 Tbar_5 Qseasnake_5 Qbar_5 Rs_5 Rl_5 rain_5 zi P zu zt zq];
yship_5=cor26(xship_5);
ysnake_5=cor26(xsnake_5);

%*************************Minute 50-60********************************************

mean_sog_6=mean(ship_sog(fifthr+1:sixthr));
if mean_sog_1 < 1.5, mean_sog_6=0;end;
std_sog_6=std(ship_sog(fifthr+1:sixthr));
min_sog_6=min(ship_sog(fifthr+1:sixthr));
max_sog_6=max(ship_sog(fifthr+1:sixthr));

mean_cog_6=atan2(mean(sin(ship_cog(fifthr+1:sixthr)/180*pi)),mean(cos(ship_cog(fifthr+1:sixthr)/180*pi)))*180/pi;
std_cog_6=std(ship_cog(fifthr+1:sixthr));
min_cog_6=min(ship_cog(fifthr+1:sixthr));
max_cog_6=max(ship_cog(fifthr+1:sixthr));
if mean_cog_6 < 0, mean_cog_6=mean_cog_6+360;end;
if min_cog_6 < 0, min_cog_6=min_cog_6+360;end;
if max_cog_6 < 0, max_cog_6=max_cog_6+360;end;

Result(1,fifthq+1:sixthq)=Result(1,fifthq+1:sixthq)+mean_sog_6.*cos(mean_cog_6/180*pi)+mean(Platvel(1,fifthq+1:sixthq));
Result(2,fifthq+1:sixthq)=Result(2,fifthq+1:sixthq)+mean_sog_6.*sin(mean_cog_6/180*pi)+mean(Platvel(1,fifthq+1:sixthq));
Result(3,fifthq+1:sixthq)=Result(3,fifthq+1:sixthq)+mean(Platvel(3,fifthq+1:sixthq));

uvw(1,fifthq+1:sixthq)=uvw(1,fifthq+1:sixthq)+mean_sog_6.*cos(mean_cog_6/180*pi)+mean(uvwplat(1,fifthq+1:sixthq));
uvw(2,fifthq+1:sixthq)=uvw(2,fifthq+1:sixthq)+mean_sog_6.*sin(mean_cog_6/180*pi)+mean(uvwplat(1,fifthq+1:sixthq));
uvw(3,fifthq+1:sixthq)=uvw(3,fifthq+1:sixthq)+mean(uvwplat(3,fifthq+1:sixthq));

[Ustream_6,Vstream_6,Wstream_6,alpha_6,beta_6]=rotcorr(Result(1,fifthq+1:sixthq),Result(2,fifthq+1:sixthq),Result(3,fifthq+1:sixthq));           %converts [u,v,w] to [Ustream, Vstream, Wstream] for ten minutes
[Uj_6,Vj_6,Wj_6,alphaj_6,betaj_6]=rotcorr(uvw(1,fifthq+1:sixthq),uvw(2,fifthq+1:sixthq),uvw(3,fifthq+1:sixthq));

Ubarj_6=mean(Uj_6);
Tbar_6=mean(ship_Tair(fifthr+1:sixthr));      %Best
Tbarv_6=mean(Tairv(fifths+1:sixths));
rho_6=P*100./(273.15+Tbar_6)./287;

wu_6=cov(Ustream_6,Wstream_6);
wv_6=cov(Vstream_6,Wstream_6);
wt_6=cov(Wstream_6,sst(fifthq+1:sixthq));
wtd_6=cov(Wstream_6,Tdry(fifthq+1:sixthq));
Ustr_6=sqrt(mean(su(fifthu+1:sixthu)).^2 + mean(sv(fifthu+1:sixthu)).^2 + mean(sw(fifthu+1:sixthu)).^2);

wuj_6=cov(Uj_6,Wj_6);
wvj_6=cov(Vj_6,Wj_6);
wtj_6=cov(Wj_6,sst(fifthq+1:sixthq));
waj_6=cov(Wj_6,abshum(fifthq+1:sixthq)./rho_6);       % m/s * g/kg

mean_head_6=atan2(mean(sin(heading(fifthq+1:sixthq))),mean(cos(heading(fifthq+1:sixthq))))*180/pi;
std_head_6=std(heading(fifthq+1:sixthq))*180/pi;
min_head_6=min(heading(fifthq+1:sixthq))*180/pi;
max_head_6=max(heading(fifthq+1:sixthq))*180/pi;
if mean_head_6<0, mean_head_6=mean_head_6+360;end;
if min_head_6<0, min_head_6=min_head_6+360;end;
if max_head_6<0, max_head_6=max_head_6+360;end;

mean_wind_6=mean(ship_wind(fifthr+1:sixthr));
std_wind_6=std(ship_wind(fifthr+1:sixthr));
max_wind_6=max(ship_wind(fifthr+1:sixthr));
min_wind_6=min(ship_wind(fifthr+1:sixthr));

mean_dir_6=mean(ship_dir(fifthr+1:sixthr));
std_dir_6=std(ship_dir(fifthr+1:sixthr));
max_dir_6=max(ship_dir(fifthr+1:sixthr));
min_dir_6=min(ship_dir(fifthr+1:sixthr));

Reldir_6=atan2(mean(sv(fifthu+1:sixthu)),-mean(su(fifthq+1:sixthu)))*180/pi;
Un_6=mean(su(fifthu+1:sixthu)).*cos(mean_head_6/180*pi)-mean(-sv(fifthu+1:sixthu)).*sin(mean_head_6/180*pi)+mean_sog_6.*cos(mean_cog_6/180*pi);   %to north
Ue_6=mean(su(fifthu+1:sixthu)).*sin(mean_head_6/180*pi)+mean(-sv(fifthu+1:sixthu)).*cos(mean_head_6/180*pi)+mean_sog_6.*sin(mean_cog_6/180*pi);   %to east
Ubar_6=sqrt(Un_6.^2 + Ue_6.^2);

%disp(['Ubar_6 = ',num2str(Ubar_6)]);
%disp(['Ubarj_6 = ',num2str(Ubarj_6)]);

Twindir_6=atan2(-Ue_6,-Un_6).*180/pi;
if Reldir_6<0,Reldir_6=Reldir_6+360;end;
if Twindir_6<0,Twindir_6=Twindir_6+360;end;

wa_6=cov(Wstream_6,abshum(fifthq+1:sixthq)./rho_6);       % m/s * g/kg
wh_6=cov(Wstream_6,H2O(fifthq+1:sixthq)./rho_6);
wc_6=cov(Wstream_6,CO2(fifthq+1:sixthq));

RH_6=mean(ship_hum(fifthr+1:sixthr));
RHv_6=mean(RHv(fifths+1:sixths));        %Best
Qbar_6=0.622*RH_6/100*6.112*exp((17.67*Tbar_6)/(Tbar_6+243.5))/1013.25*1000;
Qbarv_6=0.622*RHv_6/100*6.112*exp((17.67*Tbar_6)/(Tbar_6+243.5))/1013.25*1000;

Tseasnake_6=mean(Tsea_snake(fifths+1:sixths));
Esatsnake_6=6.112.*exp(17.67.*Tseasnake_6./(Tseasnake_6+243.5))*0.98;
Qseasnake_6=Esatsnake_6*0.622/(P-0.378*Esatsnake_6)*1000;

Tsea_6=mean(ship_Tsea(fifthr+1:sixthr));
Esat_6=6.112.*exp(17.67.*Tsea_6./(Tsea_6+243.5))*0.98;
Qsea_6=Esat_6*0.622/(P-0.378*Esat_6)*1000;

Ssin_6=std(shead(fifthq+1:sixthq));
Scos_6=std(chead(fifthq+1:sixthq));
Sp1_6=std(Platvel(1,fifthq+1:sixthq));
Sp2_6=std(Platvel(2,fifthq+1:sixthq));
Sp3_6=std(Platvel(3,fifthq+1:sixthq));

rain_6=mean(rain(fifths+1:sixths));
Rl_6=(median(LongIR1(fifths+1:sixths))+median(LongIR2(fifths+1:sixths)))/2;
Rs_6=(median(Solar1(fifths+1:sixths))+median(Solar2(fifths+1:sixths)))/2;
std_Rl_6=std((LongIR1(fifths+1:sixths) + LongIR2(fifths+1:sixths))/2);
std_Rs_6=std((Solar1(fifths+1:sixths) + Solar2(fifths+1:sixths))/2);

mean_CO2_6=mean(CO2(fifthq+1:sixthq));
mean_H2Oa_6=mean(H2O(fifthq+1:sixthq)./rho_6);
mean_H2Op_6=mean(abshum(fifthq+1:sixthq)./rho_6);       % g/kg
mean_H2Ow_6=mean(Quw(fifthq+1:sixthq));
med_CO2_6=median(CO2(fifthq+1:sixthq));
med_H2Oa_6=median(H2O(fifthq+1:sixthq)./rho_6);
med_H2Op_6=median(abshum(fifthq+1:sixthq)./rho_6);       % g/kg
med_H2Ow_6=median(Quw(fifthq+1:sixthq));
std_CO2_6=std(CO2(fifthq+1:sixthq));
std_H2Oa_6=std(H2O(fifthq+1:sixthq)./rho_6);
std_H2Op_6=std(abshum(fifthq+1:sixthq)./rho_6);       % g/kg
std_H2Ow_6=std(Quw(fifthq+1:sixthq));

std_absorb_6=std(Oabsorb(fifthq+1:sixthq));
std_clear_6=std(Oclear(fifthq+1:sixthq));
med_absorb_6=median(Oabsorb(fifthq+1:sixthq));
med_clear_6=median(Oclear(fifthq+1:sixthq));

mean_gpscog_6=atan2(mean(sin(gpshead(fiftht+1:sixtht)/180*pi)),mean(cos(gpshead(fiftht+1:sixtht)/180*pi)))*180/pi;
std_gpscog_6=std(gpshead(fiftht+1:sixtht));
min_gpscog_6=min(gpshead(fiftht+1:sixtht));
max_gpscog_6=max(gpshead(fiftht+1:sixtht));
if mean_gpscog_6<0, mean_gpscog_6=mean_gpscog_6+360;end;
if min_gpscog_6<0, min_gpscog_6=min_gpscog_6+360;end;
if max_gpscog_6<0, max_gpscog_6=max_gpscog_6+360;end;

mean_gpsspeed_6=mean(gpsspeed(fiftht+1:sixtht))/1.94;
min_gpsspeed_6=min(gpsspeed(fiftht+1:sixtht))/1.94;
max_gpsspeed_6=max(gpsspeed(fiftht+1:sixtht))/1.94;

Lat_6=gpslat(fiftht+1);
Lon_6=gpslon(fiftht+1);

mean_gyro_6=atan2(mean(sin(ship_gyro(fifthr+1:sixthr)/180*pi)),mean(cos(ship_gyro(fifthr+1:sixthr)/180*pi)))*180/pi;
std_gyro_6=std(ship_gyro(fifthr+1:sixthr));
min_gyro_6=min(ship_gyro(fifthr+1:sixthr));
max_gyro_6=max(ship_gyro(fifthr+1:sixthr));
if mean_gyro_6 < 0, mean_gyro_6=mean_gyro_6+360;end;
if min_gyro_6 < 0, min_gyro_6=min_gyro_6+360;end;
if max_gyro_6 < 0, max_gyro_6=max_gyro_6+360;end;

xship_6=[Ubar_6 Us Tsea_6 Tbar_6 Qsea_6 Qbar_6 Rs_6 Rl_6 rain_6 zi P zu zt zq];
xsnake_6=[Ubar_6 Us Tseasnake_6 Tbar_6 Qseasnake_6 Qbar_6 Rs_6 Rl_6 rain_6 zi P zu zt zq];
yship_6=cor26(xship_6);
ysnake_6=cor26(xsnake_6);

clear su sv ship_dir ship_time ship_wind
clear uvw Ustream Vstream Wstream LongIR1 LongIR2 P Platvel RHv Result Solar1 Solar2
clear LongIR1 LongIR2 Platvel RHv Result Solar1 Solar2 Tdry
clear zi zq zt zu shead chead rain ppp px py pz heading fz 
clear firstr firsts firstt firstu
clear seconr secons secont seconu
clear thirdr thirds thirdt thirdu
clear fourtr fourts fourtt fourtu
clear fifthr fifths fiftht fifthu
clear sixthr sixths sixtht sixthu
clear Oabsorb Oclear Quw fx fy fxyz dfx1 dfx2 chead bf af Zf Qbar Qbarv
clear Velrot Tsea_snake Tairv Esat ship_Tair ship_Tsea ship_cog ship_gyro ship_hum ship_sog
clear ship_cog ship_gyro RH ship_sog gpslat_sign gpslon_sign gpslat gpslon gpshour
clear gpsmin gpssec gpshead gpstime gpsspeed
clear uvwplat sw

%   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);

clear tau2 spinup2 ans

%[wu_1(1,2) wuj_1(1,2)
%wu_2(1,2) wuj_2(1,2)
%wu_3(1,2) wuj_3(1,2)
%wu_4(1,2) wuj_4(1,2)
%wu_5(1,2) wuj_5(1,2)
%wu_6(1,2) wuj_6(1,2)]

