function [xp,fp,p,mxtake,iretcd]=lnsrch(n,x,f,g,p,func,stepmx,steptl,sx)
if any(imag(x)) | any(imag(p));
  disp('x,p entering in lnsrch:')
  [x,p]
end;
mxtake=0;
iretcd=2;
sln=sqrt(sum(sx.*sx.*p.*p));
if sln>stepmx;
  scl=stepmx/sln;
  p=scl*p;
  sln=stepmx;
end;
%slp=0;
%m=rem(n,5);
%if m;
%  slp=sum(g.*p);
%end;
%if (~m | (m & n>=5) ); 
%  ind1=[m+1:5:n]';
%  ind2=[m+2:5:n]';
%  ind3=[m+3:5:n]';
%  ind4=[m+4:5:n]';
%  slp=slp+sum(g(ind1).*p(ind1)+g(ind2).*p(ind2)+g(ind3).*p(ind3)+  ..
%              g(ind4).*p(ind4));
%end;
slp=sum(g.*p);
rln=max(abs(p)./max(abs(x),ones(n,1)./sx));
rmnlmb=steptl/rln;
lambda=1.0;
%
% LOOP
% CHECK IF NEW ITERATE SATISFACTORY.  GENERATE NEW LAMBDA IF NECESSARY.
%
while (iretcd>=2);
if any(imag(x)) | any(imag(p)) | imag(lambda);
  disp('x,p,lambda in lnsrch loop:')
  [x,p]
  lambda
end;

  xp=x+lambda*p;
  eval(['fp=',func,'(xp);'])
  if fp<=f+slp*1e-4*lambda;
    iretcd=0;
    if lambda==1.0 & sln>.99*stepmx;
      mxtake=1;
    end;
  else;
    if lambda<rmnlmb;
      iretcd=1;
    else;
      if lambda==1.0;
        tlmbda=-slp/(2*(fp-f-slp));
      else;
        t1=fp-f-lambda*slp;
        t2=pfp-f-plmbda*slp;
        t3=1/(lambda-plmbda);
        a=t3*(t1/(lambda*lambda)-t2/(plmbda*plmbda));
        b=t3*(t2*lambda/(plmbda*plmbda)-t1*plmbda/(lambda*lambda));
        disc=b*b-3*a*slp;
        if disc>b*b;
          tlmbda=(-b+sign(a)*sqrt(disc))/(3*a);
        else;
          tlmbda=(-b-sign(a)*sqrt(disc))/(3*a);
        end;
        if tlmbda>.5*lambda; tlmbda=.5*lambda; end;
      end;
      plmbda=lambda;
      pfp=fp;
      if tlmbda<lambda*.1;
        lambda=lambda*.1;
      else;
        lambda=tlmbda;
      end;
    end;
  end;
end;
