disp('eval_scs_rhb_06')
%reads etl ship data system files, computes 1-min medians of selected variables, writes processed
%ascii file
%data files are expected to be in directory named 'dayddd'

%ddd is julian day it must be specified from the keyboard or in a program
%that calls read_scs
%cwf 11/03

%changed by cwf to inlcde imet PIR  12/2004


fclose('all');
st2=[];timegps=[];lat=[];latns=[];lon=[];cog=[];sog=[];tsg=[];tss=[];
rh=[];ta=[];solar=[];imrain=[];lrg=[];odec=[];imu=[];imd=[];ir=[];lons=[];

st3=[];gpstime=[];latm=[];lonm=[];cogx=[];sogx=[];tsgm=[];tssm=[];
rhm=[];tam=[];solarm=[];imrainm=[];lrgm=[];odecm=[];imum=[];imdm=[];irm=[];
secs=[];
sx=[];headr=[];
stx1=[];headr1=[];
scstime=[];

%make string out of file jul day
jd=num2str(ddd);
if ddd<100
    jd=['0' num2str(ddd)];
end;
if ddd<10
    jd=['00' num2str(ddd)];
end;

hr='00';
g=0.0098;							%adiabatic lapse rate
eps_w=0.97;							%emissivity of water
sig_sb=5.67e-8;					    %Stefan Boltzmann constant
blk=1;								%10-min index
jax=1;                              %line counter
backchk=0;                          %first time index
stx=[]; st2=[];
clear data

for jam=hhh:hhh,						%cycle thru 24 hourly stats files
    if jam<10,
        hr=['0' num2str(jam)];
    else
        hr=num2str(jam);
    end;
    e=[way_raw_data, 'day' ,jd,'\P6_' jd hr];
    disp(['Reading scs file from hour ',jd,':',hr]);
    flist=fopen(e,'r');
    if flist>0,					%if the file exists,
        while ~feof(flist),
            dmp=fgets(flist);       %read header line
            stx=fscanf(flist,'%g,%g,%g%*c,%10f%*c,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g',[17,inf]);
            %             if isempty(stx)
            %                 stx=ones(17,1683).*NaN
            %             end;
            [nr1,nl1]=size(stx);
            if nr1==17
                secs = [garbage stx(1,:)];
                timegps = [timegps stx(2,:)];
                lat =     [lat stx(3,:)./100];
                lon =     [lon stx(4,:)./100];
                cog =     [cog stx(5,:)];
                sog =     [sog stx(6,:)*0.514];
                tsg =     [tsg stx(7,:)];
                tss =     [tss stx(8,:)];
                rh =      [rh stx(9,:)];
                ta =      [ta stx(10,:)];
                solar =   [solar stx(11,:)];
                imrain =  [imrain stx(12,:)];
                lrg =     [lrg stx(13,:)];
                odec =    [odec stx(14,:)];
                imu =     [imu stx(15,:)*.514];
                imd =     [imd stx(16,:)];
                ir  =     [ir stx(17,:)];
            else
                x=fgets(flist);
            end;
        end;
    end;
    flist2=fopen(e,'r');
    if flist2>0,					%if the file exists,
        dmp2=fgetl(flist2);
        while ~feof(flist2),
            headr=fgetl(flist2);
            if length(headr)>100
                latns =   [latns headr(32)];
                lons =    [lons headr(44)];
            else
                dmp3=fgets(flist2);
            end;
        end;
    end;
end;         %0-23 loop
%    if jam<10,
%        hr=['0' num2str(jam)];
%    else
%        hr=num2str(jam);
%    end;
%    e=['d:\data\day',jd,'\P6_' jd hr];
%    %e=[way_data '\day',jd,'\P6_' jd hr];
%    disp(['Reading scs file from hour ',jd,':',hr]);
%    flist=fopen(e,'r');
%    if flist>0,
%    while ~feof(flist),
%       			%if the file exists,
%          dmp=fgets(flist);       %read header line
%          %stx=fscanf(flist,'%g,%g,%g%c,%g%c,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,',[18,inf]);
%          stx=fgets(flist);
%          stx=deblank(stx);
%          if length(stx)==137,
%             garbage = str2num(stx(1:14));
%             timegps = [timegps str2num(stx(16:21))];
%             lat =     [lat str2num(stx(23:31))/100];
%             latns =   [latns stx(32)];
%             lon =     [lon str2num(stx(34:43))/100];
%             cog =     [cog str2num(stx(46:50))];
%             sog =     [sog str2num(stx(52:56))*0.514];
%             tsg =     [tsg str2num(stx(58:63))];
%             tss =     [tss str2num(stx(65:71))];
%             rh =      [rh str2num(stx(73:79))];
%             ta =      [ta str2num(stx(81:87))];
%             solar =   [solar str2num(stx(89:95))];
%             imrain =  [imrain str2num(stx(97:102))];
%             lrg =     [lrg str2num(stx(104:108))];
%             odec =    [odec str2num(stx(110:117))];
%             imu =     [imu str2num(stx(119:124))*.514];
%             imd =     [imd str2num(stx(126:131))];
%             ir  =     [ir str2num(stx(133:137))];
%          else
%             disp('bad line')
%             stx=[];
%          end;
%          jax=jax+1;
%       end;                %end flist
%    end;                   %end while
% end;         %0-23 loop

if ~isempty(timegps)
    %     if timegps(1)>10,
    %         timegps(1)=str2num(hr).*10000;
    %     end;
    if timegps(1)>timegps(2),
        timegps(1)=0;
    end;
    for ii=2:length(timegps)
        if timegps(ii)<timegps(ii-1),
            timegps(ii)=timegps(ii-1)+2;
        end;
    end;

    gpshour=fix(timegps./10000);
    gpsmin=fix(timegps./100)-gpshour.*100;
    gpssec=timegps-gpshour.*10000-gpsmin.*100;
    gpstim=gpshour+gpsmin./60+gpssec./60./60; %Decimal time in hours

    dates = datestr(datenum(1601,1,1) + datenum(secs*10000*100e-9/86400));
    hh = str2num(dates(:,13:14));
    mm = str2num(dates(:,16:17));
    ss = str2num(dates(:,19:20));
    scstim=(hh+((mm+(ss/60))/60));
length(scstim)
    rdcon=pi/180;

    for iii=2:length(gpstim),
        if gpstim(iii)<=gpstim(1),
            gpstim(iii)=gpstim(iii-1)+gpstim(2)-gpstim(1);
        end;
        if rh(iii)==0,
            rh(iii)=rh(iii-1);
        end;
        if tsg(iii)==0,
            tsg(iii)=tsg(iii-1);
        end;
        if tss(iii)==0,
            tss(iii)=tss(iii-1);
        end;
        if ta(iii)==0,
            ta(iii)=ta(iii-1);
        end;
        if rh(iii)==0,
            solar(iii)=solar(iii-1);
        end;
        if imrain(iii)==0,
            imrain(iii)=imrain(iii-1);
        end;
        if imd(iii)==0,
            imd(iii)=imd(iii-1);
        end;
        if imu(iii)==0,
            imu(iii)=imu(iii-1);
        end;
        if cog(iii)==0,
            cog(iii)=cog(iii-1);
        end;
        if lrg(iii)==0,
            lrg(iii)=lrg(iii-1);
        end;
        if lat(iii)==0,
            lat(iii)=lat(iii-1);
        end;
        if lon(iii)==0,
            lon(iii)=lon(iii-1);
        end;
        if ir(iii)==0,
            ir(iii)=ir(iii-1);
        end;
    end;

    st2=zeros(17,length(scstim));

    %convert deg.min to dec deg
    latmin=(lat-floor(lat))*100;
    lat2=floor(lat)+latmin/60;
    slats=find(latns == 'S');
    if (length(slats) >= 1)
        lat2(slats)=-lat2(slats);
    end
    lonmin=(lon-floor(lon))*100;
    lon2=floor(lon)+lonmin/60;

    lat=lat2;lon=-lon2;

    %*********************  account for change of SW sensor on day 342
    %if ddd>341;solar=solar/1.04;end;
    %********************


    st2(1,:)=gpstim;
    st2(2,:)=lat;
    st2(3,:)=lon;
    st2(4,:)=cog;
    st2(5,:)=sog;
    st2(6,:)=tsg;
    st2(7,:)=tss;
    st2(8,:)=rh;
    st2(9,:)=ta;
    st2(10,:)=solar;
    st2(11,:)=imrain;
    st2(12,:)=lrg;
    st2(13,:)=odec;
    st2(14,:)=imu;
    st2(15,:)=imd;
    st2(16,:)=ir;
    st2(17,:)=scstim;

    st2=st2';

    st2(:,4)=sog'.*cos(cog'*rdcon);    %convert to N
    st2(:,5)=sog'.*sin(cog'*rdcon);    %convert to E
    st2(:,14)=imu'.*cos(imd'*rdcon);   %convert to N
    st2(:,15)=imu'.*sin(imd'*rdcon);   %convert to E
    lrgc=cos(lrg'*rdcon);        %convert to N
    lrgs=sin(lrg'*rdcon);        %convert to E

    jk=[];
    flagup=1;

    for iii=2:length(st2);
        if flagup,
            oldgps=st2(iii-1,1).*60;
            flagup=0;
        end;
        if oldgps+1 < st2(iii,1).*60,
            jk=[jk
                iii-1];
            flagup=1;
        end;
    end;


    st3=zeros(length(jk),17);
    lrgsm=zeros(1,length(jk));
    lrgcm=zeros(1,length(jk));

    %1-min medians of all quantities
    if jk(1)<3,
        jk(1)=3;
    end;
    st3(1,:)=median(st2(1:jk(1)-1,:));
    lrgcm(1)=median(lrgc(1:jk(1)-1,:));
    lrgsm(1)=median(lrgs(1:jk(1)-1,:));

    for iii=2:length(jk);
        st3(iii,1:17)=median(st2(jk(iii-1):jk(iii)-1,1:17));
        lrgcm(iii)=median(lrgc(jk(iii-1):jk(iii)-1));
        lrgsm(iii)=median(lrgs(jk(iii-1):jk(iii)-1));
    end;

    gpstime=st3(:,1);
    scstime=st3(:,17);
    latm=st3(:,2);%./100;
    lonm=st3(:,3);%./100;
    sogx=st3(:,4);
    sogy=st3(:,5);
    tsgm=st3(:,6);           %Tsea
    tssm=st3(:,7);          %salinity
    rhm=st3(:,8);           %relative humidity
    tam=st3(:,9);           %Tair
    solarm=st3(:,10);        %solar
    rainm=st3(:,11);         %Rain
    odecm=st3(:,13);         %ODEC
    imux=st3(:,14);
    imuy=st3(:,15);
    irm=st3(:,16);

    sogm=sqrt(sogx.*sogx+sogy.*sogy);				%1-min sog
    cogm=mod(atan2(sogy,sogx)/rdcon+360,360);
    imum=sqrt(imux.*imux+imuy.*imuy);				%1-min imet speed, m/s
    imdm=mod(atan2(imuy,imux)/rdcon+360,360);
    lrgm=mod(atan2(lrgsm,lrgcm)/rdcon+360,360);		%laser ring gyro

    for iii=2:length(gpstim),
        if gpstim(iii)<gpstim(1),
            gpstim(iii)=gpstim(iii-1)+gpstim(2)-gpstim(1);
        end;
    end;

    if length(imum)<60,
        imum=[imum
            imum(length(imum))+0.*(1:60-length(imum))'];
    end;
    if length(gpstime)<60,
        gpstime=[gpstime
            gpstime(length(gpstime))+(gpstime(2)-gpstime(1)).*(1:60-length(gpstime))'];
    end;
    if length(imdm)<60,
        imdm=[imdm
            imdm(length(imdm))+0.*(1:60-length(imdm))'];
    end;
    if length(rhm)<60,
        rhm=[rhm
            rhm(length(rhm))+0.*(1:60-length(rhm))'];
    end;
    if length(latm)<60,
        latm=[latm
            latm(length(latm))+0.*(1:60-length(latm))'];
    end;
    if length(lonm)<60,
        lonm=[lonm
            lonm(length(lonm))+0.*(1:60-length(lonm))'];
    end;
    if length(tam)<60,
        tam=[tam
            tam(length(tam))+0.*(1:60-length(tam))'];
    end;
    if length(cogm)<60,
        cogm=[cogm
            cogm(length(cogm))+0.*(1:60-length(cogm))'];
    end;
    lrgm=lrgm';
    if length(lrgm)<60,
        lrgm=[lrgm
            lrgm(length(lrgm))+0.*(1:60-length(lrgm))'];
    end;
    if length(sogm)<60,
        sogm=[sogm
            sogm(length(sogm))+0.*(1:60-length(sogm))'];
    end;
    if length(tsgm)<60,
        tsgm=[tsgm
            tsgm(length(tsgm))+0.*(1:60-length(tsgm))'];
    end;
    if length(solarm)<60,
        solarm=[solarm
            solarm(length(solarm))+0.*(1:60-length(solarm))'];
    end;
    if length(rainm)<60,
        rainm=[rainm
            rainm(length(rainm))+0.*(1:60-length(rainm))'];
    end;
    if length(irm)<60,
        irm=[irm
            irm(length(irm))+0.*(1:60-length(irm))'];
    end;
end;