anchor.diff:         case(2)
             p%energy = normalized_B_nu(nu, src%temperature)
          case(3)
  -          p%dust_id = select_dust_specific_energy_rho(icell)
  -          p%emiss_var_id = jnu_var_id(icell%ic, p%dust_id)
  -          p%emiss_var_frac = jnu_var_frac(icell%ic, p%dust_id)
  +          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
   
  
positive.code:     type(photon),intent(inout) :: p
      ! the emitted photon
  
      type(grid_cell),intent(out) :: icell
  
      ! Sample position in map
      call grid_sample_pdf_map(src%luminosity_map, icell)
      call random_position_cell(icell, p%r)
  
      ! Sample isotropic angle
      call random_sphere_angle3d(p%a)
  
negative.code:           case(2)
               call random_planck_frequency(p%nu, src%spot(ispot)%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
  