function [DATA,status]=hdfsd_read(SD,dim_lims,jset)
%
% Call [DATA,status]=hdfsd_read(SD,dim_lims,jset)
%
% Reads subset of dataset number jset from HDF file opened with SD=hdfsd_open(hdf_file);
% Subset given by cell array dim_lims ([] ==> get whole dataset);
%
% SD		= structure from SD=hdfsd_open(filename);  Assumes this has been called.
% dim_lims	= {ylims,xlims,zlims,...}, an ndim by 1 cell array of index specifications, where
%	ylims 	= [] ==> all, OR [ymin,ymax], OR [ymin,ymax,+increment], OR [ymin,ymax,-|dy|], 
%			  OR [jy1;jy2; ... jyn], and similarly for xlims, zlims, etc;  Set dim_lims=[] to 
%			  get entire dataset.
% jset		= data set number for when HDF file contains multiple data sets (default=1);
%
% OUTPUTS
% DATA			= Structure containing data, dimensions, scales, and labels;
% DATA.title    = Name of HDF file from which data was read;
% DATA.data		= Data matrix;
% DATA.dlab     = Data label;
% DATA.ndim     = Number of dimensions;
% DATA.dims		= Array containing data matrix dimensions;
% DATA.scales	= Cell array containing axis values;
% DATA.labels	= Cell array containing axis labels in (1:ndim) and data label in (ndim+1);
% status	= 0 for success, negative for failure;
%
% HDF dataset is automatically closed at then end of the routine.
%
if nargin<3, jset=1; end
if jset<1 | jset>SD.nset, error(['Illegal jset for SD: ',int2str([jset,SD.nset])]); end
if isempty(dim_lims),
	for jrank=1:SD.rank{jset}, dim_lims{jrank}=[]; end
end
%
[nr,nc]=size(dim_lims);
%if nr~=1 | nc~=1, disp('Illegals dim_lims: should be rank by 1 cell array'); end
ndim=nr*nc;
if ndim<SD.rank{jset},
%	disp('Illegals dim_lims: must be rank by 1 cell array');
%	disp([SD.rank{jset},ndim]); status=-1; return;
    for jrank=ndim+1,SD.rank{jset}, dim_lims{jrank}=[]; end
    ndim=SD.rank{jset};
elseif ndim>SD.rank{jset},
%%    disp('dim_lims cell array has too many elements.  Ignoring extra dimensions.');
    ndim=SD.rank{jset};
end
start=zeros(1,ndim); stride=zeros(1,ndim); ntake=zeros(1,ndim); ndim_max=7;
ii_pick=zeros(1,ndim_max);
sd_dims=SD.dim_sizes{jset};
DATA.title=SD.hdf_file;
DATA.dlab=SD.data_name{jset}; DATA.data=[];
DATA.ndim=ndim; DATA.dims=zeros(1,ndim); DATA.labels=cell(1,ndim); DATA.scales=cell(1,ndim);
for jdim=1:ndim,
    ldim=sd_dims(jdim); axvals=SD.scales{jset}{jdim};
	lims=dim_lims{jdim}; lmin=min(lims); lmax=max(lims); [nr,nc]=size(lims);
    if nr==1 & nc==2 & all(lims==0), lims=[]; nr=0; nc=0; end       % Make 0,0 mean take all;
	if nr==1 & nc==2,
		ii=find(axvals>=lims(1) & axvals<=lims(2));
        if isempty(ii),
            axlim_str=num2str([axvals(1),axvals(end)]);
            disp(['Dim limits outside range for dim ',int2str(jdim),': lims = ',int2str(lims(:)'), ...
                '; Axis limits = ',axlim_str]); status=-1;
            return; 
        end
		start(jdim)=ii(1)-1; ntake(jdim)=length(ii); stride(jdim)=1;
	elseif isempty(lims) | (nr==1 & nc==1 & lims(1)==0),
		start(jdim)=0; ntake(jdim)=ldim; stride(jdim)=1;
	elseif nr==1 & nc==3,
		ii=find(axvals>=lims(1) & axvals<=lims(2));
		start(jdim)=ii(1)-1; ntake(jdim)=length(ii); 
		if lims(3)>0 & round(lims(3))==lims(3),
			stride(jdim)=lims(3);
		elseif lims(3)<0,
			delta=axvals(2)-axvals(1);
			stride(jdim)=nint(abs(lims(3)/delta));
		else
			disp(['Illegal dim_lims for jdim = ',jdim,'; skip factor must be integer.']);
			return
		end
	elseif nc==1
        if ~all(round(lims)==lims) | lmin<1 | lmax>ldim,
% For dim_lims in column vector of axis values, take closest:
            diff_mat=abs(repmat(axvals(:),[1,nr])-repmat(lims(:)',[ldim,1]));
    		[dum,ii]=min(diff_mat,[],1);
            dim_lims{jdim}=unique(ii); lims=dim_lims{jdim}; lmin=min(lims); lmax=max(lims); 
            [nr,nc]=size(lims);
        end
        ii=[lmin:lmax];
        if isempty(ii),
            axlim_str=num2str([axvals(1),axvals(end)]);
            disp(['Dim limits outside range for dim ',int2str(jdim),': lims = ',int2str(lims(:)'), ...
                '; Axis limits = ',axlim_str]); status=-1; return; 
        end            
		start(jdim)=ii(1)-1; ntake(jdim)=ii(end)-ii(1)+1; stride(jdim)=1;
		if ntake(jdim)~= nr, ii_pick(jdim)=1; end
%        keyboard
	else
		disp(['Illegal dim_lims for jdim = ',int2str(jdim)]); disp(lims); nr,nc, disp(dim_lims); return
	end
	ii=[start(jdim)+1:stride(jdim):start(jdim)+ntake(jdim)*stride(jdim)];
	DATA.labels{jdim}=SD.dim_names{jset}{jdim};
%%    jdim,ii,SD.scales{jset}, keyboard
	DATA.scales{jdim}=axvals(ii);
	DATA.dims(jdim)=length(DATA.scales{jdim});
end
%
% We must reverse start,stride,ntake to be compatible with HDF's assumption of (page,col,row) 
% rather than Matlab's assumption of (row,col,page)
% NO, WE DON'T HAVE TO DO THIS AFTER ALL.
start_rev=start(ndim:-1:1); stride=stride(ndim:-1:1); ntake=ntake(ndim:-1:1); 
% start_rev(1:ndim) is zero-based index of starting value for each dimension;
% stride(1:ndim) is index increment for each dimension;
% ntake(1:ndim) is total number of indices to be read in each dimension;
sd_id = hdfsd('start',SD.hdf_file,'read');
if sd_id<0, disp('Unable to open HDF file for reading'); return; end
%
[data,status] = hdfsd('readdata',SD.data_sds_id{jset},start_rev,stride,ntake);
if status<0,
	disp('status<0 in readdata'); keyboard; return;
end
DATA.data=double(data);
status = hdfsd('end',sd_id);
if status<0,
	disp('status<0 upon hdfsd(''end'') '); keyboard; return;
end
%
% Check if data units are dB:
DATA.iidB=0;
if ~isempty(findstr('dB',DATA.dlab)), DATA.iidB=1; end
%
%%return
%
if any(ii_pick),
	ii=cell(ndim_max,1);
	for jdim=1:ndim
		if ii_pick(jdim),
			ii{jdim}=dim_lims{jdim}(:) - start(jdim);
	%        nval=length(dim_lims{jdim}(:)); ii{jdim}=zeros(nval,1); jwarn=0;
	%        for jval=1:nval, 
	%            [dum,ii{jdim}(jval)]=min(abs(DATA.scales{jdim}-dim_lims{jdim}(jval))); 
	%            if dum~=0, jwarn=jwarn + 1; end
	%        end
	%        if jwarn>0, 
	%            disp(['WARNING: one or more values not found exactly for dim ',int2str(jdim)]); 
	%        end
            DATA.scales{jdim}=DATA.scales{jdim}(ii{jdim});
            DATA.dims(jdim)=length(ii{jdim});
		else
			ii{jdim}=[1:DATA.dims(jdim)];
		end
	end
	for jdim=ndim+1:ndim_max,
		ii{jdim}=1;
	end
	%
	DATA.data=DATA.data(ii{1},ii{2},ii{3},ii{4},ii{5},ii{6},ii{7});
end
%
% Set limits for view5:
cmin=min(DATA.data,[],1); cmax=max(DATA.data,[],1);
for jdim=2:ndim,
    cmin=min(cmin,[],jdim); cmax=max(cmax,[],jdim);
end
DATA.clim=[cmin,cmax];
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
