anchor.diff:      function OMEGA(rindex,w)
       double precision :: sigma_N, GF,vesc,Omega,mu,muplus,muminus,u,w
       integer i, rindex
  +
       vesc = tab_vesc(rindex)
       u = sqrt(w**2-vesc**2)
       Omega = 0.d0
  
positive.code:     double precision :: x, dummyf
      dummyf = 1.d0
      end function dummyf
  
      subroutine captn_general(mx_in,sigma_0_in,niso_in,nq_in,nv_in,capped)
      use capmod
  
negative.code:     result = 0.d0
      ri_for_omega = ri !accessed via the module
      !call integrator
      call dsntdqagse(integrand,dummyf,1.d0,vesc_halo/3., &
      epsabs,epsrel,limit,result,abserr,neval,ier,alist,blist,rlist,elist,iord,last)
  
      u_int_res(ri) = result*sigma_0
  
      capped = capped + tab_r(ri)**2*u_int_res(ri)*tab_dr(ri)
      ! print*,"Capgen capped ", capped, tab_r(ri), u_int_res(ri), tab_dr(ri)
      end do
      capped = 4.d0*pi*Rsun**3*capped
  
  