disp('nau_correct_hiwings2013.m')

sampleRate=10;  %every data coming here are at 10Hz.
Lendata=3600*sampleRate;
Nlendata=floor(Lendata/6);

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

ratex=detrend(ratex);
ratey=detrend(ratey);
ratez=detrend(ratez);
maccx=nanmean(accx);
maccy=nanmean(accy);
maccz=nanmean(accz);

Fc=1/30;                 %60 seconds
bc=Fc/(sampleRate/2);
[bf af]=butter(1,bc);    %designs an 1st order lowpass digital Butterworth filter with a cutoff frequency = bc
tau=1/(Fc*2*pi);
dt=1/sampleRate;
om_c=tan(dt/2/tau);tau_c=dt/2/om_c;

% wp=1/120/(20/2);
% ws=1/60/(20/2);   %20/2 is the Nyquist frequency.
wp=1/40/(sampleRate/2);
ws=1/20/(sampleRate/2);   %Edson coeff...
[n,wn]=buttord(wp,ws,1,sampleRate);  %returns the order N of the lowest order digital Butterworth filter that loses no more than Rp dB in 
%the passband and has at least Rs dB of attenuation in the stopband. Wp and Ws are the passband and stopband edge frequencies, 
%normalized from 0 to 1 (where 1 corresponds to pi radians/sample).
%Rp is the maximum permissible passband loss in decibels. Rs is the number of decibels the stopband is down from the passband.

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

clear Fc wn wp ws n

sdg=['load ', drive_w, '\fxyz.prn fxyz -ascii'];
disp(sdg)
eval(sdg)

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=nanmean(heading)/180*pi;
end;

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

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

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

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

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

psislowh = (-heading - filtfilt(bhigh,ahigh,-heading));
psi = psislowh + filtfilt(bhigh,ahigh,(cumsum(ratez) -.5*ratez - .5*ratez(1))/sampleRate);

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= 0.7;
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';


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

Vel=[U              %
     V              %
     W];            %

Plat=[accx
      accy
      accz];
  
Platacc=Plat;
Velrot=Vel;

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

deg_rate=[ratex
   ratey
   ratez];

if clearit
    clear ratex ratey ratez
end;

R=[Rx
   Ry
   Rz];

[bhi,ahi,bf2,af2,tau2,spinup2]=Filtcoef(sampleRate,1/30);

[euler,dr]=angles(ahi,bhi,sampleRate,Plat,deg_rate,heading);
%     its = 5; %Iterate 5 times
%     [euler,dr] = angles_UH(ahi,bhi,sampleRate,Plat,deg_rate,heading,its);

roll=euler(1,:);  %roll=resample(roll,1,sampleRate);   roll=roll';
pitch=euler(2,:); %pitch=resample(pitch,1,sampleRate); pitch=pitch';
yaw=euler(3,:); 

tilts=[roll pitch];


[uvwplat,xyzplat]=accels(bhi,ahi,sampleRate,Plat,euler);%clear xyzplat
%     [acc,uvwplat,xyzplat] = accels_UH(bhi,ahi,sampleRate,Plat,euler);
[uvw]=sonic_f(Vel,dr,euler,uvwplat,R);

%save the motion data for future motion sensor intercomparisons or to correct Riegl...
motmatall=[tref' uvwplat(1,:)' uvwplat(2,:)' uvwplat(3,:)' xyzplat(1,:)' xyzplat(2,:)' xyzplat(3,:)' euler(1,:)' euler(2,:)' euler(3,:)' deg_rate(1,:)' deg_rate(2,:)' deg_rate(3,:)'];
sdg=['save ', drive_w '\MotionData\HiWinGS2013_motionfile_',jd hr '.prn motmatall -ascii'];
disp(sdg)
eval(sdg)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


AAAA=uvw;
% % SpecCorr;


% % for ll=0:5
% % [uvw(1,6000*ll+1:6000*(ll+1)),uvwplat(1,6000*ll+1:6000*(ll+1)),Xlagflag] = lagcorrection(uvw(1,6000*ll+1:6000*(ll+1)),uvwplat(1,6000*ll+1:6000*(ll+1)),16);
% % [uvw(2,6000*ll+1:6000*(ll+1)),uvwplat(2,6000*ll+1:6000*(ll+1)),Xlagflag] = lagcorrection(uvw(2,6000*ll+1:6000*(ll+1)),uvwplat(2,6000*ll+1:6000*(ll+1)),16);
% % [uvw(3,6000*ll+1:6000*(ll+1)),uvwplat(3,6000*ll+1:6000*(ll+1)),Xlagflag] = lagcorrection(uvw(3,6000*ll+1:6000*(ll+1)),uvwplat(3,6000*ll+1:6000*(ll+1)),16);
% % % Xlagflag
% % % if Xlagflag > 0,
% % %     uvw(1,1:6000-Xlagflag) = uvw(1,1:6000-Xlagflag) +uvwplat(1,Xlagflag +1:6000);
% % %     uvw(2,1:6000-Xlagflag) = uvw(2,1:6000-Xlagflag) +uvwplat(2,Xlagflag +1:6000);
% % % else
% % %     uvw(1,-Xlagflag +1:6000) = uvw(1,-Xlagflag +1:6000) + uvwplat(1,1:6000+ Xlagflag);
% % %     uvw(2,-Xlagflag +1:6000) = uvw(2,-Xlagflag +1:6000) + uvwplat(2,1:6000+ Xlagflag);
% % % end
% % end    
% Compute motion corrected winds by checking lag time...
% [uvw(1,:),uvwplat(1,:),Xlagflag] = lagcorrection(uvw(1,:),uvwplat(1,:),16);
% [uvw(2,:),uvwplat(2,:),Xlagflag] = lagcorrection(uvw(2,:),uvwplat(2,:),16);
[uvw(3,:),uvwplat(3,:),Xlagflag] = lagcorrection(uvw(3,:),uvwplat(3,:),16);
Xlagflag
if Xlagflag > 0,
    uvw(1,1:length(uvw(3,:))-Xlagflag) = uvw(1,1:length(uvw(3,:))-Xlagflag) +uvwplat(1,Xlagflag +1:length(uvw(3,:)));
    uvw(2,1:length(uvw(3,:))-Xlagflag) = uvw(2,1:length(uvw(3,:))-Xlagflag) +uvwplat(2,Xlagflag +1:length(uvw(3,:)));
else
    uvw(1,-Xlagflag +1:length(uvw(3,:))) = uvw(1,-Xlagflag +1:length(uvw(3,:))) + uvwplat(1,1:length(uvw(3,:))+ Xlagflag);
    uvw(2,-Xlagflag +1:length(uvw(3,:))) = uvw(2,-Xlagflag +1:length(uvw(3,:))) + uvwplat(2,1:length(uvw(3,:))+ Xlagflag);
end

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

if decorr==1
    UVWdecorJE;
end;
% 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;   
% lagflag
DDDD=uvw;
DDDD=DDDD';

lagflagJE=Xlagflag;


% % %  Remove edge effects (due to filters) - shall remove that on u, v, w but
% % %  do it only on w here
% % npts=30*fs;              %30s will be removed on either end seconds for now
% % uvw(:,1:npts)=NaN;  uvw(:,end-npts:end)=NaN;


% % specsize_1=6000;
% % %some more tests
% % i=1;
% % [Cuwj_1,Fson_uwj]=csd2(detrend(uvw(3,00001:Nlendata)),detrend(uvw(i,00001:Nlendata)),specsize_1,sampleRate,hamming(specsize_1));[Cuwj_1,Fson_uwj]=specsmoo(Cuwj_1,sampleRate);
% % [Cuwj_2,Fson_uwj]=csd2(detrend(uvw(3,Nlendata+1:Nlendata*2)),detrend(uvw(i,Nlendata+1:Nlendata*2)),specsize_1,sampleRate,hamming(specsize_1));[Cuwj_2,Fson_uwj]=specsmoo(Cuwj_2,sampleRate);
% % [Cuwj_3,Fson_uwj]=csd2(detrend(uvw(3,Nlendata*2+1:Nlendata*3)),detrend(uvw(i,Nlendata*2+1:Nlendata*3)),specsize_1,sampleRate,hamming(specsize_1));[Cuwj_3,Fson_uwj]=specsmoo(Cuwj_3,sampleRate);
% % [Cuwj_4,Fson_uwj]=csd2(detrend(uvw(3,Nlendata*3+1:Nlendata*4)),detrend(uvw(i,Nlendata*3+1:Nlendata*4)),specsize_1,sampleRate,hamming(specsize_1));[Cuwj_4,Fson_uwj]=specsmoo(Cuwj_4,sampleRate);
% % [Cuwj_5,Fson_uwj]=csd2(detrend(uvw(3,Nlendata*4+1:Nlendata*5)),detrend(uvw(i,Nlendata*4+1:Nlendata*5)),specsize_1,sampleRate,hamming(specsize_1));[Cuwj_5,Fson_uwj]=specsmoo(Cuwj_5,sampleRate);
% % [Cuwj_6,Fson_uwj]=csd2(detrend(uvw(3,Nlendata*5+1:Nlendata*6)),detrend(uvw(i,Nlendata*5+1:Nlendata*6)),specsize_1,sampleRate,hamming(specsize_1));[Cuwj_6,Fson_uwj]=specsmoo(Cuwj_6,sampleRate);
% % yuj=mean([Cuwj_1' Cuwj_2' Cuwj_3' Cuwj_4' Cuwj_5' Cuwj_6']');
% % figure(7);semilogx(Fson_uwj,Fson_uwj.*yuj,'--')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%*****************************************************************
%            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

%matrix for a right-handed coordinate system
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;

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

if clearit
    clear Plat T OM ssx ssy ssz accy accz;
end;

maccx=nanmean(Platacc(1,:));
maccy=nanmean(Platacc(2,:));
maccz=nanmean(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;                 %60 seconds

tau=1/(Fc*2*pi);bc=Fc/(sampleRate/2);[bf af]=butter(1,bc);dt=1/sampleRate;
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/(sampleRate/2);[bf af]=butter(1,bc);dt=1/sampleRate;
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];
ii=find(isnan(fxyz));
if length(ii)>0 
    fxyz=[0 0 0 0 0 0];
end;
sdg=['save ', drive_w, '\fxyz.prn fxyz -ascii'];
disp(sdg)
eval(sdg)

if clearit
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;
end;

Result=Velrot;

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

for lag=1:16;
   Jcovtestlag=cov(Velrot(3,1:length(Velrot)-lag),Platvel(3,lag+1:length(Velrot)));
   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:-16;
   Jcovtestlag=cov(Velrot(3,-lag+1:length(Velrot)),Platvel(3,1:length(Velrot)+lag));
   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;

if clearit
clear Jcovtest Jcovtestlag lag
end;

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;
lagflag
lagflagJH=lagflag;
if clearit
clear lagflag
end;

if decorr==1
    UVWdecorJH;
end

lagmat=[ddd hhh lagflagJE lagflagJH];
sdg=['save ', drive_w '\HiWinGS2013_motcorrlagdata.prn lagmat -ascii -append'];
eval(sdg)

% [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 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

Us=0;
zi=600;

if isfinite(press)
    P=nanmean(press);
else
    P=1013;
end;

%******************************************************************

zu=16;zu_ship=15;    % ~ sensors height
zt=14;zt_ship=15;
zq=14;zq_ship=15;


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

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

mean_sog_1=nanmean(sogm(01:10));if mean_sog_1 < 0.2, mean_sog_1=0;end;std_sog_1=std(sogm(01:10));min_sog_1=min(sogm(01:10));max_sog_1=max(sogm(01:10));
mean_sog_2=nanmean(sogm(11:20));if mean_sog_2 < 0.2, mean_sog_2=0;end;std_sog_2=std(sogm(11:20));min_sog_2=min(sogm(11:20));max_sog_2=max(sogm(11:20));
mean_sog_3=nanmean(sogm(21:30));if mean_sog_3 < 0.2, mean_sog_3=0;end;std_sog_3=std(sogm(21:30));min_sog_3=min(sogm(21:30));max_sog_3=max(sogm(21:30));
mean_sog_4=nanmean(sogm(31:40));if mean_sog_4 < 0.2, mean_sog_4=0;end;std_sog_4=std(sogm(31:40));min_sog_4=min(sogm(31:40));max_sog_4=max(sogm(31:40));
mean_sog_5=nanmean(sogm(41:50));if mean_sog_5 < 0.2, mean_sog_5=0;end;std_sog_5=std(sogm(41:50));min_sog_5=min(sogm(41:50));max_sog_5=max(sogm(41:50));
mean_sog_6=nanmean(sogm(51:end));if mean_sog_6 < 0.2, mean_sog_6=0;end;std_sog_6=std(sogm(51:end));min_sog_6=min(sogm(51:end));max_sog_6=max(sogm(51:end));

mean_cog_1=mod(360+180.*atan2(mean(sin(cogm(01:10)./180.*pi)),mean(cos(cogm(01:10)./180.*pi)))./pi,360);
mean_cog_2=mod(360+180.*atan2(mean(sin(cogm(11:20)./180.*pi)),mean(cos(cogm(11:20)./180.*pi)))./pi,360);
mean_cog_3=mod(360+180.*atan2(mean(sin(cogm(21:30)./180.*pi)),mean(cos(cogm(21:30)./180.*pi)))./pi,360);
mean_cog_4=mod(360+180.*atan2(mean(sin(cogm(31:40)./180.*pi)),mean(cos(cogm(31:40)./180.*pi)))./pi,360);
mean_cog_5=mod(360+180.*atan2(mean(sin(cogm(41:50)./180.*pi)),mean(cos(cogm(41:50)./180.*pi)))./pi,360);
mean_cog_6=mod(360+180.*atan2(mean(sin(cogm(51:end)./180.*pi)),mean(cos(cogm(51:end)./180.*pi)))./pi,360);

std_cog_1=std(unwrap(cog(0001:0300)./180.*pi)*180/pi);
std_cog_2=std(unwrap(cog(0301:0600)./180.*pi)*180/pi);
std_cog_3=std(unwrap(cog(0601:0900)./180.*pi)*180/pi);
std_cog_4=std(unwrap(cog(0901:1200)./180.*pi)*180/pi);
std_cog_5=std(unwrap(cog(1201:1500)./180.*pi)*180/pi);
std_cog_6=std(unwrap(cog(1501:1800)./180.*pi)*180/pi);

max_cog_1=max(cog(0001:0300));
max_cog_2=max(cog(0301:0600));
max_cog_3=max(cog(0601:0900));
max_cog_4=max(cog(0901:1200));
max_cog_5=max(cog(1201:1500));
max_cog_6=max(cog(1501:1800));

min_cog_1=min(cog(0001:0300));
min_cog_2=min(cog(0301:0600));
min_cog_3=min(cog(0601:0900));
min_cog_4=min(cog(0901:1200));
min_cog_5=min(cog(1201:1500));
min_cog_6=min(cog(1501:1800));

mean_sinsog_1=nanmean(sogm(01:10).*sin(cogm(01:10)./180.*pi));std_sinsog_1=std(sogm(01:10).*sin(cogm(01:10)./180.*pi));min_sinsog_1=min(sogm(01:10).*sin(cogm(01:10)./180.*pi));max_sinsog_1=max(sogm(01:10).*sin(cogm(01:10)./180.*pi));
mean_sinsog_2=nanmean(sogm(11:20).*sin(cogm(11:20)./180.*pi));std_sinsog_2=std(sogm(11:20).*sin(cogm(11:20)./180.*pi));min_sinsog_2=min(sogm(11:20).*sin(cogm(11:20)./180.*pi));max_sinsog_2=max(sogm(11:20).*sin(cogm(11:20)./180.*pi));
mean_sinsog_3=nanmean(sogm(21:30).*sin(cogm(21:30)./180.*pi));std_sinsog_3=std(sogm(21:30).*sin(cogm(21:30)./180.*pi));min_sinsog_3=min(sogm(21:30).*sin(cogm(21:30)./180.*pi));max_sinsog_3=max(sogm(21:30).*sin(cogm(21:30)./180.*pi));
mean_sinsog_4=nanmean(sogm(31:40).*sin(cogm(31:40)./180.*pi));std_sinsog_4=std(sogm(31:40).*sin(cogm(31:40)./180.*pi));min_sinsog_4=min(sogm(31:40).*sin(cogm(31:40)./180.*pi));max_sinsog_4=max(sogm(31:40).*sin(cogm(31:40)./180.*pi));
mean_sinsog_5=nanmean(sogm(41:50).*sin(cogm(41:50)./180.*pi));std_sinsog_5=std(sogm(41:50).*sin(cogm(41:50)./180.*pi));min_sinsog_5=min(sogm(41:50).*sin(cogm(41:50)./180.*pi));max_sinsog_5=max(sogm(41:50).*sin(cogm(41:50)./180.*pi));
mean_sinsog_6=nanmean(sogm(51:end).*sin(cogm(51:end)./180.*pi));std_sinsog_6=std(sogm(51:end).*sin(cogm(51:end)./180.*pi));min_sinsog_6=min(sogm(51:end).*sin(cogm(51:end)./180.*pi));max_sinsog_6=max(sogm(51:end).*sin(cogm(51:end)./180.*pi));

mean_cossog_1=nanmean(sogm(01:10).*cos(cogm(01:10)./180.*pi));std_cossog_1=std(sogm(01:10).*cos(cogm(01:10)./180.*pi));min_cossog_1=min(sogm(01:10).*cos(cogm(01:10)./180.*pi));max_cossog_1=max(sogm(01:10).*cos(cogm(01:10)./180.*pi));
mean_cossog_2=nanmean(sogm(11:20).*cos(cogm(11:20)./180.*pi));std_cossog_2=std(sogm(11:20).*cos(cogm(11:20)./180.*pi));min_cossog_2=min(sogm(11:20).*cos(cogm(11:20)./180.*pi));max_cossog_2=max(sogm(11:20).*cos(cogm(11:20)./180.*pi));
mean_cossog_3=nanmean(sogm(21:30).*cos(cogm(21:30)./180.*pi));std_cossog_3=std(sogm(21:30).*cos(cogm(21:30)./180.*pi));min_cossog_3=min(sogm(21:30).*cos(cogm(21:30)./180.*pi));max_cossog_3=max(sogm(21:30).*cos(cogm(21:30)./180.*pi));
mean_cossog_4=nanmean(sogm(31:40).*cos(cogm(31:40)./180.*pi));std_cossog_4=std(sogm(31:40).*cos(cogm(31:40)./180.*pi));min_cossog_4=min(sogm(31:40).*cos(cogm(31:40)./180.*pi));max_cossog_4=max(sogm(31:40).*cos(cogm(31:40)./180.*pi));
mean_cossog_5=nanmean(sogm(41:50).*cos(cogm(41:50)./180.*pi));std_cossog_5=std(sogm(41:50).*cos(cogm(41:50)./180.*pi));min_cossog_5=min(sogm(41:50).*cos(cogm(41:50)./180.*pi));max_cossog_5=max(sogm(41:50).*cos(cogm(41:50)./180.*pi));
mean_cossog_6=nanmean(sogm(51:end).*cos(cogm(51:end)./180.*pi));std_cossog_6=std(sogm(51:end).*cos(cogm(51:end)./180.*pi));min_cossog_6=min(sogm(51:end).*cos(cogm(51:end)./180.*pi));max_cossog_6=max(sogm(51:end).*cos(cogm(51:end)./180.*pi));

%remove the ship's speed at higher resolution rather than just the mean...
Ngps=sog.*cos(cog./180.*pi)./1.946;
Egps=sog.*sin(cog./180.*pi)./1.946;
% shipn=sog.*cos(cog./180.*pi)./1.946;
% shipe=sog.*sin(cog./180.*pi)./1.946;
% tgyro=(0:2:3599);
% % Ngps=interp1(tgyro,shipn,tref,'spline','extrap');
% % Egps=interp1(tgyro,shipe,tref,'spline','extrap');
% Ngps=interp1(tgyro,shipn,tref);
% Egps=interp1(tgyro,shipe,tref);
%correct for some bad points due to spline+extrap at the end.
% Ngps(tgyro(end)*10:36000)=Ngps(tgyro(end)*10);
% Egps(tgyro(end)*10:36000)=Egps(tgyro(end)*10);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fs=10;         %Sonic are output at 10 Hz
fc=1/200;       %Cutoff frequency for angle calculation

%    DEFINE A HIGH PASS FILTER WHICH RETAINS 
%    REAL ACCELERATION BUT REMOVES DRIFT

% nfreq = fs/2.0;                   %Nyquist frequency   
% wp = fc/nfreq;          
% ws = .7*wp;   
% [n,wn] = buttord(wp,ws,10,25);        %Filter coefficients
% [bhi,ahi] = butter(n,wn,'high');  %High order Butterworth filter
% Ngps = Ngps - filtfilt(bhi,ahi,Ngps);
% Egps = Egps - filtfilt(bhi,ahi,Egps);

wp=1/25/(sampleRate/2);
ws=0.8*wp;
% [n,wn]=buttord(wp,ws,10,25);
[n,wn]=buttord(wp,ws,3,7);
% wp=1/120/(sampleRate/2);
% ws=.8*wp;
% [n,wn]=buttord(wp,ws,3,5);
% wp=1/240/(sampleRate/2);
% ws=1/120/(sampleRate/2);
% [n,wn]=buttord(wp,ws,1,14);
% [n,wn]=buttord(wp,ws,9,30);
% % %disp([num2str(ws),'  ',num2str(wn),'  ',num2str(wp)])
[bhi,ahi]=butter(n,wn,'low');

% switch leg(1:end-1)
%     case {'leg1','leg2','leg3','leg4'}
%         %uses WRT water here for legs2and3
%         CN=Cspd.*cos(Cdir./180.*pi);  %'from'...(-Cn) !!!!!!!!!!!!!!!!!!!
%         CE=Cspd.*sin(Cdir./180.*pi);  %'from'... (-Ce)      !!!!!!!!!!!!!!!!!!!
%         Ngps = filtfilt(bhi,ahi,Ngps);
%         Egps = filtfilt(bhi,ahi,Egps);
%         
%         if ddd==310 | ddd==311 | ddd==342 %| length(CN)<6  %no nav data yet
%             CN=zeros(1,6);CE=zeros(1,6);
%         end;
%         
%         Ngps(1,00001:Nlendata)=Ngps(1,00001:Nlendata)+CN(1);
%         Ngps(1,Nlendata+1:Nlendata*2)=Ngps(1,Nlendata+1:Nlendata*2)+CN(2);
%         Ngps(1,Nlendata*2+1:Nlendata*3)=Ngps(1,Nlendata*2+1:Nlendata*3)+CN(3);
%         Ngps(1,Nlendata*3+1:Nlendata*4)=Ngps(1,Nlendata*3+1:Nlendata*4)+CN(4);
%         Ngps(1,Nlendata*4+1:Nlendata*5)=Ngps(1,Nlendata*4+1:Nlendata*5)+CN(5);
%         Ngps(1,Nlendata*5+1:Nlendata*6)=Ngps(1,Nlendata*5+1:Nlendata*6)+CN(6);   
%       
%         Egps(1,00001:Nlendata)=Egps(1,00001:Nlendata)+CE(1);
%         Egps(1,Nlendata+1:Nlendata*2)=Egps(1,Nlendata+1:Nlendata*2)+CE(2);
%         Egps(1,Nlendata*2+1:Nlendata*3)=Egps(1,Nlendata*2+1:Nlendata*3)+CE(3);
%         Egps(1,Nlendata*3+1:Nlendata*4)=Egps(1,Nlendata*3+1:Nlendata*4)+CE(4);
%         Egps(1,Nlendata*4+1:Nlendata*5)=Egps(1,Nlendata*4+1:Nlendata*5)+CE(5);
%         Egps(1,Nlendata*5+1:Nlendata*6)=Egps(1,Nlendata*5+1:Nlendata*6)+CE(6);   

%     case {'leg1'}   
        Ngps = filtfilt(bhi,ahi,Ngps);
        Egps = filtfilt(bhi,ahi,Egps);      
% end
%LB- 10/22/2012 - watch out sign convention 'from' and 'to'...
% Result_u=Result_u+Ngps=Result_u+Ngps-Cn
% Result_v=Result_v-Egps=Result_v-Egps+Ce

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [Sc1,Fc1]=psd2((shipn),1800,0.5,hamming(1800));[Sc1,Fc1]=specsmoo(Sc1,0.5);
% [Sc2,Fc2]=psd2((Ngps),36000,10,hamming(36000));[Sc2,Fc2]=specsmoo(Sc2,10);
% [Se1,Fe1]=psd2((shipe),1800,0.5,hamming(1800));[Se1,Fe1]=specsmoo(Se1,0.5);
% [Se2,Fe2]=psd2((Egps),36000,10,hamming(36000));[Se2,Fe2]=specsmoo(Se2,10);

Result(1,:)=Result(1,:)+Ngps;
Result(2,:)=Result(2,:)-Egps;

Result(1,00001:Nlendata)=Result(1,00001:Nlendata)-mean(Platvel(1,00001:Nlendata));
Result(1,Nlendata+1:Nlendata*2)=Result(1,Nlendata+1:Nlendata*2)-mean(Platvel(1,Nlendata+1:Nlendata*2));
Result(1,Nlendata*2+1:Nlendata*3)=Result(1,Nlendata*2+1:Nlendata*3)-mean(Platvel(1,Nlendata*2+1:Nlendata*3));
Result(1,Nlendata*3+1:Nlendata*4)=Result(1,Nlendata*3+1:Nlendata*4)-mean(Platvel(1,Nlendata*3+1:Nlendata*4));
Result(1,Nlendata*4+1:Nlendata*5)=Result(1,Nlendata*4+1:Nlendata*5)-mean(Platvel(1,Nlendata*4+1:Nlendata*5));
Result(1,Nlendata*5+1:Nlendata*6)=Result(1,Nlendata*5+1:Nlendata*6)-mean(Platvel(1,Nlendata*5+1:Nlendata*6));   %replace + by - 18 October 2005
Result(2,00001:Nlendata)=Result(2,00001:Nlendata)-mean(Platvel(2,00001:Nlendata));
Result(2,Nlendata+1:Nlendata*2)=Result(2,Nlendata+1:Nlendata*2)-mean(Platvel(2,Nlendata+1:Nlendata*2));
Result(2,Nlendata*2+1:Nlendata*3)=Result(2,Nlendata*2+1:Nlendata*3)-mean(Platvel(2,Nlendata*2+1:Nlendata*3));
Result(2,Nlendata*3+1:Nlendata*4)=Result(2,Nlendata*3+1:Nlendata*4)-mean(Platvel(2,Nlendata*3+1:Nlendata*4));
Result(2,Nlendata*4+1:Nlendata*5)=Result(2,Nlendata*4+1:Nlendata*5)-mean(Platvel(2,Nlendata*4+1:Nlendata*5));
Result(2,Nlendata*5+1:Nlendata*6)=Result(2,Nlendata*5+1:Nlendata*6)-mean(Platvel(2,Nlendata*5+1:Nlendata*6));
Result(3,00001:Nlendata)=Result(3,00001:Nlendata)-mean(Platvel(3,00001:Nlendata));
Result(3,Nlendata+1:Nlendata*2)=Result(3,Nlendata+1:Nlendata*2)-mean(Platvel(3,Nlendata+1:Nlendata*2));
Result(3,Nlendata*2+1:Nlendata*3)=Result(3,Nlendata*2+1:Nlendata*3)-mean(Platvel(3,Nlendata*2+1:Nlendata*3));
Result(3,Nlendata*3+1:Nlendata*4)=Result(3,Nlendata*3+1:Nlendata*4)-mean(Platvel(3,Nlendata*3+1:Nlendata*4));
Result(3,Nlendata*4+1:Nlendata*5)=Result(3,Nlendata*4+1:Nlendata*5)-mean(Platvel(3,Nlendata*4+1:Nlendata*5));
Result(3,Nlendata*5+1:Nlendata*6)=Result(3,Nlendata*5+1:Nlendata*6)-mean(Platvel(3,Nlendata*5+1:Nlendata*6));





%############ Motion correction view #################
if hhh==12 | hhh==00
    figure(1);clf;plot(tmotion(1:length(tmotion)),Result(1,1:length(tmotion)),'g',tmotion(1:length(tmotion)),Platvel(1,1:length(tmotion)),'r',tmotion(1:length(tmotion)),Velrot(1,1:length(tmotion)))
    title([num2str(ddd),' ', num2str(hhh),' blue= u raw, green= u corrected, red= u correction, Hare']);axis([0 100 -25 25])
    if printit==1
        print(graphdevice,[drive_w '\u_Harecorrection_' num2str(ddd) num2str(hhh) graphformat]);
    end;
    figure(2);clf;plot(tmotion(1:length(tmotion)),Result(2,1:length(tmotion)),'g',tmotion(1:length(tmotion)),Platvel(2,1:length(tmotion)),'r',tmotion(1:length(tmotion)),Velrot(2,1:length(tmotion)))
    title([num2str(ddd),' ', num2str(hhh),' blue= v raw, green= v corrected, red= v correction, Hare']);axis([0 100 -25 25])
    if printit==1
        print(graphdevice,[drive_w '\v_Harecorrection_' num2str(ddd) num2str(hhh) graphformat]);
    end;
    figure(3);clf;plot(tmotion(1:length(tmotion)),Result(3,1:length(tmotion))-mean(Result(3,1:length(tmotion))),'g',tmotion(1:length(tmotion)),Platvel(3,1:length(tmotion)),'r',tmotion(1:length(tmotion)),Velrot(3,1:length(tmotion))-mean(Velrot(3,1:length(tmotion))))
    hold on;plot([1 3600],[0 0],'k--');grid
    title([num2str(ddd),' ', num2str(hhh),' blue= w raw, green= w corrected, red= w correction, Hare']);axis([0 100 -6 6])
    if printit==1
        print(graphdevice,[drive_w '\w_Harecorrection_' num2str(ddd) num2str(hhh) graphformat]);
    end;
end
%#####################################################


uvw(1,:)=uvw(1,:)+Ngps;
uvw(2,:)=uvw(2,:)-Egps;

uvw(1,00001:Nlendata)=uvw(1,00001:Nlendata)-mean(uvwplat(1,00001:Nlendata));
uvw(1,Nlendata+1:Nlendata*2)=uvw(1,Nlendata+1:Nlendata*2)-mean(uvwplat(1,Nlendata+1:Nlendata*2));
uvw(1,Nlendata*2+1:Nlendata*3)=uvw(1,Nlendata*2+1:Nlendata*3)-mean(uvwplat(1,Nlendata*2+1:Nlendata*3));
uvw(1,Nlendata*3+1:Nlendata*4)=uvw(1,Nlendata*3+1:Nlendata*4)-mean(uvwplat(1,Nlendata*3+1:Nlendata*4));
uvw(1,Nlendata*4+1:Nlendata*5)=uvw(1,Nlendata*4+1:Nlendata*5)-mean(uvwplat(1,Nlendata*4+1:Nlendata*5));
uvw(1,Nlendata*5+1:Nlendata*6)=uvw(1,Nlendata*5+1:Nlendata*6)-mean(uvwplat(1,Nlendata*5+1:Nlendata*6));   %replace + by - 18 October 2005
uvw(2,00001:Nlendata)=uvw(2,00001:Nlendata)-mean(uvwplat(2,00001:Nlendata));
uvw(2,Nlendata+1:Nlendata*2)=uvw(2,Nlendata+1:Nlendata*2)-mean(uvwplat(2,Nlendata+1:Nlendata*2));
uvw(2,Nlendata*2+1:Nlendata*3)=uvw(2,Nlendata*2+1:Nlendata*3)-mean(uvwplat(2,Nlendata*2+1:Nlendata*3));
uvw(2,Nlendata*3+1:Nlendata*4)=uvw(2,Nlendata*3+1:Nlendata*4)-mean(uvwplat(2,Nlendata*3+1:Nlendata*4));
uvw(2,Nlendata*4+1:Nlendata*5)=uvw(2,Nlendata*4+1:Nlendata*5)-mean(uvwplat(2,Nlendata*4+1:Nlendata*5));
uvw(2,Nlendata*5+1:Nlendata*6)=uvw(2,Nlendata*5+1:Nlendata*6)-mean(uvwplat(2,Nlendata*5+1:Nlendata*6));
uvw(3,00001:Nlendata)=uvw(3,00001:Nlendata)-mean(uvwplat(3,00001:Nlendata));
uvw(3,Nlendata+1:Nlendata*2)=uvw(3,Nlendata+1:Nlendata*2)-mean(uvwplat(3,Nlendata+1:Nlendata*2));
uvw(3,Nlendata*2+1:Nlendata*3)=uvw(3,Nlendata*2+1:Nlendata*3)-mean(uvwplat(3,Nlendata*2+1:Nlendata*3));
uvw(3,Nlendata*3+1:Nlendata*4)=uvw(3,Nlendata*3+1:Nlendata*4)-mean(uvwplat(3,Nlendata*3+1:Nlendata*4));
uvw(3,Nlendata*4+1:Nlendata*5)=uvw(3,Nlendata*4+1:Nlendata*5)-mean(uvwplat(3,Nlendata*4+1:Nlendata*5));
uvw(3,Nlendata*5+1:Nlendata*6)=uvw(3,Nlendata*5+1:Nlendata*6)-mean(uvwplat(3,Nlendata*5+1:Nlendata*6));

% %     save motion data for other flux computation
% Reldirby=atan2(V,-U)*180/pi;ii=find(Reldirby<0);Reldirby(ii)=Reldirby(ii)+360;
% [Uj,Vj,Wj,alphaj,betaj]=rotcorr(uvw(1,:),uvw(2,:),uvw(3,:));
% uvwmat=[tmotion' Uj' Vj' Wj' Reldirby'];
% sdg=['save ', drive_w '\MotionData\DYNAMO2011_UVWstream_',jd hr,'.prn uvwmat -ascii'];
% eval(sdg)


%############ Motion correction view #################
if hhh==12 | hhh==00
    figure(4);clf;plot(tmotion,uvw(1,:),'g',tmotion,uvwplat(1,:),'r',tmotion,AAAA(1,:))
    title([num2str(ddd),' ', num2str(hhh),' blue= u raw, green= u corrected, red= u correction, Edson']);axis([0 100 -25 25])
    if printit==1
        print(graphdevice,[drive_w '\u_Edsoncorrection_' num2str(ddd) num2str(hhh) graphformat]);
    end;    
    figure(5);clf;plot(tmotion,uvw(2,:),'g',tmotion,uvwplat(2,:),'r',tmotion,AAAA(2,:))
    title([num2str(ddd),' ', num2str(hhh),' blue= v raw, green= v corrected, red= v correction, Edson']);axis([0 100 -25 25])
    if printit==1
        print(graphdevice,[drive_w '\v_Edsoncorrection_' num2str(ddd) num2str(hhh) graphformat]);
    end;  
    figure(6);clf;plot(tmotion,uvw(3,:)-mean(uvw(3,:)),'g',tmotion,uvwplat(3,:),'r',tmotion,AAAA(3,:)-mean(AAAA(3,:)))
     hold on;plot([1 3600],[0 0],'k--');grid
    title([num2str(ddd),' ', num2str(hhh),' blue= w raw, green= w corrected, red= w correction, Edson']);axis([0 100 -6 6])
    if printit==1
        print(graphdevice,[drive_w '\w_Edsoncorrection_' num2str(ddd) num2str(hhh) graphformat]);
    end;  
end;
%######################################################

[Ustream_1,Vstream_1,Wstream_1,alpha_1,beta_1]=rotcorr(Result(1,00001:Nlendata),Result(2,00001:Nlendata),Result(3,00001:Nlendata));           %converts [u,v,w] to [Ustream, Vstream, Wstream] for ten minutes
[Ustream_2,Vstream_2,Wstream_2,alpha_2,beta_2]=rotcorr(Result(1,Nlendata+1:Nlendata*2),Result(2,Nlendata+1:Nlendata*2),Result(3,Nlendata+1:Nlendata*2));
[Ustream_3,Vstream_3,Wstream_3,alpha_3,beta_3]=rotcorr(Result(1,Nlendata*2+1:Nlendata*3),Result(2,Nlendata*2+1:Nlendata*3),Result(3,Nlendata*2+1:Nlendata*3));
[Ustream_4,Vstream_4,Wstream_4,alpha_4,beta_4]=rotcorr(Result(1,Nlendata*3+1:Nlendata*4),Result(2,Nlendata*3+1:Nlendata*4),Result(3,Nlendata*3+1:Nlendata*4));
[Ustream_5,Vstream_5,Wstream_5,alpha_5,beta_5]=rotcorr(Result(1,Nlendata*4+1:Nlendata*5),Result(2,Nlendata*4+1:Nlendata*5),Result(3,Nlendata*4+1:Nlendata*5));
[Ustream_6,Vstream_6,Wstream_6,alpha_6,beta_6]=rotcorr(Result(1,Nlendata*5+1:Nlendata*6),Result(2,Nlendata*5+1:Nlendata*6),Result(3,Nlendata*5+1:Nlendata*6));

Tsonic_1=Tsonic(00001:Nlendata);
Tsonic_2=Tsonic(Nlendata+1:Nlendata*2);
Tsonic_3=Tsonic(Nlendata*2+1:Nlendata*3);
Tsonic_4=Tsonic(Nlendata*3+1:Nlendata*4);
Tsonic_5=Tsonic(Nlendata*4+1:Nlendata*5);
Tsonic_6=Tsonic(Nlendata*5+1:Nlendata*6);

Licor_CO2_1=Licor_CO2(0*Nlenlic+1:1*Nlenlic);
Licor_CO2_2=Licor_CO2(1*Nlenlic+1:2*Nlenlic);
Licor_CO2_3=Licor_CO2(2*Nlenlic+1:3*Nlenlic);
Licor_CO2_4=Licor_CO2(3*Nlenlic+1:4*Nlenlic);
Licor_CO2_5=Licor_CO2(4*Nlenlic+1:5*Nlenlic);
Licor_CO2_6=Licor_CO2(5*Nlenlic+1:length(Licor_CO2));

Licor_H2O_1=Licor_H2O(0*Nlenlic+1:1*Nlenlic);
Licor_H2O_2=Licor_H2O(1*Nlenlic+1:2*Nlenlic);
Licor_H2O_3=Licor_H2O(2*Nlenlic+1:3*Nlenlic);
Licor_H2O_4=Licor_H2O(3*Nlenlic+1:4*Nlenlic);
Licor_H2O_5=Licor_H2O(4*Nlenlic+1:5*Nlenlic);
Licor_H2O_6=Licor_H2O(5*Nlenlic+1:length(Licor_CO2));

if Licor_Temp(1)>40,
   Licor_Temp(1)=25;
end;
if Licor_Temp(1)<10,
   Licor_Temp(1)=25;
end;

if Licor_Press(1)>120,
   Licor_Press(1)=100.5;
end;
if Licor_Press(1)<90,
   Licor_Press(1)=100.5;
end;

for ii=2:length(Licor_Temp);
    if Licor_Temp(ii)>40,
        Licor_Temp(ii)=Licor_Temp(ii-1);
    end;
    if Licor_Temp(ii)<10,
        Licor_Temp(ii)=Licor_Temp(ii-1);
    end;
end;

for ii=2:length(Licor_Press);
    if Licor_Press(ii)>40,
        Licor_Press(ii)=Licor_Press(ii-1);
    end;
    if Licor_Press(ii)<10,
        Licor_Press(ii)=Licor_Press(ii-1);
    end;
end;
Licor_CO2_1=Licor_CO2(0*Nlenlic+1:1*Nlenlic);

Licor_Temp_1=Licor_Temp(0*Nlenlic+1:1*Nlenlic);
Licor_Temp_2=Licor_Temp(1*Nlenlic+1:2*Nlenlic);
Licor_Temp_3=Licor_Temp(2*Nlenlic+1:3*Nlenlic);
Licor_Temp_4=Licor_Temp(3*Nlenlic+1:4*Nlenlic);
Licor_Temp_5=Licor_Temp(4*Nlenlic+1:5*Nlenlic);
Licor_Temp_6=Licor_Temp(5*Nlenlic+1:length(Licor_Temp));

Licor_Press_1=Licor_Press(0*Nlenlic+1:1*Nlenlic);
Licor_Press_2=Licor_Press(1*Nlenlic+1:2*Nlenlic);
Licor_Press_3=Licor_Press(2*Nlenlic+1:3*Nlenlic);
Licor_Press_4=Licor_Press(3*Nlenlic+1:4*Nlenlic);
Licor_Press_5=Licor_Press(4*Nlenlic+1:5*Nlenlic);
Licor_Press_6=Licor_Press(5*Nlenlic+1:length(Licor_Press));

[Uj_1,Vj_1,Wj_1,alphaj_1,betaj_1]=rotcorr(uvw(1,00001:Nlendata),uvw(2,00001:Nlendata),uvw(3,00001:Nlendata));
[Uj_2,Vj_2,Wj_2,alphaj_2,betaj_2]=rotcorr(uvw(1,Nlendata+1:Nlendata*2),uvw(2,Nlendata+1:Nlendata*2),uvw(3,Nlendata+1:Nlendata*2));
[Uj_3,Vj_3,Wj_3,alphaj_3,betaj_3]=rotcorr(uvw(1,Nlendata*2+1:Nlendata*3),uvw(2,Nlendata*2+1:Nlendata*3),uvw(3,Nlendata*2+1:Nlendata*3));
[Uj_4,Vj_4,Wj_4,alphaj_4,betaj_4]=rotcorr(uvw(1,Nlendata*3+1:Nlendata*4),uvw(2,Nlendata*3+1:Nlendata*4),uvw(3,Nlendata*3+1:Nlendata*4));
[Uj_5,Vj_5,Wj_5,alphaj_5,betaj_5]=rotcorr(uvw(1,Nlendata*4+1:Nlendata*5),uvw(2,Nlendata*4+1:Nlendata*5),uvw(3,Nlendata*4+1:Nlendata*5));
[Uj_6,Vj_6,Wj_6,alphaj_6,betaj_6]=rotcorr(uvw(1,Nlendata*5+1:Nlendata*6),uvw(2,Nlendata*5+1:Nlendata*6),uvw(3,Nlendata*5+1:Nlendata*6));

Ubarj_1=nanmean(Uj_1);
Ubarj_2=nanmean(Uj_2);
Ubarj_3=nanmean(Uj_3);
Ubarj_4=nanmean(Uj_4);
Ubarj_5=nanmean(Uj_5);
Ubarj_6=nanmean(Uj_6);

press_1=nanmean(press(1:10));
press_2=nanmean(press(11:20));
press_3=nanmean(press(21:30));
press_4=nanmean(press(31:40));
press_5=nanmean(press(41:50));
press_6=nanmean(press(51:end));

Tbar_1=nanmean(Tvais(1:10));
Tbar_2=nanmean(Tvais(11:20));
Tbar_3=nanmean(Tvais(21:30));
Tbar_4=nanmean(Tvais(31:40));
Tbar_5=nanmean(Tvais(41:50));
Tbar_6=nanmean(Tvais(51:end));

Tship_1=nanmean(imta(1:10));
Tship_2=nanmean(imta(11:20));
Tship_3=nanmean(imta(21:30));
Tship_4=nanmean(imta(31:40));
Tship_5=nanmean(imta(41:50));
Tship_6=nanmean(imta(51:end));

rho_1=P*100./(273.15+Tbar_1)./287;
rho_2=P*100./(273.15+Tbar_2)./287;
rho_3=P*100./(273.15+Tbar_3)./287;
rho_4=P*100./(273.15+Tbar_4)./287;
rho_5=P*100./(273.15+Tbar_5)./287;
rho_6=P*100./(273.15+Tbar_6)./287;

wu_1=cov(Ustream_1,Wstream_1);
wu_2=cov(Ustream_2,Wstream_2);
wu_3=cov(Ustream_3,Wstream_3);
wu_4=cov(Ustream_4,Wstream_4);
wu_5=cov(Ustream_5,Wstream_5);
wu_6=cov(Ustream_6,Wstream_6);

wv_1=cov(Vstream_1,Wstream_1);
wv_2=cov(Vstream_2,Wstream_2);
wv_3=cov(Vstream_3,Wstream_3);
wv_4=cov(Vstream_4,Wstream_4);
wv_5=cov(Vstream_5,Wstream_5);
wv_6=cov(Vstream_6,Wstream_6);

wt_1=cov(Wstream_1,Tsonic_1);
wt_2=cov(Wstream_2,Tsonic_2);
wt_3=cov(Wstream_3,Tsonic_3);
wt_4=cov(Wstream_4,Tsonic_4);
wt_5=cov(Wstream_5,Tsonic_5);
wt_6=cov(Wstream_6,Tsonic_6);

Ustr_1=sqrt(nanmean(U(00001:Nlendata)).^2 + nanmean(V(00001:Nlendata)).^2 + nanmean(W(00001:Nlendata)).^2);
Ustr_2=sqrt(nanmean(U(Nlendata+1:Nlendata*2)).^2 + nanmean(V(Nlendata+1:Nlendata*2)).^2 + nanmean(W(Nlendata+1:Nlendata*2)).^2);
Ustr_3=sqrt(nanmean(U(Nlendata*2+1:Nlendata*3)).^2 + nanmean(V(Nlendata*2+1:Nlendata*3)).^2 + nanmean(W(Nlendata*2+1:Nlendata*3)).^2);
Ustr_4=sqrt(nanmean(U(Nlendata*3+1:Nlendata*4)).^2 + nanmean(V(Nlendata*3+1:Nlendata*4)).^2 + nanmean(W(Nlendata*3+1:Nlendata*4)).^2);
Ustr_5=sqrt(nanmean(U(Nlendata*4+1:Nlendata*5)).^2 + nanmean(V(Nlendata*4+1:Nlendata*5)).^2 + nanmean(W(Nlendata*4+1:Nlendata*5)).^2);
Ustr_6=sqrt(nanmean(U(Nlendata*5+1:Nlendata*6)).^2 + nanmean(V(Nlendata*5+1:Nlendata*6)).^2 + nanmean(W(Nlendata*5+1:Nlendata*6)).^2);

wuj_1=cov(Uj_1,Wj_1);
wuj_2=cov(Uj_2,Wj_2);
wuj_3=cov(Uj_3,Wj_3);
wuj_4=cov(Uj_4,Wj_4);
wuj_5=cov(Uj_5,Wj_5);
wuj_6=cov(Uj_6,Wj_6);

wvj_1=cov(Vj_1,Wj_1);
wvj_2=cov(Vj_2,Wj_2);
wvj_3=cov(Vj_3,Wj_3);
wvj_4=cov(Vj_4,Wj_4);
wvj_5=cov(Vj_5,Wj_5);
wvj_6=cov(Vj_6,Wj_6);

wtj_1=cov(Wj_1,Tsonic_1);
wtj_2=cov(Wj_2,Tsonic_2);
wtj_3=cov(Wj_3,Tsonic_3);
wtj_4=cov(Wj_4,Tsonic_4);
wtj_5=cov(Wj_5,Tsonic_5);
wtj_6=cov(Wj_6,Tsonic_6);

    
Licor_w_1=Wj_1;
Licor_w_2=Wj_2;
Licor_w_3=Wj_3;
Licor_w_4=Wj_4;
Licor_w_5=Wj_5;
Licor_w_6=Wj_6;

I_1=Licor_CO2_1.*0;
I_2=Licor_CO2_2.*0;
I_3=Licor_CO2_3.*0;
I_4=Licor_CO2_4.*0;
I_5=Licor_CO2_5.*0;
I_6=Licor_CO2_6.*0;

for ii=1:length(Licor_CO2_1);
    if Licor_CO2_1(ii)>90000;
        I_1(ii)=1;
    end;
end;
for ii=1:length(Licor_CO2_2);
    if Licor_CO2_2(ii)>90000;
        I_2(ii)=1;
    end;
end;
for ii=1:length(Licor_CO2_3);
    if Licor_CO2_3(ii)>90000;
        I_3(ii)=1;
    end;
end;
for ii=1:length(Licor_CO2_4);
    if Licor_CO2_4(ii)>90000;
        I_4(ii)=1;
    end;
end;
for ii=1:length(Licor_CO2_5);
    if Licor_CO2_5(ii)>90000;
        I_5(ii)=1;
    end;
end;
for ii=1:length(Licor_CO2_6);
    if Licor_CO2_6(ii)>90000;
        I_6(ii)=1;
    end;
end;

Licor_CO2_1=Licor_CO2_1.*8.314.*(nanmean(Licor_Temp_1)+273.15)./nanmean(Licor_Press_1);   %convert from mmol/m3 to umol/mol
Licor_CO2_2=Licor_CO2_2.*8.314.*(nanmean(Licor_Temp_2)+273.15)./nanmean(Licor_Press_2);   %convert from mmol/m3 to umol/mol
Licor_CO2_3=Licor_CO2_3.*8.314.*(nanmean(Licor_Temp_3)+273.15)./nanmean(Licor_Press_3);   %convert from mmol/m3 to umol/mol
Licor_CO2_4=Licor_CO2_4.*8.314.*(nanmean(Licor_Temp_4)+273.15)./nanmean(Licor_Press_4);   %convert from mmol/m3 to umol/mol
Licor_CO2_5=Licor_CO2_5.*8.314.*(nanmean(Licor_Temp_5)+273.15)./nanmean(Licor_Press_5);   %convert from mmol/m3 to umol/mol
Licor_CO2_6=Licor_CO2_6.*8.314.*(nanmean(Licor_Temp_6)+273.15)./nanmean(Licor_Press_6);   %convert from mmol/m3 to umol/mol

Licor_H2O_1=Licor_H2O_1.*18./1000./rho_1;         %convert from mmol/m3 to g/kg
Licor_H2O_2=Licor_H2O_2.*18./1000./rho_2;         %convert from mmol/m3 to g/kg
Licor_H2O_3=Licor_H2O_3.*18./1000./rho_3;         %convert from mmol/m3 to g/kg
Licor_H2O_4=Licor_H2O_4.*18./1000./rho_4;         %convert from mmol/m3 to g/kg
Licor_H2O_5=Licor_H2O_5.*18./1000./rho_5;         %convert from mmol/m3 to g/kg
Licor_H2O_6=Licor_H2O_6.*18./1000./rho_6;         %convert from mmol/m3 to g/kg

%*  changed to /, cwf 12/04
Licor_CO2_1=Licor_CO2_1./(1-Licor_H2O_1./18.*1000.*rho_1.*8.314.*(nanmean(Licor_Temp_1)+273.15)./1000./nanmean(Licor_Press_1)./1000);   %dilution correction
Licor_CO2_2=Licor_CO2_2./(1-Licor_H2O_2./18.*1000.*rho_2.*8.314.*(nanmean(Licor_Temp_2)+273.15)./1000./nanmean(Licor_Press_2)./1000);   %dilution correction
Licor_CO2_3=Licor_CO2_3./(1-Licor_H2O_3./18.*1000.*rho_3.*8.314.*(nanmean(Licor_Temp_3)+273.15)./1000./nanmean(Licor_Press_3)./1000);   %dilution correction
Licor_CO2_4=Licor_CO2_4./(1-Licor_H2O_4./18.*1000.*rho_4.*8.314.*(nanmean(Licor_Temp_4)+273.15)./1000./nanmean(Licor_Press_4)./1000);   %dilution correction
Licor_CO2_5=Licor_CO2_5./(1-Licor_H2O_5./18.*1000.*rho_5.*8.314.*(nanmean(Licor_Temp_5)+273.15)./1000./nanmean(Licor_Press_5)./1000);   %dilution correction
Licor_CO2_6=Licor_CO2_6./(1-Licor_H2O_6./18.*1000.*rho_6.*8.314.*(nanmean(Licor_Temp_6)+273.15)./1000./nanmean(Licor_Press_6)./1000);   %dilution correction

LicCO2_1=[];LicH2O_1=[];WCO2_1=[];
LicCO2_2=[];LicH2O_2=[];WCO2_2=[];
LicCO2_3=[];LicH2O_3=[];WCO2_3=[];
LicCO2_4=[];LicH2O_4=[];WCO2_4=[];
LicCO2_5=[];LicH2O_5=[];WCO2_5=[];
LicCO2_6=[];LicH2O_6=[];WCO2_6=[];

for ii=1:length(Licor_CO2_1);
    if I_1(ii)==0,
        LicCO2_1=[LicCO2_1 Licor_CO2_1(ii)];
        LicH2O_1=[LicH2O_1 Licor_H2O_1(ii)];
        WCO2_1=[WCO2_1 Licor_w_1(ii)];
    end;
end;
for ii=1:length(Licor_CO2_2);
    if I_2(ii)==0,
        LicCO2_2=[LicCO2_2 Licor_CO2_2(ii)];
        LicH2O_2=[LicH2O_2 Licor_H2O_2(ii)];
        WCO2_2=[WCO2_2 Licor_w_2(ii)];
    end;
end;
for ii=1:length(Licor_CO2_3);
    if I_3(ii)==0,
        LicCO2_3=[LicCO2_3 Licor_CO2_3(ii)];
        LicH2O_3=[LicH2O_3 Licor_H2O_3(ii)];
        WCO2_3=[WCO2_3 Licor_w_3(ii)];
    end;
end;
for ii=1:length(Licor_CO2_4);
    if I_4(ii)==0,
        LicCO2_4=[LicCO2_4 Licor_CO2_4(ii)];
        LicH2O_4=[LicH2O_4 Licor_H2O_4(ii)];
        WCO2_4=[WCO2_4 Licor_w_4(ii)];
    end;
end;
for ii=1:length(Licor_CO2_5);
    if I_5(ii)==0,
        LicCO2_5=[LicCO2_5 Licor_CO2_5(ii)];
        LicH2O_5=[LicH2O_5 Licor_H2O_5(ii)];
        WCO2_5=[WCO2_5 Licor_w_5(ii)];
    end;
end;
for ii=1:length(Licor_CO2_6);
    if I_6(ii)==0,
        LicCO2_6=[LicCO2_6 Licor_CO2_6(ii)];
        LicH2O_6=[LicH2O_6 Licor_H2O_6(ii)];
        WCO2_6=[WCO2_6 Licor_w_6(ii)];
    end;
end;

%LicTemp_1=[LicTemp_1 Licor_Temp_1(ii)];
%LicPress_1=[LicPress_1 Licor_Press_1(ii)];

if length(LicCO2_1)==0,
    LicCO2_1=[1 1];
    LicH2O_1=[1 1];
    WCO2_1=[1 1];
end;
if length(LicCO2_2)==0,
    LicCO2_2=[1 1];
    LicH2O_2=[1 1];
    WCO2_2=[1 1];
end;
if length(LicCO2_3)==0,
    LicCO2_3=[1 1];
    LicH2O_3=[1 1];
    WCO2_3=[1 1];
end;
if length(LicCO2_4)==0,
    LicCO2_4=[1 1];
    LicH2O_4=[1 1];
    WCO2_4=[1 1];
end;
if length(LicCO2_5)==0,
    LicCO2_5=[1 1];
    LicH2O_5=[1 1];
    WCO2_5=[1 1];
end;
if length(LicCO2_6)==0,
    LicCO2_6=[1 1];
    LicH2O_6=[1 1];
    WCO2_6=[1 1];
end;

wlcj_1=cov(WCO2_1,LicCO2_1);wlcj_1=wlcj_1(1,2);
wlcj_2=cov(WCO2_2,LicCO2_2);wlcj_2=wlcj_2(1,2);
wlcj_3=cov(WCO2_3,LicCO2_3);wlcj_3=wlcj_3(1,2);
wlcj_4=cov(WCO2_4,LicCO2_4);wlcj_4=wlcj_4(1,2);
wlcj_5=cov(WCO2_5,LicCO2_5);wlcj_5=wlcj_5(1,2);
wlcj_6=cov(WCO2_6,LicCO2_6);wlcj_6=wlcj_6(1,2);

wlhj_1=cov(WCO2_1,LicH2O_1);wlhj_1=wlhj_1(1,2);
wlhj_2=cov(WCO2_2,LicH2O_2);wlhj_2=wlhj_2(1,2);
wlhj_3=cov(WCO2_3,LicH2O_3);wlhj_3=wlhj_3(1,2);
wlhj_4=cov(WCO2_4,LicH2O_4);wlhj_4=wlhj_4(1,2);
wlhj_5=cov(WCO2_5,LicH2O_5);wlhj_5=wlhj_5(1,2);
wlhj_6=cov(WCO2_6,LicH2O_6);wlhj_6=wlhj_6(1,2);

disp('Licor H2O +')
for iii=0:20;
   A1=cov(WCO2_1(1:length(WCO2_1)-iii),LicH2O_1(1+iii:length(LicH2O_1)));
   A2=cov(WCO2_2(1:length(WCO2_2)-iii),LicH2O_2(1+iii:length(LicH2O_2)));
   A3=cov(WCO2_3(1:length(WCO2_3)-iii),LicH2O_3(1+iii:length(LicH2O_3)));
   A4=cov(WCO2_4(1:length(WCO2_4)-iii),LicH2O_4(1+iii:length(LicH2O_4)));
   A5=cov(WCO2_5(1:length(WCO2_5)-iii),LicH2O_5(1+iii:length(LicH2O_5)));
   A6=cov(WCO2_6(1:length(WCO2_6)-iii),LicH2O_6(1+iii:length(LicH2O_6)));
   if length(A1)<2,
       A1=[1 1
           1 1];
   end;
   if length(A2)<2,
       A2=[1 1
           1 1];
   end;
   if length(A3)<2,
       A3=[1 1
           1 1];
   end;
   if length(A4)<2,
       A4=[1 1
           1 1];
   end;
   if length(A5)<2,
       A5=[1 1
           1 1];
   end;
   if length(A6)<2,
       A6=[1 1
           1 1];
   end;
%   [wlhj_1 wlhj_2 wlhj_3 wlhj_4 wlhj_5 wlhj_6]
   B1=A1(1,2);
   B2=A2(1,2);
   B3=A3(1,2);
   B4=A4(1,2);
   B5=A5(1,2);
   B6=A6(1,2);
   if B1>wlhj_1,
      wlhj_1=B1;
   end;
   if B2>wlhj_2,
      wlhj_2=B2;
   end;
   if B3>wlhj_3,
      wlhj_3=B3;
   end;
   if B4>wlhj_4,
      wlhj_4=B4;
  end;
   if B5>wlhj_5,
      wlhj_5=B5;
   end;
   if B6>wlhj_6,
      wlhj_6=B6;
   end;
end;

disp('Licor H2O -')
for iii=-20:0;
   A1=cov(WCO2_1(1-iii:length(WCO2_1)),LicH2O_1(1:length(LicH2O_1)+iii));
   A2=cov(WCO2_2(1-iii:length(WCO2_2)),LicH2O_2(1:length(LicH2O_2)+iii));
   A3=cov(WCO2_3(1-iii:length(WCO2_3)),LicH2O_3(1:length(LicH2O_3)+iii));
   A4=cov(WCO2_4(1-iii:length(WCO2_4)),LicH2O_4(1:length(LicH2O_4)+iii));
   A5=cov(WCO2_5(1-iii:length(WCO2_5)),LicH2O_5(1:length(LicH2O_5)+iii));
   A6=cov(WCO2_6(1-iii:length(WCO2_6)),LicH2O_6(1:length(LicH2O_6)+iii));
   if length(A1)<2,
       A1=[1 1
           1 1];
   end;
   if length(A2)<2,
       A2=[1 1
           1 1];
   end;
   if length(A3)<2,
       A3=[1 1
           1 1];
   end;
   if length(A4)<2,
       A4=[1 1
           1 1];
   end;
   if length(A5)<2,
       A5=[1 1
           1 1];
   end;
   if length(A6)<2,
       A6=[1 1
           1 1];
   end;
   B1=A1(1,2);
   B2=A2(1,2);
   B3=A3(1,2);
   B4=A4(1,2);
   B5=A5(1,2);
   B6=A6(1,2);
%   [wlhj_1 wlhj_2 wlhj_3 wlhj_4 wlhj_5 wlhj_6]
   if B1>wlhj_1,
      wlhj_1=B1;
   end;
   if B2>wlhj_2,
      wlhj_2=B2;
   end;
   if B3>wlhj_3,
      wlhj_3=B3;
   end;
   if B4>wlhj_4,
      wlhj_4=B4;
  end;
   if B5>wlhj_5,
      wlhj_5=B5;
   end;
   if B6>wlhj_6,
      wlhj_6=B6;
   end;
end;

wlcj_1=0;
wlcj_2=0;
wlcj_3=0;
wlcj_4=0;
wlcj_5=0;
wlcj_6=0;

disp('Licor CO2 +')
for iii=0:20;
   A1=cov(WCO2_1(1:length(WCO2_1)-iii),LicCO2_1(1+iii:length(LicCO2_1)));
   A2=cov(WCO2_2(1:length(WCO2_2)-iii),LicCO2_2(1+iii:length(LicCO2_2)));
   A3=cov(WCO2_3(1:length(WCO2_3)-iii),LicCO2_3(1+iii:length(LicCO2_3)));
   A4=cov(WCO2_4(1:length(WCO2_4)-iii),LicCO2_4(1+iii:length(LicCO2_4)));
   A5=cov(WCO2_5(1:length(WCO2_5)-iii),LicCO2_5(1+iii:length(LicCO2_5)));
   A6=cov(WCO2_6(1:length(WCO2_6)-iii),LicCO2_6(1+iii:length(LicCO2_6)));
   if length(A1)<2,
       A1=[1 1
           1 1];
   end;
   if length(A2)<2,
       A2=[1 1
           1 1];
   end;
   if length(A3)<2,
       A3=[1 1
           1 1];
   end;
   if length(A4)<2,
       A4=[1 1
           1 1];
   end;
   if length(A5)<2,
       A5=[1 1
           1 1];
   end;
   if length(A6)<2,
       A6=[1 1
           1 1];
   end;
   B1=A1(1,2);
   B2=A2(1,2);
   B3=A3(1,2);
   B4=A4(1,2);
   B5=A5(1,2);
   B6=A6(1,2);
%   [wlcj_1 wlcj_2 wlcj_3 wlcj_4 wlcj_5 wlcj_6]
   if B1>wlcj_1,
      wlcj_1=B1;
   end;
   if B2>wlcj_2,
      wlcj_2=B2;
   end;
   if B3>wlcj_3,
      wlcj_3=B3;
   end;
   if B4>wlcj_4,
      wlcj_4=B4;
   end;
   if B5>wlcj_5,
      wlcj_5=B5;
   end;
   if B6>wlcj_6,
      wlcj_6=B6;
   end;
end;

disp('Licor CO2 -')
for iii=-20:0;
   A1=cov(WCO2_1(1-iii:length(WCO2_1)),LicCO2_1(1:length(LicCO2_1)+iii));
   A2=cov(WCO2_2(1-iii:length(WCO2_2)),LicCO2_2(1:length(LicCO2_2)+iii));
   A3=cov(WCO2_3(1-iii:length(WCO2_3)),LicCO2_3(1:length(LicCO2_3)+iii));
   A4=cov(WCO2_4(1-iii:length(WCO2_4)),LicCO2_4(1:length(LicCO2_4)+iii));
   A5=cov(WCO2_5(1-iii:length(WCO2_5)),LicCO2_5(1:length(LicCO2_5)+iii));
   A6=cov(WCO2_6(1-iii:length(WCO2_6)),LicCO2_6(1:length(LicCO2_6)+iii));
   if length(A1)<2,
       A1=[1 1
           1 1];
   end;
   if length(A2)<2,
       A2=[1 1
           1 1];
   end;
   if length(A3)<2,
       A3=[1 1
           1 1];
   end;
   if length(A4)<2,
       A4=[1 1
           1 1];
   end;
   if length(A5)<2,
       A5=[1 1
           1 1];
   end;
   if length(A6)<2,
       A6=[1 1
           1 1];
   end;
   B1=A1(1,2);
   B2=A2(1,2);
   B3=A3(1,2);
   B4=A4(1,2);
   B5=A5(1,2);
   B6=A6(1,2);
%   [wlcj_1 wlcj_2 wlcj_3 wlcj_4 wlcj_5 wlcj_6]
   if abs(B1)>abs(wlcj_1),
      wlcj_1=B1;
   end;
   if abs(B2)>abs(wlcj_2),
      wlcj_2=B2;
   end;
   if abs(B3)>abs(wlcj_3),
      wlcj_3=B3;
   end;
   if abs(B4)>abs(wlcj_4),
      wlcj_4=B4;
   end;
   if abs(B5)>abs(wlcj_5),
      wlcj_5=B5;
   end;
   if abs(B6)>abs(wlcj_6),
      wlcj_6=B6;
   end;
end;

wohj_1=0;
wohj_2=0;
wohj_3=0;
wohj_4=0;
wohj_5=0;
wohj_6=0;

wlcj_1=wlcj_1.*28.96.*365.24.*24.*60.*60.*1e-6;   %units are mol/m2/year
wlcj_2=wlcj_2.*28.96.*365.24.*24.*60.*60.*1e-6;   %units are mol/m2/year
wlcj_3=wlcj_3.*28.96.*365.24.*24.*60.*60.*1e-6;   %units are mol/m2/year
wlcj_4=wlcj_4.*28.96.*365.24.*24.*60.*60.*1e-6;   %units are mol/m2/year
wlcj_5=wlcj_5.*28.96.*365.24.*24.*60.*60.*1e-6;   %units are mol/m2/year
wlcj_6=wlcj_6.*28.96.*365.24.*24.*60.*60.*1e-6;   %units are mol/m2/year


disp('COVARIANCE')
disp('  licor')

[wlhj_1.*2.5e3.*rho_1
 wlhj_2.*2.5e3.*rho_2
 wlhj_3.*2.5e3.*rho_3
 wlhj_4.*2.5e3.*rho_4
 wlhj_5.*2.5e3.*rho_5
 wlhj_6.*2.5e3.*rho_6]
   
disp('MEAN')
disp(' licor')

[mean(LicH2O_1)
 mean(LicH2O_2)
 mean(LicH2O_3)
 mean(LicH2O_4)
 mean(LicH2O_5)
 mean(LicH2O_6)]
   
%disp('     oak   licor')

%[wocj_1*29*3600*24*365/1e6 wlcj_1*29*3600*24*365/1e6
%   wocj_2*29*3600*24*365/1e6 wlcj_2*29*3600*24*365/1e6
%   wocj_3*29*3600*24*365/1e6 wlcj_3*29*3600*24*365/1e6
%   wocj_4*29*3600*24*365/1e6 wlcj_4*29*3600*24*365/1e6
%   wocj_5*29*3600*24*365/1e6 wlcj_5*29*3600*24*365/1e6
%   wocj_6*29*3600*24*365/1e6 wlcj_6*29*3600*24*365/1e6]

mean_wind_1=nanmean(imum(01:10));std_wind_1=std(imum(01:10));max_wind_1=max(imum(01:10));min_wind_1=min(imum(01:10));
mean_wind_2=nanmean(imum(11:20));std_wind_2=std(imum(11:20));max_wind_2=max(imum(11:20));min_wind_2=min(imum(11:20));
mean_wind_3=nanmean(imum(21:30));std_wind_3=std(imum(21:30));max_wind_3=max(imum(21:30));min_wind_3=min(imum(21:30));
mean_wind_4=nanmean(imum(31:40));std_wind_4=std(imum(31:40));max_wind_4=max(imum(31:40));min_wind_4=min(imum(31:40));
mean_wind_5=nanmean(imum(41:50));std_wind_5=std(imum(41:50));max_wind_5=max(imum(41:50));min_wind_5=min(imum(41:50));
mean_wind_6=nanmean(imum(51:end));std_wind_6=std(imum(51:end));max_wind_6=max(imum(51:end));min_wind_6=min(imum(51:end));

mean_dirC_1=nanmedian(cos(imdm(01:10)./180.*pi));
mean_dirC_2=nanmedian(cos(imdm(11:20)./180.*pi));
mean_dirC_3=nanmedian(cos(imdm(21:30)./180.*pi));
mean_dirC_4=nanmedian(cos(imdm(31:40)./180.*pi));
mean_dirC_5=nanmedian(cos(imdm(41:50)./180.*pi));
mean_dirC_6=median(cos(imdm(51:end)./180.*pi));

mean_dirS_1=nanmedian(sin(imdm(01:10)./180.*pi));
mean_dirS_2=nanmedian(sin(imdm(11:20)./180.*pi));
mean_dirS_3=nanmedian(sin(imdm(21:30)./180.*pi));
mean_dirS_4=nanmedian(sin(imdm(31:40)./180.*pi));
mean_dirS_5=nanmedian(sin(imdm(41:50)./180.*pi));
mean_dirS_6=nanmedian(sin(imdm(51:end)./180.*pi));

mean_dir_1=mod(atan2(mean_dirS_1,mean_dirC_1).*180./pi+360,360);		%
mean_dir_2=mod(atan2(mean_dirS_2,mean_dirC_2).*180./pi+360,360);		%
mean_dir_3=mod(atan2(mean_dirS_3,mean_dirC_3).*180./pi+360,360);		%
mean_dir_4=mod(atan2(mean_dirS_4,mean_dirC_4).*180./pi+360,360);		%
mean_dir_5=mod(atan2(mean_dirS_5,mean_dirC_5).*180./pi+360,360);		%reported IMET true wind direction
mean_dir_6=mod(atan2(mean_dirS_6,mean_dirC_6).*180./pi+360,360);		%

clear mean_dirC_1 mean_dirC_2 mean_dirC_3 mean_dirC_4 mean_dirC_5 mean_dirC_6 
clear mean_dirS_1 mean_dirS_2 mean_dirS_3 mean_dirS_4 mean_dirS_5 mean_dirS_6 

std_dir_1=std(unwrap(imdm(01:10)));max_dir_1=max(unwrap(imdm(01:10)));min_dir_1=min(unwrap(imdm(01:10)));
std_dir_2=std(unwrap(imdm(11:20)));max_dir_2=max(unwrap(imdm(11:20)));min_dir_2=min(unwrap(imdm(11:20)));
std_dir_3=std(unwrap(imdm(21:30)));max_dir_3=max(unwrap(imdm(21:30)));min_dir_3=min(unwrap(imdm(21:30)));
std_dir_4=std(unwrap(imdm(31:40)));max_dir_4=max(unwrap(imdm(31:40)));min_dir_4=min(unwrap(imdm(31:40)));
std_dir_5=std(unwrap(imdm(41:50)));max_dir_5=max(unwrap(imdm(41:50)));min_dir_5=min(unwrap(imdm(41:50)));
std_dir_6=std(unwrap(imdm(51:end)));max_dir_6=max(unwrap(imdm(51:end)));min_dir_6=min(unwrap(imdm(51:end)));

Reldir_1=atan2(nanmean(V(00001:Nlendata)),-nanmean(U(00001:Nlendata)))*180/pi;
Reldir_2=atan2(nanmean(V(Nlendata+1:Nlendata*2)),-nanmean(U(Nlendata+1:Nlendata*2)))*180/pi;
Reldir_3=atan2(nanmean(V(Nlendata*2+1:Nlendata*3)),-nanmean(U(Nlendata*2+1:Nlendata*3)))*180/pi;
Reldir_4=atan2(nanmean(V(Nlendata*3+1:Nlendata*4)),-nanmean(U(Nlendata*3+1:Nlendata*4)))*180/pi;
Reldir_5=atan2(nanmean(V(Nlendata*4+1:Nlendata*5)),-nanmean(U(Nlendata*4+1:Nlendata*5)))*180/pi;
Reldir_6=atan2(nanmean(V(Nlendata*5+1:Nlendata*6)),-nanmean(U(Nlendata*5+1:Nlendata*6)))*180/pi;


%     save motion data for other flux computation
Reldirby=[atan2((V(00001:Nlendata)),-(U(00001:Nlendata)))*180/pi, atan2((V(Nlendata+1:Nlendata*2)),-(U(Nlendata+1:Nlendata*2)))*180/pi, atan2((V(Nlendata*2+1:Nlendata*3)),-(U(Nlendata*2+1:Nlendata*3)))*180/pi, atan2((V(Nlendata*3+1:Nlendata*4)),-(U(Nlendata*3+1:Nlendata*4)))*180/pi, atan2((V(Nlendata*4+1:Nlendata*5)),-(U(Nlendata*4+1:Nlendata*5)))*180/pi, atan2((V(Nlendata*5+1:Nlendata*6)),-(U(Nlendata*5+1:Nlendata*6)))*180/pi];
ii=find(Reldirby<0);Reldirby(ii)=Reldirby(ii)+360;
Uj=[Uj_1, Uj_2, Uj_3, Uj_4, Uj_5, Uj_6];
Vj=[Vj_1, Vj_2, Vj_3, Vj_4, Vj_5, Vj_6];
Wj=[Wj_1, Wj_2, Wj_3, Wj_4, Wj_5, Wj_6];
try
    uvwmat=[tmotion' Uj' Vj' Wj' Reldirby'];
catch
    tmotion=tmotion';
    uvwmat=[tmotion' Uj' Vj' Wj' Reldirby'];
end
sdg=['save ', drive_w '\MotionData\HiWinGS2013_UVWstream_',jd hr, dcname '.prn uvwmat -ascii'];
eval(sdg)


% % % % % % % %save various data for big spectrum project (Simon/Jeff request 031309)
% %     Reldirby=atan2(V,-U)*180/pi;    ii=find(Reldirby<0);Reldirby(ii)=Reldirby(ii)+360;
% %     [Uj,Vj,Wj,alphaj,betaj]=rotcorr(uvw(1,:),uvw(2,:),uvw(3,:));
% %     MCO2=[Licor_CO2_1 Licor_CO2_2 Licor_CO2_3 Licor_CO2_4 Licor_CO2_5 Licor_CO2_6]; %umol/mol
% %     MH2O=[Licor_H2O_1 Licor_H2O_2 Licor_H2O_3 Licor_H2O_4 Licor_H2O_5 Licor_H2O_6]; %g/kg
% %     
% %     specmatall=[tref' Uj' Vj' Wj' Tsonic' MH2O' MCO2' Reldirby' heading'];
% %     sdg=['save ', 'C:\Data\CRUISES\VOCALS\data\CorrectedData4spec\specmat_',jd hr,'.prn specmatall -ascii'];
% %     eval(sdg)
% %     
% %     specmatall2=[press' org'];
% %     sdg=['save ', 'C:\Data\CRUISES\VOCALS\data\CorrectedData4spec\filtermat_',jd hr,'.prn specmatall2 -ascii'];
% %     eval(sdg)
% %     
% %     odecmat=odec';
% %     sdg=['save ', 'C:\Data\CRUISES\VOCALS\data\CorrectedData4spec\odecmat_',jd hr,'.prn odecmat -ascii'];
% %     eval(sdg)
% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% %
% % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % %


U=U./0.95;              %Flow Distortion correction
V=V./1.15;

Un_1=nanmean(U(00001:Nlendata).*cos(heading(00001:Nlendata))+V(00001:Nlendata).*sin(heading(00001:Nlendata)))+nanmean(sog(00001:Nlendata)./1.94.*cos(cog(00001:Nlendata)./180.*pi));
Un_2=nanmean(U(Nlendata+1:Nlendata*2).*cos(heading(Nlendata+1:Nlendata*2))+V(Nlendata+1:Nlendata*2).*sin(heading(Nlendata+1:Nlendata*2)))+nanmean(sog(Nlendata+1:Nlendata*2)./1.94.*cos(cog(Nlendata+1:Nlendata*2)./180.*pi));
Un_3=nanmean(U(Nlendata*2+1:Nlendata*3).*cos(heading(Nlendata*2+1:Nlendata*3))+V(Nlendata*2+1:Nlendata*3).*sin(heading(Nlendata*2+1:Nlendata*3)))+nanmean(sog(Nlendata*2+1:Nlendata*3)./1.94.*cos(cog(Nlendata*2+1:Nlendata*3)./180.*pi));
Un_4=nanmean(U(Nlendata*3+1:Nlendata*4).*cos(heading(Nlendata*3+1:Nlendata*4))+V(Nlendata*3+1:Nlendata*4).*sin(heading(Nlendata*3+1:Nlendata*4)))+nanmean(sog(Nlendata*3+1:Nlendata*4)./1.94.*cos(cog(Nlendata*3+1:Nlendata*4)./180.*pi));
Un_5=nanmean(U(Nlendata*4+1:Nlendata*5).*cos(heading(Nlendata*4+1:Nlendata*5))+V(Nlendata*4+1:Nlendata*5).*sin(heading(Nlendata*4+1:Nlendata*5)))+nanmean(sog(Nlendata*4+1:Nlendata*5)./1.94.*cos(cog(Nlendata*4+1:Nlendata*5)./180.*pi));
Un_6=nanmean(U(Nlendata*5+1:Nlendata*6).*cos(heading(Nlendata*5+1:Nlendata*6))+V(Nlendata*5+1:Nlendata*6).*sin(heading(Nlendata*5+1:Nlendata*6)))+nanmean(sog(Nlendata*5+1:Nlendata*6)./1.94.*cos(cog(Nlendata*5+1:Nlendata*6)./180.*pi));

Ue_1=nanmean(U(00001:Nlendata).*sin(heading(00001:Nlendata))-V(00001:Nlendata).*cos(heading(00001:Nlendata)))+nanmean(sog(00001:Nlendata)./1.94.*sin(cog(00001:Nlendata)./180.*pi));
Ue_2=nanmean(U(Nlendata+1:Nlendata*2).*sin(heading(Nlendata+1:Nlendata*2))-V(Nlendata+1:Nlendata*2).*cos(heading(Nlendata+1:Nlendata*2)))+nanmean(sog(Nlendata+1:Nlendata*2)./1.94.*sin(cog(Nlendata+1:Nlendata*2)./180.*pi));
Ue_3=nanmean(U(Nlendata*2+1:Nlendata*3).*sin(heading(Nlendata*2+1:Nlendata*3))-V(Nlendata*2+1:Nlendata*3).*cos(heading(Nlendata*2+1:Nlendata*3)))+nanmean(sog(Nlendata*2+1:Nlendata*3)./1.94.*sin(cog(Nlendata*2+1:Nlendata*3)./180.*pi));
Ue_4=nanmean(U(Nlendata*3+1:Nlendata*4).*sin(heading(Nlendata*3+1:Nlendata*4))-V(Nlendata*3+1:Nlendata*4).*cos(heading(Nlendata*3+1:Nlendata*4)))+nanmean(sog(Nlendata*3+1:Nlendata*4)./1.94.*sin(cog(Nlendata*3+1:Nlendata*4)./180.*pi));
Ue_5=nanmean(U(Nlendata*4+1:Nlendata*5).*sin(heading(Nlendata*4+1:Nlendata*5))-V(Nlendata*4+1:Nlendata*5).*cos(heading(Nlendata*4+1:Nlendata*5)))+nanmean(sog(Nlendata*4+1:Nlendata*5)./1.94.*sin(cog(Nlendata*4+1:Nlendata*5)./180.*pi));
Ue_6=nanmean(U(Nlendata*5+1:Nlendata*6).*sin(heading(Nlendata*5+1:Nlendata*6))-V(Nlendata*5+1:Nlendata*6).*cos(heading(Nlendata*5+1:Nlendata*6)))+nanmean(sog(Nlendata*5+1:Nlendata*6)./1.94.*sin(cog(Nlendata*5+1:Nlendata*6)./180.*pi));

%  Do that in DA codes...
% switch leg(1:end-1)
%     case {'leg2','leg3','leg4'}
%         %uses WRT water here for legs2and3
%     Un_1=Un_1+CN(1);
%     Un_2=Un_2+CN(2);
%     Un_3=Un_3+CN(3);
%     Un_4=Un_4+CN(4);
%     Un_5=Un_5+CN(5);
%     Un_6=Un_6+CN(6);
% 
%     Ue_1=Ue_1+CE(1);
%     Ue_2=Ue_2+CE(2);
%     Ue_3=Ue_3+CE(3);
%     Ue_4=Ue_4+CE(4);
%     Ue_5=Ue_5+CE(5);
%     Ue_6=Ue_6+CE(6);
%    
% end       

Ubar_1=sqrt(Un_1.^2 + Ue_1.^2);
Ubar_2=sqrt(Un_2.^2 + Ue_2.^2);
Ubar_3=sqrt(Un_3.^2 + Ue_3.^2);
Ubar_4=sqrt(Un_4.^2 + Ue_4.^2);
Ubar_5=sqrt(Un_5.^2 + Ue_5.^2);
Ubar_6=sqrt(Un_6.^2 + Ue_6.^2);

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


wocj_1=cov(Licor_w_1,Licor_CO2_1);   % m/s * mic atm
wocj_2=cov(Licor_w_2,Licor_CO2_2);   % m/s * mic atm
wocj_3=cov(Licor_w_3,Licor_CO2_3);   % m/s * mic atm
wocj_4=cov(Licor_w_4,Licor_CO2_4);   % m/s * mic atm
wocj_5=cov(Licor_w_5,Licor_CO2_5);   % m/s * mic atm
wocj_6=cov(Licor_w_6,Licor_CO2_6);   % m/s * mic atm

wocj_1=wocj_1(1,2);
wocj_2=wocj_2(1,2);
wocj_3=wocj_3(1,2);
wocj_4=wocj_4(1,2);
wocj_5=wocj_5(1,2);
wocj_6=wocj_6(1,2);

ship_RH_1=nanmean(imrh(1:10));
ship_RH_2=nanmean(imrh(11:20));
ship_RH_3=nanmean(imrh(21:30));
ship_RH_4=nanmean(imrh(31:40));
ship_RH_5=nanmean(imrh(41:50));
ship_RH_6=nanmean(imrh(51:end));

RHv_1=nanmean(Rhvais(1:10));
RHv_2=nanmean(Rhvais(11:20));
RHv_3=nanmean(Rhvais(21:30));
RHv_4=nanmean(Rhvais(31:40));
RHv_5=nanmean(Rhvais(41:50));
RHv_6=nanmean(Rhvais(51:end));

Qship_1=qair_p([Tship_1 ship_RH_1 press_1]);
Qship_2=qair_p([Tship_2 ship_RH_2 press_2]);
Qship_3=qair_p([Tship_3 ship_RH_3 press_3]);
Qship_4=qair_p([Tship_4 ship_RH_4 press_4]);
Qship_5=qair_p([Tship_5 ship_RH_5 press_5]);
Qship_6=qair_p([Tship_6 ship_RH_6 press_6]);

Qbarv_1=qair_p([Tbar_1 RHv_1 press_1]);
Qbarv_2=qair_p([Tbar_2 RHv_2 press_2]);
Qbarv_3=qair_p([Tbar_3 RHv_3 press_3]);
Qbarv_4=qair_p([Tbar_4 RHv_4 press_4]);
Qbarv_5=qair_p([Tbar_5 RHv_5 press_5]);
Qbarv_6=qair_p([Tbar_6 RHv_6 press_6]);

Tseasnake_1=nanmean(Tsea(1:10));
Tseasnake_2=nanmean(Tsea(11:20));
Tseasnake_3=nanmean(Tsea(21:30));
Tseasnake_4=nanmean(Tsea(31:40));
Tseasnake_5=nanmean(Tsea(41:50));
Tseasnake_6=nanmean(Tsea(51:end));

Qseasnake_1=qsea_p([Tseasnake_1 press_1]);
Qseasnake_2=qsea_p([Tseasnake_2 press_2]);
Qseasnake_3=qsea_p([Tseasnake_3 press_3]);
Qseasnake_4=qsea_p([Tseasnake_4 press_4]);
Qseasnake_5=qsea_p([Tseasnake_5 press_5]);
Qseasnake_6=qsea_p([Tseasnake_6 press_6]);

Tsea_1=nanmean(tsg(1:10));
Tsea_2=nanmean(tsg(11:20));
Tsea_3=nanmean(tsg(21:30));
Tsea_4=nanmean(tsg(31:40));
Tsea_5=nanmean(tsg(41:50));
Tsea_6=nanmean(tsg(51:end));

Qsea_1=qsea_p([Tsea_1 press_1]);
Qsea_2=qsea_p([Tsea_2 press_2]);
Qsea_3=qsea_p([Tsea_3 press_3]);
Qsea_4=qsea_p([Tsea_4 press_4]);
Qsea_5=qsea_p([Tsea_5 press_5]);
Qsea_6=qsea_p([Tsea_6 press_6]);

Ssin_1=std(shead(1:Nlendata));Scos_1=std(chead(1:Nlendata));
Ssin_2=std(shead(Nlendata+1:Nlendata*2));Scos_2=std(chead(Nlendata+1:Nlendata*2));
Ssin_3=std(shead(Nlendata*2+1:Nlendata*3));Scos_3=std(chead(Nlendata*2+1:Nlendata*3));
Ssin_4=std(shead(Nlendata*3+1:Nlendata*4));Scos_4=std(chead(Nlendata*3+1:Nlendata*4));
Ssin_5=std(shead(Nlendata*4+1:Nlendata*5));Scos_5=std(chead(Nlendata*4+1:Nlendata*5));
Ssin_6=std(shead(Nlendata*5+1:Nlendata*6));Scos_6=std(chead(Nlendata*5+1:Nlendata*6));

% Sp1_1=std(Platvel(1,1:Nlendata));Sp2_1=std(Platvel(2,1:Nlendata));Sp3_1=std(Platvel(3,1:Nlendata));
% Sp1_2=std(Platvel(1,Nlendata+1:Nlendata*2));Sp2_2=std(Platvel(2,Nlendata+1:Nlendata*2));Sp3_2=std(Platvel(3,Nlendata+1:Nlendata*2));
% Sp1_3=std(Platvel(1,Nlendata*2+1:Nlendata*3));Sp2_3=std(Platvel(2,Nlendata*2+1:Nlendata*3));Sp3_3=std(Platvel(3,Nlendata*2+1:Nlendata*3));
% Sp1_4=std(Platvel(1,Nlendata*3+1:Nlendata*4));Sp2_4=std(Platvel(2,Nlendata*3+1:Nlendata*4));Sp3_4=std(Platvel(3,Nlendata*3+1:Nlendata*4));
% Sp1_5=std(Platvel(1,Nlendata*4+1:Nlendata*5));Sp2_5=std(Platvel(2,Nlendata*4+1:Nlendata*5));Sp3_5=std(Platvel(3,Nlendata*4+1:Nlendata*5));
% Sp1_6=std(Platvel(1,Nlendata*5+1:Nlendata*6));Sp2_6=std(Platvel(2,Nlendata*5+1:Nlendata*6));Sp3_6=std(Platvel(3,Nlendata*5+1:Nlendata*6));
% 
Sp1_1=std(uvwplat(1,1:Nlendata));Sp2_1=std(uvwplat(2,1:Nlendata));Sp3_1=std(uvwplat(3,1:Nlendata));
Sp1_2=std(uvwplat(1,Nlendata+1:Nlendata*2));Sp2_2=std(uvwplat(2,Nlendata+1:Nlendata*2));Sp3_2=std(uvwplat(3,Nlendata+1:Nlendata*2));
Sp1_3=std(uvwplat(1,Nlendata*2+1:Nlendata*3));Sp2_3=std(uvwplat(2,Nlendata*2+1:Nlendata*3));Sp3_3=std(uvwplat(3,Nlendata*2+1:Nlendata*3));
Sp1_4=std(uvwplat(1,Nlendata*3+1:Nlendata*4));Sp2_4=std(uvwplat(2,Nlendata*3+1:Nlendata*4));Sp3_4=std(uvwplat(3,Nlendata*3+1:Nlendata*4));
Sp1_5=std(uvwplat(1,Nlendata*4+1:Nlendata*5));Sp2_5=std(uvwplat(2,Nlendata*4+1:Nlendata*5));Sp3_5=std(uvwplat(3,Nlendata*4+1:Nlendata*5));
Sp1_6=std(uvwplat(1,Nlendata*5+1:Nlendata*6));Sp2_6=std(uvwplat(2,Nlendata*5+1:Nlendata*6));Sp3_6=std(uvwplat(3,Nlendata*5+1:Nlendata*6));

rain_1=nanmean(org(1:10));
rain_2=nanmean(org(11:20));
rain_3=nanmean(org(21:30));
rain_4=nanmean(org(31:40));
rain_5=nanmean(org(41:50));
rain_6=nanmean(org(51:end));

ship_rain_1=nanmean(imrain(1:10));
ship_rain_2=nanmean(imrain(11:20));
ship_rain_3=nanmean(imrain(21:30));
ship_rain_4=nanmean(imrain(31:40));
ship_rain_5=nanmean(imrain(41:50));
ship_rain_6=nanmean(imrain(51:end));

Rl_1=nanmean(pir(1:10));Rs_1=nanmean(psp(1:10));
Rl_2=nanmean(pir(11:20));Rs_2=nanmean(psp(11:20));
Rl_3=nanmean(pir(21:30));Rs_3=nanmean(psp(21:30));
Rl_4=nanmean(pir(31:40));Rs_4=nanmean(psp(31:40));
Rl_5=nanmean(pir(41:50));Rs_5=nanmean(psp(41:50));
Rl_6=nanmean(pir(51:end));Rs_6=nanmean(psp(51:end));

std_Rl_1=std(pir(1:10));std_Rs_1=std(psp(1:10));
std_Rl_2=std(pir(11:20));std_Rs_2=std(psp(11:20));
std_Rl_3=std(pir(21:30));std_Rs_3=std(psp(21:30));
std_Rl_4=std(pir(31:40));std_Rs_4=std(psp(31:40));
std_Rl_5=std(pir(41:50));std_Rs_5=std(psp(41:50));
std_Rl_6=std(pir(51:end));std_Rs_6=std(psp(51:end));

Rs_ship_1=nanmean(imsol(1:10));
Rs_ship_2=nanmean(imsol(11:20));
Rs_ship_3=nanmean(imsol(21:30));
Rs_ship_4=nanmean(imsol(31:40));
Rs_ship_5=nanmean(imsol(41:50));
Rs_ship_6=nanmean(imsol(51:end));

Rl_ship_1=nanmean(imir(1:10));
Rl_ship_2=nanmean(imir(11:20));
Rl_ship_3=nanmean(imir(21:30));
Rl_ship_4=nanmean(imir(31:40));
Rl_ship_5=nanmean(imir(41:50));
Rl_ship_6=nanmean(imir(51:end));


mean_gpscog_1=mod(360+180.*atan2(nanmean(sin(gpshead(0001:0600)./180.*pi)),nanmean(cos(gpshead(0001:0600)./180.*pi)))./pi,360);
mean_gpscog_2=mod(360+180.*atan2(nanmean(sin(gpshead(0601:1200)./180.*pi)),nanmean(cos(gpshead(0601:1200)./180.*pi)))./pi,360);
mean_gpscog_3=mod(360+180.*atan2(nanmean(sin(gpshead(1201:1800)./180.*pi)),nanmean(cos(gpshead(1201:1800)./180.*pi)))./pi,360);
mean_gpscog_4=mod(360+180.*atan2(nanmean(sin(gpshead(1801:2400)./180.*pi)),nanmean(cos(gpshead(1801:2400)./180.*pi)))./pi,360);
mean_gpscog_5=mod(360+180.*atan2(nanmean(sin(gpshead(2401:3000)./180.*pi)),nanmean(cos(gpshead(2401:3000)./180.*pi)))./pi,360);
mean_gpscog_6=mod(360+180.*atan2(nanmean(sin(gpshead(3001:3600)./180.*pi)),nanmean(cos(gpshead(3001:3600)./180.*pi)))./pi,360);

std_gpscog_1=std(unwrap(gpshead(0001:0600)./180.*pi)*180/pi);
std_gpscog_2=std(unwrap(gpshead(0601:1200)./180.*pi)*180/pi);
std_gpscog_3=std(unwrap(gpshead(1201:1800)./180.*pi)*180/pi);
std_gpscog_4=std(unwrap(gpshead(1801:2400)./180.*pi)*180/pi);
std_gpscog_5=std(unwrap(gpshead(2401:3000)./180.*pi)*180/pi);
std_gpscog_6=std(unwrap(gpshead(3001:3600)./180.*pi)*180/pi);

max_gpscog_1=max(gpshead(0001:0600));
max_gpscog_2=max(gpshead(0601:1200));
max_gpscog_3=max(gpshead(1201:1800));
max_gpscog_4=max(gpshead(1801:2400));
max_gpscog_5=max(gpshead(2401:3000));
max_gpscog_6=max(gpshead(3001:3600));

min_gpscog_1=min(gpshead(0001:0600));
min_gpscog_2=min(gpshead(0601:1200));
min_gpscog_3=min(gpshead(1201:1800));
min_gpscog_4=min(gpshead(1801:2400));
min_gpscog_5=min(gpshead(2401:3000));
min_gpscog_6=min(gpshead(3001:3600));

gpsspeed=gpsspeed./1.94;

mean_gpsspeed_1=nanmean(gpsspeed(0001:0600));min_gpsspeed_1=min(gpsspeed(0001:0600));max_gpsspeed_1=max(gpsspeed(0001:0600));std_gpsspeed_1=std(gpsspeed(0001:0600));
mean_gpsspeed_2=nanmean(gpsspeed(0601:1200));min_gpsspeed_2=min(gpsspeed(0601:1200));max_gpsspeed_2=max(gpsspeed(0601:1200));std_gpsspeed_2=std(gpsspeed(0601:1200));
mean_gpsspeed_3=nanmean(gpsspeed(1201:1800));min_gpsspeed_3=min(gpsspeed(1201:1800));max_gpsspeed_3=max(gpsspeed(1201:1800));std_gpsspeed_3=std(gpsspeed(1201:1800));
mean_gpsspeed_4=nanmean(gpsspeed(1801:2400));min_gpsspeed_4=min(gpsspeed(1801:2400));max_gpsspeed_4=max(gpsspeed(1801:2400));std_gpsspeed_4=std(gpsspeed(1801:2400));
mean_gpsspeed_5=nanmean(gpsspeed(2401:3000));min_gpsspeed_5=min(gpsspeed(2401:3000));max_gpsspeed_5=max(gpsspeed(2401:3000));std_gpsspeed_5=std(gpsspeed(2401:3000));
mean_gpsspeed_6=nanmean(gpsspeed(3001:3600));min_gpsspeed_6=min(gpsspeed(3001:3600));max_gpsspeed_6=max(gpsspeed(3001:3600));std_gpsspeed_6=std(gpsspeed(3001:3600));

Lat_1=gpslat(01);Lon_1=gpslon(01);
Lat_2=gpslat(601);Lon_2=gpslon(601);
Lat_3=gpslat(1201);Lon_3=gpslon(1201);
Lat_4=gpslat(1801);Lon_4=gpslon(1801);
Lat_5=gpslat(2401);Lon_5=gpslon(2401);
Lat_6=gpslat(3001);Lon_6=gpslon(3001);


mean_gyro_1=mod(360+180.*atan2(nanmean(sin(heading(00001:Nlendata))),nanmean(cos(heading(00001:Nlendata))))./pi,360);
mean_gyro_2=mod(360+180.*atan2(nanmean(sin(heading(Nlendata+1:Nlendata*2))),nanmean(cos(heading(Nlendata+1:Nlendata*2))))./pi,360);
mean_gyro_3=mod(360+180.*atan2(nanmean(sin(heading(Nlendata*2+1:Nlendata*3))),nanmean(cos(heading(Nlendata*2+1:Nlendata*3))))./pi,360);
mean_gyro_4=mod(360+180.*atan2(nanmean(sin(heading(Nlendata*3+1:Nlendata*4))),nanmean(cos(heading(Nlendata*3+1:Nlendata*4))))./pi,360);
mean_gyro_5=mod(360+180.*atan2(nanmean(sin(heading(Nlendata*4+1:Nlendata*5))),nanmean(cos(heading(Nlendata*4+1:Nlendata*5))))./pi,360);
mean_gyro_6=mod(360+180.*atan2(nanmean(sin(heading(Nlendata*5+1:Nlendata*6))),nanmean(cos(heading(Nlendata*5+1:Nlendata*6))))./pi,360);

std_gyro_1=std(unwrap(heading(00001:Nlendata)).*180./pi);
std_gyro_2=std(unwrap(heading(Nlendata+1:Nlendata*2)).*180./pi);
std_gyro_3=std(unwrap(heading(Nlendata*2+1:Nlendata*3)).*180./pi);
std_gyro_4=std(unwrap(heading(Nlendata*3+1:Nlendata*4)).*180./pi);
std_gyro_5=std(unwrap(heading(Nlendata*4+1:Nlendata*5)).*180./pi);
std_gyro_6=std(unwrap(heading(Nlendata*5+1:Nlendata*6)).*180./pi);

min_gyro_1=min(heading(00001:Nlendata)).*180./pi;
min_gyro_2=min(heading(Nlendata+1:Nlendata*2)).*180./pi;
min_gyro_3=min(heading(Nlendata*2+1:Nlendata*3)).*180./pi;
min_gyro_4=min(heading(Nlendata*3+1:Nlendata*4)).*180./pi;
min_gyro_5=min(heading(Nlendata*4+1:Nlendata*5)).*180./pi;
min_gyro_6=min(heading(Nlendata*5+1:Nlendata*6)).*180./pi;

max_gyro_1=max(heading(00001:Nlendata)).*180./pi;
max_gyro_2=max(heading(Nlendata+1:Nlendata*2)).*180./pi;
max_gyro_3=max(heading(Nlendata*2+1:Nlendata*3)).*180./pi;
max_gyro_4=max(heading(Nlendata*3+1:Nlendata*4)).*180./pi;
max_gyro_5=max(heading(Nlendata*4+1:Nlendata*5)).*180./pi;
max_gyro_6=max(heading(Nlendata*5+1:Nlendata*6)).*180./pi;

mean_Tsonic_1=nanmean(Tsonic_1);
mean_Tsonic_2=nanmean(Tsonic_2);
mean_Tsonic_3=nanmean(Tsonic_3);
mean_Tsonic_4=nanmean(Tsonic_4);
mean_Tsonic_5=nanmean(Tsonic_5);
mean_Tsonic_6=nanmean(Tsonic_6);

xship_1=[mean_wind_1 Us Tsea_1 Tship_1 Qsea_1 Qship_1 Rs_ship_1 Rl_ship_1 ship_rain_1 zi P zu_ship zt_ship zq_ship];
xship_2=[mean_wind_2 Us Tsea_2 Tship_2 Qsea_2 Qship_2 Rs_ship_2 Rl_ship_2 ship_rain_2 zi P zu_ship zt_ship zq_ship];
xship_3=[mean_wind_3 Us Tsea_3 Tship_3 Qsea_3 Qship_3 Rs_ship_3 Rl_ship_3 ship_rain_3 zi P zu_ship zt_ship zq_ship];
xship_4=[mean_wind_4 Us Tsea_4 Tship_4 Qsea_4 Qship_4 Rs_ship_4 Rl_ship_4 ship_rain_4 zi P zu_ship zt_ship zq_ship];
xship_5=[mean_wind_5 Us Tsea_5 Tship_5 Qsea_5 Qship_5 Rs_ship_5 Rl_ship_5 ship_rain_5 zi P zu_ship zt_ship zq_ship];
xship_6=[mean_wind_6 Us Tsea_6 Tship_6 Qsea_6 Qship_6 Rs_ship_6 Rl_ship_6 ship_rain_6 zi P zu_ship zt_ship zq_ship];

yship_1=coul_lkb(xship_1);
yship_2=coul_lkb(xship_2);
yship_3=coul_lkb(xship_3);
yship_4=coul_lkb(xship_4);
yship_5=coul_lkb(xship_5);
yship_6=coul_lkb(xship_6);

xsnake_1=[Ubar_1 Us Tseasnake_1 Tbar_1 Qseasnake_1 Qbarv_1 Rs_1 Rl_1 rain_1 zi P zu-0.8 zt zq];
xsnake_2=[Ubar_2 Us Tseasnake_2 Tbar_2 Qseasnake_2 Qbarv_2 Rs_2 Rl_2 rain_2 zi P zu-0.8 zt zq];
xsnake_3=[Ubar_3 Us Tseasnake_3 Tbar_3 Qseasnake_3 Qbarv_3 Rs_3 Rl_3 rain_3 zi P zu-0.8 zt zq];
xsnake_4=[Ubar_4 Us Tseasnake_4 Tbar_4 Qseasnake_4 Qbarv_4 Rs_4 Rl_4 rain_4 zi P zu-0.8 zt zq];
xsnake_5=[Ubar_5 Us Tseasnake_5 Tbar_5 Qseasnake_5 Qbarv_5 Rs_5 Rl_5 rain_5 zi P zu-0.8 zt zq];
xsnake_6=[Ubar_6 Us Tseasnake_6 Tbar_6 Qseasnake_6 Qbarv_6 Rs_6 Rl_6 rain_6 zi P zu-0.8 zt zq];

ysnake_1=coul_lkb(xsnake_1);
ysnake_2=coul_lkb(xsnake_2);
ysnake_3=coul_lkb(xsnake_3);
ysnake_4=coul_lkb(xsnake_4);
ysnake_5=coul_lkb(xsnake_5);
ysnake_6=coul_lkb(xsnake_6);


% % % % % % % % % %save data for Alan  %%%%%%%%%%%%%%%%
% % if ~isempty(Uj_1(~isnan(Uj_1))) && ~isempty(Wj_1(~isnan(Wj_1)))
% %     Ujskew_1=skewness(Uj_1(~isnan(Uj_1)));Vjskew_1=skewness(Vj_1(~isnan(Vj_1)));Wjskew_1=skewness(Wj_1(~isnan(Wj_1)));
% %     Ujskew_2=skewness(Uj_2(~isnan(Uj_2)));Vjskew_2=skewness(Vj_2(~isnan(Vj_2)));Wjskew_2=skewness(Wj_2(~isnan(Wj_2)));
% %     Ujskew_3=skewness(Uj_3(~isnan(Uj_3)));Vjskew_3=skewness(Vj_3(~isnan(Vj_3)));Wjskew_3=skewness(Wj_3(~isnan(Wj_3)));
% %     Ujskew_4=skewness(Uj_4(~isnan(Uj_4)));Vjskew_4=skewness(Vj_4(~isnan(Vj_4)));Wjskew_4=skewness(Wj_4(~isnan(Wj_4)));
% %     Ujskew_5=skewness(Uj_5(~isnan(Uj_5)));Vjskew_5=skewness(Vj_5(~isnan(Vj_5)));Wjskew_5=skewness(Wj_5(~isnan(Wj_5)));
% %     Ujskew_6=skewness(Uj_6(~isnan(Uj_6)));Vjskew_6=skewness(Vj_6(~isnan(Vj_6)));Wjskew_6=skewness(Wj_6(~isnan(Wj_6))); 
% % else
% %     Ujskew_1=NaN;Vjskew_1=NaN;Wjskew_1=NaN;
% %     Ujskew_2=NaN;Vjskew_2=NaN;Wjskew_2=NaN;
% %     Ujskew_3=NaN;Vjskew_3=NaN;Wjskew_3=NaN;
% %     Ujskew_4=NaN;Vjskew_4=NaN;Wjskew_4=NaN;
% %     Ujskew_5=NaN;Vjskew_5=NaN;Wjskew_5=NaN;
% %     Ujskew_6=NaN;Vjskew_6=NaN;Wjskew_6=NaN; 
% % end
% % 
% %     Alanview=[ddd+hhh/24+00/24/60 nanmean(Uj_1) nanmean(Vj_1) nanmean(Wj_1) nanvar(Uj_1) nanvar(Vj_1) nanvar(Wj_1) alphaj_1 betaj_1 wuj_1(1,2) wvj_1(1,2) wtj_1(1,2) Ujskew_1 Vjskew_1 Wjskew_1
% %           ddd+hhh/24+10/24/60 nanmean(Uj_2) nanmean(Vj_2) nanmean(Wj_2) nanvar(Uj_2) nanvar(Vj_2) nanvar(Wj_2) alphaj_2 betaj_2 wuj_2(1,2) wvj_2(1,2) wtj_2(1,2) Ujskew_2 Vjskew_2 Wjskew_2
% %           ddd+hhh/24+20/24/60 nanmean(Uj_3) nanmean(Vj_3) nanmean(Wj_3) nanvar(Uj_3) nanvar(Vj_3) nanvar(Wj_3) alphaj_3 betaj_3 wuj_3(1,2) wvj_3(1,2) wtj_3(1,2) Ujskew_3 Vjskew_3 Wjskew_3
% %           ddd+hhh/24+30/24/60 nanmean(Uj_4) nanmean(Vj_4) nanmean(Wj_4) nanvar(Uj_4) nanvar(Vj_4) nanvar(Wj_4) alphaj_4 betaj_4 wuj_4(1,2) wvj_4(1,2) wtj_4(1,2) Ujskew_4 Vjskew_4 Wjskew_4
% %           ddd+hhh/24+40/24/60 nanmean(Uj_5) nanmean(Vj_5) nanmean(Wj_5) nanvar(Uj_5) nanvar(Vj_5) nanvar(Wj_5) alphaj_5 betaj_5 wuj_5(1,2) wvj_5(1,2) wtj_5(1,2) Ujskew_5 Vjskew_5 Wjskew_5
% %           ddd+hhh/24+50/24/60 nanmean(Uj_6) nanmean(Vj_6) nanmean(Wj_6) nanvar(Uj_6) nanvar(Vj_6) nanvar(Wj_6) alphaj_6 betaj_6 wuj_6(1,2) wvj_6(1,2) wtj_6(1,2) Ujskew_6 Vjskew_6 Wjskew_6];
% %       
% % 
% %     sdg=['save ', drive_w '\AlanBrewerRequest\Alanmat_',jd hr,'.prn Alanview -ascii'];
% %     eval(sdg)
    
% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% %
% % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % % % %% % %


if clearit
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 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 Tseasnake 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
end 
%   xship_1=[mean_wind_1 Us Tsea_1 Tship_1 Qsea_1 Qship_1 Rs_ship_1 Rl_1 ship_rain_1 zi P zu_ship zt_ship zq_ship];

%   1)   mean_wind_1 
%   2)   Us 
%   3)   Tsea_1 
%   4)   Tship_1 
%   5)   Qsea_1 
%   6)   Qship_1 
%   7)   Rs_ship_1 
%   8)   Rl_1 
%   9)   ship_rain_1 
%   10)  zi 
%   11)  P 
%   12)  zu_ship 
%   13)  zt_ship 
%   14)  zq_ship

%   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);
if clearit
clear tau2 spinup2 ans Rhvais Tsea Tvais  cog cogm gpstim lat lon lrgm %U V W
clear imdm imum imta imsol imrh imrain sog sogm tsg tss org pir psp
end;
%clear Uj_1 Uj_2 Uj_3 Uj_4 Uj_5 Uj_6 
%clear Vj_1 Vj_2 Vj_3 Vj_4 Vj_5 Vj_6 
%clear Wj_1 Wj_2 Wj_3 Wj_4 Wj_5 Wj_6 
%clear Ustream_1 Ustream_2 Ustream_3 Ustream_4 Ustream_5 Ustream_6 
%clear Vstream_1 Vstream_2 Vstream_3 Vstream_4 Vstream_5 Vstream_6 
%clear Wstream_1 Wstream_2 Wstream_3 Wstream_4 Wstream_5 Wstream_6 
%clear Tsonic Tsonic_1 Tsonic_2 Tsonic_3 Tsonic_4 Tsonic_5 Tsonic_6 
%clear abshum abshum_1 abshum_2 abshum_3 abshum_4 abshum_5 abshum_6