anchor.diff:        end subroutine synchemisnoabs
   
         subroutine sympolemisth(nu,n,b,t,theta,e)
  +        use bessel, only: besselk0,besselk1,besselk
           use phys_constants, ec=>e
   ! Pandya+2016 fitting functions for thermal pol synchrotron
           real(kind=8), intent(in), dimension(:) :: nu,n,t,b,theta
           real(kind=8), intent(out), dimension(size(n),11) :: e
           real(kind=8), dimension(size(n)) :: thetae,nuc,xm,ji, &
                jq,jv,ju,bnutnu
  +        real(kind=8), dimension(size(n)) :: xarg,eps11m22,eps12, &
  +             tm,tp,targ,omega0,wp2,rhou,rhoq,rhov
           real(kind=8), dimension(size(n)) :: fac,jis,jqs,jvs,x
           real(kind=8) :: nucminval = 1d0
           real(kind=8) :: thetaemin = 1d-10
  
positive.code: 
        e=reshape((/ji,jq,ju,jv,ai,aq,au,av, &
        rhoq,rhou,rhov/),(/size(ji),11/),order=(/1,2/))
  !      testindx=maxloc(B,1)
  !      write(6,*) 'emis testindx: ',e(testindx,:)
  !      write(6,*) 'emis ang nu: ',theta(testindx),nu(testindx)
  !      write(6,*) 'emis n b t: ',n(testindx),thetae(testindx),b(testindx)
  !      write(6,*) 'ai: ',ai, e(:,5), e(:,1), ji
   !     write(6,*) 'reshape: ',maxval(emis(:,1)), maxval(emis(:,5))
        if(any(isnan(jq))) then
           write(6,*) 'NaN in polsynchemis.f90'
           !write(6,*) 'xm: ',xm
  
negative.code:          ivx(xm)
        ju=0d0
  
        bnutnu=bnu(T,nu)
  
        ai=ji/bnutnu
  