anchor.diff:       ! emissivity of the dust types, binned to the image spectral resolution
   
        real(dp),allocatable :: tmp_spectrum(:)
  +     real(dp),allocatable :: tmp_stokes(:)
   
        character(len=10) :: track_origin
        logical :: uncertainties
  
positive.code:        io = 1
      end if
  
      stokes_vector = [p%s%i,p%s%q,p%s%u,p%s%v]
  
      if(inu >= 1 .and. inu <= img%n_nu) then
         if(img%compute_image) then
            call find_image_bin(img,x_image,y_image,ix,iy)
            if(ix >= 1 .and. ix <= img%n_x) then
               if(iy >= 1 .and. iy <= img%n_y) then
                  img%img(inu,ix,iy,im,io,:) = img%img(inu,ix,iy,im,io,:) + stokes_vector * p%energy
                  if(img%uncertainties) then
                     img%img2(inu,ix,iy,im,io,:) = img%img2(inu,ix,iy,im,io,:) + stokes_vector**2._dp * p%energy**2._dp
                     img%imgn(inu,ix,iy,im,io,:) = img%imgn(inu,ix,iy,im,io,:) + 1._dp
                  end if
               end if
  
negative.code:        if(img%compute_sed) then
            call find_sed_bin(img,x_image,y_image,ir)
            if(ir >= 1 .and. ir <= img%n_ap) then
               img%sed(inu,ir,im,io,:) = img%sed(inu,ir,im,io,:) + img%tmp_stokes * p%energy
               if(img%uncertainties) then
                  img%sed2(inu,ir,im,io,:) = img%sed2(inu,ir,im,io,:) + img%tmp_stokes**2._dp * p%energy**2._dp
                  img%sedn(inu,ir,im,io,:) = img%sedn(inu,ir,im,io,:) + 1._dp
               end if
            end if
  