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

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

    %disp('read_altimeter')

    % ddd=85;jam=7;
    % read_parameters_GASEXIII;

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


    for jam=00:06:23,					    %cycle thru 24 hourly stats files
        hr=sprintf('%02i',jam);
        e=[way_raw_data_flux, yyyy(3:4),jd,'\lsr0' yyyy(3:4) jd hr '_raw.txt'];
        disp(['Reading altimeter file ',e]);
        flist=fopen(e,'r');
        stx=[];temp=[];

        if flist>0
            while ~feof(flist) % end-of-file not reached
                try
                    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)'];
                catch
                    % if one line is cut short, blank it out totally
                    [l6,c6]=size(temp{1,6});
                    if l6==0 || c6==0; continue; end % read next batch of lines
                    for ll=1:6
                        if length(temp{1,ll})~=l6
                            temp{1,ll}(l6+1)=[];
                        end;
                    end;
                    if isempty(temp{1,1});
                        temp{1,6}=temp{1,6}';
                    end;
                    stx=[stx cell2mat(temp)'];
                end; %end try/catch
            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 + (jam+(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 ' yyyy]);ylabel('Wave height non motion corrected  (m)')
            %     title(['Wave ht ', cruise, yyyy]);
            %  if prtit; print('-djpeg ', [way_images_flux 'WaveHt_' num2str(ddd),'_',hr, '.png']); end

            Owave=wave;
            if max(Owave) > 0;   % DEW 16 Oct 2011 Days with glassy seas have no signal and crash program
                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.*abs(Swaves-2.5e-4),'--');xlabel('Frequency (Hz)');ylabel('f*Sw(f)  (m^2)')
                title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD altimeter,  spectrum: hour %02i',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,jam),'FontWeight','Bold','Interpreter','none')
                xlabel('Frequency (Hz)');
                ylabel('f*Sw(f)  (m^2)');grid
                annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
                if prtit; print(graphdevice,[way_images_flux, cruise  '_WaveSpectra_hour' hr '_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); 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)])
                title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD altimeter,  spectrum: hour %02i',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,jam),'FontWeight','Bold','Interpreter','none')
                xlabel('Time (s)');legend('Wave height raw (m)','Wave height despiked (m)','Location','Best')
                ylabel('f*Sw(f)  (m^2)');grid;xlim([0 3600])
                annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');

                if prtit; print(graphdevice,[way_images_flux, cruise  '_WavHeight_hour' hr '_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end
            end  % test for data in hour DEW 16 Oct 2011
        else
            disp(sprintf('\b<--File not found.'))
        end;         %end if flist
    end;

    fclose all;
