% read_wave
% Reads Riegl 10Hz laser range meter for wave height. (slow)
% NTAS 2009 :: 2009-06-12 :: LB

% clear all;fclose all;hhh=9;read_Riegl_wave_VOCALS

disp('read_altimeter')

% ddd=85;hhh=7;
% read_parameters_GASEXIII;

fclose('all');
jd=sprintf('%03i',ddd);
clear temp stx;

stx=[];temp=[];

for hhh=00:06:23,					    %cycle thru 24 hourly stats files
    hr=sprintf('%02i',hhh);
    e=[way_raw_data_flux, year(3:4),jd,'\lsr0' year(3:4) jd hr '_raw.txt'];
    disp(['Reading altimeter file from hour ',num2str(hhh)]);
    flist=fopen(e,'r');
    
    temp=textscan(flist,'%2f%2f%3f %*f%f %*f%f %*f%f %*[^\n]','headerLines',1,'delimiter','; ','emptyValue',NaN,'treatAsEmpty',{'r','a','q'});
    
%     stx=[stx cell2mat(temp)'];
    stx=cell2mat(temp)';
% end;

ii=find(~isnan(stx(4,:)));
minu=stx(1,ii);
sec=stx(2,ii);
msec=stx(3,ii);
wave=stx(4,ii);
amplitude=stx(5,ii);

jd_wave=ddd + (hhh+(minu+(sec+msec/1000)/60)/60)/24;
timer=minu*60+sec+msec/1000;

mwave=median(wave);
wave=wave-mwave;

figure;plot(jd_wave,wave,'.-');
xlabel(['JD ' year]);ylabel('Wave height non motion corrected  (m)')
title(['Wave ht ', cruise, year]);
%  if prtit; print('-djpeg ', [way_images_flux 'WaveHt_' num2str(ddd),'_',hr, '.png']); end

Owave=detrend(wave);
med_Owave=median(Owave);
sort_Owave=sort(Owave);
jk=floor(0.5*length(sort_Owave));
jkp=floor(.667*length(sort_Owave));
jkm=floor(.333*length(sort_Owave));
sig_Owave=0.5*(sort_Owave(jkp)-sort_Owave(jkm));
Otimer=timer;

if abs(Owave(1)-med_Owave) > 7*sig_Owave,
   Owave(1)=med_Owave;
end;
for kkk=2:length(Owave);
   if abs(Owave(kkk)-med_Owave) > 7*sig_Owave,
      Owave(kkk)=Owave(kkk-1);
%       Otimer(kkk)=Otimer(kkk-1)+0.5;
   end;
end;

[Swave,Fwave]=psd2(detrend(wave),length(wave),10);[Swaves,Fwaves]=specsmoo(Swave,10);
[SOwave,FOwave]=psd2(detrend(Owave),length(Owave),10);[SOwaves,FOwaves]=specsmoo(SOwave,10);

figure;loglog(Fwaves,Fwaves.*Swaves,Fwaves,Fwaves.*Swaves,'--');xlabel('Frequency (Hz)');ylabel('f*Sw(f)  (m^2)')
 if prtit; print('-dpng ', [way_images_flux 'WaveSpec_' num2str(ddd),'_',hr, '.png']); end
 
sig_wave_ht=4.*sqrt(var(wave));
sig_Owave_ht=4.*sqrt(var(Owave));
disp(['Non-compensated significant wave height = ',num2str(sig_wave_ht)])

figure;plot(timer,wave,Otimer,Owave);title(['Sig W Ht ',num2str(sig_wave_ht),'  ',num2str(sig_Owave_ht)])
xlabel(['Time (sec)']);legend('Wave height raw (m)','Wave height despiked (m)')
 if prtit; print('-djpeg ', [way_images_flux 'WaveHt_' num2str(ddd),'_',hr, '.png']); end
 
 
end;
