function E=cell3(the,phi,aaa,bbb)
k=2*pi;
a=sin(the).*cos(phi);
b=sin(the).*sin(phi);
% if ((a==0)&&(b==0))
%     E=(1+cos(the)).*(4/k/k)*6/pi*6/pi;
% elseif a==0
%     E=(1+cos(the)).*(4/k/k)*6/pi.*sin(pi/6*b)./b;
% elseif b==0
%         E=(1+cos(the)).*(4/k/k)*6/pi.*sin(pi/6*a)./a;
% else
Fac=[1,aaa,1;aaa,bbb,aaa;1,aaa,1];
%Fac=[1,2,1;2,4,2;1,2,1];
PF=exp(i/3*pi*(-a-b))*Fac(1,1)+...
    exp(i/3*pi*(-b))*Fac(1,2)+...
    exp(i/3*pi*(a-b))*Fac(1,3)+...
    exp(i/3*pi*(-a))*Fac(2,1)+...
    ones(size(a))*Fac(2,2)+...
    exp(i/3*pi*(a))*Fac(2,3)+...
    exp(i/3*pi*(-a+b))*Fac(3,1)+...
    exp(i/3*pi*(b))*Fac(3,2)+...
    exp(i/3*pi*(a+b))*Fac(3,3);
    E=(1+cos(the))/2.*sin(pi/6*a)./a.*sin(pi/6*b)./b.*PF;
%    E=(1+cos(x))/(2).*sin(pi/6*sin(x).*cos(y))./(sin(x).*cos(y)).*sin(pi/6*sin(x).*sin(y))./sin(x)./sin(y);
    q=(abs(E.^2)).*sin(the);