function testderivative(varargin);
% testderivative(S, [options]);
% Function does not work anymore well
% Computes derivatives of basic limit functions in \RR^2 in two ways and compares the outcome.
% 1) Via the difference scheme
% 2) Via differences of the basic limit function
%
% Input: 
%   S           a subdivision scheme (passed to getS)  
%
% Options:
%   'it',val            iterations. Default = 6
%   'res',val           Resolution. Default = 0.1
%   'derivative',val    1x2 vector. Order of the partial derivative. Default=[0 1];
%   
% Output:
%   Nice plots
%
% See also: getS, blf
%
% Written by: tommsch, 2016

S=getS(varargin{1});
it=parsem('it',varargin,6);
STEP=parsem('res',varargin,.1);
dirvec=parsem('derivative',varargin,[0 1]);

dirnum=2;

subplot(3,1,1);
[c,PM,xyzv]=blf(S,'iteration',it); %blf


%interpolate to uniform grid. Then I can compute the real derivative
[xq,yq] = ndgrid(min(xyzv{1}(1,:)):STEP:max(xyzv{1}(1,:)), min(xyzv{1}(2,:)):STEP:max(xyzv{1}(2,:))); %We must use ndgrid. meshgrid has different x- and y-axis
vq1 = griddata(xyzv{1}(1,:),xyzv{1}(2,:),xyzv{1}(3,:),xq,yq);
%mesh(xq,yq,vq1);
d1vq=diffsequence(vq1,dirvec,'samesize')/STEP;
subplot(3,1,2);
    mesh(xq,yq,d1vq,'FaceAlpha',.4); %Derivative of blf
    xlabel('x')
    ylabel('y')
    zlabel('z')
    title(['Real derivative and Differences'])

%Try to compute derivative from differences of sequence c
dc=diffsequ(c,1);
dc=PM'*dc;
dcxyzv=compresscoordinates(dc.M{dirnum}.c, PM,'idx',dc.M{dirnum}.idx);
dcvq = griddata(dcxyzv(1,:),dcxyzv(2,:),dcxyzv(3,:),xq,yq);
mesh(xq,yq,dcvq);

err=d1vq-dcvq;
val=err(:);
val(isnan(val))=[];
norm(val,2)


subplot(3,1,3);
mesh(xq,yq,d1vq-dcvq);
xlabel('x')
ylabel('y')
zlabel('z')
title(['Error'])
axis equal

end

function dummy; end %#ok<DEFNU> %Generates an error, if the 'end' of a function is missing. 