function A=readL2GP(file, swath, READATT, GETFIELDS)
%
%function A=readL2GP(FILE, SWATH, READATT, GETFIELDS))
%function A=readL2GP(DIRECTORY, SWATH, READATT, GETFIELDS))
%i
% eg
%FILE  = '/data/emls/l2gp/v02.20/2005/028/MLS-Aura_L2GP-O3_v02-20-c06_2005d243.he5'
%SWATH = 'O3'
%  or
%FILE  = '/data/emls/l2gp/v02.20/2005/028/MLS-Aura_L2GP-DGG_v02-020-c06_2005d243.he5'
%SWATH = 'O3-640'   %(a diagnostic product stored in the DGG file)
%  or
%DIRECTORY = '/data/emls/l2gp/v02.20/2005/028/'
%SWATH     = 'O3-640'  % (a diagnostic product stored in the DGG file)
%
%If the first argument is a directory, the program attempts to guess the
%correct file, first looking for files ending with '.he5' that have SWATH in their
%names and, if none are found, then for those that have 'DGG' in their
%names.  In either case, if multiple files are found, the last (alphabetically) is
%used.  This selects the latest "cycle" (eg 'c06', above) if level-2 has been run repeatedly.
%
%READATT = true [false]
%  true causes attributes to be returned in a
%  substructure, "A.att.(attribute_name)"
%
%the structure field 'L2gpValue' is the product and 'L2gpPrecision' is the precision
%
%
%omitted or empty swath brings up list of swaths in the specified file
%
%Note: some fields are single or int32 and may need to be caste to double
% (eg   double(A.Status) ) for compatibility with certain matlab operations
%
%
%for version 7.4.0.336 (R2007a) compatibility, added
% 'V71Dimensions', true
%to hdf5read
%
%
% Michael J. Schwartz  04-Oct-2006
%             updated  03-Aug-2016
%
%Copyright 2006, 2016 by the California Institute of Technology. ALL RIGHTS RESERVED.
%United States Government Sponsorship acknowledged.
%Any commercial use must be negotiated with the Office of Technology Transfer
%at the California Institute of Technology.
%
%This software is subject to U.S. export control laws and regulations and has been
%classified as 4D994.  By accepting this software, the user agrees to comply with all
%applicable U.S. export laws and regulations.  User has the responsibility to obtain
%export licenses, or other export authority as may be required before exporting such
%information to foreign countries or providing access to foreign persons.

if exist('GETFIELDS')~=1, GETFIELDS=false; end

if exist('READATT') ~=1, READATT=false; end

if (~exist('swath', 'var') || isempty(swath) || strcmp(swath, '*')) && exist(file, 'file')==2,
    fi=h5info(file, '/HDFEOS/SWATHS/');
    swaths={fi.Groups.Name}';
    for j=1:numel(swaths),
        r=find(swaths{j}=='/');
        swaths{j}=swaths{j}((r(3)+1):end);
    end
    colvec(catcellstr(sprintfcell('%2d ', 1:numel(swaths)), swaths(:)))
    result=input('enter number of desired swath or return to exit:  ');

    if ~isempty(result) && result>0 && result<=numel(swaths),
        swath=swaths{result};
    else
        return
    end
end
    
    

%if iscell(file), file=[file{:}]; end
if iscell(file), file=file{1}; end
if isstruct(file), 
    if iscell(file.l2gpdir), file=file.l2gpdir{1}; else file=file.l2gpdir; end
end
if ~strncmp(swath, 'WMOTPPressure-StdProd', 21) && strncmp(swath, 'WMOTPPressure', 13)
    swath_='Temperature'; 
else
    swath_=swath;
end
    
if exist(file, 'dir'),
    Path=file;
    foo=dir(file);
    r=zeros(numel(foo));
    for i=1:length(foo);
        [a, Name, Ext]=fileparts(foo(i).name);
        if strcmpi(Ext, '.he5'), r(i)=1; end
        if (length(foo(i).name)>9 && foo(i).name(end-8)=='t' &&  ~(file(end-8)== 't')), r(i)=0; end %funky ones I want to avoid
    end
    foo=foo(find(r)); %only files ending .he5  and not some funky ones that I want to avoid
    if length(foo)==1,
        file=foo.name;
    else
        a=strfind({foo.name}, swath_);
        for i=1:length(a), aa(i)=~isempty(a{i}); end
        aa=find(aa);
        if isempty(aa),
            a=strfind({foo.name}, 'DGG');
            for i=1:length(a), aa(i)=~isempty(a{i}); end
            aa=find(aa);
            if ~isempty(aa),
                file=foo(aa(end)).name;
            else
                a=strfind({foo.name}, 'DMP');
                for i=1:length(a), aa(i)=~isempty(a{i}); end
                aa=find(aa);
                if ~isempty(aa),
                    file=foo(aa(end)).name;
                else   
                    warning(['file ' file ' not found.']);
                    return;
                end
            end
        else
            file=foo(aa(end)).name;
        end
    end;
elseif exist(file, 'file')
    [Path, file, ext]=fileparts(file);
    file=[file ext];
else
    error('''file'' is not a directory or a file');
end
if Path(end)~='/',
    Path(end+1)='/';
end


b=file;
for i=1:3,
    [a,b]=strtok(b, '_');
end

r=find(a=='-');
try
    A.version=a(1:(r(end)-1));  %parse 'version' from the file name
catch
    A.version='';
end
try
    A.l2cycle=a((r(end)+1):end);  %parse 'l2cycle' from the file name
catch
    A.l2cycle='';
end
[a,b]=strtok(b, '_');
A.date=strtok(a, '.');



file=fullfile(Path, file);
if strcmp(file(end+(-3:0)), '.he5'),
    
    A.file=file;
    A.swath=swath;
    
    if iscell(GETFIELDS) | (GETFIELDS)
        %fi=hdf5info(file, 'ReadAttributes', false);
%         for i=1:length(fi.GroupHierarchy.Groups(1).Groups(2).Groups),
%             SWATH_{i}=fi.GroupHierarchy.Groups(1).Groups(2).Groups(i).Name;
%             SWATH{i}=strtokn(SWATH_{i}, 'last', '/');
%         end
         SWATH=getSwathNames(file);
        N=find(strcmp(SWATH, swath));
        %       for i=1:length(fi.GroupHierarchy.Groups(1).Groups(2).Datasets),
        %           GEO_{i}=fi.GroupHierarchy.Groups(1).Groups(2).Datasets(i).Name;
        %           GEO{i}=strtokn(GEO_{i}, 'last', '/');
        %       end
        
        N=find(strcmp(SWATH, swath), 1);
        if isempty(N), error(['Can not find swath ' swath ' in ' file '.']); end
        
        if iscellstr(GETFIELDS), 
            datafields=GETFIELDS;
        else
            datafields=    getGroupNames(file, ['/HDFEOS/SWATHS/' swath '/Data Fields']);
        end
        for i=1:numel(datafields), datafields_{i}=    ['/HDFEOS/SWATHS/' swath '/Data Fields/' datafields{i}]; end
        
        geofields=    getGroupNames(file, ['/HDFEOS/SWATHS/' swath '/Geolocation Fields']);
         for i=1:numel(geofields), geofields_{i}=    ['/HDFEOS/SWATHS/' swath '/Geolocation Fields/' geofields{i}]; end
      % geofields_=    {fi.GroupHierarchy.Groups(1).Groups(2).Groups(N).Groups(2).Datasets.Name};
        
        for i=1:length(datafields_),
            datafields{i}=strrep(strrep(datafields{i}, ' ', '_'), '.', 'p');
            %if isempty(field), continue; end
            if ~isstrprop(datafields{i}, 'alpha'), datafields{i}=['A' datafields{i}]; end
            A.(fixfieldname(datafields{i}))=hdf5read(A.file, datafields_{i}, 'V71Dimensions', true);
        end
        for i=1:length(geofields_),
            geofields{i}=strrep(strrep(geofields{i}, ' ', '_'), '.', 'p');
            %field=strrep(strrep(strtokn(geofields_{i}, 'last', '/'), ' ', '_'), '.', 'p');
            % if isempty(field), continue; end
            if ~isstrprop(geofields{i}, 'alpha'), geofields{i}=['A' geofields{i}]; end
           
            A.(fixfieldname(geofields{i}))=hdf5read(A.file, geofields_{i}, 'V71Dimensions', true);
        end
        
    else
 %       datafields={'L2gpValue' 'L2gpPrecision' 'Quality' 'Status', 'Convergence'};
        datafields={'L2gpValue' 'L2gpPrecision' 'Quality' 'Status', 'Convergence'};
        geofields= {'ChunkNumber' 'Latitude' 'LineOfSightAngle' 'LocalSolarTime' 'Longitude' 'OrbitGeodeticAngle' 'Pressure' 'SolarZenithAngle' 'Time'};
        
        for i=1:length(datafields),
            try
                A.(datafields{i})=hdf5read(A.file, ['/HDFEOS/SWATHS/' swath '/Data Fields/' datafields{i}], 'V71Dimensions', true);
            catch
                A.(datafields{i})=[];
            end
        end
        for i=1:length(geofields),
            try
                A.(geofields{i})=hdf5read(A.file, ['/HDFEOS/SWATHS/' swath '/Geolocation Fields/' geofields{i}], 'V71Dimensions', true);
            catch
                A.(geofields{i})=[];
            end
        end
        if READATT, fi=hdf5info(file, 'ReadAttributes', false); end
    end
    
    
    
    
    if READATT,
        fi=hdf5info(file, 'ReadAttributes', false);
        Att=fi.GroupHierarchy.Groups(1).Groups(1).Groups.Attributes;
        r=find(Att(1).Name=='/');
        for i=1:length(Att),
            name=strrep(Att(i).Name((r(end)+1):end), ' ', '_');
            name=strrep(name , '(', '_');
            name=strrep(name , ')', '_');
            %    A.(name)=hdf5read(file, Att(i).Name);
            %    if strcmp(class(A.(name)), 'hdf5.h5string'), A.(name)=A.(name).Data; end
            lastslash=find(Att(i).Name=='/',1,'last');
            if ~isempty(lastslash)
            A.att.(name)=hdf5read(file, Att(i).Name(1:lastslash), Att(i).Name((lastslash+1):end));
            else
            A.att.(name)=hdf5read(file, Att(i).Name);
            end 
            if strcmp(class(A.att.(name)), 'hdf5.h5string'), A.att.(name)=A.att.(name).Data; end
        end
    else
        A.att=struct([]);  %included so that structure arrays may be created including elements with and without
    end
    
    if strncmp(A.swath, 'WMOTPPressure', 13)
        A.L2gpValue=A.L2gpValue';
        A.L2gpPrecision=A.L2gpPrecision';
        if strfind(A.swath, 'MLS')
            T=readL2GP(strrep(A.file, 'DGG', 'Temperature'), 'Temperature');
            A.Status=T.Status;
            A.Convergence=T.Convergence;
            A.Quality=T.Quality;
        end
    end
    
    
    
    
else
    warning([file '  does not end .he5'])
end

%if isfield(A, 'PCF1'), A=rmfield(A, 'PCF1'); end
%if isfield(A, 'PCF2'), A=rmfield(A, 'PCF2'); end

% function datafields=fixfieldname(datafields)
% if ischar(datafields),
%     notcell=1;
%     datafields={datafields};
% end
% for i=1:numel(datafields), 
%     if ~ischar(datafields{i}), error('not a character string'); end
%     if isempty(datafields{i}), error('empty'); end
%     if ~((datafields{i}(1)>=97 && datafields{i}(1)<=122) || (datafields{i}(1)>=65 && datafields{i}(1)<=90)),
%         datafields{i}=['X' datafields{i}];
%     end
% end
% if exist('notcell', 'var'), datafields=datafields{1}; end

 
function S=fixfieldname(S, delim)
%function S=fixfieldname(S, delim)
%  convert string into valid field name-  only alpha numeric or '_'.  Leading char is alpha.
if ~exist('delim', 'var'), delim=','; end
ISCHAR=0;
if ischar(S),
    ISCHAR=1;
    foo=S;
    clear('S');
    i=0;
    while numel(foo)>=1,
        i=i+1;
        [S{i},foo]=strtok(foo, delim);
        while(S{i}(1)==' ' && numel(S{i})>1), S{i}=S{i}(2:end); end
    end
    
end
if ~iscell(S); error('Input must be a string or cell array of strings.'); end


for j=1:length(S),
    
%S=S(((abs(S)>='0' & abs(S)<='9') | (abs(S)>='a' & abs(S)<='z') | (abs(S)>='0' & abs(S)<='Z')));
S{j}(~(((abs(S{j})>='0' & abs(S{j})<='9') | (abs(S{j})>='a' & abs(S{j})<='z') | (abs(S{j})>='A' & abs(S{j})<='Z'))))='_';


N=length(S{j});
for i=1:N,
    if upper(S{j}(1))<'A' || upper(S{j}(1))>'Z',
        S{j}=['x_' S{j}];
        %S{j}=S{j}(2:end); 
    else break; 
    end
       
end
%if i==N, warning('TruncToZero', 'string truncated to zero length');end
end

if ISCHAR && numel(S)==1, S=S{1}; end
%if length(S)==1, S=S{1}; end

function x=colvec(x)
x=x(:);


function C=catcellstr(varargin)
%function C=catcellstr(varargin)
N=0;
for i=1:nargin,
  if ischar(varargin{i}), varargin{i}={varargin{i}}; end
  if ~iscellstr(varargin{i}), error('inputs are not cell strings or char'); end
  N(i)=numel(varargin{i});
end



for i=find(N==1)
    varargin{i}=varargin{i}(ones(max(N),1));
end
C=varargin{1};
for i=1:max(N),
  for j=2:nargin,
    C{i}=cat(2, C{i}, varargin{j}{i});
  end
end


function C=sprintfcell(S, varargin)
%function C=sprintfcell(S, varargin)
%like sprintf but vectorized so that there must be one
%input for each token in the format string and they may
%be vectors, scalars, or cell arrays for string inputs
for i=1:length(varargin),
L(i)=length(varargin{i}(:));
end

N=max(L);
r=find(L ~= max(L));
if any(L(r) ~= 1), error('inputs must be same length or scalar'); end


for i=1:N,
    for j=1:length(varargin),
       args{j}=varargin{j}(1 + (i-1)*(L(j)~=1));
       if iscell(args{j}), args{j}=args{j}{1}; end
    end
    C{i}=sprintf(S, args{:});
end

function [S,R]=strtokn(R, N, D)
%function [S,R]=strtokn(R, N, D)
%
% N may be string 'last' for last or -1 (negative numbers count from end)
%STRTOKN Nth token in the supplied string.
%        STRTOK(S,D) returns the Nth token delimited by the specified
%        delimiters D in the supplied string S.  It also optionally 
%        returns the remainder of the original string.  D is an array
%        of character delimiters. If no delimiters are supplied, 
%        "white space" is assumed.
% 
%        See also strtok, ISSPACE.
if isnumeric(N) && N<0,
    rev=1;
    R=reverse(R);
    N=abs(N);
else
    rev=0;
end

if ~exist('D', 'var'), D = char([9:13 32]); end% White space characters
if nargin==1,
  for i=1:N,
    [S,R]=strtok(R);
  end;
  return;
elseif strncmpi(N, 'last', 1)
 
   
   S=reverse(strtok(reverse(R), D));
   R='';

elseif strncmpi(N, 'all', 1)
  i=1;  
  while ~isempty(R),  
    [S{i},R]=strtok(R, D);
    i=i+1;
  end;
   
elseif ischar(N),
    error('bad string input N');
else
  for i=1:N,
    [S,R]=strtok(R, D);
  end;
end;
if isempty(R) && isempty(S), error('Ran out of string.'); end;

if rev, S=reverse(S); end



function S=ls_(varargin)
%function S=ls_(regexp, n)
%function S=ls_(-lsargs, regexp, strtoknargs)
%returns a single string that can be used as an input
%by default, it is the last token returned matching the wildcard
% 
lsargs=0;
for i=1:numel(varargin)
    if strncmp(varargin{i}, '-', 1), lsargs=i; else break; end
end
if numel(varargin)<lsargs+1, varargin{lsargs+1}=''; end

if numel(varargin)<lsargs+2, varargin{lsargs+2}='last'; end
    
try
    if isempty(varargin{lsargs+1}),
    S=strtokn(ls(varargin{1:lsargs}), varargin{lsargs+2:end});
    else
    S=strtokn(ls(varargin{1:(lsargs+1)}), varargin{lsargs+2:end});
    end
catch
    S=''; 
end

