anchor.diff:    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,
  
positive.code: 
      end do
              Tw_out = sum(Tw_out_num)/sum(Tw_out_denom);
              ! print*,Tw_out, Tw
              dT = abs(Tw - Tw_out)/Tw;
              Tw = Tw_out;
              ! print*, "dT ", dT
    end do
    Tw = Tw/TcGeV;
  
  
  
  
  
negative.code:   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,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.617d-5, GN =6.674d-8, kB=1.3806d-16
  
    integer :: i, j, k
    mdmg = mdm*1.78d-24
  