% read_scs_day
% VOCALS 2008 :: 2008-10-17 :: Simon de Szoeke
%reads ship data system files
disp(['read_scs_day: ', cruise, ' ', year]);

jd=sprintf('%03i',ddd);
% plotit,prtit,saveit inherited from flux_eval
% plotit=true;  % set true to enable plots
% prtit=true    % set true to output plots to files
% saveit=true;  % set true 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=[];
stx=[];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,'\scs0', year(3:4) jd hr '_raw.txt'];
    disp(['Reading scs file from hour ',jd,':',hr]);
    flist=fopen(e,'r');
    temp=[];stx=[];
    if flist>0,                                 %if the file exists,
        while ~feof(flist)
            try
                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);
                temp{1,4}=double(temp{1,4});
                temp{1,6}=double(temp{1,6});
                % reading gps time should be more like this, but hard to
                % handle different data types in matlab structures.
                %temp=textscan(flist,'%f %02d%02d%05f %f%c %f%c %f %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]','delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN);
                stx=[stx cell2mat(temp)'];
            catch
                [l19,c19]=size(temp{1,19});
                for ll=1:19
                     if length(temp{1,ll})~=l19
                         temp{1,ll}(l19+1)=[];
                     end;
                end;
                if isempty(temp{1,1}); continue; end
                if isempty(temp{1,1});
                    disp('Should never do this!')
                    temp{1,7}=temp{1,7}';
                end;
             stx=[stx cell2mat(temp)'];
            end;
        end;  %end while
        
        [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,:)*.514];
        imu =     [imu stx(17,:)*.514];   %converted from knot to m/s
        imd =     [imd stx(18,:)];
        ir  =     [ir stx(19,:)];
    end;
    %     end;
end;         %0-23 hour 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 true %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;
lat2(char(latns)=='S')=-lat2(char(latns)=='S');
lonmin=(lon-floor(lon))*100;
lon2=floor(lon)+lonmin/60;
lon2(char(lons)=='W')=-lon2(char(lons)=='W');
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 in hours

%%%%%%%%
% 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 ~isempty(dates)
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;
end;

% why can't rh, rain, or course be zero?
%ii=find(imrain==0);imrain(ii)=imrain(ii-1);
ii=find(rh==0);rh(ii)=rh(ii-1);
ii=find(tsg==0);tsg(ii)=tsg(ii-1);
ii=find(tss==0);tss(ii)=tss(ii-1);
ii=find(ta==0);ta(ii)=ta(ii-1);
ii=find(rh==0);rh(ii)=rh(ii-1);
ii=find(imd==0);imd(ii)=imd(ii-1);
ii=find(imu==0);imu(ii)=imu(ii-1);
%ii=find(cog==0);cog(ii)=cog(ii-1);
ii=find(lrg==0);lrg(ii)=lrg(ii-1);
%ii=find(lat==0);lat(ii)=lat(ii-1);
ii=find(lon==0);lon(ii)=lon(ii-1);
ii=find(ir==0);ir(ii)=ir(ii-1);

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

%Learn to use a function.

% find close to integral minutes--doesn't .always. work right
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;
% Remove empty intervals, or median works on the wrong dimension. SPdeS
jk(find(diff(jk)<2)+1)=[];

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,:),1);
lrgcm(1)=median(lrgc(1:jk(1)-1,:),1);
lrgsm(1)=median(lrgs(1:jk(1)-1,:),1);

for iii=2:length(jk);
    if jk(iii)-1-jk(iii-1)>0 % if there are data in the interval
        st3(iii,1:16)=median(st2(jk(iii-1):jk(iii)-1,1:16),1);
        lrgcm(iii)=median(lrgc(jk(iii-1):jk(iii)-1),1);
        lrgsm(iii)=median(lrgs(jk(iii-1):jk(iii)-1),1);
    else
        st3(iii)=NaN;
        lrgcm(iii)=NaN;
        lrgsm(iii)=NaN;
    end
end;

gpstime=st3(:,1);     % you DON'T use the median of time.
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;

ii=find(lrg>500);lrg(ii)=NaN;
ii=find(abs(diff(lon))>1);lon(ii)=lon(ii-1);
ii=find(abs(diff(lat))>1);lat(ii)=lat(ii-1);
ii=find(solar>3000);solar(ii)=NaN;
ii=find(ir>3000);ir(ii)=NaN;
ii=find(tss>1000);tss(ii)=NaN;
ii=find(tsg>100);tsg(ii)=NaN;
ii=find(cog>500);cog(ii)=NaN;
ii=find(sog>100);sog(ii)=NaN;
ii=find(ta>100);ta(ii)=NaN;
ii=find(imrain>1000);imrain(ii)=NaN;

if plotit
    figure;plot(gpstim./24+ddd,cog,gpstim./24+ddd,lrg);title('SCS course and gyro');xlim([ddd ddd+1])
    figure;plot(gpstim./24+ddd,imu);title('SCS wind speed');xlim([ddd ddd+1])
    figure;plot(gpstim./24+ddd,imrain);title('SCS rain');xlim([ddd ddd+1])
    figure;plot(gpstim./24+ddd,imd);title('SCS wind direction');xlim([ddd ddd+1])
    figure;plot(gpstim./24+ddd,sog);title('SCS speed over ground');xlim([ddd ddd+1])
%     figure;plot(lon,lat,'.');title('SCS Lat vs Lon');
    figure;plot(gpstim./24+ddd,tsg,gpstim./24+ddd,ta);title('SCS Tsea   Tair');xlim([ddd ddd+1])
    figure;plot(gpstim./24+ddd,tss);title('SCS salinity');xlim([ddd ddd+1])
    figure;plot(gpstim./24+ddd,rh);title('SCS RH');xlim([ddd ddd+1])
    figure;plot(gpstim./24+ddd,solar,gpstim./24+ddd,ir);title('SCS solar & IR');xlim([ddd ddd+1])
end

lond=lonm;
latd=latm;
if plotit && true
    figure;
    map_epicz(lond,latd,[-90 -60 -25 10]); % maps track
    if prtit
        print('-depsc2','-tiff',[way_raw_images_flux 'track' cruise year '_' jd '.eps']);
        print('-djpeg',[way_raw_images_flux 'track' cruise year '_' jd '.jpg']);
    end
elseif plotit
    figure;
    plot(lonm,latm,'.');title('SCS lat vs lon');
    if prtit; print('-djpeg90 ',[way_images_flux '\RHB_track_' num2str(ddd) '.jpg']); end
end

jd_scs=gpstime./24+ddd;

ii=find(lrgm>400);lrgm(ii)=NaN;
ii=find(lrgm<0);lrgm(ii)=NaN;

% % Already plotted above...
% 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');

%Plots using HP time conversion jd_scs DEW 5/20/2007
% jd_scs not congruent with variables...
% 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')

np=length(sogm);

if saveit
    prt_jas_scs;
    
    % save scs data for comparisons
    [r,c] = size(st3);
    st3(:,c+1) = jd_scs;
    scs = st3;
    save H:\RHB\flux\Raw\scs scs
end;