function SP=data_avg(DATA,kdim,dim_lims,iiavg,iiview,iipr);
% Call SP=data_avg(DATA,kdim,dim_lims,iiavg,iiview,iipr);
%
% kdim		= dimension(s) over which to operate;
% dim_lims	= cell array of dimension limits to process: {ylims xlims zlims z2lims ...},
%       where ylims = ([] for all, [ymin,ymax] or [jy1;jy2; ...jyn]), etc;
%       Set dim_lims=[] to process all indices of all dimensions;
% iiavg     = 1=average, 2=integrate, 3=max, 4=median, 5=max-median, 6=X-spectra, -P=P'th percentile (0<-P<100)
%             (can be vector of same length as kdim);
% iiview	= 1 to call view5 on the output SP (default is 0);
% iipr      = 1 to print output labels;
%
%
% 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.
if nargin==1,
    str=sprintf('    %s   ',DATA.labels{1:DATA.ndim}); disp(['Axis labels: ',str]);
    str=['Axis limits: '];
    for jdim=1:DATA.ndim, str=[str,'   ',sprintf('[%1.2f %1.2f]',DATA.scales{jdim}([1,end]))]; end
    disp(str)
    str=sprintf('       %i       ',DATA.dims); disp(['Axis lengths: ',str]);
    kdim=    str2num(input(['                                Enter dimension numbers over which to collapse: '],'s'));
    iiavg=   str2num(input(['For each collapse dimension, enter type of operation (1=avg,2=int,3=max,4=med): '],'s'));
    if isempty(iiavg), iiavg=1; end
    dim_lims=str2num(input(['For all dims, enter limits in form "{[min,max] [min,max] ...}" ([] for all)   : '],'s'));
    iiview=  str2num(input(['                             Automatically view output? (1=yes,0=no,[]=yes)   : '],'s'));
    if isempty(iiview), iiview=1; end
else
	if nargin<3 dim_lims=[]; end
    if nargin<4, iiavg=1; end
    if nargin<5 | isempty(iiview), iiview=0; end
end
if nargin<6, iipr=1; end
if isempty(dim_lims), 
    for jdim=1:DATA.ndim, dim_lims{jdim}=[]; end
end
%
if isfield(DATA,'iidB'), iidB=DATA.iidB; else, iidB=1; end
%
if any(kdim<1) | any(kdim>DATA.ndim), 
    error(['kdim must be between 1 and ndim: ',int2str(kdim(:)',DATA.ndim)]); 
end
%
SP=structure_copy(DATA,{'data'});
% loaded_fields=fieldnames(DATA);
% for j=1:length(loaded_fields),
% 	if ~(strcmp(loaded_fields{j},{'data'})),
% 		expr=['SP.',loaded_fields{j},'=DATA.',loaded_fields{j},';']; eval(expr);
%     else,
%         SP.data=[];
%     end
% end
%
KD=cell(5,1); dims_out=zeros(1,5);
if length(dim_lims)~=DATA.ndim,
    error('INPUT ERROR: dim_lims must be a cell array with DATA.ndim'); 
end
for jdim=1:5,
    if jdim>DATA.ndim,
        KD{jdim}=1; dims_out(jdim)=1;
    else
		[KD{jdim},lims]=find_indices(DATA.scales{jdim},DATA.dims(jdim), ...
            dim_lims{jdim},DATA.labels{jdim});
        SP.scales{jdim}=DATA.scales{jdim}(KD{jdim});
        if isfield(DATA,'scales2'), SP.scales2{jdim}=DATA.scales2{jdim}(KD{jdim}); end
        dims_out(jdim)=length(SP.scales{jdim});
    end
end
if iidB==1,
    SP.data=zeros(dims_out); %%size(SP.data), size(DATA.data),
    for j=1:length(KD{1}),
    	SP.data(j,:,:,:,:)=10.^(DATA.data(KD{1}(j),KD{2},KD{3},KD{4},KD{5})/10);
    end
% 	SP.data=10.^(DATA.data(KD{1},KD{2},KD{3},KD{4},KD{5})/10);
elseif iidB==0,
%	SP.data=zmagsq(DATA.data(KD{1},KD{2},KD{3},KD{4},KD{5}));
%	SP.dlab=[SP.dlab,'=>dB']; iidB=1;
	SP.data=DATA.data(KD{1},KD{2},KD{3},KD{4},KD{5});
	iidB=0;
elseif iidB==-1,
    if all(iiavg==3) | any(iiavg==6) | all(iiavg==1.5),
        SP.data=DATA.data(KD{1},KD{2},KD{3},KD{4},KD{5});
    else,
		SP.data=zmagsq(DATA.data(KD{1},KD{2},KD{3},KD{4},KD{5}));
        SP.dlab=[SP.dlab,'=>dB']; iidB=1;
    end
end
dims_cur=size(SP.data);
% This for trailing ones in size(SP.data):
if length(dims_cur)<DATA.ndim, dims_cur=[dims_cur,ones(1,DATA.ndim-length(dims_cur))]; end
ndim_proc=length(kdim);
if length(iiavg)==1,
    iiavg=repmat(iiavg,[1,ndim_proc]);
elseif length(iiavg)~=ndim_proc,
    error('iiavg must be length 1 or same length as kdim');
end
%
if ~isfield(SP,'ind_max'),
    P.ind_max=repmat({[]},[1,length(dims_cur)]); 
elseif length(SP.ind_max)~=length(dims_cur),
    P.ind_max=repmat({[]},[1,length(dims_cur)]); 
end
for jdim_proc=1:ndim_proc,
    mdim=kdim(jdim_proc);
    lims=DATA.scales{mdim}(KD{mdim}([1,end])); iiav=iiavg(jdim_proc);
% Still need to do to set SP.ind_max:
%    if dims_cur(mdim)<=1,       % Do nothing if nothing to do;
%        if iipr, disp(['Doing nothing since dims(mdim)=1 for mdim = ',int2str(mdim)]); end 
    if iiav==1 | iiav==1.5,
        SP.labels{mdim}=str_make('Avg over',SP.labels{mdim},lims);
        SP.data=mean(SP.data,mdim);
    elseif iiav==2, 
    	delta=diff(DATA.scales{mdim}([1:2]));
        SP.labels{mdim}=str_make('Int over',SP.labels{mdim},lims);
        SP.data=sum(SP.data,mdim)*delta; 
    elseif iiav==3,
        SP.labels{mdim}=str_make('Max over',SP.labels{mdim},lims);
        [SP.data,SP.ind_max{mdim}]=max(SP.data,[],mdim);
        if isempty(SP.ind_max{mdim}), 'empty', end
%    elseif mod(iiav,10)==3 & floor(iiav/10)>=1 & floor(iiav/10)<=ndim,
%        SP.labels{mdim}=str_make('Col max over',SP.labels{mdim},lims);
%        [SP.data,SP.ind_max{mdim}]=max(SP.data,[],mdim);
    elseif iiav==4,
        SP.labels{mdim}=str_make('Med over',SP.labels{mdim},lims);
        SP.data=median(SP.data,mdim);
    elseif iiav==5,
        SP.labels{mdim}=str_make('Max-Med over',SP.labels{mdim},lims);
        [data,SP.ind_max{mdim}]=max(SP.data,[],mdim);
        SP.data=data./median(SP.data,mdim);
    elseif iiav==6,
        if length(iiavg)>1 & jdim_proc~=ndim_proc, error('X-spectra must be last operation'); end
        for jdim=SP.ndim:-1:mdim+1,
            SP.labels{jdim+1}=SP.labels{jdim};
            SP.scales{jdim+1}=SP.scales{jdim};
            if isfield(SP,'scales2'), SP.scales2{jdim+1}=SP.scales2{jdim}; SP.labels2{jdim+1}=SP.labels2{jdim}; end
        end
        SP.labels{mdim+1}=SP.labels{mdim}; SP.scales{mdim+1}=SP.scales{mdim};
        if isfield(SP,'scales2'), SP.labels2{mdim+1}=SP.labels2{mdim}; SP.scales2{mdim+1}=SP.scales2{mdim}; end
        if mdim>1, SP.data=permute(SP.data,[mdim,1:mdim-1,mdim+1:SP.ndim]); end                       % Permute so mdim 1st dim;
        dims_old=size(SP.data);
        if length(dims_old)>2, SP.data=reshape(SP.data,[dims_old(1),prod(dims_old(2:end))]); end      % Reshape to 2-D;
        [nr,nc]=size(SP.data); data2=zeros([nr,nr,nc]);
        for jcol=1:nc,
            data2(:,:,jcol)=SP.data(:,jcol) * SP.data(:,jcol)';
        end
        if length(dims_old)>2, data2=reshape(data2,[nr,nr,dims_old(2:end)]); end                % Reshape to original;
        if mdim>1,
            SP.data=permute(data2,[3:mdim+1,1,2,mdim+2:SP.ndim+1]);                              % Permute to original
        else    
            SP.data=data2;
        end
        SP.dims=size(SP.data); SP.ndim=length(SP.dims); SP.dlab=['X-Spectra of ',SP.dlab];
    elseif iiav>-99 & iiav<-1,
% Find percentile value:
        SP.labels{mdim}=str_make([int2str(-iiav),'%ile over '],SP.labels{mdim},lims);
        [data,SP.ind_max{mdim}]=sort(SP.data,mdim);
        pct_indx=round(-.01*iiav*(SP.dims(mdim)-1)) + 1;
        if mdim==1, SP.data=data(pct_indx,:,:,:,:,:,:);
        elseif mdim==2, SP.data=data(:,pct_indx,:,:,:,:,:);
        elseif mdim==3, SP.data=data(:,:,pct_indx,:,:,:,:);
        elseif mdim==4, SP.data=data(:,:,:,pct_indx,:,:,:);
        elseif mdim==5, SP.data=data(:,:,:,:,pct_indx,:,:);
        elseif mdim==6, SP.data=data(:,:,:,:,:,pct_indx,:);
        elseif mdim==7, SP.data=data(:,:,:,:,:,:,pct_indx);
        end
    elseif iiav>-199 & iiav<-101,
% Find percentile value:
        SP.labels{mdim}=str_make(['%ile over '],SP.labels{mdim},lims);
        [data,SP.ind_max{mdim}]=sort(SP.data,mdim);
        pct_ile1=-(iiav+100); pct_iles=[pct_ile1,50,100-pct_ile1];
        pct_indx=round(.01*pct_iles*(SP.dims(mdim)-1)) + 1;
        if mdim==1, SP.data=data(pct_indx,:,:,:,:,:,:);
        elseif mdim==2, SP.data=data(:,pct_indx,:,:,:,:,:);
        elseif mdim==3, SP.data=data(:,:,pct_indx,:,:,:,:);
        elseif mdim==4, SP.data=data(:,:,:,pct_indx,:,:,:);
        elseif mdim==5, SP.data=data(:,:,:,:,pct_indx,:,:);
        elseif mdim==6, SP.data=data(:,:,:,:,:,pct_indx,:);
        elseif mdim==7, SP.data=data(:,:,:,:,:,:,pct_indx);
        end
    elseif ~any(iiav==[1:5]),
		disp(['jdim_proc,mdim,iiav = ',int2str([jdim_proc,mdim,iiav])]);
		error(['Illegal iiavg in data_avg (1=average, 2=integrate, 3=max): ',int2str(iiav)]); %SP=[]; return
    end
    if dims_cur(mdim)>1 & iiav~=6,
        if iiav<-100,
            SP.scales{mdim}=pct_iles;
            if isfield(SP,'scales2'), SP.scales2{mdim}=SP.scales{mdim}; end
        else,
            SP.scales{mdim}=mean(lims);
            if isfield(SP,'scales2'), SP.scales2{mdim}=mean(DATA.scales2{mdim}(KD{mdim}([1,end]))); end
        end
% Reshape SP.data to have averaged dimension have length 1:
        dims_cur(mdim)=size(SP.data,mdim);
        SP.data=reshape(SP.data,dims_cur); %%size(SP.data), dims_cur,
    end
end
if iiav~=6, SP.dims=dims_cur(1:DATA.ndim); end
SP.iidB=iidB;
%
if iipr, 
    for mdim=1:ndim_proc,
        disp(SP.labels{kdim(mdim)}); 
    end
end
%
if iidB==1,
	SP.data=10*log10(SP.data);
	cmax=ceil(max(SP.data(:))/5)*5; cmin=max(floor(min(SP.data(:))/5)*5,cmax-60);
elseif iidB==0,
%	SP.data=sqrt(SP.data);
	cmax=max(SP.data(:)); cmin=min(SP.data(:));
elseif iidB==-1,
    cmax=max(dB(SP.data(:))); cmin=min(dB(SP.data(:)));
end
SP.clim=[cmin,cmax];
%
if iiview, view5(SP); end
%
%%%
%
function [iiz,zlims]=find_indices(z,nz,zlims,zlab)
%
[nr,nc]=size(zlims);
if nr==0,
    iiz=[1:nz];
elseif nr==1 & nc==2,
	dz=(z(nz) - z(1))/(max(1,nz-1));
	iiz=find(z>=zlims(1)-dz/2 & z<=zlims(2)+dz/2);
	if isempty(iiz), 
		zlims,z(1),z(nz), 
		error(['USER ERROR: no overlap in ',zlab])
	end
elseif nc==1,
    iiz=zlims;
%    iiz=zeros(length(zlims),1);
%    for jz=1:length(zlims);
%        [dum,iimin]=min(abs(z-zlims(jz))); iiz(jz)=iimin;
%    end
else
    disp(z),disp(nz),disp(zlims),disp(zlab),
	error('zlims must be [zmin,zmax] or [z1;z2; ... zn]');
end
%%z([1,end]),nz,length(z),zlims,length(iiz)
zlims=z(iiz([1,end])); zlims=zlims(:)';
%
function str=str_make(avg_lab,dim_lab,lims)
%
str=[avg_lab,' ',dim_lab,'(',num2str(lims(1)),',',num2str(lims(2)),')'];

