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:   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 :: mdmg,mnucg, Tc, rhoc,mn, Tw,nin,escFrac(nlines),vescc(nlines),nxIso(nlines),mfp(nlines),scatrate(nlines)
    double precision, parameter :: kBeV=8.617e-5
  
    integer :: i, j, k
    mdmg = mdm*1.78d-24
  
negative.code:   ! else if ((nq .eq. )) !q, v dependence goes here
    end if
  
    rchi = (3.*(kB*Tc)/(2.*pi*GN*rhoc*mdmg))**.5;
    Knud = mfp(1)/rchi;
    ! print*,"Knud", Knud
    if (Knud .lt. 0.1) then
      print*, "WARNING, K = ", knud, " is < 0.1. Approximate evaporation scheme is likely very wrong."
    end if
  
    suparg = trapz(tab_r,rsun/mfp,nlines)
    Earg = escFrac*scatrate*exp(-suparg)*c0
  
    open(55,file = "sv.dat")
    do j=1,nlines
    write(55,*) tab_r(j), Earg(j),suparg,nxIso(j),escFrac(j),scatrate(j),exp(-suparg)
    end do
    close(55)
  
  
    EvapRate = Nwimps*4.*pi*trapz(tab_r,tab_r**2*nxIso*Earg,nlines)
  
    if (isnan(EvapRate)) then
      stop "NaN evap rate, check it"
    end if
  
  end subroutine fastevap
  
  