Fix spot spectra for stars with spots

hyperion/model/tests/test_spot_source.py

@@ -0,0 +1,70 @@
+from __future__ import print_function, division
+
+import os
+import tempfile
+
+import numpy as np
+import pytest
+
+from .. import Model
+from ...util.functions import random_id
+
+
+@pytest.mark.requires_hyperion_binaries
+def test_spot_uses_its_own_spectrum():
+    # Regression test: photons emitted from a spot must use the spot's own
+    # spectrum, not the parent sphere's. We give the sphere and the spot
+    # disjoint emission bands (no dust, so photons escape directly), so any
+    # flux in the spot's band can only come from the spot having used its own
+    # spectrum. Before the fix, source_emit overwrote the spot-sampled
+    # frequency with the sphere spectrum, so the spot band got zero flux.
+
+    m = Model()
+    m.set_cartesian_grid([-1e12, 1e12], [-1e12, 1e12], [-1e12, 1e12])
+
+    nu = np.logspace(np.log10(3e12), np.log10(1e15), 300)
+    fnu_sphere = np.where((nu > 1e13) & (nu < 2e13), 1., 0.)   # ~15-30 micron
+    fnu_spot = np.where((nu > 3e14) & (nu < 6e14), 1., 0.)     # ~0.5-1 micron (disjoint)
+
+    s = m.add_spherical_source()
+    s.radius = 1e11
+    s.position = (0., 0., 0.)
+    s.luminosity = 1.
+    s.spectrum = (nu, fnu_sphere)
+
+    spot = s.add_spot()
+    spot.longitude = 0.
+    spot.latitude = 0.
+    spot.radius = s.radius
+    spot.luminosity = 1.
+    spot.spectrum = (nu, fnu_spot)
+
+    sed = m.add_peeled_images(sed=True, image=False)
+    sed.set_viewing_angles([45.], [45.])
+    sed.set_wavelength_range(60, 0.1, 100.)
+    sed.set_aperture_radii(1, 1e12, 1e12)
+
+    m.set_n_initial_iterations(0)
+    m.set_n_photons(imaging=20000)
+
+    tmpdir = tempfile.mkdtemp()
+    m.write(os.path.join(tmpdir, random_id()))
+    mo = m.run()
+
+    wav, nufnu = mo.get_sed()
+    wav = np.array(wav)
+    nufnu = np.squeeze(np.array(nufnu))
+
+    sphere_band = (wav > 10.) & (wav < 40.)    # ~15-30 micron
+    spot_band = (wav > 0.4) & (wav < 1.2)      # ~0.5-1 micron
+
+    sphere_flux = np.nansum(nufnu[sphere_band])
+    spot_flux = np.nansum(nufnu[spot_band])
+
+    # the model actually ran and the sphere emitted in its band
+    assert sphere_flux > 0.
+    # the spot emitted in its OWN band (exactly zero before the fix)
+    assert spot_flux > 0.
+    assert spot_flux > 0.2 * sphere_flux
+    # and essentially nothing leaked between the two disjoint bands
+    assert spot_flux + sphere_flux > 0.9 * np.nansum(nufnu)

src/sources/source_type.f90

@@ -177,7 +177,7 @@ subroutine source_read(group, s)
           s%spot(i)%cost = cos(spot_size * deg2rad)
 
           g_spot = mp_open_group(group, spot_names(i))
-          call set_spectrum(group, s%spot(i)%freq_type, s%spot(i)%spectrum, s%spot(i)%temperature)
+          call set_spectrum(g_spot, s%spot(i)%freq_type, s%spot(i)%spectrum, s%spot(i)%temperature)
           call mp_close_group(g_spot)
 
           if(s%freq_type == 3) call error("source_read", "Spot cannot have LTE spectrum")
@@ -456,19 +456,22 @@ subroutine source_emit(src,p,nu)
              p%emiss_var_frac = jnu_var_frac(p%icell%ic, p%dust_id)
              call dust_sample_emit_probability(d(p%dust_id),p%emiss_var_id,p%emiss_var_frac,nu,p%energy)
           end select
-       end if
 
-       select case(src%freq_type)
-       case(1)
-          p%energy = interpolate_pdf(src%spectrum, nu, bounds_error=.false., fill_value=0._dp)
-       case(2)
-          p%energy = normalized_B_nu(nu, src%temperature)
-       case(3)
-          p%dust_id = select_dust_specific_energy_rho(p%icell)
-          p%emiss_var_id = jnu_var_id(p%icell%ic, p%dust_id)
-          p%emiss_var_frac = jnu_var_frac(p%icell%ic, p%dust_id)
-          call dust_sample_emit_probability(d(p%dust_id),p%emiss_var_id,p%emiss_var_frac,nu,p%energy)
-       end select
+       else
+
+          select case(src%freq_type)
+          case(1)
+             p%energy = interpolate_pdf(src%spectrum, nu, bounds_error=.false., fill_value=0._dp)
+          case(2)
+             p%energy = normalized_B_nu(nu, src%temperature)
+          case(3)
+             p%dust_id = select_dust_specific_energy_rho(p%icell)
+             p%emiss_var_id = jnu_var_id(p%icell%ic, p%dust_id)
+             p%emiss_var_frac = jnu_var_frac(p%icell%ic, p%dust_id)
+             call dust_sample_emit_probability(d(p%dust_id),p%emiss_var_id,p%emiss_var_frac,nu,p%energy)
+          end select
+
+       end if
 
     else
 
@@ -486,19 +489,22 @@ subroutine source_emit(src,p,nu)
              p%emiss_var_frac = jnu_var_frac(p%icell%ic, p%dust_id)
              call dust_sample_j_nu(d(p%dust_id),p%emiss_var_id,p%emiss_var_frac,p%nu)
           end select
-       end if
 
-       select case(src%freq_type)
-       case(1)
-          p%nu = sample_pdf(src%spectrum)
-       case(2)
-          call random_planck_frequency(p%nu, src%temperature)
-       case(3)
-          p%dust_id = select_dust_specific_energy_rho(p%icell)
-          p%emiss_var_id = jnu_var_id(p%icell%ic, p%dust_id)
-          p%emiss_var_frac = jnu_var_frac(p%icell%ic, p%dust_id)
-          call dust_sample_j_nu(d(p%dust_id),p%emiss_var_id,p%emiss_var_frac,p%nu)
-       end select
+       else
+
+          select case(src%freq_type)
+          case(1)
+             p%nu = sample_pdf(src%spectrum)
+          case(2)
+             call random_planck_frequency(p%nu, src%temperature)
+          case(3)
+             p%dust_id = select_dust_specific_energy_rho(p%icell)
+             p%emiss_var_id = jnu_var_id(p%icell%ic, p%dust_id)
+             p%emiss_var_frac = jnu_var_frac(p%icell%ic, p%dust_id)
+             call dust_sample_j_nu(d(p%dust_id),p%emiss_var_id,p%emiss_var_frac,p%nu)
+          end select
+
+       end if
 
     end if