function [alfa,x] = StepSize(fun, x, d, alfa, params)

x0 = x.p;
Dphi0 = x.g'*d;
if ( (alfa <= 0) | (Dphi0 > 0) )
  error('Initialization of step incorrect');
end;

phi0 = x.f;
alfaL = 0;
alfaR = inf;
iparams = struct('fl',phi0,'gl',Dphi0,'safetol',5e-2,'prev',alfaR-alfaL);

for iter=1:params.maxit
  x.p = x0 + alfa*d;
  x.f = feval(fun,x.p,1);
  
  if (x.f >= phi0 + alfa*params.c1*Dphi0) 
    alfaR = alfa;
    iparams.prev2 = iparams.prev; iparams.prev = alfaR-alfaL; 
    iparams.fr = x.f;
    alfa = Interp(alfaL, alfaR, iparams);
  else
    x.g = feval(fun,x.p,2);
    DphiAlfa = x.g'*d;
    
    if (DphiAlfa >= params.c2*Dphi0)
      if (DphiAlfa <= params.c2*abs(Dphi0))
        return;
      else
        alfaR = alfa;
        iparams.fr = x.f;
      end
    else
      alfaL = alfa;
      iparams.fl = x.f;
      iparams.gl = DphiAlfa;
    end
    iparams.prev2 = iparams.prev; iparams.prev = alfaR-alfaL;
    if isinf(alfaR)
      alfa = Extrap(alfaL);
    else
      alfa = Interp( alfaL, alfaR, iparams);
    end	
  end
end
fprintf('step size criteria were not met\n');
fprintf('after %d step size iterations.\n', iter);
error('STEP SIZE FAILURE');
return;

function mid = Interp(left, right, iparams)
diff = right-left;
a = (iparams.fr-iparams.fl - iparams.gl*diff)/(diff^2);
if (a <= 1e-6)
  mid = (left + right)/2.0;
  return;
  error('non convex quadratic: impossible');
end
adjust = -iparams.gl/(2*a);
safetol = iparams.safetol*diff;
if (adjust > safetol)
  mid = left + adjust;
  if (mid > right - safetol || max(right-mid,mid-left) > 0.66*iparams.prev2 )
    mid = (left + right)/2.0;
  end
else
  mid = (left + right)/2.0;
end
return;

function mid = Extrap(left)
mid = 2*left;
return;