function [theta1res,theta2res,delta,Xi,Q] = lyngby_lzit_thopt(T, n, ...
    delta, betahat, theta1, theta2, Yi, xi, fhatvek);

% lyngby_lzit_thopt    - Lange-Zeger, Optimize theta parameters
%
%	function [theta1res,theta2res,delta,Xi,Q] = ...
%		lyngby_lzit_thopt(T, n, delta, betahat, theta1,
%		theta2, Yi, xi, fhatvek);
%
%	Input:	
%	
%	Output:	
%	
%	The function calculates the first order derivatives 
%	and perform a soft linesearch.
%
%       See also: lyngby_lzit_main

% cvs : $Id: lyngby_lzit_thopt.m,v 1.2 1997/10/28 10:31:41 fnielsen Exp $
%       $Revision: 1.2 $
	
    opt = 1;	% 1: first order	
                % 2: second order

    theta1res = theta1;
    theta2res = theta2;
    
    dLi = lyngby_lzit_ftlamb(T,n,theta1,theta2);
    Xi = dLi.*xi;
    Qlast = sum((Yi-Xi*betahat).*conj(Yi-Xi*betahat)./fhatvek); 

    w = [1:(T/(2*n))]'*n*2*pi/T;
    
    % Gradient
    dt1 = (-log(1+j*w/theta2)).*Xi;
    dt2 = ( (j*w*theta1)./ (theta2^2 + j*w*theta2) ).*Xi;
    gradt1 = -2*real(sum(conj(betahat)*conj(dt1).*(Yi-Xi*betahat)./fhatvek));
    gradt2 = -2*real(sum(conj(betahat)*conj(dt2).*(Yi-Xi*betahat)./fhatvek));
    

    if opt == 1
      Dtheta1 = -gradt1;
      Dtheta2 = -gradt2;
      
    elseif opt == 2
      % Hessian
      dt1t1 = (log((theta2+j)/theta2)^2).*Xi;
      dt1t2 = ((j*w.*(1-theta1*(log((theta2+j*w)/theta2))))./(theta2*(j*w+theta2))).*Xi;
      dt2t1 = dt1t2;
      dt2t2 = ((j*theta1*w.*(j*w*(theta1-1)-2*theta2))./((theta2*(j*w+theta2)).*(theta2*(j*w+theta2)))).*Xi;
      
      t11 = -2*real(sum(conj(betahat)*(-betahat.*dt1.*conj(dt1)+(Yi-Xi*betahat).*conj(dt1t1))./fhatvek));
      t12 = -2*real(sum(conj(betahat)*(-betahat.*dt2.*conj(dt1)+(Yi-Xi*betahat).*conj(dt1t2))./fhatvek));
      t21 = -2*real(sum(conj(betahat)*(-betahat.*dt1.*conj(dt2)+(Yi-Xi*betahat).*conj(dt2t1))./fhatvek));
      t22 = -2*real(sum(conj(betahat)*(-betahat.*dt2.*conj(dt2)+(Yi-Xi*betahat).*conj(dt2t2))./fhatvek));
      
      H = [t11 t12 ; t21 t22];
      
      % Avoid negative definit
      if min(eig(H))<0
	alfa = sqrt(-min(eig(H)));
	H = [t11+alfa,t12;t21 t22+alfa];
      end;
      
      stheta = - H \ [gradt1;gradt2];
      
      Dtheta1 = stheta(1,1);
      Dtheta2 = stheta(2,1);
    end
    
    theta1next = max([0 theta1 + delta*Dtheta1]);
    theta2next = max([0 theta2 + delta*Dtheta2]);
    dLi = lyngby_lzit_ftlamb(T,n,theta1next,theta2next);
    Xi = dLi.*xi;
    
    Q = sum((Yi-Xi*betahat).*conj(Yi-Xi*betahat)./fhatvek);
    diffQ = Q - Qlast;
    
    % Soft linesearch
    while ((diffQ>0) & (delta>10^(-16)))
      delta = delta*0.5;
      
      theta1next = max([0 theta1 + delta*Dtheta1]);
      theta2next = max([0 theta2 + delta*Dtheta2]);
      
      dLi = lyngby_lzit_ftlamb(T,n,theta1next,theta2next);
      Xi = dLi.*xi;
      Q = sum((Yi-Xi*betahat).*conj(Yi-Xi*betahat)./fhatvek);
      diffQ = Q-Qlast;
    end;
    
    % If costfunction improve, then update parametres else use old parametres
    if (diffQ<=0)
      theta1res = theta1next;
      theta2res = theta2next;
    end
    









