disp('read_scs_day_rhb_07')
%reads ship data system files
%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


jd=num2str(ddd);
if ddd<100,
    jd=['0' num2str(ddd)];
end;
if ddd<10,
    jd=['00' num2str(ddd)];
end;

plotit=1;   %set to 1 to enable plots
prtit=1;    %set to 1 to call the routine that writes results in ascii files

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=0;                              %line counter
backchk=0;                          %first time index
timegps=[];lat=[];latns=[];lon=[];lons=[];cog=[];sog=[];tsg=[];tss=[];
rh=[];ta=[];solar=[];imrain=[];lrg=[];odec=[];imu=[];imd=[];ir=[];garbage=[];
secs=[];  %DEW 5/21/2007
sx=[];headr=[];
stx1=[];headr1=[];
clear st3;

for jam=00:23,                                         %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]);
            [nr1,nl1]=size(stx);
            if nr1==17
                %secs = [garbage stx(1,:)];  %DEW 5/21/2007
                secs = [secs stx(1,:)];   %DEW 5/21/2007
                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;
            %length(secs) % DEW 5/23/2007

        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
                % find extra S's
                %if headr(33) ~= ','
                %   headr(32:35)
                %end
                latns =   [latns headr(32)];
                lons =    [lons headr(44)];
            else
                dmp3=fgetl(flist2);
            end;
        end;
    end;
end;         %0-23 loop

if isempty(timegps)
    timegps = [timegps 0];
    lat =     [lat 0];
    latns =   [latns 0];
    lon =     [lon 0];
    lons =     [lons 0];
    cog =     [cog 0];
    sog =     [sog 0];
    tsg =     [tsg 0];
    tss =     [tss 0];
    rh =      [rh 0];
    ta =      [ta 0];
    solar =   [solar 0];
    imrain =  [imrain 0];
    lrg =     [lrg 0];
    odec =    [odec 0];
    imu =     [imu 0];
    imd =     [imd 0];
    ir =      [ir 0];
else

    %%%%%%%%%%%%%%%%%%%%%%%%%  in case first line is from previous day

    %%DEW 5/21/2007 This section of code messes up timegps for AMMA07 JD
    %%140 when the gps time in the P6 file shits by almost an hour.
    %% should be using secs
    if 0 %DEW 5/21/2007
        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;
    end;  % DEW 5/21/2007
end;

%convert deg.min to dec deg
latmin=(lat-floor(lat))*100;
lat2=floor(lat)+latmin/60;
slats=find(latns == 'S');
if 0  %The following code was taken out because of JD287 STRATUS06 corrupt P6 files
    if (length(slats) >= 1)
        lat2(slats)=-lat2(slats);
    end
end
if 1
    % The following code is to fix corrupt file error
    % It assumes that all P6 data are South of the equator
    % JD 287 STRATUS06 the ltns file is larger than the lat2 file
    % This is not possilbe unless there are some extra "S" in the SCS files
    %DEW Oct 15, 2006 fix assumes if latns > lat2 all LATs need to be made negative
    if length(slats) < length(lat2) & (length(slats) >= 1)  %DEW
        if (length(slats) >= 1)  %DEW
            lat2(slats)=-lat2(slats);
        else   %DEW
            lat2=-lat2;   %DEW
        end
    end
end
lonmin=(lon-floor(lon))*100;
lon2=floor(lon)+lonmin/60;

lat=lat2;
lon=-lon2;

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

%%%%%%%%
% DEW Oct 18, 2006
% This code isn't used for calculating the Julian decimal day
% GPS time sent with SCS is used and should be kept as a check on SCS time
% Time in flux data files is the number of 100-nano second intervals since
% Jan 1, 1601 divided by 10000. It was chosen to maintain compatibility
% with an HP-UNIX version of the code from long ago
%  0             'dd-mmm-yyyy HH:MM:SS'   01-Mar-2000 15:45:17
dates = datestr(datenum(1601,1,1) + datenum(secs*10000*100e-9/86400));
if length(dates) > 0
hh = str2num(dates(:,13:14));
mm = str2num(dates(:,16:17));
ss = str2num(dates(:,19:20));
jdscs = ddd + (hh+((mm+(ss/60))/60))/24;
end

lgpstim=length(gpstim)
ljdscs=length(jdscs)
%%%%%%%DEW 05/21/2007
%gpstim = jdscs;
%%%%%%%%%%%%%%%%%%%


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;

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


st2=zeros(16,length(gpstim));

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=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;
[n m]=size(st2);
if n>1
    for iii=2:n;
        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;
else
    jk=2;
end;

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

%1-min medians of all quantities
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:16)=median(st2(jk(iii-1):jk(iii)-1,1:16));
    lrgcm(iii)=median(lrgc(jk(iii-1):jk(iii)-1));
    lrgsm(iii)=median(lrgs(jk(iii-1):jk(iii)-1));
end;

gpstime=st3(:,1);
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   1 %% DEW 5/21/2007 skip plotting using gps time
    figure;plot(gpstim./24+ddd,cog,gpstim./24+ddd,lrg);title('SCS course and gyro')
    figure;plot(gpstim./24+ddd,imu);title('SCS wind speed')
    figure;plot(gpstim./24+ddd,imrain);title('SCS rain')
    figure;plot(gpstim./24+ddd,imd);title('SCS wind direction')
    figure;plot(gpstim./24+ddd,sog);title('SCS speed over ground')
    figure;plot(lon,lat,'.');title('SCS Lat vs Lon')
    figure;plot(gpstim./24+ddd,tsg,gpstim./24+ddd,ta);title('SCS Tsea   Tair')
    figure;plot(gpstim./24+ddd,tss);title('SCS salinity')
    figure;plot(gpstim./24+ddd,rh);title('SCS RH')
    figure;plot(gpstim./24+ddd,solar,gpstim./24+ddd,ir);title('SCS solar & IR')

end %% DEW 5/21/2007 skip plotting using gps time

save d:\data\SCSrad jdscs solar ir gpstim  % save data for comparison with SCS using rad_stat_06.m
rad_stat_06
lond=lonm;
latd=latm;
figure;map_amma2007;


figure;plot(lonm,latm,'.');title('SCS lat vs lon')
figure;plot(gpstime./24+ddd,sogm,gpstime./24+ddd,odecm);title('SCS speed and ODEC')
figure;plot(gpstime./24+ddd,cogm);title('SCS course')
figure;plot(gpstime./24+ddd,imum);title('SCS wind speed')
figure;plot(gpstime./24+ddd,imdm);title('SCS wind direction')
figure;plot(gpstime./24+ddd,lrgm);title('SCS gyro')

if 1  % %Plots using HP time conversion jdscs DEW 5/20/2007

    %figure;plot(jdscs,cog,jdscs,lrg);title('SCS course and gyro')
    %figure;plot(jdscs,imu);title('SCS IMET wind speed')
    %figure;plot(jdscs,imrain);title('SCS IMET rain')
    %figure;plot(jdscs,imd);title('SCS IMET wind direction')
    %figure;plot(jdscs,sog);title('SCS speed over ground')
    %figure;plot(lon,lat,'.');title('SCS Lat vs Lon')
    %figure;plot(jdscs,tsg,jdscs,ta);title('SCS Tsea   Tair')
    %figure;plot(jdscs,tss);title('SCS salinity')
    %figure;plot(jdscs,rh);title('SCS RH')
    figure;plot(jdscs,solar,jdscs,ir);title('SCS solar & IR')

end

np=length(sogm);

length(jdscs)  % DEW 5/23/2007

if prtit==1,
    prt_jas_scs_rhb_07;
end;