function F = flux2(nth,nphi)
% Flux from a limb darkened nonspherical star.
%
% D. Holmgren holmgren@phobos.astro.uwo.ca
b = zeros(3,11);
b(1,1) = 1.0001956; b(3,2) = -1.183952e-09; b(2,3) = 4.2098469e-05;
b(1,10) = -0.080843068; b(2,10) = 2.5709925e-04; b(3,10) = -5.598946e-08;
b(1,11) = 0.080509197; b(2,11) = -2.545361e-04; b(3,11) = 5.9308193e-08;
incl = input('Enter incl: ');
incl = incl*pi/180.;
w = input('Enter w: ');
theta = linspace(0,pi,nth)';
phi = linspace(-pi/2,pi/2,nphi)';
dth = theta(2)-theta(1); dphi = phi(2)-phi(1);
sum = 0.; u = 0.6; ci = cos(incl); si = sin(incl);
for k = 1:nth
  sth = sin(theta(k)); cth = cos(theta(k));
  for l = 1:nphi
    cph = cos(phi(l));
    cth1=cth*si+sth*cph*ci;
    theta1 = acos(cth1);
    xr = 0.;
    for i = 1:3
     for j = 1:11
      xr=xr+b(i,j)*(w^j)*theta(k)^(2*i);
     end
    end
    if w == 0.
     xr = 1.;
    end
%    tph1=sth*sin(phi(k))/(sth*cph*si-cth*ci);
%    phi1 = atan(tph1);
    a = 8*w^2*xr^3*cth*sth/(27-8*w^2*xr^3*sth*sth);
    mgr = (sth*si*cph+cth*ci)-a*(si*cph*cth-sth*si); mgr = abs(mgr);
    sum = sum + (1 - u + u*mgr)*xr^2*mgr*sth;
  end
end
F = sum*dth*dphi;