anchor.diff:                      end if
                   end do !q_pow
               end do !w_pow
  -            u_int_res(ri) = u_int_res(ri) * (2*mnuc*a)/(2*J+1)
  -            u_int_res(ri) = u_int_res(ri) * NAvo*tab_starrho(ri)*tab_mfr(ri,eli)/(mnuc*a)
  -            u_int_res(ri) = u_int_res(ri) * (hbar*c0)**2
  -            capped = capped + tab_r(ri)**2*u_int_res(ri)*tab_dr(ri)
  -
  -            ! capped = capped + u_int_res(ri) * tab_r(ri)**2*tab_dr(ri) * (2*mnuc*a)/(2*J+1) * NAvo*tab_starrho(ri)*tab_mfr(ri,eli)/(mnuc*a) * (hbar*c0)**2
  +            ! u_int_res(ri) = u_int_res(ri) * (2*mnuc*a)/(2*J+1)
  +            ! u_int_res(ri) = u_int_res(ri) * NAvo*tab_starrho(ri)*tab_mfr(ri,eli)/(mnuc*a)
  +            ! u_int_res(ri) = u_int_res(ri) * (hbar*c0)**2
  +            ! capped = capped + tab_r(ri)**2*u_int_res(ri)*tab_dr(ri)
  +
  +            factor_final = (2*mnuc*a)/(2*J+1) *NAvo*tab_starrho(ri)*tab_mfr(ri,eli)/(mnuc*a) *tab_r(ri)**2*tab_dr(ri) *(hbar*c0)**2
  +            capped = capped + u_int_res(ri) * factor_final
           end do !eli
       end do !ri
   
  
positive.code:     ! populate_array will place the non-zero value into a chosen slot at runtime
      coupling_Array = reshape((/0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, &
                                  0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0/), (/14, 2/))
  end subroutine captn_init_oper
  
  !   this is the integral over R in eqn 2.3 in 1501.03729
  
negative.code:     ! specific to captn_oper
      integer :: funcType, tau, taup, term_R, term_W, q_pow, w_pow ! loop indicies
      integer :: q_functype, q_index
      double precision :: J, j_chi, RFuncConst, WFuncConst, yConverse, mu_T, prefactor_functype, factor_final, prefactor_current
      double precision :: RD, RM, RMP2, RP1, RP2, RS1, RS1D, RS2 !R functions stored in their own source files
      double precision :: prefactor_array(niso,11,2)
  
  