function [xout, vout, sout, hout, fout] = cont(varargin)
%
% CONTINUE(cds.curve, x0, v0, options)
%
% Continues the curve from x0 with optional directional vector v0
% options is a option-vector created with CONTSET
% The first two parameters are mandatory.

% WM: Rearranged many things throughout the code to put it
% all in a more logical order
global cds driver_window MC
%    Interactive plotting
stop;
drawnow;
if isfield(MC,'mainwindow')
    status = MC.mainwindow.mstatus;
    if isfield(MC.mainwindow,'options_window')
        plotpoints = get(MC.mainwindow.options_window,'Userdata');
    else 
        plotpoints = 1;
    end
else MC.D2 = [];MC.D3 = [];MC.PRC=[];MC.dPRC=[];MC.numeric_fig = []; 
    status = [];plotpoints = inf; 
    MC.mainwindow.duration = [];
    MC.mainwindow.mstatus = [];
end
set(status,'String','computing');

if (nargin < 6)% case not extend
    cds = [];
    cds.curve = '';
    warning off;
    lastwarn('');
    
    % Parse given options
    [cds.curve, x0, v0, opt] = ParseCommandLine(varargin{:});
    
    % Do some "stupid user" checks
    checkstupid(x0);
    cds.ndim = length(x0);
    curvehandles = feval(cds.curve);
    cds.curve_func = curvehandles{1};
    cds.curve_defaultprocessor = curvehandles{2};
    cds.curve_options = curvehandles{3};
    cds.curve_jacobian = curvehandles{4};
    cds.curve_hessians = curvehandles{5};
    cds.curve_testf = curvehandles{6};
    cds.curve_userf = curvehandles{7};
    cds.curve_process = curvehandles{8};
    cds.curve_singmat = curvehandles{9};
    cds.curve_locate = curvehandles{10};
    cds.curve_init = curvehandles{11};
    cds.curve_done = curvehandles{12};
    cds.curve_adapt = curvehandles{13};
    % Read out all options or set default values
    %
    % Get options from curve file
    options = feval(cds.curve_options);
%     options = feval(cds.curve, 'options');
    % Merge options from curve with cmdline, cmdline overrides curve
    cds.options = contmrg(opt, options);
    
    cds.options.MoorePenrose      = contget(cds.options, 'MoorePenrose', 1);
    cds.options.SymDerivative     = contget(cds.options, 'SymDerivative', 0);
    cds.options.SymDerivativeP    = contget(cds.options, 'SymDerivativeP', 0);
    cds.options.AutDerivative     = contget(cds.options, 'AutDerivative', 1);
    cds.options.AutDerivativeIte     = contget(cds.options, 'AutDerivativeIte',24);
    cds.options.Increment         = contget(cds.options, 'Increment', 1e-5);
    cds.options.MaxNewtonIters    = contget(cds.options, 'MaxNewtonIters', 3);
    cds.options.MaxCorrIters      = contget(cds.options, 'MaxCorrIters', 10);
    cds.options.MaxTestIters      = contget(cds.options, 'MaxTestIters', 10);
    cds.options.FunTolerance      = contget(cds.options, 'FunTolerance', 1e-6);
    cds.options.VarTolerance      = contget(cds.options, 'VarTolerance', 1e-6);
    cds.options.TestTolerance     = contget(cds.options, 'TestTolerance', 1e-5);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    Userfunctions         = contget(cds.options, 'Userfunctions',0);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Singularities         = contget(cds.options, 'Singularities', 0);
    WorkSpace             = contget(cds.options, 'WorkSpace', 0);
    Backward              = contget(cds.options, 'Backward',0);
    CheckClosed           = contget(cds.options, 'CheckClosed', 50);
    npoints               = contget(cds.options, 'MaxNumPoints', 300);
    Adapt                 = contget(cds.options, 'Adapt',3);
    
    Locators              = contget(cds.options, 'Locators', []);
    IgnoreSings           = contget(cds.options, 'IgnoreSingularity', []);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    UserInfo             = contget(cds.options, 'UserfunctionsInfo', []);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cds.h                 = contget(cds.options, 'InitStepsize' , 0.01);
    cds.h_max             = contget(cds.options, 'MaxStepsize'  , 0.1);
    cds.h_min             = contget(cds.options, 'MinStepsize'  , 1e-5); 
    
    cds.pJac = [];
    cds.pJacX = [];
    
    % Set algorithm specific variables & declarations
    %
    cds.h_inc_fac = 1.3;
    cds.h_dec_fac = 0.5;
    
    x1   = zeros(cds.ndim,1);           %
    x2   = zeros(cds.ndim,1);           % convenient to have these 4 vectors,
    v1   = zeros(cds.ndim,1);           % x1 and x2 are last two points
    v2   = zeros(cds.ndim,1);           % v1 and v2 the corresponding tangent vectors
    
    i    = 0;                       % our iteration parameter
    ind  = [];                       % our iteration parameter to plot
    it   = 0;                       % number of newton steps
    
    %interactive plotting
    numsing = sing_numeric(IgnoreSings);
    % Get info about testfunctions and singularities
    %
    
    cds.nActTest = 0;
    if Singularities
        [cds.S,SingLables] = feval(cds.curve_singmat);
        [cds.nSing, cds.nTest] = size(cds.S);
        Singularities = cds.nSing>0 & cds.nTest>0;
        % setup testfunction variables and stuff
        %
        
        if Singularities
            % Ignore Singularities
            cds.S(IgnoreSings,:) = 8;
            ActSing = setdiff(1:cds.nSing, IgnoreSings);
            cds.ActSing = ActSing;
            cds.nActSing = length(ActSing);
            
            % Which test functions must vanish somewhere in S?
            ActTest = find( sum((cds.S==0),1) > 0 );
            cds.ActTest = ActTest;
            cds.nActTest = length(ActTest);
            % WM: Build matrix with indices of testfunctions which
            % should be zero for each singularity            
            cds.SZ = zeros(cds.nTest+1,cds.nActSing+1);
            ml = 2;
            for si=ActSing(1:cds.nActSing)
                t = find( cds.S(si,:) == 0 )';
                l = size(t,1);
                cds.SZ(1:l,si) = t;
                ml = max(ml,l);
            end
            cds.SZ = cds.SZ(1:ml,:);
             cds.atv = 1;    cds.testvals = zeros(2, cds.nActTest);  
            % 1st row: testvals at x1, 2nd: testvals at x2
            cds.testzero = zeros(cds.ndim, cds.nActTest);  % coords where testf is zero
            cds.testvzero = zeros(cds.ndim, cds.nActTest);  % v where testf is zero
            
            if isempty(ActTest)
                Singularities = 0;
            end
        end
    end
    if Userfunctions
        cds.nUserf = size(UserInfo,2);cds.utv = 1;
        cds.uservals = zeros(2,cds.nUserf);%1st row testvals at x1,2nd testvals at x2
    end
    
    
    xout = zeros(cds.ndim,npoints);     % each point on curve
    vout = zeros(cds.ndim,npoints);     % all tangent vectors
    %hout = zeros(2+cds.nActTest+size(UserInfo,2),npoints); % keep track of all stepsizes & testfunctions
    sout = [];                          % all occured singularities
%     fout = [];
    
    %interactive plotting
    if size(MC.D2)>0|size(MC.D3)>0|~isempty(MC.numeric_fig)|~isempty(MC.PRC)|~isempty(MC.dPRC)
        output(numsing,xout,[],[],[],1);
    end
    
    % Algorithm starts here
    %

    StartTime = clock;
        
    % restarting point...
    
    % Init user space
    %
    if WorkSpace
        if feval(cds.curve_init, x0, v0)~=0
            delete(driver_window);
            set(status,'String','error');
            errordlg('Initializer failed.');
            return;
        end
    end
    % if x0 and v0 known: assume point on curve
    % if v0 unknown: correct x0

    if isempty(v0)
        [x0, v0] = CorrectStartPoint(x0, v0);
        if isempty(x0)
            delete(driver_window);
            set(status,'String','No convergence at x0');
            debug('No convergence at x0.\n') 
            debug('elapsed time  = %.1f secs\n', etime(clock, StartTime));
            debug('0 npoints');            
            errordlg('No convergence at x0.');   
            xout = []; vout = []; sout = []; hout = []; fout = [];
            return;            
        end
    end
    if Backward
        v0 = -v0;
    end
    
    % WM: added call to the default processor
    s.index = 1;
    s.label = '00';
    s.data  = [];
    s.msg   = 'This is the first point of the curve';
    
    [failed,f,s] = DefaultProcessor(x0,v0,s);
    debug('first point found\n');%printv(x0);
    debug('tangent vector to first point found\n'); %printv(v0);
    x1 = x0;
    v1 = v0;
    cds.testvals = [];
    cds.uservals = [];
    if Singularities
        % WM: calculate all testfunctions at once
        [tfvals,failed] = EvalTestFunc(ActTest,x0,v0);
        cds.testvals(2,:) = tfvals(ActTest);
        if ~isempty(failed)
            delete(driver_window);
            set(status,'String','error');
            errordlg('Evaluation of test functions failed at starting point.');
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if Userfunctions
        [ufvals,failed]=feval(cds.curve_userf, UserInfo,1:cds.nUserf, x0, v0);
        cds.uservals(2,:)=ufvals;
        if ~isempty(failed)
            delete(driver_window);
            set(status,'String','error');
            errordlg('Evaluation of userfunctions failed at starting point.'); 
        end        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    xout(:,1) = x0;
    vout(:,1) = v0;
    sout = [s]; 
    if Singularities & ~Userfunctions
        hout(:,1) = [0;0;cds.testvals(2,:)'];
    elseif Userfunctions & ~Singularities
        hout(:,1) = [0;0;cds.uservals(2,:)'];
    elseif Userfunctions & Singularities
        hout(:,1) = [0;0;cds.uservals(2,:)';cds.testvals(2,:)'];
    else
        hout(:,1) = [0;0];
    end
    sf = size(f,1);
    if sf > 0
        fout = zeros(size(f,1),npoints);
    else
        fout = zeros(1,npoints);
    end
    fout(:,1) = f;
    i   = 1;
    ind = [1];
else
    xout = varargin{1};  i  = size(xout,2);
    x0 = xout(:,1);   x1 = xout(:,end);
    vout = varargin{2};
    v1 = vout(:,end);    ind = [1];
    npoints = contget(cds.options, 'MaxNumPoints', 300);
    xout = [ xout zeros(cds.ndim,npoints)];
    vout = [ vout zeros(cds.ndim,npoints)]; 
    sout =  varargin{3};
    hout =  varargin{4};
    fout = [ varargin{5} zeros(size(varargin{5},1),npoints)];
    cds  = varargin{6};
    npoints = npoints + i;
    sout(end) = [];
    s = sout(end);
    
  
    Userfunctions         = contget(cds.options, 'Userfunctions',0);
    Singularities         = contget(cds.options, 'Singularities', 0);
    WorkSpace             = contget(cds.options, 'WorkSpace', 0);
    CheckClosed           = contget(cds.options, 'CheckClosed', 50);
    Adapt                 = contget(cds.options, 'Adapt', 1);
    
    Locators              = contget(cds.options, 'Locators', [])
    IgnoreSings           = contget(cds.options, 'IgnoreSingularity', []);
    UserInfo              = contget(cds.options, 'UserfunctionsInfo', []);
  
    if Singularities
        ActTest = cds.ActTest;
        ActSing = cds.ActSing;
        [cds.S,SingLables] = feval(cds.curve_singmat);
        
        if isempty(ActTest)
            Singularities = 0;
        end
    end
    if Userfunctions
        cds.nUserf = size(UserInfo,2);cds.utv = 1;
        cds.uservals = zeros(2,cds.nUserf);%1st row testvals at x1,2nd testvals at x2
    end
    if Userfunctions
        [ufvals,failed]=feval(cds.curve_userf, UserInfo,1:cds.nUserf, x1, v1);
        cds.uservals(2,:)=ufvals;
        if ~isempty(failed)
            delete(driver_window);
            set(status,'String','error');
            errordlg('Evaluation of userfunctions failed at starting point.'); 
        end        
    end


    
     %interactive plotting
    numsing = sing_numeric(IgnoreSings);

    if size(MC.D2)>0|size(MC.D3)>0|~isempty(MC.numeric_fig)|~isempty(MC.PRC)|~isempty(MC.dPRC)
        output(numsing,xout,[],[],[],1);
    end
    
    StartTime = clock;
    debug('start computing extended curve\n');
end
while i < npoints
    % correction
    %
    corrections = 1;
    while 1
      % predict next point  
      xpre = x1 + cds.h * v1;
      
      [x2, v2, it] = newtcorr(xpre, v1);
        
      if ~isempty(x2) & v1'*v2 > 0.9
            % we may have found a new point, call default processor