function Fnu = flux3(nth,nphi,Inu)
% Flux from non-spherical star.
%
% Preliminary version.  Appears to work.
% This version to interpolate in model fluxes.
% To use: Fnu = flux3(nth,nphi,Inu);
%
% by D. Holmgren, Jan 18 95 holmgren@phobos.astro.uwo.ca

% stellar radius for Collins' example...
Rp = 2.7851e09;
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); u = 0.6;
sum = 0.; ci = cos(incl); si = sin(incl); w2 = w*w;
for k = 1:nth
 cth = cos(theta(k)); sth = sin(theta(k));
 for l = 1:nphi
  cph = cos(phi(l)); sph = sin(phi(l));
  cths = cth*si + sth*cph*ci;
  ths = acos(cths); sths =sin(ths); thd = ths*180/pi;
  xr = frot2(w,thd); x3 = xr^3;
% note that [g] = cgs units...
  g = grot(w,xr,ths); logg = log10(g); g1= g/100;
  Teff = tempr(g1,w);
%  disp(Teff)
  Ilocal = table2(Inu,logg,Teff); Ilocal = 10^Ilocal;
  tphs = sth*sph/(sth*cph*ci-cth*ci);
  phs = atan(tphs); cphs = cos(phs);
  a = 8*w2*x3*sths*cths/(27-8*w2*x3*sths*sths);
  mgr = sths*ci*cphs+cths*ci - a*(si*cphs*cths-sths*si);
  sum=sum + Ilocal*xr*xr*abs(mgr)*sth;
 end
end
% compute flux...
Fnu = Rp^2*sum*dth*dphi;