anchor.diff:  
                                       ! calculates the total number of q^2
                                       q_index = 1 + q_functype + term_W - 1 + floor((term_R-1.)/2.)
  +                                    prefactor_current = prefactor_functype*RFuncConst*WFuncConst*yConverse_array(eli)**(term_W-1)
   
                                       ! check if term_R is even (in my index convention this corresponds to it having a v^2 in the term)
                                       ! v^2 = w^2 - q^2/(2mu_T)^2
                                       if ( mod(term_R,2).eq.0 ) then
                                           ! this is the -q^2/(2mu_T)^2 contribution
  -                                        prefactor_array(eli,q_index+1,1) = prefactor_array(eli,q_index+1,1) - &
  -                                            prefactor_functype*RFuncConst*WFuncConst*yConverse**(term_W-1) * (c0**2/(4.*mu_T**2)) ! 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)
  +                                        prefactor_array(eli,q_index+1,1) = prefactor_array(eli,q_index+1,1) - prefactor_current * &
  +                                            (c0**2/(4.*mu_T**2)) ! 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)
                                           ! this is the +w^2 contribution
  -                                        prefactor_array(eli,q_index,2) = prefactor_array(eli,q_index,2) + &
  -                                            prefactor_functype*RFuncConst*WFuncConst*yConverse**(term_W-1)
  +                                        prefactor_array(eli,q_index,2) = prefactor_array(eli,q_index,2) + prefactor_current
                                           
                                       else
  -                                        prefactor_array(eli,q_index,1) = prefactor_array(eli,q_index,1) + &
  -                                            prefactor_functype*RFuncConst*WFuncConst*yConverse**(term_W-1)
  +                                        prefactor_array(eli,q_index,1) = prefactor_array(eli,q_index,1) + prefactor_current
   
                                       end if
                                   end if
  
positive.code:                     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
          end do !eli
      end do !ri
  
  
negative.code:                                                         0., 0., 0., 0., 0., 0./) !spins pulled from https://physics.nist.gov/PhysRefData/Handbook/element_name.htm
      double precision :: coupling_Array(14,2)
      double precision :: W_array(8,16,2,2,7)
      double precision :: yConverse_array(16)
  
      integer :: q_shared
      logical :: w_shared
  