%MXCURVEC_CIRCLE_TEST simple test for mxcurvec
%
% determine 2 arrows in a strictly circular flow field
% wehre all arrows end up at their own tail
%
%See also: mxcurvec, mxcurvec_linear_test

%% define scenario: stream function: 
%  * perfect circle:        wherever you start,
%  * identical "rondetijd": wherever you start.

    C = 1e-6;
    
    x = [-1e3:100:1e3];
    y = [-1e3:100:1e3];
    
   [XX,YY] = meshgrid(x,y);
    X    = corner2center(XX);
    Y    = corner2center(YY);
   [TH,R]= cart2pol(XX,YY);
    PHI  = -C.*R.^2;
   [V,U] = gradient2(XX,YY,PHI);V=-V;
  % U    = -2*C.*R.*sin(TH);
  % V    = +2*C.*R.*cos(TH);
     
    OPT.dt               = 1*pi/C; % still to calculate analytical "rondetijd" f(C)
    OPT.nrvertices       = 100;%floor(OPT.dt./1e5);
    OPT.headthickness    = 0.15;
    OPT.arrowthickness   = 0.05;
    OPT.nhead            = 3;
    OPT.iopt             = 0;% 0: (x,y), 1: (lat,lon)
    
    x0 = [400 800];
    y0 = [  0   0];
    relwdt=ones(size(x0)).*0;
   
%% plot scenario

    pcolorcorcen(XX,YY,PHI);
    colorbarwithtitle('\Phi')
    hold on
    quiver(X,Y,U,V,'k')
    plot(x0,y0,'r*')
    axis equal

%% call mxcurvec

    [xar,yar,xax,yax,len] = mxcurvec(x0(:),y0,... % should be [n x 1], so use x(:), and can not have length==1
                                 X,Y,U,V,U,V,...  % does not have to be [n x 1]
                                 OPT.dt,...
                                 OPT.nrvertices,...
                                 OPT.headthickness,...
                                 OPT.arrowthickness,...
                                 OPT.nhead,...
                                 relwdt,...       % does not have to be [n x 1]
                                 0);
%% curvec.m code:

%% Set 999.999 arrow separators to nan values
    
    xar(xar<1000.0 & xar>999.998)=NaN;
    yar(yar<1000.0 & yar>999.998)=NaN;
    xax(xax<1000.0 & xax>999.998)=NaN;
    yax(yax<1000.0 & yax>999.998)=NaN;
    
%% Count number of points per arrow
    
    ic=1;while ~isnan(xar(ic,1));ic=ic+1;end % NOT: ic =  (OPT.nrvertices - OPT.nhead)*2+5 due to nhead cut-off in mxcurvec

%% Put all arrows in 2D matrix
    
    polx=reshape(xar,[ic               length(x0)]);
    poly=reshape(yar,[ic               length(x0)]);
    xax =reshape(xax,[OPT.nrvertices+1 length(x0)]);
    yax =reshape(yax,[OPT.nrvertices+1 length(x0)]);

%% add mxcurvec result
    
    plot  (xax         ,yax         ,'k.-')
    plot  (xax(end-1,:),yax(end-1,:),'ko')
    
%% check whether particles are back to square one after time required for one one full circle

    OK = all(abs(xax(1,:) - xax(end-1,:)) < 1e-6)
    
    [xax(1,:) xax(end-1,:)]