function SP=cube_avg(LG,kkavg,kkint,kkmax,ylims,xlims,zlims,iiview,iipr);
% Call SP=cube_avg(LG,kkavg,kkint,kkmax,ylims,xlims,zlims,iiview,iipr);
%
% kkavg		= list of dimensions over which to average;
% kkint		= list of dimenstions over which to integrate;
% kkmax		= list of dimensions over which to max;
% ylims		= y-axis interval [ymin,ymax] or y-axis indices [jy1;jy2; ... jyn] ([] for all);
% xlims		= x-axis interval [xmin,xmax] or x-axis indices [jx1;jx2; ... jxn] ([] for all);
% zlims		= z-axis interval [zmin,zmax] or z-axis indices [jz1;jz2; ... jzn] ([] for all);
%	NOTE: for LG.data with dimension>3, zlims=[zmin zmax; z2min z2max; z3min z3max];
% iiview	= 1 to call view3 on the output SP (default is 0);
%
% Use data_avg to average data structures in DATA.scales{xax yax zax ...} format.
%
[ny,nx,nz,nz2,nz3]=size(LG.data);
if ~isfield(LG,'z2'),LG.z2=0; end
if ~isfield(LG,'z3'),LG.z3=0; end
if nargin<2, kkavg=[1 2 3]; end
if nargin<3, kkint=[]; end
if nargin<4, kkmax=[]; end
if nargin<5 | isempty(ylims), ylims=[LG.y(1) LG.y(ny)]; end
if nargin<6 | isempty(xlims), xlims=[LG.x(1) LG.x(nx)]; end
if nargin<7 | isempty(zlims), 
	zlims=[LG.z(1) LG.z(nz)];
	z2lims=[LG.z2(1) LG.z2(nz2)];
	z3lims=[LG.z3(1) LG.z3(nz3)]; 
end
if nargin<8 | isempty(iiview), iiview=0; end
if nargin<9, iipr=1; end
%
if isfield(LG,'iidB'), iidB=LG.iidB; else, iidB=1; end
%
[nr,nc]=size(zlims);
if nr==3,
	z3lims=zlims(3,:); z2lims=zlims(2,:); zlims=zlims(1,:);
elseif nr==2,
	if nz3>1, z3lims=[LG.z3(1) LG.z3(nz3)]; else, z3lims=[1 1]; end
	z2lims=zlims(2,:); zlims=zlims(1,:);
elseif nr==1,
	if nz2>1, z2lims=[LG.z2(1) LG.z2(nz2)]; else, z2lims=[1 1]; end
	if nz3>1, z3lims=[LG.z3(1) LG.z3(nz3)]; else, z3lims=[1 1]; end
end
%
if any(kkavg<1 | kkavg>5) | any(kkint<1 | kkint>5) | any(kkmax<1 | kkmax>5),
	error('USER ERROR: all elements of kkavg,kkint,kkmax must be between 1 and 5');
end
%
[iiy,ylims]=find_indices(LG.y,ny,ylims,'y-axis');
[iix,xlims]=find_indices(LG.x,nx,xlims,'x-axis');
[iiz,zlims]=find_indices(LG.z,nz,zlims,'z-axis');
if nz2>1, [iiz2,z2lims]=find_indices(LG.z2,nz2,z2lims,'z2-axis'); else, iiz2=[1]; end
if nz3>1, [iiz3,z3lims]=find_indices(LG.z3,nz3,z3lims,'z3-axis'); else, iiz3=[1]; end
%
nxout=length(iix); nyout=length(iiy); 
nzout=length(iiz); nz2out=length(iiz2); nz3out=length(iiz3);
%
loaded_fields=fieldnames(LG);
for j=1:length(loaded_fields),
	if ~any(strcmp(loaded_fields{j},{'data'})),
		expr=['SP.',loaded_fields{j},'=LG.',loaded_fields{j},';']; eval(expr);
	end
end
%%SP.fft_root=LG.fft_root;
%%if isfield(LG,'title'), SP.title=LG.title; else, SP.title=SP.fft_root; end
SP.y=LG.y(iiy); ylims=SP.y([1,end]); ylims=ylims(:)';
SP.x=LG.x(iix); xlims=SP.x([1,end]); xlims=xlims(:)';
SP.z=LG.z(iiz); zlims=SP.z([1,end]); zlims=zlims(:)';
if isfield(LG,'z2'),
	SP.z2=LG.z2(iiz2); z2lims=SP.z2([1,end]); z2lims=z2lims(:)';
end
if isfield(LG,'z3');
	SP.z3=LG.z3(iiz3); z3lims=SP.z3([1,end]); z3lims=z3lims(:)'; 
end
%
if ny>1, dy=LG.y(2)-LG.y(1); else, dy=1; end
if nx>1, dx=LG.x(2)-LG.x(1); else, dx=1; end
if nz>1, dz=LG.z(2)-LG.z(1); else, dz=1; end
if nz2>1, dz2=LG.z2(2)-LG.z2(1); else, dz2=1; end
if nz3>1, dz3=LG.z3(2)-LG.z3(1); else, dz3=1; end
%
if iidB==1,
	data=10.^(LG.data(iiy,iix,iiz,iiz2,iiz3)/10);
else
	data=zmagsq(LG.data(iiy,iix,iiz,iiz2,iiz3));
	SP.iidB=1;
	SP.dlab=[SP.dlab,'=>dB'];
end
SP.iidB=iidB;
for javg=1:length(kkavg),
	data=mean(data,kkavg(javg));
	switch kkavg(javg)
	case 1
		SP.ylab=['Avg over ',LG.ylab,' (',num2str(ylims),')']; plab=SP.ylab;
		SP.y=mean(LG.y(iiy));
	case 2,
		SP.xlab=['Avg over ',LG.xlab,' (',num2str(xlims),')']; plab=SP.xlab;;
		SP.x=mean(LG.x(iix));
%%		size(data), size(LG.data), size(LG.data(iiy,iix,iiz,iiz2,iiz3)), keyboard
	case 3
		SP.zlab=['Avg over ',LG.zlab,' (',num2str(zlims),')']; plab=SP.zlab;;
		SP.z=mean(LG.z(iiz));
	case 4
		SP.z2lab=['Avg over ',LG.z2lab,' (',num2str(z2lims),')']; plab=SP.z2lab;;
		SP.z2=mean(LG.z2(iiz2));
	case 5
		SP.z3lab=['Avg over ',LG.z3lab,' (',num2str(z3lims),')']; plab=SP.z3lab;;
		SP.z3=mean(LG.z3(iiz3));
	end
	if iipr, disp(plab); end
end
for jint=1:length(kkint),
	data=sum(data,kkint(jint));
	switch kkint(jint)
	case 1
		SP.ylab=['Integ over ',LG.ylab,' (',num2str(ylims),')']; plab=SP.ylab;
		SP.y=mean(LG.y(iiy));
		data=data*dy;
	case 2,
		SP.xlab=['Integ over ',LG.xlab,' (',num2str(xlims),')']; plab=SP.xlab;
		SP.x=mean(LG.x(iix));
		data=data*dx;
	case 3
		SP.zlab=['Integ over ',LG.zlab,' (',num2str(zlims),')']; plab=SP.zlab;
		SP.z=mean(LG.z(iiz));
		data=data*dz;
	case 4
		SP.z2lab=['Integ over ',LG.z2lab,' (',num2str(z2lims),')']; plab=SP.z2lab;
		SP.z2=mean(LG.z2(iiz2));
		data=data*dz2;
	case 5
		SP.z3lab=['Integ over ',LG.z3lab,' (',num2str(z3lims),')']; plab=SP.z3lab;
		SP.z3=mean(LG.z3(iiz3));
		data=data*dz3;
	end
	if iipr, disp(plab); end
end
for jmax=1:length(kkmax),
%%	size(data),
	[data,jjmax]=max(data,[],kkmax(jmax));
%%	size(data),
	switch kkmax(jmax)
	case 1
		SP.ylab=sprintf(['Max over ',LG.ylab,' (%4.2g,%4.2g)'],ylims); plab=SP.ylab;
		SP.y=mean(LG.y(iiy));
%		size(iiy),size(iiy(jjmax)),[nxout,nzout,nz2out,nz3out], keyboard
		SP.iymax=reshape(iiy(jjmax),[nxout,nzout,nz2out,nz3out]);
		SP.ymax=reshape(LG.y(iiy(jjmax)),[nxout,nzout,nz2out,nz3out]);
		SP.jjmax=reshape(jjmax,[nxout,nzout,nz2out,nz3out]);
	case 2,
		SP.xlab=sprintf(['Max over ',LG.xlab,' (%4.2g,%4.2g)'],xlims); plab=SP.xlab;
		SP.x=mean(LG.x(iix)); 
		SP.ixmax=reshape(iix(jjmax),[nyout,nzout,nz2out,nz3out]);
		SP.xmax=reshape(LG.x(iix(jjmax)),[nyout,nzout,nz2out,nz3out]);
		SP.jjmax=reshape(jjmax,[nyout,nzout,nz2out,nz3out]);
	case 3,
		SP.zlab=sprintf(['Max over ',LG.zlab,' (%4.2g,%4.2g)'],zlims); plab=SP.zlab;
		SP.z=mean(LG.z(iiz)); 
		SP.izmax=reshape(iiz(jjmax),[nyout,nxout,nz2out,nz3out]);
		SP.zmax=reshape(LG.z(iiz(jjmax)),[nyout,nxout,nz2out,nz3out]);
	case 4,
		SP.z2lab=sprintf(['Max over ',LG.z2lab,' (%4.2g,%4.2g)'],z2lims); plab=SP.z2lab;
		SP.z2=mean(LG.z2(iiz2)); 
		SP.z2max=reshape(LG.z2(iiz2(jjmax)),[nyout,nxout,nzout,nz3out]);
	case 5,
		SP.z3lab=sprintf(['Max over ',LG.z3lab,' (%4.2g,%4.2g)'],z3lims); plab=SP.z3lab;
		SP.z3=mean(LG.z3(iiz3)); 
		SP.z3max=reshape(LG.z3(iiz3(jjmax)),[nyout,nxout,nzout,nz2out]);
	end
	if iipr, disp(plab); end
end
%
%%size(data), SP, keyboard
if iidB,
	SP.data=10*log10(data);
	cmax=ceil(max(SP.data(:))/5)*5; cmin=max(floor(min(SP.data(:))/5)*5,cmax-60);
else
	SP.data=sqrt(data);
	cmax=max(SP.data(:)); cmin=min(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==1 & nc==2,
%	iiz1=find(z<=zlims(1)); if isempty(iiz1), iiz1=1; end
%	iiz2=find(z>=zlims(2)); if isempty(iiz2), iiz2=nz; end
%	iiz=iiz1(length(iiz1)):iiz2(1);
	dz=(z(nz) - z(1))/(max(1,nz-1));
	iiz=find(z>=zlims(1)-dz/2 & z<=zlims(2)+dz/2);
%%	length(iiz),length(z), keyboard,
	if isempty(iiz), 
		zlims,z(1),z(nz), 
		error(['USER ERROR: no overlap in ',zlab])
	end
elseif nc==1,
	iiz=zlims;
else
	disp('zlims must be [zmin,zmax] or [jz1;jz2; ... jzn]');
	error
end
%%z([1,end]),nz,length(z),zlims,length(iiz)
zlims=z(iiz([1,end])); zlims=zlims(:)';

	
%
% 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.
