anchor.diff:      
       contains
   
  +    ! having removed the scaling momentum, are the units off here? I'm looking at the p/c0 in particular
       function GFFI_H_oper(w,vesc,mq)
           double precision :: p, mu,w,vesc,u,muplus,GFFI_H_oper,G
           integer mq
  
positive.code:                                             !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
                                          end if
                                          !!!!!!!!!!!!!!!!!!!!!!!!!!! the rest of these are multiplicative factors pulled out of the integral (and ideally canceled out to be minimally impactful on performance)
                                          ! 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) !!!!!!!!!!!!!!!!!! we have a rogue 'w' here too, flag it for the integrand!
                                          !!!!!!!!!!!!!!!!!!!!!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/mx_in)**2
                                          ! u_int_res(ri) = u_int_res(ri)*NAvo*tab_starrho(ri)*tab_mfr(ri,eli)/(mnuc*a)
  
                                          ! !!!! trying to collect all multiplicative factors here
                                          ! u_int_res(ri) = u_int_res(ri) * WFuncConst*yConverse**(term_W-1) &
                                          !     * RFuncConst*(hbar*c0)**2 * 2/(2*J+1) &
                                          !     * NAvo*tab_starrho(ri)*tab_mfr(ri,eli)
                                          ! 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
  
                                          capped = capped + tab_r(ri)**2*u_int_res(ri)*tab_dr(ri)
  
  
negative.code:                                     ! 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
  