anchor.diff:           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"
  
positive.code:          ivx(xm)
        ju=0d0
  
        !write(6,*) 'min/max B: ', minval(B),maxval(B)
        !write(6,*) 'min/max thetae: ', minval(thetae),maxval(thetae)
        !write(6,*) 'min/max sin(theta): ', minval(sin(theta)),maxval(sin(theta))
        !write(6,*) 'min/max n: ', minval(n),maxval(n)
        !write(6,*) 'min/max nu: ', minval(nu),maxval(nu)
        !write(6,*) 'min/max xm: ', minval(xm),maxval(xm)
  
  
  !      write(6,*) 'vars'
  !      write(6,'(3E9.4)') n,B,T
  !      write(*,*) 'nu: ',minval(nu),maxval(nu),minval(nuc),&
  !        maxval(nuc),minval(b),maxval(b),minval(theta),maxval(theta)
  !      write(*,*) 'xm: '
  !      write(6,'(4E9.4)') xm,iix(xm),iqx(xm),ivx(xm)
  !      write(6,*) 'j: ', jq/ji
        bnutnu=bnu(T,nu)
  
        ai=ji/bnutnu
  
negative.code:       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
  