anchor.diff:  end function integrand_oper
   
   
  -!   this is the integral over R in eqn 2.3 in 1501.03729
  -!THIS IS THE IMPORTANT FUNCTION: the integrand for the integral over u
  -function integrand_oper_extrawterm(u,foveru)
  -    use opermod
  -    double precision :: u, w, integrand_oper_extrawterm, foveru
  -    external foveru
  -
  -    w = sqrt(u**2+vesc_shared_arr(rindex_shared)**2)
  -
  -    !Switch depending on whether we are capturing on Hydrogen or not
  -    if (a_shared .gt. 2.d0) then
  -        integrand_oper_extrawterm = foveru(u)*GFFI_A_oper(w,vesc_shared_arr(rindex_shared),a_shared,q_shared+1)
  -    else
  -        integrand_oper_extrawterm = foveru(u)*GFFI_H_oper(w,vesc_shared_arr(rindex_shared),q_shared+1)
  -    end if
  -end function integrand_oper_extrawterm
  -
  -
   ! call captn_oper to run capt'n with the effective operator method
   subroutine captn_oper(mx_in, jx_in, niso, capped)
       use opermod
  
positive.code:     end if
      allocate(u_int_res(nlines))
  
      !As far as I can tell, the second argument (fofuoveru) does nothing in this integrator.
      !I've sent it to an inert dummy just in case.
      
      
      capped = 0.d0
      ! Bottom part of the integral is always zero -- happy little slow DM particles can always be captured.
      umin = 0.d0
      ! the current response function type in order: M, S2, S1, P2, MP2, P1, D, S1D
      do funcType = 1,8
          ! the first index on each response function
          do tau=1,2
              ! the second index on each response function
              do taup=1,2
                  ! the possible terms for each R function in order: c, v2, q2, v2q2, q4, v2q4
                  do term_R = 1,6
                      ! pick out the appropriate term's constant from a given R function of tau, taup, and term_R
                      ! currently passes mnuc, and c0 - these are constants that could be shared to it through the shared module?
                      select case (funcType)
                      case (1)
                          RFuncConst = RM(mnuc,c0,tau,taup,term_R-1,j_chi,coupling_Array) !!!!!!!!!!!!!!! in the R functions the R term starts at zero, should change it to start at 1 like other Fortran things do for consistency
                      case (2)
                          RFuncConst = RS2(mnuc,c0,tau,taup,term_R-1,j_chi,coupling_Array)
                      case (3)
                          RFuncConst = RS1(mnuc,c0,tau,taup,term_R-1,j_chi,coupling_Array)
                      case (4)
                          RFuncConst = RP2(mnuc,tau,taup,term_R-1,j_chi,coupling_Array)
                      case (5)
                          RFuncConst = RMP2(mnuc,tau,taup,term_R-1,j_chi,coupling_Array)
                      case (6)
                          RFuncConst = RP1(mnuc,tau,taup,term_R-1,j_chi,coupling_Array)
                      case (7)
                          RFuncConst = RD(mnuc,tau,taup,term_R-1,j_chi,coupling_Array)
                      case (8)
                          RFuncConst = RS1D(tau,taup,term_R-1,j_chi,coupling_Array)
                      case default
                          RFuncConst = 0.
                          print*, "Um, I ran out of R functions to choose from?"
                      end select
  
                      !!!!!!!!!!!!!!!!!!!!!! if I check here to see if RFuncConst is zero, I could skip the integrator call, no?
                      !!!!!!!!!!!!!!!!!!!!!! what ever pops out of the integrator is just multiplied by this term, so if it's zero it'll kill the integrator's result off (and move on to the next iteration of the loop)
                      if (RFuncConst.ne.0.) then
                          qOffset = 0 ! qOffset is used to tell the GFFI how many extra q^2 terms come from the current R function term
                          if (funcType.gt.3) then ! note that the functions following SigmaPrime in the P_tot sum are all multiplied by an extra factor of q^2/M_N 
                              qOffset = qOffset + 1
                          end if
                          if ((term_R.eq.3) .or. (term_R.eq.4)) then ! add on the single q^2 term
                              qOffset = qOffset + 1
                          else if ((term_R.eq.5) .or. (term_R.eq.6)) then ! add on the two q^2 terms
                              qOffset = qOffset + 2
                          end if
                          
                          w_shared = .false.
                          if((term_R.eq.2).or.(term_R.eq.4).or.(term_R.eq.6)) then ! if the current R term is one with a v_T^2 in it, this accounts for the v^2/c^2 upon expansion (we deal with the negative just after this do loop)
                              w_shared = .true.
                          end if
  
                          !Loop over the different elements
                          do eli = 1, niso
                              a = AtomicNumber_oper(eli)
                              a_shared = a !make accessible via the module
  
                              mu = mdm/(mnuc*a)
                              muplus = (1.+mu)/2.
                              muminus = (mu-1.d0)/2.
  
                              J = AtomicSpin_oper(eli)
                              !!!!!!!!!!!!!!!!!!!!!!!!!simplify yConverse mathematical opperation to help loop speed?
                              ! yConverse = 2/3.*((0.91*(mnuc*a)**(1./3)+0.3)*10**-13)**2/(2*hbar*c0)**2 !the conversion factor between q and y: y = yConverse * q^2
                              yConverse = (1.38*10.**-27) * ( ((mnuc*a)**(1./3) + 0.33) / (hbar*c0) )**2 !should be less computationally demanding I think, unless the 10**-27 makes floats do weird stuff - could combine the 10**-27 with the (hbar*c0)**2 to get a less extreme float prefactor
                              ! get the target nucleus-dm reduced mass mu_T
                              mu_T = (mnuc*a*mdm)/(mnuc*a+mdm)
  
                              ! for all terms, do this base calculation
                              ! the possible y-terms for each W function in order: y^0, y^1, y^2, y^3, y^4, y^5, y^6
                              do term_W = 1,7
                                  WFuncConst = W_array(funcType,eli,tau,taup,term_W)
                                  ! skip if the result gets multiplied by zero in the WFunction
                                  if (WFuncConst.ne.0.) then
                                      q_shared = term_W-1+qOffset ! the GFFI call in the integrand needs to know this
  
                                      !Loop over the shells of constant radius in the star
                                      ! use OMP on this loop: shares vesc_shared(needs to be an array of length ri to be shared with thread safety), umax(depends on vesc, not shared with other functions - make private), and the arrays over the index ri: u_int_res(ri), tab_starrho(ri), tab_mfr(ri,eli), tab_r(ri), tab_dr(ri)
                                      !$OMP parallel 
                                      !$OMP default(private) !(none) 
                                      !$OMP shared(vesc_shared_arr)  ! w_shared is sent to the integrand oper to turn on multiplycation by w^2 there!
                                      do ri = 1, nlines
                                          vesc = tab_vesc(ri)
                                          rindex_shared = ri !make accessible via the module
                                          vesc_shared_arr(ri) = vesc !make accessible via the module
                                          ! Chop the top of the integral off at the smaller of the halo escape velocity or the minimum velocity required for capture.
                                          umax = min(vesc * sqrt(mu)/abs(muminus), vesc_halo)
  
                                          u_int_res(ri) = 0.
                                          !Call integrator
                                          call dsntdqagse(integrand_oper,vdist_over_u,umin,umax, &
                                          epsabs,epsrel,limit,u_int_res(ri),abserr,neval,ier,alist,blist,rlist,elist,iord,last)
  
                                          ! now we deal with the extra negative term mentioned just earlier
                                          if(w_shared) then 
                                              result = 0.
                                              !Call integrator
                                              call dsntdqagse(integrand_oper_extrawterm,vdist_over_u,umin,umax, &
                                              epsabs,epsrel,limit,result,abserr,neval,ier,alist,blist,rlist,elist,iord,last)
                                              u_int_res(ri) = u_int_res(ri) - c0**2/(4.*mu_T**2) * result  ! The Rfunctions are programmed with the 1/c0^2 in their v_perp^2 term (so I need to un-correct it for the- q^2/(2*mu_T)^2, and leave it be for the w^2/c^2)
                                              ! when w_shared is true, u_int_res(ri) contains the result of the integral performed with the w^2/c^2 term, here it is done on _extrawterm with the -q^2/(2*mu_T)^2 term (extra q^2 accounted for by _extrawterm)
                                          end if
                                          !!!!!!!!!!!!!!!!!!!!!!!!!!! the rest of these are multiplicative factors pulled out of the integral (and ideally canceled out to be minimally impactful on performance)
                                          u_int_res(ri)=u_int_res(ri)*(2*mnuc*a)/(2*J+1) *RFuncConst*WFuncConst*yConverse**(term_W-1)
                                          ! extra terms from sumW
                                          ! u_int_res(ri) = u_int_res(ri) * WFuncConst * yConverse**(term_W-1)
                                          ! extra terms from ptot
                                          if (funcType.gt.3) then ! note that the functions following Sigma` in the P_tot sum are all multiplied by an extra factor of q^2/mnuc^2 
                                              u_int_res(ri) = u_int_res(ri) * 1/mnuc**2
                                          end if
                                          ! u_int_res(ri) = u_int_res(ri) * RFuncConst * (hbar*c0)**2
                                          ! extra terms from omega_oper
                                          ! u_int_res(ri) = u_int_res(ri) * (2*mnuc*a)/(2*J+1)
                                          !!!!!!!!!!!!!!!!!!!!!extra terms from captngeneral - check these, it seems Pat has different terms than what I get from my math, specifically this 2.d0*(muplus/mu)**2
                                          ! u_int_res(ri) = u_int_res(ri)*2.d0*NAvo*tab_starrho(ri)*tab_mfr(ri,eli)*(muplus/mdm)**2
                                          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)
  
                                          if (isnan(capped)) then
                                              capped = 0.d0
                                              stop 'NaN encountered whilst trying compute capture rate.'
                                          end if
                                      end do !ri
                                      !$OMP end parallel
                                  end if !WFuncConst>0
                              end do !term_W
                          end do !eli
                      end if !RFuncConst>0
                  end do !term_R
              end do !taup
          end do !tau
      end do !funcType
  
  
      ! ! completes integral (2.3) in paper 1501.03729 (gives dC/dV as fn of radius)
      ! do ri=1,nlines !loop over the star
      !     result = 0.d0
      !     ri_for_omega = ri !accessed via the module
      !     !call integrator
      !     call dsntdqagse(integrand_oper,dummyf,1.d0,vesc_halo, &
      !         epsabs,epsrel,limit,result,abserr,neval,ier,alist,blist,rlist,elist,iord,last)
      !     u_int_res(ri) = result
      !     capped = capped + tab_r(ri)**2*u_int_res(ri)*tab_dr(ri)
      ! end do
  
      !completes integral (2.4) of paper 1501.03729
      capped = 4.d0*pi*Rsun**3*capped
  
      if (capped .gt. 1.d100) then
  
negative.code:     double precision, allocatable :: u_int_res(:)
  
      ! 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
      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)
  
      dimension alist(1000),blist(1000),elist(1000),iord(1000),rlist(1000)!for integrator
      !external gausstest !this is just for testing
  