function SP=data_avg2(DATA,kdim,iiavg,parms,dim_lims,iiview,iipr);
% Call SP=data_avg2(DATA,kdim,iiavg,parms,dim_lims,iiview,iipr);
%
% kdim		= dimension(s) over which to operate;
% iiavg     = 1=average, 2=integrate, 3=max, 4=median, 5=var, -P=P'th percentile (0<-P<100)
%             (can be vector of same length as kdim);
% parms     = [npt_avg,npt_step,pct_cens] or [-scale_avg,-scale_step,pct_sens] where
%             npt_avg,npt_step          = the number of points to avg over and to step between output points;
%             -scale_avg,-scale_step    = time (or other scale) to avg over and to step between output points;
%             pct_cens                  = percentage of top values to censor beforehand (0 to include all);
% 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;
% iiview	= 1 to call view3 on the output SP (default is 0);
% iipr      = 1 to print output label;
%
%
% 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,
    disp(['Dimension labels: ']); disp(sprintf('   %s   ',DATA.labels{1:DATA.ndim}));
    disp(['Dimension limits: ']);
    str=[];
    for jdim=1:DATA.ndim,
        str=[str,'   ',sprintf('[%1.2f %1.2f]',DATA.scales{jdim}([1,end]))];
    end
    disp(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'));
    parms=   str2num(input(['For all dims, enter npt_avg,npt_step,pct_cens in form: "[na,ns;%; na,ns,% ...]": '],'s'));
    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, iiavg=1; end
    if nargin<4, parms=[1,1,0]; end
	if nargin<5 dim_lims=[]; end
    if nargin<6 | isempty(iiview), iiview=0; end
end
if isempty(dim_lims), 
    for jdim=1:DATA.ndim, dim_lims{jdim}=[]; end
end
ndim_proc=length(kdim);
%
[nr,nc]=size(parms);
if nc==2, parms=[parms zeros(nr,1)]; nc=3; end
if nc~=3 | nr~=ndim_proc,
    disp('Bad dimensions on parms: '); eval('help data_avg2'); disp(parms); return
end
if nargin<6, iipr=1; 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); LD=KD; 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(KD{jdim});
    end
    dims_cur(jdim)=length(KD{jdim});
    LD{jdim}=[1:dims_cur(jdim)];
end
if iidB==1,
    SP.data=zeros(dims_out); %%size(SP.data), size(DATA.data),
    for j=1:length(KD{1}),
    	SP.data(j,LD{2},LD{3},LD{4},LD{5})=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});
elseif iidB==-1,
    if all(iiavg==3),
        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
% 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
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
%
navg=zeros(1,ndim_proc); j1j2=zeros(3,ndim_proc,50);
SP.ind_max=repmat({[]},[1,length(dims_cur)]);
for jdim_proc=1:ndim_proc,
    mdim=kdim(jdim_proc); nax=dims_cur(mdim);
    lims=SP.scales{mdim}([1,end]);
    delta=abs(diff(lims)/max(1,nax-1));
    npt_avg(jdim_proc)=parms(jdim_proc,1);      % Number of points to average
    if npt_avg(jdim_proc)==0, 
        npt_avg(jdim_proc)=nax;
    elseif npt_avg(jdim_proc)<0, 
        npt_avg(jdim_proc)=max(1,round(-npt_avg(jdim_proc)/delta)); 
    end
    npt_step(jdim_proc)=parms(jdim_proc,2);     % Number of points to step
    if npt_step(jdim_proc)==0,
        npt_step(jdim_proc)=4*nax; 
    elseif npt_step(jdim_proc)<0, 
        npt_step(jdim_proc)=max(1,round(-npt_step(jdim_proc)/delta));
    end
    pct_cens=parms(jdim_proc,3);
    if pct_cens<0 | pct_cens>100,
        error(['Illegal pct_cens: ',pct_cens]);
    end
    npt_ovl=max(0,npt_avg(jdim_proc) - npt_step(jdim_proc));        % Number of overlapped points
% Compute number of averages to take:
    navg(jdim_proc)=round((nax-npt_ovl)/(npt_avg(jdim_proc) - npt_ovl));
    if navg(jdim_proc)<=0, 
        disp(['Setting navg=1 since navg calculated to be zero for dim ',int2str(jdim_proc)]); navg(jdim_proc)=1; 
    end
% Adjust npt_avg so that they fit perfectly:
    npt_avg(jdim_proc)=round((nax - npt_ovl)/navg(jdim_proc)) + npt_ovl;
%%    navg(jdim_proc)=max(1,round(nax/npt_step(jdim_proc)));
    j1=1;
    for javg=1:navg(jdim_proc),
        j2=min(nax,j1+npt_avg(jdim_proc));
        if j2==nax, j1=max(1,j2-npt_avg(jdim_proc)); end
        jmid=round((j1+j2)/2);
        j1j2(1:3,jdim_proc,javg)=[j1 j2 jmid];
        j1=j1 + npt_step(jdim_proc);
    end
%%    navg(jdim_proc)=1 + floor((nax - npt_avg(jdim_proc))/npt_step(jdim_proc));
    dims_out(mdim)=navg(jdim_proc);     % Reduce output dimension from input dimension;
end
%dims_out,navg,npt_avg,npt_step,
for jdim_proc=1:ndim_proc,
    mdim=kdim(jdim_proc); nax=DATA.dims(mdim);
    lims=SP.scales{mdim}([1,end]); iiav=iiavg(jdim_proc);
    delta=abs(diff(lims)/max(1,nax-1));
% 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
    jpt0=0;
    sc0=SP.scales{mdim}; SP.scales{mdim}=zeros(1,navg(jdim_proc));
    if isfield(SP,'scales2'),sc02=SP.scales2{mdim}; SP.scales2{mdim}=zeros(1,navg(jdim_proc)); end
    dims_next=dims_cur; dims_next(mdim)=navg(jdim_proc); data_out=zeros(dims_next);
    pct_cens=parms(jdim_proc,3); 
    for javg=1:navg(jdim_proc),
%        LD{mdim}=jpt0 + [1:npt_avg(jdim_proc)];
        LD{mdim}=[j1j2(1,jdim_proc,javg):j1j2(2,jdim_proc,javg)];
        data0=SP.data(LD{1},LD{2},LD{3},LD{4},LD{5});
        if pct_cens>0,
            npt_avg=round((1-.01*pct_cens)*size(data0,mdim));
            data0=sort(data0,mdim); 
            ND=LD; ND{mdim}=[1:npt_avg];        % Take first npt_avg points of dimension mdim;
            data0=data0(ND{1},ND{2},ND{3},ND{4},ND{5});
        end
        jmid=j1j2(3,jdim_proc,javg); SP.scales{mdim}(javg)=sc0(jmid);
        if isfield(SP,'scales2'), SP.scales2{mdim}(javg)=sc02(jmid);; end
        if iiav==1,
            if javg==1, SP.labels{mdim}=str_make('Avg over',SP.labels{mdim},lims); end
            data=mean(data0,mdim);
		elseif iiav==2, 
        	delta=diff(DATA.scales{mdim}([1:2]));
            if javg==1, SP.labels{mdim}=str_make('Int over',SP.labels{mdim},lims); end
            data=sum(data0,mdim)*delta; 
		elseif iiav==3,
            if javg==1, SP.labels{mdim}=str_make('Max over',SP.labels{mdim},lims); end
            [data,SP.ind_max{mdim}]=max(data0,[],mdim);
            if isempty(SP.ind_max{mdim}),
                'empty'
            end
        elseif iiav==4,
            if javg==1, SP.labels{mdim}=str_make('Med over',SP.labels{mdim},lims); end
            data=median(data0,mdim);
        elseif iiav==5,
            if javg==1, SP.labels{mdim}=str_make('Var over',SP.labels{mdim},lims); end
            data=std(data0,0,mdim);
%             data=std(sqrt(data0),0,mdim); data=data.^2;
% %         elseif iiav==5,
% %             if javg==1, SP.labels{mdim}=str_make('Max-Med over',SP.labels{mdim},lims); end
% %             [data2,SP.ind_max{mdim}]=max(data0,[],mdim);
% %             data=data2./median(data0,mdim);
        elseif iiav<0 & iiav>-100,
	% Find percentile value:
            if javg==1, SP.labels{mdim}=str_make([int2str(-iiav),'%ile over '],SP.labels{mdim},lims); end
            [data2,SP.ind_max{mdim}]=sort(data0,mdim);
            pct_indx=round(-.01*iiav*(SP.dims(mdim)-1)) + 1;
            if mdim==1, data=data2(pct_indx,:,:,:,:,:,:);
            elseif mdim==2, data=data2(:,pct_indx,:,:,:,:,:);
            elseif mdim==3, data=data2(:,:,pct_indx,:,:,:,:);
            elseif mdim==4, data=data2(:,:,:,pct_indx,:,:,:);
            elseif mdim==5, data=data2(:,:,:,:,pct_indx,:,:);
            elseif mdim==6, data=data2(:,:,:,:,:,pct_indx,:);
            elseif mdim==7, data=data2(:,:,:,:,:,:,pct_indx);
            end
        else
			disp(['jdim_proc,mdim,iiav = ',int2str([jdim_proc,mdim,iiav])]);
			error('Illegal iiavg in data_avg (1=average, 2=integrate, 3=max)'); %SP=[]; return
        end
        LD{mdim}=javg;
        data_out(LD{1},LD{2},LD{3},LD{4},LD{5})=data;
        jpt0=jpt0 + npt_step(jdim_proc);
    end
    LD{mdim}=[1:navg(jdim_proc)];
    SP.data=data_out;
    if dims_cur(mdim)>1,	
% Reshape SP.data to have averaged dimension have length 1:
        dims_cur(mdim)=navg(jdim_proc);
        SP.data=reshape(SP.data,dims_cur); %%size(SP.data), dims_cur,
    end
end
SP.dims=dims_cur(1:DATA.ndim);
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,
	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)),')'];

