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:       wp2=4d0*pi*n*ec**2/m
        omega0=ec*B/m/c
        xarg=thetae*sqrt(sqrt(2d0)*sin(theta)*(1d3*omega0/2d0/pi/nu))
  ! my slightly modified versions
        !eps11m22=jffunc(xarg)*wp2*omega0**2/(2d0*pi*nu)**4* &
        !(beselk(1d0/thetae,1)/beselk(1d0/thetae,2)+6d0*thetae)* &
        !sin(theta)**2
        !eps12=wp2*omega0/(2d0*pi*nu)**3* &
        !(beselk(1d0/thetae,0)-shgmfunc(xarg))/beselK(1d0/thetae,2)*cos(theta)
  
  ! with updated bessel functions
  !      eps11m22=jffunc(xarg)*wp2*omega0**2/(2d0*pi*nu)**4* &
  !     (besselk1(1d0/thetae)/besselk(2,1d0/thetae)+6d0*thetae)* &
  !      sin(theta)**2
  !      eps12=wp2*omega0/(2d0*pi*nu)**3* &
  !      (besselk0(1d0/thetae)-(sign(1d0,thetae+1d0)-1d0)/2d0*shgmfunc(xarg))/besselk(2,1d0/thetae)*cos(theta)
  
  ! UPDATING to fix NaNs for small temperatures thetae <~ 10^-2 by setting bessel ratios there = 1
        where(thetae.gt.1d-2)
  
negative.code:          ivx(xm)
        ju=0d0
  
        bnutnu=bnu(T,nu)
  
        ai=ji/bnutnu
  