clear
[theta,phi]=meshgrid(0:0.02*pi:2*pi,0:pi/12:1.1*pi);
the=pi/2-theta;
q=zeros(201);
for kk=1:201
    aaa=-20+(kk-1)/5;
    bbb=-3.0926*aaa-2.7998;
F=@(the,phi)cell5(the,phi,aaa,bbb).^2.*sin(the);
A=dblquad(F,0,pi,0,2*pi);
B=max(max(cell5(the,phi,aaa,bbb)));
q(kk)=4*pi/A*B*B;
log10(q)*10;
kk
end