function [Opttheta1,Optslope,Optbeta] = lyngby_lzgs_search(P, y, T,...
    mintheta1, maxtheta1, Notheta1,...
    minslope, maxslope, Noslope,...
    alpha, zoomlevel)

% lyngby_lzgs_search   - Lange Zeger main algorithm for grid search
%
%      function [Opttheta1,Optslope,Optbeta] =
%          lyngby_lzgs_search(P, y, T, ...
%          mintheta1, maxtheta1, Notheta1,...
%          minslope, maxslope, Noslope,...
%          alpha,zoomlevel)
%
%      See also LYNGBY, LYNGBY_LZGS_MAIN
% 
% $Id: lyngby_lzgs_search.m,v 1.3 2003/04/28 15:12:28 fnielsen Exp $


    if sum(abs(T)) == 0 
      Opttheta1 = 0;
      Optslope = 0;
      OptBeta = 0;
      return
    end

    Deltatheta1 = (maxtheta1-mintheta1) / Notheta1;
    theta1 = mintheta1 + Deltatheta1*(0:Notheta1-1);

    Deltaslope = (maxslope-minslope)/Noslope;
    slope = minslope + Deltaslope*(0:Noslope-1);

    sumy = sum(y);
    err = zeros(Notheta1,Noslope);
    for i = 1:Notheta1,
      for j = 1:Noslope,
	theta2 = (theta1(i)-1)/slope(j);
	Lambda = lyngby_lz_lambda(theta1(i),theta2,T,1);
	Zun = conv(Lambda,P);
	Zun = Zun(1:length(y));
	Beta(i,j)=(Zun'*y)/(Zun'*Zun);    
	err(i,j)=sum((y-Beta(i,j)*Zun).^2)+alpha*Beta(i,j)^2;
      end
    end


    % Minimum and its index
    minerr = min(min(err));
    [ithetaar, islopear] = find(err==minerr);
    
    %iindex = find(err==minerr);
    %ithetaar = mod(iindex(1),size(err,1));
    %islopear = floor(iindex(1)/size(err,1))+1;
    
    
    % Matlab 5 would do
    %[itheta1ar,islopear]=ind2sub(size(err),find(err==minerr));
    %disp(sprintf('Found %i Minima',length(itheta1ar)));
    
    itheta1 = ithetaar(1);
    islope = islopear(1);
    
    Opttheta1 = theta1(itheta1);
    Optslope = slope(islope);
    Optbeta = Beta(itheta1,islope);
    
    if zoomlevel>1
      [Opttheta1,Optslope,Optbeta] = lyngby_lzgs_search(P,y,T,...
	  max(1.5, Opttheta1-Deltatheta1/Notheta1/2),...
	  max(1.5+Deltatheta1, Opttheta1+Deltatheta1/Notheta1/2),...
	  Notheta1,...
	  max(0.5, Optslope-Deltaslope/Noslope/2),...
	  max(0.5, Optslope+Deltaslope/Noslope/2),...
	  Noslope,...
	  alpha, zoomlevel-1);
    end
    












