anchor.diff:      S%U = - sin_2_beta * RLS2 + cos_2_beta * RLS3
       S%V = RLS4
   
  -    a = a_final
  -
     end subroutine scatter_stokes
   
     elemental real(dp) function B_nu(nu,T)
  
positive.code:   !
    !#############################################################################
  
    subroutine scatter_stokes(s,a,a_scat,P1,P2,P3,P4)
  
      implicit none
  
      type(angle3d_dp),intent(inout) :: a           ! The photon direction angle
      type(angle3d_dp),intent(in)    :: a_scat      ! The photon scattering angle
      type(stokes_dp),intent(inout)  :: s           ! The Stokes parameters of the photon
      real(dp)                       :: P1,P2,P3,P4 ! 4-element matrix elements
  
      ! Local
      type(angle3d_dp) :: a_original ! Original direction
      type(angle3d_dp) :: a_final    ! Final direction
      real(dp) :: cos_i2,cos_2_i2
      real(dp) :: sin_i2,sin_2_i2
      real(dp) :: cos_2_alpha,cos_2_beta
      real(dp) :: sin_2_alpha,sin_2_beta
      real(dp) :: RLS1,RLS2,RLS3,RLS4
  
      a_original = a
  
      call rotate_angle3d(a_scat,a_original,a_final,cos_big_a=cos_i2,sin_big_a=sin_i2)
  
      cos_2_alpha = 1._dp - 2._dp * a_scat%sinp * a_scat%sinp
      sin_2_alpha =       - 2._dp * a_scat%sinp * a_scat%cosp
  
      cos_2_i2    = 1._dp - 2._dp * sin_i2 * sin_i2
      sin_2_i2    =         2._dp * sin_i2 * cos_i2
  
      if(a_scat%sinp < 0.) then
         cos_2_beta =  cos_2_i2
  
negative.code:     real(dp) :: nu
  
      type(angle3d_dp) :: a_scat
      type(angle3d_dp) :: a_final
  
      real(dp) :: P1,P2,P3,P4,norm
  
  