Verify the solid angle ratio is calculated correctly for star and disk

[gully]

Right now star_BB.py computes the relative flux ratio of the star to disk like so:

F_bol1 = self.F_bol_interp.interp(p.grid)
F_bol2 = stef_boltz * p.T_BB**4.0 ## Radiance from black body
self.qq = F_bol2/F_bol1[0]

and then computes the net model like so:

model1 = self.Omega * (self.chebyshevSpectrum.k * self.flux_mean + X.dot(self.mus))
# Black body
model2 = self.Omega2 * self.qq * self.BB_lam.value * self.chebyshevSpectrum.k

# Return a "metadata blob" of both spectrum models
raw_models = [model1, model2]

net_model = model1 + model2

The absolute values of the solid angles have no meaning since we do not have absolutely flux-calibrated spectra. However the relative value of the solid angles is physically meaningful. How much bigger is the disk than the star (since they share the same distance).

However, the posterior values I'm getting for logOmega2 are much less than those for logOmega, suggesting that something is wrong in the co-factors I'm assuming for the relative flux ratio of the Phoenix models and Black Body.

TODO's:

• Assess the units of Phoenix and astropy's Black Body.
• Revisit the F_bol assumptions

[gully]

In order to make sure the solid angles for star and disk are accurate relative to each other, the flux units on the black body and stellar photosphere model grid have to be the same.

This is a bit of a subtle point, and ties back to Starfish Issue #38. In Issue 38 I computed the ratios of the bolometric flux of two stellar models, interpolating in-between grid points. I'm doing the same thing here, except now the second spectrum is simply a black body.

In grid_tools.py around line 258:

if self.norm:
    f *= 1e-8 #convert from erg/cm^2/s/cm to erg/cm^2/s/A
    F_bol = trapz(f, self.wl_full)
    f = f * (C.F_sun / F_bol) #bolometric luminosity is always 1 L_sun

If the models were absolutely flux calibrated, we would simply do:

# Stellar photosphere
model1 = self.Omega * (self.chebyshevSpectrum.k * self.flux_mean + X.dot(self.mus))
# Black body
model2 = self.Omega2 * self.BB_lam.value * self.chebyshevSpectrum.k

But they're not, so to get back to absolute units, we should adapt the models like this:

# Stellar photosphere
model1 = self.Omega * (self.chebyshevSpectrum.k * self.flux_mean + X.dot(self.mus)) / (C.F_sun / F_bol)
# Black body
model2 = self.Omega2 * self.BB_lam.value * self.chebyshevSpectrum.k * np.pi

The extra factor 1/(C.F_sun / F_bol) on model1 comes from de-normalizing per above, and the factor of pi on model2 comes from integrating over half the sky on the stellar surface.
Finally, we can multiply both models by the extra factor:

# Stellar photosphere
model1 = self.Omega * (self.chebyshevSpectrum.k * self.flux_mean + X.dot(self.mus))
# Black body
model2 = self.Omega2 * (self.BB_lam.value * np.pi) * self.chebyshevSpectrum.k  * (C.F_sun / F_bol)

So that in the limit Omega2 -> 0, we get the same fitting behavior as Starfish Classic. (In the actual implementation I multiply in the factor of pi earlier.)

[gully]

This has not been settled!
The emulator has to be updated.