function tsout = mtimes(ts1, ts2)
%MTIMES Overload matrix multiplication
%
% MTIMES successively performs matrix multiplication of each sample of ts1 
% by the corresponding sample in ts2. For example, if the size of ts1 is
% 5-by-2-by-3 with GridFirst true and the size of ts2 is 5-by-3-by-4 with
% GridFirst true, the output timeseries will have size of 5-by-2-by-4 with
% GridFirst true. 

%   Author(s): James G. Owen, Rong Chen
%   Copyright 1986-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.1 $ $Date: 2004/12/26 21:35:11 $

if strcmp(class(ts2),'tsdata.timeseries')
    % Check sizes match
    if ~isequal(getGridSize(ts1),getGridSize(ts2))
        error('timeseries:mtimes:badsizes',...
            'Time series have differing lengths')
    end

    % Merge time vectors onto a common basis
    [ts1timevec, ts2timevec, outprops] = ...
        timemerge(ts1.timeInfo, ts2.timeInfo,ts1.time,ts2.time);

    % Relative time vectors - remove initial values
    if isempty(outprops.ref) 
        ts1timevec = ts1timevec-ts1timevec(1);
        ts2timevec = ts2timevec-ts2timevec(1);
    end

    % Check that the time vectors match
    intervalLen = ((ts1timevec(end)-ts1timevec(1))/length(ts1timevec));
    if norm(ts1timevec-ts2timevec)/intervalLen>1e-6
        error('timeseries:mtimes:badtimes','Time vectors do not match')
    end

    % Load the data so its not extracted for each pass of the loop and check
    % the dimensions
    data1 = ts1.Data;
    s1 = hdsGetSize(data1);
    data2 = ts2.Data;
    s2 = hdsGetSize(data2);
    if length(s1)>3
        error('timeseries:mtimes:non3dimmax', ...
            'Cannot multiply time series with data of dimension > 3')
    end

    % carry out data alignment to
    [data1,data2,s1,s2] = tsAlignSizes(data1, ...
        ts1.DataInfo.GridFirst,data2,ts2.DataInfo.GridFirst);
    if length(s1)~=length(s2)
        error('timeseries:mtimes:nonconfdim', ...
            'Cannot multiply time series where data has differing dimensions')
    end

    % For mtimes if both input data is 2d then arithmatic can be performed
    % directly on the ord data. Otherwise we must step through one sample at
    % a time and perform the mtimes on a one sample slice of ord data

    if length(s1)==2
        if s1(2)==s2(2)
            s = [s1(1) 1]; % Output array size, in which each sample is a scalor
            try  
                dataout = hdsNewArray(data1,s); 
                for k=1:s(1)
                    dataout(k) = data1(k)*data2(k)';
                end
            catch
                % Try casting to double in case ordinate data type does not support
                % mtimes arithmatic 
                dataout = hdsNewArray(1,s); 
                try
                    data1 = double(data1);
                    data2 = double(data2);
                    for k=1:s(1)
                        dataout(k) = data1(k)*data2(k)';
                    end
                catch
                   errstr = sprintf('%s\n%s','mrdivide failed for this type of ordinate data object.', ...
                             'Try implementing a double cast method or overloaded addition for this object');
                   error('timeseries:mtimes:badcast',errstr)
                end
            end
        % Check conformability
        else
            msg = ['Sample sizes of ', num2str(s1) , ' and ',...
                num2str(s2), ' are not conformable'];
            error('timeseries:mrdivide:nonconfdim', msg)
        end
    else % Perform mtimes on each sample for 3d arrays -> slow	
        if s1(3)==s2(2)
            s = [s1(1) s1(2) s2(3)]; % Output array size 
            try  
               dataout = hdsNewArray(data1,s); 
               for k=1:s(1)
                  r = [{k} repmat({':'},[1 length(s(2:end))])];
                  dataout = hdsSetSlice(dataout,r,feval('mtimes',hdsReshapeArray(hdsGetSlice(data1,r),s1(2:end)), ...
                      hdsReshapeArray(hdsGetSlice(data2,r),s2(2:end))));
               end
            catch
               % Try casting to double in case ordinate data type does not support
               % mtimes arithmatic 
               dataout = hdsNewArray(1,s); 
               try
                   data1 = double(data1);
                   data2 = double(data2);
                   for k=1:s(1)
                       r = [{k} repmat({':'},[1 length(s(2:end))])];
                       dataout = hdsSetSlice(dataout,r,feval('mtimes',hdsReshapeArray(hdsGetSlice(data1,r),s1(2:end)), ...
                          hdsReshapeArray(hdsGetSlice(data2,r),s2(2:end))));  
                   end
               catch
                   errstr = sprintf('%s\n%s','mtimes failed for this type of ordinate data object.', ...
                             'Try implementing a double cast method or overloaded addition for this object');
                   error('timeseries:mtimes:badcast',errstr)
               end
            end
        else
            msg = ['Sample sizes of ', num2str(s1(2:end)) , ' and ',...
                num2str(s2(2:end)), ' are not conformable'];
            error('timeseries:mtimes:nonconfdim', msg)
        end
    end    

    % Build output time series. If both time series are subclasses 
    % try to make the output the same classs
    if classhandle(ts1) == classhandle(ts2)
        c = classhandle(ts1);
        tsout = eval(sprintf('%s.%s',c.Package.Name,c.Name));
        set(tsout,'Data', dataout,'Time',ts1timevec)
    else
        tsout = tsdata.timeseries;
        init(tsout,dataout,ts1timevec);
    end

    set(tsout,'Name',['Mtimes_' ts1.Name '_' ts2.Name]);
    set(tsout.timeInfo,'Startdate',outprops.ref,'Units', ...
        outprops.outunits,'Format',outprops.outformat);

    % Quality arithmatic - merge quality info and combine quality codes
    if ~isempty(ts1.qualityInfo) && ~isempty(ts2.qualityInfo)
        tsout.qualityInfo = merge(ts1.qualityInfo,ts2.qualityInfo);
    end
    if ~isempty(get(get(tsout,'qualityInfo'),'Code')) && ~isempty(ts1.quality) && ~isempty(ts2.quality)
        tsout.quality = min(ts1.quality,ts2.quality);
    end

    % Merge ordinate metadata
    tsout.dataInfo = mtimes(ts1.dataInfo,ts2.dataInfo);

    % Merge events
    addevent(tsout,horzcat(ts1.Events,ts2.Events));
elseif isnumeric(ts2)
    if isscalar(ts2)
        tsout = ts1;
        tsout.data=ts1.data.*ts2;
        return
    end

    if ~isequal(getGridSize(ts1),[size(ts2,1) 1])
        error('timeseries:mrdicide:badsizes',...
            'The first dimension of the numeric array has to agree with the length of the time series.  ')
    end
    % the second input argument is a double and 
    data1 = ts1.Data;
    s1 = hdsGetSize(data1);
    if ts1.datainfo.gridfirst==false
        data1=reshape(data1,[s1(end) s1(1:end-1)]);
    end
    s1 = hdsGetSize(data1);
    s2=size(ts2);
    if length(s1)>3
        error('timeseries:mrdicide:non3dimmax', ...
            'Cannot multiply time series with data of dimension > 3')
    end
    % carry out data alignment to
    if length(s1)~=length(s2)
        error('timeseries:mrdicide:sizemis','Size of the data arrays are mismatching')
    end
    
    % For mrdivide if both input data is 2d then arithmatic can be performed
    % directly on the ord data. Otherwise we must step through one sample at
    % a time and perform the mtimes on a one sample slice of ord data

    if length(s1)==2
        if s1(2)==s2(2)
            s = [s1(1) 1]; % Output array size, in which each sample is a scalor
            try  
                dataout = hdsNewArray(data1,s); 
                for k=1:s(1)
                    dataout(k) = data1(k,:)*ts2(k,:)';
                end
            catch
                % Try casting to double in case ordinate data type does not support
                % mtimes arithmatic 
                dataout = hdsNewArray(1,s); 
                try
                    data1 = double(data1);
                    for k=1:s(1)
                        dataout(k) = data1(k,:)*ts2(k,:)';
                    end
                catch
                   errstr = sprintf('%s\n%s','mrdivide failed for this type of ordinate data object.', ...
                             'Try implementing a double cast method or overloaded addition for this object');
                   error('timeseries:mtimes:badcast',errstr)
                end
            end
        % Check conformability
        else
            msg = ['Sample sizes of ', num2str(s1) , ' and ',...
                num2str(s2), ' are not conformable'];
            error('timeseries:mrdivide:nonconfdim', msg)
        end
    else 
        if s1(3)==s2(2)
            s = [s1(1) s1(2) s2(3)]; % Output array size 
            tmp=repmat({':'},[1 length(s(2:end))]);
            try  
               dataout = hdsNewArray(data1,s); 
               for k=1:s(1)
                  r = [{k} tmp];
                  dataout = hdsSetSlice(dataout,r,feval('mtimes',hdsReshapeArray(hdsGetSlice(data1,r),s1(2:end)), ...
                      reshape(ts2(k,:,:),s2(2:end))));
               end
            catch
               % Try casting to double in case ordinate data type does not support
               % mtimes arithmatic 
               dataout = hdsNewArray(1,s); 
               try
                   data1 = double(data1);
                   data2 = double(data2);
                   for k=1:s(1)
                       r = [{k} tmp];
                       dataout = hdsSetSlice(dataout,r,feval('mtimes',hdsReshapeArray(hdsGetSlice(data1,r),s1(2:end)), ...
                           reshape(ts2(k,:,:),s2(2:end))));
                   end
               catch
                   errstr = sprintf('%s\n%s','mtimes failed for this type of ordinate data object.', ...
                             'Try implementing a double cast method or overloaded addition for this object');
                   error('timeseries:mtimes:badcast',errstr)
               end
            end
        else
            msg = ['Sample sizes of ', num2str(s1(2:end)) , ' and ',...
                num2str(s2(2:end)), ' are not conformable'];
            error('timeseries:mtimes:nonconfdim', msg)
        end
    end    

    tsout = ts1;
    if ts1.datainfo.gridfirst==true
        tsout.data=dataout;
        tsout.datainfo.gridfirst=true;
    else
        s_out = hdsGetSize(dataout);
        tsout.data=reshape(dataout,[s_out(2:end) s_out(1)]);
        tsout.datainfo.gridfirst=false;
    end
else
    error('timeseries:mtimes:sizemis','Types of the data arrays are mismatching')
end