function estimated_location = newton_optimize( ...
    init_t, max_interval_width, min_fct_val, max_fct_val, ...
    two_sided, fdobj, deriv_degree, error_level )

    % Based on an initial estimate, looks for a zero crossing with
    % arbitrary precision
    %
    % Input parameters
    % ----------------
    % init_t              Initial estimate of the zero
    % max_interval_width  The zero has to be within a certain range. 
    % [min_fct_val, max_fct_val]
    %                     Range over which the function is defined
    % two_sided           Depending on two_sided, the desired zero must be
    %                     within: 
    %                     two_sided = true: current_xvalue +/- max_interval_width/2
    %                     two_sided = false: current_xvalue - max_interval_width
    %                     if the estimate falls outside this range,
    %                     we got too far, and the system will restart
    % fdobj               Functional data object that represents the function
    % deriv_degree        Derivative number of fdobj that where we are looking
    %                     for zeros (e.g., if deriv_degree = 2, we are looking
    %                     for zero crossings in the accelertion of dfobj
    % error_level         Required precision of the estimate. Success is reached
    %                     when the y value is within 0 +/- error_level
    %
    % Output value
    % ------------
    % estimated_location  Estimate of the zero crossing location within
    global_constants;


    pause_on_warnings = false; % to investigate problems, set to true so that the 
                               % program pauses after showing a warning 
                               % warnings are usually not problematic -- they
                               % occur during normal processing. 
                               % so as default, leave this value to false
                               
    display_warnings = true;

    learning_rate = 0.5;  % indicates how fast we allow the system to converge
                          % large learning rates can be faster, but might
                          % also make the system unstable. If an unstable
                          % system is dectected, learning rate is
                          % automatically decreased 
                          % Give a value between 0 and 1 (default is 0.5)

    restarts      = 0;  % number of times optimization failed (due to instability) 
                        % and had to be restarted

    iter          = 0;  % number of iterations in the loop. Used to detect infinite loops

    max_iter      = 25;  % maximum of iterations allowed
    max_restarts  = 50;  % maximum number of restarts allowed

    % We expect the real zero to be in the [min_interval_value, max_interval_value]
    % interval (the initial estimate should be in this range or at its boundary)
    if (two_sided == true)
        % the real zero can be on either side of the initial estimate
        max_interval_value = init_t + max_interval_width/2;
        min_interval_value = init_t - max_interval_width/2;
    else
        % the real zero must be lower than the current estimate
        max_interval_value = init_t;
        min_interval_value = init_t - max_interval_width;
    end

    % now make sure we don't exceed the range defined for the function
    min_interval_value = max(min_interval_value, min_fct_val);
    max_interval_value = min(max_interval_value, max_fct_val);
    
    estimated_min_val = eval_fd(min_interval_value, fdobj, deriv_degree);
    estimated_max_val = eval_fd(max_interval_value, fdobj, deriv_degree);
    

    % Initializing the estimate of location and function value
    estimated_location = init_t;  % initialize current estimate of x location of the zero crossing
    estimated_y        = eval_fd(estimated_location, fdobj, deriv_degree);

    if (verbose)
        disp(['new optimization -- Looking for a zero crossing (deriv = ', ...
              num2str(deriv_degree), ') between (', ...
               num2str(min_interval_value), ', ' num2str(estimated_min_val), ') and (',...
               num2str(max_interval_value), ', ' num2str(estimated_max_val), ')'])
    end
    
    if (sign(estimated_min_val)*sign(estimated_max_val) ~= -1)
        error('Function does not change sign!!!')
    end
    
%     disp(['   Initial estimate: X = ', ...
%            num2str(estimated_location), '    Y = ', ...
%            num2str(estimated_y)])

    if (init_t < min_interval_value || init_t > max_interval_value)
        error('Initial x value is inconsistent with the interval where the zero crossing is expected')
    end


    % Loop until precision criterion is reached (i.e., error is low enough)
    while (abs(estimated_y) > error_level)
        iter = iter+1;

        % calculate the slope at the current estimated location
        slope = eval_fd(estimated_location, fdobj, deriv_degree+1);

        % if the slope is in the right direction (it usually is when we are 
        % near enough to the actual zero crossing), we can compute a better
        % estimation as follows:
        new_estimated_location = estimated_location - learning_rate * estimated_y / slope;

        % But if the function is highly non-linear near the zero crossing, 
        % the slope sends us in the wrong direction, or very far off.
        % We test that the new estimate is within the range allowed
        if ((new_estimated_location < min_interval_value) || ...
            (new_estimated_location > max_interval_value))

            if (display_warnings)
                disp(['*** Restart : ' num2str(restarts), ...
                      '   Reason: New estimate (X = ', ...
                      num2str(new_estimated_location), ') fell outside the allowed range'])
            end

            % If the value is too far off, we restart with a random
            % value within the interval allowed
            new_estimated_location = min_interval_value + rand() * max_interval_width;

            % make sure this new value is within range
            new_estimated_location = max(new_estimated_location, min_fct_val);
            new_estimated_location = min(new_estimated_location, max_fct_val);

            % Reduce learning rate if it is not already low
            % this helps stabilize learning
            if (learning_rate > 0.2)
                learning_rate = 0.99 * learning_rate;
            end

            % Finally, it could also be that the interval defined is too small
            % so enlarge it slightly
            min_interval_value = min_interval_value - 0.01 * max_interval_width;
            max_interval_value = max_interval_value + 0.01 * max_interval_width;

            min_interval_value = max(min_interval_value, min_fct_val);
            max_interval_value = min(max_interval_value, max_fct_val);

            % Update the interval width
            max_interval_width = max_interval_value - min_interval_value;

            % ok, so let's restart with these new parameters
            if (display_warnings)
                disp(['    New estimated location = ' num2str(new_estimated_location), ...
                     ' with range [', num2str(min_interval_value), ...
                     ', ', num2str(max_interval_value), ']'])
            end
            restarts = restarts + 1;

            if(pause_on_warnings)
                pause 
            end
        end

        % if the max number of iterations was reached, also restart
        if (iter == max_iter)
            new_estimated_location = min_interval_value + rand() * max_interval_width;

            % make sure this new value is within range
            new_estimated_location = max(new_estimated_location, min_fct_val);
            new_estimated_location = min(new_estimated_location, max_fct_val);

            if (display_warnings)
                disp(['*** Restart : ' num2str(restarts), ...
                      '   Reason: Could not find solution in ', ...
                      num2str(max_iter), ' iterations'])
            end
            iter = 0;
            restarts = restarts+1;

            if(pause_on_warnings)
                pause 
            end
        end

        % exit when optimization has failed too many times
        if (restarts >= max_restarts)
            if (die_on_error)
                % report an error and stop
                error('Optimization failed');
            else
                % return the initial estimate - no optimization done
                % but let the program continue
                disp('Optimization failed - No optimization performed');
                estimated_location = init_t;
                return;
            end
        end

        % update x and y values for the next round
        estimated_location = new_estimated_location;
        estimated_y        = eval_fd(estimated_location, fdobj, deriv_degree);
    end

    if (verbose)
        disp('Optimization completed -- Solution found');
    end

    % at this point, estimated_location gives a y value within error level
return;