%This program takes a 1hr wband data set that has been read in and
%formatted by the WbandPlot_2 program.
%It finds a mean velocity for the cloud layer and separates lowpass and
%highpass signals.  The highpass signal is approximately the vertical
%motion of the ship.
%At the end it reads the motion file and compares

h=nanmedian(htt);
%h=cldtop(:,2);
h=nanmedian(cldtop(:,2));
kk=find(ht<h);
topi=max(kk);
cldthk=nanmean(cldtop(:,8))/12600*25/1e3;
kk=find(ht<h-cldthk);
boti=max(kk);
dkcld=min(25,topi-boti+2);


j=1;
dbb=-35:0;
jj=find(Ref(topi-2,:)<-20);
clear Pi ygg;
for dk=3:dkcld;
    jj=find(Ref(topi-dk,:)>-35);

%[P,S] = polyfit(Ref(topi-dk,jj),DV(topi-dk,jj)-DV(topi-2,jj),2);
[P,S] = polyfit(Ref(topi-dk,jj),DV(topi-dk,jj),2);

yg=polyval(P,dbb);
ygg(j,:)=yg;
Pi(j,:)=P;
j=j+1;
end
figure;plot(dbb,mean(ygg),dbb,polyval(mean(Pi),dbb),'o');

j=1;
velk=[];
for dk=1:dkcld;
vel = polyval(mean(Pi),Ref(topi-dk,:));
velk(j,:)=vel;
j=j+1;
end


vc=median(DV(topi-dkcld+1:topi,:)-velk);

%dk=dkcld;
Pg=[];
[Pg,Sg] = polyfit(mean(Ref(topi-dkcld:topi-floor(dkcld/2),:)),mean(DV(topi-dkcld:topi-floor(dkcld/2),:))-vc,2);
dbt=[ones(size(dbb')) exp(.06*dbb')];
tg=nanmean(Ref(topi-dkcld:topi-floor(dkcld/2),:))';vg=nanmean(DV(topi-dkcld:topi-floor(dkcld/2),:))'-vc';xp=[ones(size(tg)) exp(.06*tg)];aa=xp\vg;ygb=dbt*aa;
figure;plot(dbb,mean(ygg),dbb,polyval(mean(Pi),dbb),'o',dbb,polyval(Pg,dbb),'--',dbb,ygb,'-x',tg,vg,'.');

velk2=[];velk3=[];
j=1;
for dk=1:dkcld;
vel = polyval(Pg,Ref(topi-dk,:));
velk2(j,:)=vel;
tgg=Ref(topi-dk,:)';xpg=[ones(size(tgg)) exp(.06*tgg)];
velg=xpg*aa;
velk3(j,:)=velg';
j=j+1;
end

dk=5;vc2=mean(DV(topi-dkcld:topi-dk,:)-velk2(dk:dkcld,:));vc3=mean(DV(topi-dkcld:topi-dk,:)-velk3(dk:dkcld,:));

wp=[.02];ws= [.01];rp=3;rs=15;[n,wmn]=cheb1ord(wp,ws,rp,rs);
[b,a]=cheby1(n,rp,wmn);
vcfL=filter(b,a,vc-mean(vc));%\
vcfH=vc-vcfL-mean(vc);
vcfL3=filter(b,a,vc3-mean(vc3));%\
vcfH3=vc3-vcfL-mean(vc3);



km=floor(dkcld/2);figure;ii=find(Ref(topi-km,:)>-35);plot(hh(ii),DV(topi-km,ii)-velk2(km,ii)-vcfH(ii));

%     vkonsa=[];%ones(10,36000);
%     vkons=[];
%     fclose('all');
%     fm = [way_raw_motion_wband '\' year ddds hr '00Kongsberg.txt'];
%     FID2 = fopen(fm,'r');
%     dumx=fgetl(FID2)
%     dumx=fgetl(FID2);
%     dumx=fgetl(FID2)
%     dumx=fgetl(FID2);
%     for i=1:90%For some reason a fscanf(fid,format,[10,inf]); terminates after 410 lines.  No idea why, but this works. I hate Bill Gates.
%     vkonsa=fscanf(FID2,'%2d %*c %2d %*c %4d %2d %*c %2d %*c %g %g %g %g %g',[10,400]);
%     vkons=[vkons vkonsa];
%     vkonsa=[];
%     end;
%     hhk=vkons(4,:)+vkons(5,:)/60+vkons(6,:)/3600;
%      
read_wband_motion
    figure;plot(hh,vcfH,hhk,vkons(9,:));xlabel('Time (hr)');ylabel('Vertical Motion (m/s)');