anchor.diff:  
   function Tx_integral(T_x, sigma_N, Nwimps, niso)
   implicit none
  -! Calculates the Tx defining integral 
  +! Calculates the Tx defining integral
   
   integer, intent(in) :: niso
   double precision, intent(in) :: T_x, Nwimps
  
positive.code: double precision :: n_x(nlines), species_indep(nlines), species_dep(nlines), sigma_nuc(niso)
  double precision :: Etrans_sp(nlines)
  integer :: i, j
  ! T_x in K, sigma_N in cm^2, 
  
  R = tab_r*Rsun ! R in cm
  phi = -tab_vesc**2/2.d0 ! phi in erg/g
  
negative.code: r = tab_r*Rsun ! cm
  phi = -tab_vesc**2/2.d0 ! erg/g
  mxg = mdm*1.782662d-24
  ! open(95,"randphi.dat")
  ! write(95,*) r , phi
  ! close(95)
  ! WIMP number density in isothermal approximation
  nx_isothermal = exp(-mxg*(phi-phi(1))/kB/T_x)
  n_0 = Nwimps/trapz(r, 4.d0*pi*r**2.d0*nx_isothermal, nlines) ! Normalize so that integral(nx) = Nwimps
  nx_isothermal = n_0*nx_isothermal
  
  