anchor.diff:                              a = AtomicNumber_oper(eli)
                               a_shared = a !make accessible via the module
   
  -                            mu = mx_in/(mnuc*a)
  +                            mu = mdm/(mnuc*a)
                               muplus = (1.+mu)/2.
                               muminus = (mu-1.d0)/2.
   
  
positive.code:     
      contains
  
      function GFFI_H_oper(w,vesc,mq)
          double precision :: p, mu,w,vesc,u,muplus,GFFI_H_oper,G
          integer mq
  
negative.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  ! 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)
  
  