function out = read_superdat(filename,varargin)
% Reads an SMS Super ASCII data (*.dat) file
%
% written by Alex Sanchez, USACE

out.time = []; 
out.data = [];
out.tjulian = [];

%Set Defaults
tread = [];
indread = [];
ncol = 1;

%Get properties
for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'times'
            tread = varargin{i+1};
        case 'index'
            indread = varargin{i+1};
        otherwise
            warning(['Unidentified Property: ',varargin{i}])
    end
end

% filename = '..\simulations\Flow_Scalar1D_sal.dat';
fid = fopen(filename,'r');
for i=1:1000000
    tline = fgetl(fid);
    nline = length(tline);
    if ~ischar(tline), break, end
    if strcmp(tline(1:min(nline,9)),'RT_JULIAN')
        out.tjulian = sscanf(tline,'%*s %f');
    end
    if strcmp(tline(1:min(nline,2)),'NC')
        nlen = sscanf(tline,'%*s %d');
        if isempty(indread)
           indread=1:nlen; 
        end
    end
    if strcmp(tline(1:min(nline,6)),'BEGSCL')
        ncol = 1;
    end
    if strcmp(tline(1:min(nline,9)),'TIMEUNITS')
        out.tunits = sscanf(tline,'%*s %s');
    end    
    if strcmp(tline(1:min(nline,4)),'NAME')
        out.name = sscanf(tline,'%*s %s');
    end  
    if strcmp(tline(1:min(nline,6)),'BEGVEC')
        ncol = 2;
    end
    if(strcmp(tline(1:min(nline,2)),'TS'))
        try
            ts = sscanf(tline,'%*s %*d %f');
        catch
            warning(['Problem reading file record: ',num2str(i)])
            break
        end
        if isempty(out.tjulian)
            disp(['Reading Time: ',num2str(ts),' ',tunits])
        else
            ts = datenum(gregorian(out.tjulian+ts/86400));
            disp(['Reading Time: ',datestr(ts,0)])
        end
        dat = [];
        if isempty(tread)
            out.time(end+1,1) = ts;
            try
                dat = fscanf(fid,'%f',[ncol,nlen]);
                if nlen~=size(dat,2)
                    warning(['Incomplete record: ',num2str(i)])
                    out.time(end)=[];
                    break
                end
            catch
                warning(['Problem reading file record: ',num2str(i)])
                break
            end
        else
            if any(ts==tread);
                out.time(end+1,1) = ts;
                try
                    dat = fscanf(fid,'%f',[ncol,nlen]);
                    if nlen~=size(dat,2)
                        warning(['Incomplete record: ',num2str(i)])
                        out.time(end)=[];
                        break
                    end
                catch
                    warning(['Problem reading file record: ',num2str(i)])
                    break
                end
            end
        end
        if ~isempty(dat)   
            out.data(end+1,:,1:ncol)=permute(dat(:,indread),[3,2,1]);
        end
    end
end
fclose(fid);
return
%**************************************************************************
function gcal=juliandate2gregorian(jd)    
% Converts a Julian Date to Gregorian

% Description:  
%  Converts decimal Julian days to Gregorian dates using the
%   astronomical convension with time zero starting at noon.  
%
% Input:   
%    jd  = decimal Julian days 
%
% Output:  gcal = six column Gregorian time [yyyy mmm dd HH MM SS]
jd=jd+0.5;
secs=rem(jd,1)*24*3600;
k=4*(floor(jd)-1721119)-1;
y=floor(k/146097);
j=k-146097*y;
k=4*floor(j/4)+3;
j=floor(k/1461);
d=floor(((k-1461*j)+4)/4);
k=5*d-3;
m=floor(k/153);
day=floor(((k-153*m)+5)/5);
y=y*100+j;
month=m-9;
year=y+1;
i=(m<10);
month(i)=m(i)+3;
year(i)=y(i);
hour=floor(secs/3600);
minute=floor(rem(secs,3600)/60);
sec=round(rem(secs,60));
gcal=[year(:),month(:),day(:),hour(:),minute(:),sec(:)];
return
