function SP=data_join_mult(DD,kdim,kdim_scale,kdim_lab,iiview);
% Call SP=data_join_mult(DD,kdim,kdim_scale,kdim_lab,iiview);
%
% DD            = structure containing data structures to join;
% kdim		    = dimension over which to jlin datasets;
% kdim_scale    = axis scale to assign to enlarged kdim dimension;
% kdim_lab      = axis label to assign to enlarged kdim dimension;
% iiview	    = 1 to call view5 on the output SP (default is 0);
%
% See also data_join, data_interp
%
%
% 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, DD='*'; end
if ischar(DD),
    mat_files=dir([DD,'*.mat']); nfile=length(mat_files);
    if nfile==0, disp('No files found that match input string'); SP=[]; return; end
    for jfile=1:nfile,
        disp([int2str(jfile),': ',mat_files(jfile).name])
    end
    iif=str2num(input('Enter file numbers to join ([]=all): ','s'));
    if isempty(iif), iif=[1:nfile]; end
    DD=cell(length(iif),1);
    for jfile=1:length(iif);
        str=['xx=load(''',mat_files(iif(jfile)).name,''');']; eval(str);
        fields=fieldnames(xx); str=['DD{jfile}=xx.',fields{1},';']; eval(str);
    end
end
ndata=length(DD); ndim=DD{1}.ndim; 
if nargin<2,
    for jdata=1:ndata,
        str=sprintf('    %s   ',DD{jdata}.labels{1:DD{jdata}.ndim}); 
        disp(['DATA ',int2str(jdata),' axis labels: ',str]);
        for jdim=1:DD{jdata}.ndim,
            lims=DD{jdata}.scales{jdim}([1,end]);
            str=['   Axis limits = ',sprintf('[%1.2f %1.2f] nval = %i',[lims(:)',DD{jdata}.dims(jdim)])]; 
            disp(str);
        end
    end
    kdim=       str2num(input(['Enter dimension number over which to join ([]=make new dim): '],'s'));
end
if isempty(kdim),kdim=ndim+1; end
if kdim>ndim, 
    dim_new=ndata; 
else, 
    dim_new=0; for jdata=1:ndata, dim_new=dim_new + DD{jdata}.dims(kdim); end
end
if nargin<3,
    kdim_scale= str2num(input(['             Enter scale for new dimension (',int2str(dim_new),' long, []=1:n): '],'s'));
end
if nargin<4,
    kdim_lab=           input(['Enter label for new dimension (no quotes) ([] to keep same): '],'s');
end
if nargin<5,
    iiview=     str2num(input(['          Automatically view output? (1=yes,0=no,[]=yes)   : '],'s'));
end
if isempty(kdim_scale), kdim_scale=[1:dim_new]; end
if isempty(iiview), iiview=1; end
% Interpolate if not same size:
dims=zeros(ndim,ndata); ax_lims=zeros(2,ndim,ndata);
for jdata=1:ndata,
    if DD{jdata}.ndim~=ndim, error('All data structures must have same number of dimensions'); end
    dims(:,jdata)=DD{jdata}.dims(:);
    for jdim=1:ndim,
        ax_lims(1:2,jdim,jdata)=DD{jdata}.scales{jdim}([1,end]);
    end
end
ii_interp=zeros(1,ndim); jdata_take=zeros(1,ndim);
for jdim=1:ndim,
    if jdim~=kdim & ~all(dims(jdim,:)==dims(jdim,1)),
        [nax,jdata_take(jdim)]=max(dims(jdim,:));
        new_scales{jdim}=DD{jdata_take(jdim)}.scales{jdim};
        ii_interp(jdim)=1;
    end
end
%
for jdim=1:ndim,
    if ii_interp(jdim),
		for jdata=1:ndata,
            if jdata==jdata_take(jdim),
                % Do nothing
            elseif no_ovl(new_scales{jdim},DD{jdata}.scales{jdim}),
                new_scales2=new_scales;
                new_scales2{jdim}=linspace(DD{jdata}.scales{jdim}(1), ...
                    DD{jdata}.scales{jdim}(end),DD{jdata_take(jdim)}.dims(jdim));
                DD{jdata}=data_interp(DD{jdata},jdim,1,new_scales2,0);            
            else
                DD{jdata}=data_interp(DD{jdata},jdim,1,new_scales,0);            
            end
        end
    end
end
%
if kdim>ndim,
    if kdim~=ndim+1, error('Cannot expand more than one dimension!!'); end
    for jdata=1:ndata, 
        DD{jdata}.ndim=DD{jdata}.ndim+ 1;
        DD{jdata}.dims=[DD{jdata}.dims 1]; 
        DD{jdata}.scales{kdim}=0; DD{jdata}.labels{kdim}='Dataset #';
        if isfield(DD{jdata},'scales2'), DD{jdata}.scales2{kdim}=jdata; DD{jdata}.labels{kdim}='Dataset #'; end
    end
end
SP=structure_copy(DD{1},{'data','INFO'});
SP.dims(kdim)=0; SP.scales{kdim}=[]; 
for jdata=1:ndata;
    SP.dims(kdim)=SP.dims(kdim) + DD{jdata}.dims(kdim);
    SP.scales{kdim}=[SP.scales{kdim}, DD{jdata}.scales{kdim}(:)'];
	if isfield(DD{jdata},'scales2'),
        if jdata==1, SP.scales2{kdim}=[];  end
        SP.scales2{kdim}=[SP.scales2{kdim}, DD{jdata}.scales2{kdim}(:)'];
    end
    SP.INFO{jdata}=structure_copy(DD{jdata},{'title','ndim','dims','scales','scales2','labels','labels2','data','dlab', ...
            'iidB','iiq','vdims','size','col','sty','back_col','iigrid','ii_img','clim','iidB'},0);
	if ndata<=2, SP.title=[SP.title,'&',DD{jdata}.title]; end
end
if ndata>2, SP.title=[DD{jdata}.title,'&&']; end
if ~isempty(kdim_scale),
    if length(kdim_scale)~=SP.dims(kdim), 
        error(['Length of kdim_scale not equal to sum of lengths of DATA.scales{kdim}']); 
    end
    SP.scales{kdim}=kdim_scale(:)';
end
if ~isempty(kdim_lab),
    SP.labels{kdim}=kdim_lab;
    if isfield(SP,'labels2'), SP.labels2{kdim}=kdim_lab; end
end
SP.data=zeros(SP.dims);
%
KD=cell(1,SP.ndim);
for jdim=1:SP.ndim,
    KD{jdim}=[1:SP.dims(jdim)];
end
%
indx0=0;
for jdata=1:ndata,
    KD{kdim}=[indx0+1:indx0+DD{jdata}.dims(kdim)];
    SP.data(KD{:})=DD{jdata}.data;
    indx0=indx0 + DD{jdata}.dims(kdim);
end
%
if iiview, view5(SP); end
%
%%%
%
function no_ovl=no_ovl(scales1,scales2)
%
no_ovl=0;
if (scales1(1) - scales2(1))*(scales1(1) - scales2(end))>0 & ...
        (scales1(end) - scales2(1))*(scales1(end) - scales2(end))>0 & ...
        (scales2(1) - scales1(1))*(scales2(1) - scales1(end))>0 & ...
        (scales2(end) - scales1(1))*(scales2(end) - scales1(end))>0,
    no_ovl=1;
end
