anchor.diff:  cumint(i) = cumint(i-1) + integrand*tab_dr(i)*Rsun
   end if
   nx(i) = (tab_T(i)/Tc)**(3./2.)*exp(-cumint(i))
  -! print*,nx(i)
   cumNx = cumNx + 4.*pi*tab_r(i)**2*tab_dr(i)*nx(i)*Rsun**3.
   
   nxIso(i) = Nwimps*exp(-Rsun**2*tab_r(i)**2/rchi**2)/(pi**(3./2.)*rchi**3) !normalized correctly
  -! print*,exp(-Rsun**2*tab_r(i)**2/rchi**2)
   end do
   nx = nx/cumNx*nwimps !normalize density
  -! print*, "niso 1 ", NxIso(1), tab_r(1), Nwimps, 1./(pi**(3./2.)*rchi**3)
   fgoth = 1./(1.+(K/.4)**2)
   hgoth = ((tab_r*Rsun - rchi)/rchi)**3 +1.
   hgoth(1) = 0.d0 !some floating point shenanigans.
  
positive.code: ! dLdr(nlines) = 0.d0
  ! dLdr = dLdr/Rsun
  
  print *, "r=", tab_r(int(nlines/2)), "epso=", epso, "rho=", tab_starrho(int(nlines/2)), "Rsun=", Rsun, "dLdR=", dLdR(int(nlines/2))
  
  Etrans = 1./(4.*pi*(tab_r+epso)**2*tab_starrho)*dLdR/Rsun**2;
  
  EtransTot = trapz(tab_r,abs(dLdR),nlines)
  
  open(10, file="/home/luke/summer_2020/mesa/captngen/Etrans_gr.txt")
  do i=1,nlines
  	write(10, *) tab_r(i), Etrans(i), tab_starrho(i), nx(i), tab_T(i)
  enddo
  close(10)
  
  ! Some testing bits:
   open(55,file = "/home/luke/summer_2020/mesa/captngen/captranstest.dat")
   do i=1,nlines
   write(55,*) tab_r(i), nx(i), tab_T(i), Ltrans(i), Etrans(i), dTdR(i), dLdR(i), tab_starrho(i), tab_g(i), dphidr(i)
   end do
   close(55)
  
  print *, "Nwimps = ", Nwimps
  
  ! Check input units
  !open(3, file="/home/luke/summer_2020/mesa/captngen/unit_check.txt")
  !write(3,*) "--------tab_r---------"
  !do i=1, nlines
  !	write(3,*) i, tab_r(i)
  !enddo
  
  !write(3,*) "--------tab_T---------"
  !do i=1, nlines
  !	write(3,*) i, tab_T(i)
  !enddo
  
  !write(3,*) "--------phi---------"
  !do i=1, nlines
  !	write(3,*) i, phi(i)
  !enddo
  
  !write(3,*) "--------tab_starrho---------"
  !do i=1, nlines
  !	write(3,*) i, tab_starrho(i)
  !enddo
  
  !do i=1, nlines
  !	write(3,*) tab_r(i), Etrans(i)
  !enddo
  !close(3)
  
  return
  
  ! print*,Ltrans(1),Etrans(1), dLdR(1),tab_r(1)
  
  !
  ! open(55,file = "smallarrays.dat")
  
negative.code:     alphaofR = alphaofR/(sigma_N*sum(nabund,1))
  
  
  if (nonlocal .eqv. .false.) then ! if nonlocal=false, use Gould & Raffelt regime to calculate transport
  
  !compute mean free path
  