function SP=data_join(D1,D2,kdim,iiview);
% Call SP=data_join(D1,D2,kdim,iiview);
%
% D1     = first data structure;
% D2     = second data structure;
% kdim		= dimension over which to join datasets;
% iiview	= 1 to call view5 on the output SP (default is 0);
%
% See also data_join_mult, 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==2,
    str=sprintf('    %s   ',D1.labels{1:D1.ndim}); disp(['DATA 1 Axis labels: ',str]);
    str=sprintf('    %s   ',D2.labels{1:D2.ndim}); disp(['DATA 2 Axis labels: ',str]);
    str=['DATA 1/2 Axis limits: '];
    for jdim=1:D1.ndim, 
        str=[str,'   DATA 1: ',sprintf('[%1.2f %1.2f]',D1.scales{jdim}([1,end]))]; 
        str=[str,'   DATA 2: ',sprintf('[%1.2f %1.2f]\n',D2.scales{jdim}([1,end]))]; 
    end
    disp(str)
    str=sprintf('       %i       ',D1.dims); disp(['DATA 1 Axis lengths: ',str]);
    str=sprintf('       %i       ',D2.dims); disp(['DATA 2 Axis lengths: ',str]);
    kdim=    str2num(input(['                 Enter dimension number over which to join: '],'s'));
    iiview=  str2num(input(['          Automatically view output? (1=yes,0=no,[]=yes)   : '],'s'));
    if isempty(iiview), iiview=1; end
else
    if nargin<4 | isempty(iiview), iiview=0; end
end
%
% Interpolate all non-joining dimensions if they are not the same size:
ndim_min=min(D1.ndim,D2.ndim);
for jdim=1:ndim_min,
    if jdim~=kdim,
        if D1.dims(jdim)>D2.dims(jdim),
            new_scales=D2.scales;
            if no_ovl(new_scales{jdim},D1.scales{jdim}),
                new_scales{jdim}=linspace(D2.scales{jdim}(1),D2.scales{jdim}(end),D1.dims(jdim)),
            else
                new_scales{jdim}=D1.scales{jdim};
            end
            D2=data_interp(D2,jdim,1,new_scales,0);
        elseif D1.dims(jdim)<D2.dims(jdim),
            new_scales=D1.scales; 
            if no_ovl(new_scales{jdim},D2.scales{jdim}),
                new_scales{jdim}=linspace(D1.scales{jdim}(1),D1.scales{jdim}(end),D2.dims(jdim)),
            else
                new_scales{jdim}=D2.scales{jdim};
            end
            D1=data_interp(D1,jdim,1,new_scales,0);
        end
    end
end
%
SP=structure_copy(D1,{'data'}),
SP.dims(kdim)=D1.dims(kdim) + D2.dims(kdim);
SP.scales{kdim}=[D1.scales{kdim}(:)', D2.scales{kdim}(:)'];
if isfield(D1,'scales2'),
	SP.scales2{kdim}=[D1.scales2{kdim}(:)', D2.scales2{kdim}(:)'];
end
SP.title=[D1.title,'&',D2.title];
SP.data=zeros(SP.dims);
%
KD=cell(1,SP.ndim);
for jdim=1:SP.ndim,
    KD{jdim}=[1:SP.dims(jdim)];
end
%
KD{kdim}=[1:D1.dims(kdim)];
SP.data(KD{:})=D1.data;
KD{kdim}=[D1.dims(kdim) + [1:D2.dims(kdim)]];
SP.data(KD{:})=D2.data;
%
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
