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:          eps12=wp2*omega0/(2d0*pi*nu)**3*cos(theta)
        end where
  
  ! s08 versions
  !      eps11m22=shffunc(xarg)*wp2*omega0**2/(2d0*pi*nu)**4* &
  !     (beselk(1d0/thetae,1)/beselk(1d0/thetae,2)+6d0*thetae)* &
  !      sin(theta)**2
  !      eps12=shgfunc(xarg)*wp2*omega0/(2d0*pi*nu)**3* &
  !      beselk(1d0/thetae,0)/beselK(1d0/thetae,2)*cos(theta)
  !      write(*,*) 'eps: ',eps11m22, eps12
  ! high-frequency limits
  !      eps11m22=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)/beselK(1d0/thetae,2)*cos(theta)
  !      targ=sqrt(4d0*eps12**2+eps11m22**2)
  !      tp=-(eps11m22-targ)/2d0/eps12
  !      tm=-(eps11m22+targ)/2d0/eps12
  
        rhov=2d0*pi*nu/c*eps12
        rhoq=2d0*pi*nu/2d0/c*eps11m22
  
  
  ! confusion about sign of rhoq. Shcherbakov (2008) disagrees with Huang & Shcherbakov (2011). Statement is different basis vectors. Basis looks identical between SH10, which I follow, and HS11. So I did signs as in Huang & Shcherbakov.
  ! limiting case for rhoq for small \nu / \nu_T where S08 approximations break down:
  !      rhoqlim=2d0**(2d0/3d0)/3d0**(7d0/2d0)*pi*n*ec**2d0/m/c/thetae**3/nuc* &
  !        xm**(-5d0/3d0)
  ! this is implemented directly in shffunc now, see 12/9/2014 notes or idl polsynchemis.pro
  ! test unpolarized
  !AC 
  !       jq=0.; jv=0.; aq=0.; av=0.; rhoq=0.; rhov=0.
  ! w/o Faraday:
  !      rhoq=0.!; rhov=0
  !      rhov=0.
  ! w/o total absorption:
  !      ai=0.
  ! w/o coupled absorption:
  !      aq=0.; av=0.
  
        if(any(isnan(jq))) then 
           write(6,*) "nan in jq"
  
negative.code:         shgmfunc=0.43793091*log(1d0+0.00185777*x**1.50316886)
          end function shgmfunc
  
  !        function shgfunc(x)
  !        ! Fitting function G(X) from Shcherbakov (2008)
  !        ! JAD 1/29/2010
  !        real(kind=8), intent(in), dimension(:) :: x
  !        real(kind=8), dimension(size(x)) :: shgfunc
  !        shgfunc=1d0-0.11d0*log(1d0+.035d0*x)
  !        end function shgfunc
  
          function iqx(x)
          ! Fitting function for polarized synchrotron emissivity function
  