anchor.diff:    double precision, intent(in) :: Nwimps
     double precision, intent(out) :: EvapRate
     double precision :: Earg(nlines), muarray(niso_in),nabund(niso_in,nlines),sigma_N(niso_in)
  -  double precision :: suparg,trapz
  +  double precision :: suparg,trapz,Knud,rchi
     double precision :: mdmg,mnucg, Tc, rhoc,mn, Tw,nin,escFrac(nlines),vescc(nlines),nxIso(nlines),mfp(nlines),scatrate(nlines)
  -  double precision, parameter :: kBeV=8.617e-5
  +  double precision, parameter :: kBeV=8.617d-5, GN =6.674d-8, kB=1.3806d-16
   
     integer :: i, j, k
     mdmg = mdm*1.78d-24
  
positive.code: end if
  fofn = f
  return
  end
  
  !Fast trapezoidal integral
  double precision function trapz(x,y,flen)
    implicit none
    integer, intent(in) :: flen
    double precision, intent (in) :: x(flen), y(flen)
    integer i
  
    trapz = y(1)*(x(2)-x(1))/2. + y(flen)*(x(i)-x(i-1))/2.
    do i = 2,flen-1
      trapz = trapz +y(i)*(x(i)-x(i-1))
    end do
  
  
    return
  end
  
negative.code:   nxIso = exp(mdm*vescc**2/2./Tw);
    nin = 4.d0*pi*trapz(tab_r,tab_r**2.*nxIso,nlines) !%niso norm
    nxIso = nxIso/nin
  
    ! print*,"norm guy ", nin ! "one: ", 4.d0*pi*trapz(tab_r,tab_r**2.*nxIso,nlines)
    !Fraction of the kinetic distribution above the local escape velocity
    escFrac = sqrt(2.d0/pi)*vescc*sqrt(mdm/Tw)*exp(-mdm*vescc**2/Tw/2.d0) - derf(sqrt(mdm/Tw/2.d0)*vescc) + 1.d0;
    ! print*,"EscFrac = ", escFrac,
  