interpontorq.m

MDPSAS工具箱是马里兰大学开发的

function B = interpontorq(R,A)

% interpontorq.m Interpolate a scalarfield object to a different
%                physical domain with the same number of dimensions.
%                Called as
%
%          B = interpontorq(relQG,A)
%
%	INPUT parameters:
%           relQG : relquad object defining new physical domain
%               A : scalarfield object
%
%	OUTPUT:
%	        B : a scalarfield object. NOTE: B contains
%                   a relquad object. Use rq2qg(B) to convert
%                   to a quadgrid object
%
% Note: Because it uses dtint (so is very accurate), for the
% interpolation to work, A must be completely defined (i.e., no NaN).

Aval  = get(A,'val');
Apd   = get(A,'pd');

I = find( isnan(Aval) ); % returns indices of NaN elements
if prod(length(I)) > 0, error('A must be completely defined (no NaN)'), end

geomOld  = get(Apd,'geom');
qpOld    = get(Apd,'qp');
nameOld  = get(Apd,'name');
Qinv     = get(Apd,'Qinv');

geomNew  = get(R,'geom');
qpNew    = get(R,'qp');
nameNew  = get(R,'name');

if length(nameOld) >= 3
    error('Not yet set up for higher-dimensional scalarfield objects')
end

% Check if any quadgrid directions are identical

samepd = compare(Apd,R);
val = [];

%%%%%%%%%%%%%%% 1D to 1D interpolation %%%%%%%%%%%%%%%

if length(nameOld) == 1

    if length(nameNew) ~= 1
        error('scalarfield/new quadgrid dimension mismatch')
    end

    val = interp(A,qpNew{1}+R.offset(1));
end

%%%%%%%%%%%%%%% 2D to 2D interpolation %%%%%%%%%%%%%%%

if length(nameOld) == 2

    if length(nameNew) ~= 2
        error('scalarfield/new quadgrid dimension mismatch')
    end

    %%% Cartesian (old) to Cartesian (new) interpolation
    %%% Cylindrical (old) to Cylindrical (new) interpolation

    if ( strcmp(geomOld{1},'slab') & strcmp(geomNew{1},'slab') ...
            & strcmp(geomOld{2},'slab') & strcmp(geomNew{2},'slab') ) ...
            | ( strcmp(geomOld{1},'cyln') & strcmp(geomNew{1},'cyln') ...
            & strcmp(geomOld{2},'slab') & strcmp(geomNew{2},'slab') )

        x = qpNew{1}+R.offset(1);
        y = qpNew{2}+R.offset(2);

        val = interp(A,x,y);

    end % Cartesian-to-Cartesian interpolation
    % Cylindrical-to-Cylindrical interpolation

    %%% Cartesian (old) to Polar (new) interpolation

    if strcmp(geomOld{1},'slab') & strcmp(geomNew{1},'cyln') ...
            & strcmp(geomOld{2},'slab') & strcmp(geomNew{2},'peri')

        % convert new polar coordinates to Cartesian

        nr = size(qpNew{1},1);
        nt = size(qpNew{2},1);

        xNew = qpNew{1}*cos(qpNew{2}');
        yNew = qpNew{1}*sin(qpNew{2}');

%         x = reshape(xNew,nr*nt,1)+R.offset(1);
%         y = reshape(yNew,nr*nt,1)+R.offset(2);
% 
        x = xNew+R.offset(1);
        y = yNew+R.offset(2);

        for i = 1:size(x,1)
            for j = 1:size(x,2)
                val(i,j) = interp(A,x(i,j),y(i,j));
            end
        end

    end % Cartesian-to-Polar interpolation

    %%% Polar (old) to Cartesian (new) interpolation

    if strcmp(geomOld{1},'cyln') & strcmp(geomNew{1},'slab') ...
            & strcmp(geomOld{2},'peri') & strcmp(geomNew{2},'slab')

        nx = size(qpNew{1},1);
        ny = size(qpNew{2},1);

        % first shift the (new) Cartesian grid by the offset specified
        % in the (base) polar coordinates, then transform the shifted
        % coordinates to polar in order to do the interpolation

        xOffset = R.offset(1)*cos(R.offset(2));
        yOffset = R.offset(1)*sin(R.offset(2));

        [xNew, yNew] = ndgrid(qpNew{1}+xOffset,qpNew{2}+yOffset);

        r = sqrt( xNew.^2 + yNew.^2 );
        t = atan2(yNew, xNew);

        Ineg = find( t < 0 );
        t(Ineg) = t(Ineg) + 2*pi;

        for i = 1:size(r,1)
            for j = 1:size(r,2)
                val(i,j) = interp(A,r(i,j),t(i,j));
            end
        end

    end % Polar-to-Cartesian interpolation

    %%% Polar (old) to Polar (new) interpolation

    if strcmp(geomOld{1},'cyln') & strcmp(geomNew{1},'cyln') ...
            & strcmp(geomOld{2},'peri') & strcmp(geomNew{2},'peri')

        % first convert the (new) polar grid to Cartesian coordinates,
        % then shift the (new) Cartesian grid by the offset specified
        % in the (base) polar coordinates, then transform the shifted
        % coordinates to polar in order to do the interpolation

        if samepd(1) & sum(abs(R.offset)) == 0

            % r coordinates are identical and there is not offset, so this
            % is a spacial case of rotation only

            rNew = qpNew{1};
            tNew = qpNew{2};

            for i = 1:length(rNew)
                [val(:,i),extrapflag] = ...
                    dtint(qpOld{2},Aval(i,:)',tNew,Qinv{2},'peri');
                if extrapflag, PrintExtrapolationWarning = 1; disp('dim 2'), end
            end

            val = val';

        else % most general polar-to-polar interpolation

            nr = size(qpNew{1},1);
            nt = size(qpNew{2},1);

            xNew = qpNew{1}*cos(qpNew{2}');
            yNew = qpNew{1}*sin(qpNew{2}');

            xOffset = R.offset(1)*cos(R.offset(2));
            yOffset = R.offset(1)*sin(R.offset(2));

            xNew = xNew+xOffset;
            yNew = yNew+yOffset;

            r = sqrt( xNew.^2 + yNew.^2 );
            t = atan2(yNew, xNew);

            Ineg = find( t < 0 );
            t(Ineg) = t(Ineg) + 2*pi;
            
            for i = 1:size(r,1)
                for j = 1:size(r,2)
                    val(i,j) = interp(A,r(i,j),t(i,j));
                end
            end
        end

    end % Polar-to-Polar interpolation

    %     %%% Periodic (old) to Slab (new) interpolation
    %
    %     if strcmp(geomOld{1},'slab') & strcmp(geomNew{1},'slab') ...
    %             & strcmp(geomOld{2},'peri') & strcmp(geomNew{2},'slab')
    %
    %         for i = 1:size(Aval,2)
    %             [val_int(:,i),extrapflag] = ...
    %                 dtint(qpOld{1},Aval(:,i),qpNew{1}+R.offset(1),Qinv{1},'slab');
    %             if extrapflag, PrintExtrapolationWarning = 1; end
    %         end
    %
    %         % Note! Because the limits of the periodic quadrature points do not
    %         % necessarily have to be defined in [0 2pi], we must scale everything
    %         % relative to the base quadgrid scaled to [0 2pi]
    %
    %         scale2pi = 2*pi/(qpOld{2}(end)-qpOld{2}(1));
    %
    %         for i = 1:size(val_int,1)
    %             [val_temp,extrapflag] = ...
    %                 dtint(qpOld{2}*scale2pi,val_int(i,:)', ...
    %                 qpNew{2}*scale2pi+R.offset(2)*scale2pi,Qinv{2},'peri');
    %             val(i,:) = val_temp';
    %             if extrapflag, PrintExtrapolationWarning = 1; end
    %         end
    %
    %     end % Periodic-to-Slab interpolation

end

if isempty(val), warning('requested interpolation not yet supported'), end

B = scalarfield(R,val);