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 elementsif prod(length(I)) > 0, error('A must be completely defined (no NaN)'), endgeomOld  = 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 identicalsamepd = 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 interpolationendif isempty(val), warning('requested interpolation not yet supported'), endB = scalarfield(R,val);