function [u,jac,hes,cub,du,djac,dhes]=hsu4(x,p)
d=p(1); h1=p(2); h2=p(3); theta=p(4); aa=p(5);  
%
a1=-aa*log(1-h1);
a2=-aa*log(1-h1-h2);
n1=x(4); n2=x(5); 
%
[f,fj,fh,fc,df,dfj,dfh]=hsf4(x,p(1:4));
df=df(:);

u=log(f)-a1*(n1-n2)-a2*n2;
if nargout > 1;
 jac=fj/f -[0 0 0 a1 a2-a1 0]';
end;
if nargout > 2;
  hes=fh/f-fj*fj'/(f*f);
end;

if nargout > 3;

  % cub=  vec( dhes(x,p)/dx )
 
  cb=zeros(6);
  for j=1:6;
    cb(:)=fc(:,j);
    cb=cb/f-fh*fj(j)/(f*f)-(fh(:,j)*fj'+fj*fh(j,:))/(f*f)+2*fj*fj'*fj(j)/f^3;
    cub(:,j)=cb(:);
  end;

end;

if nargout > 4;

  % du = du(x,p)/dp

  du=1/f*[df;0]-[0;aa*(n1-n2)/(1-h1)+aa*n2/(1-h1-h2);aa*n2/(1-h1-h2); ..
     0;(a1*(n1-n2)+a2*n2)/aa];
  
end;
  
if nargout > 5;

  % djac = djac(x,p)/dp

  tem1=zeros(6,5);
  tem1(4,2)=-aa/(1-h1); tem1(4,5)=log(1-h1);
  tem1(5,2)=-aa/(1-h1-h2)+aa/(1-h1); tem1(5,3)=-aa/(1-h1-h2); 
  tem1(5,5)=(a1-a2)/aa;
  djac=[dfj,zeros(6,1)]/f-fj*[df;0]'/(f*f)+tem1;
 
end;

if nargout > 6;

  % dhes = dhes(x,p)/dp

  df=[df;0]; 
  dfj=[dfj,zeros(6,1)];
  dfh=[dfh,zeros(36,1)];
  dh=zeros(6);
  dhes=zeros(36,5);
  for j=1:4;
    dh(:)=dfh(:,j);
    dh=dh/f-fh*df(j)/(f*f)-(dfj(:,j)*fj'+fj*dfj(:,j)')/(f*f)+2*fj*fj'*df(j)/f^3;
    dhes(:,j)=dh(:);
  end;
end;
