disp(['read_scs_day_', cruise, '_', year]);
%reads ship data system files

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_flux, 'DAY', jd,'\scs' year(2:4) jd hr '_raw.txt']
    %disp(['Reading scs file from hour ',jd,':',hr]);
    flist=fopen(e,'r');
    if flist>0,                                 %if the file exists,
        while ~feof(flist),
            temp =textscan(flist,'%f %f %f%c %f%c %f %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]','delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', 'Auxiliary Ship');
%             temp =textscan(flist,'%2f%2f%3f %f %f%c %f%c %f %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]','delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', 'Auxiliary Ship');
            temp{1,4}=double(temp{1,4});
            temp{1,6}=double(temp{1,6});
            
%             for ii=1:19
%                 if length(temp{1,ii})~=1800
%                     temp{1,ii}=[temp{1,ii};ones(1800-length(temp{1,ii}),1)*temp{1,ii}(length(temp{1,ii})-1)];
%                 end;
%             end;    
            stx=cell2mat(temp)';
            stx(2:19,any(isnan(stx),1))=NaN;      %put NaN into any rows containing NaNs (to remove bad lines, gap...)
            
            [nr1,nl1]=size(stx);
                garbage = [garbage stx(1,:)];
                timegps = [timegps stx(2,:)];
                lat =     [lat stx(3,:)./100];
                latns =   [latns stx(4,:)];        %in ASCII (use char to convert in characters) / 78=N, 87=W, 83=S, 69=E
                lon =     [lon stx(5,:)./100];
                lons =    [lons stx(6,:)];
                cog =     [cog stx(7,:)];
                sog =     [sog stx(8,:)*0.514];   %converted from knot to m/s
                tsg =     [tsg stx(9,:)];
                tss =     [tss stx(10,:)];
                rh =      [rh stx(11,:)];
                ta =      [ta stx(12,:)];
                solar =   [solar stx(13,:)];
                imrain =  [imrain stx(14,:)];
                lrg =     [lrg stx(15,:)];
                odec =    [odec stx(16,:)];
                imu =     [imu stx(17,:)*.514];   %converted from knot to m/s
                imd =     [imd stx(18,:)];
                ir  =     [ir stx(19,:)];

        end;
    end;
end;         %0-23 loop


    %%%%%%%%%%%%%%%%%%%%%%%%%  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


%convet from degs.minutes to dec degrees
latmin=(lat-floor(lat))*100;
lat2=floor(lat)+latmin/60;
slats=find(latns == 83);       %ASCII: 78=N, 87=W, 83=S, 69=E
if (length(slats) >= 1)
   lat2(slats)=-lat2(slats);
end
lonmin=(lon-floor(lon))*100;
lon2=floor(lon)+lonmin/60;
slons=find(lons == 87);        %ASCII: 78=N, 87=W, 83=S, 69=E
if (length(slons) >= 1)
   lon2(slons)=-lon2(slons);
end
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));
jd_scs = ddd + (hh+((mm+(ss/60))/60))/24;
end

% lgpstim=length(gpstim)
% ljd_scs=length(jd_scs)
%%%%%%%DEW 05/21/2007
%gpstim = jd_scs;
%%%%%%%%%%%%%%%%%%%


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(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 f:\AMMA_2007\RHB\flux\raw\SCSrad jd_scs solar ir gpstim  % save data for comparison with SCS using rad_stat_06.m
% rad_stat_06
lond=lonm;
latd=latm;
figure;map_Stratus2007;

jd_scs=gpstime./24+ddd;

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 0  % %Plots using HP time conversion jd_scs DEW 5/20/2007

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

np=length(sogm);

if prtit==1,
    prt_jas_scs_Stratus_07;
end;