c  Program to calculate some parameters of DM halos for different virial masses and redshifts.

c  1.04: Real M_s, V_s values are added (corresponding to enclosed mass at r_s, and circular speed at r_s).
c  1.05: NFW/Burkert input parameter

      program halo
      implicit NONE
      real*8 z, m_vir,H0,Omega_m,Omega_mz,Delta,x,r_vir,c_vir,V_vir,
     #       H,sgm,T_vir,zriddr,A,rh,tau,vs,r_s,
     #       Ms,rho_s,M0,V0
      integer i
      logical NFW
      character sprof*20
      external funcn,funcb
      common A
      include 'const.f'

      write (*,*) 'NFW or Burkert (n/b):'
      read (*,*) sprof
      write (*,*)  sprof
      write (*,*) 'm_vir (M_Sun):'
      read (*,*) m_vir
      write (*,*)  m_vir
      write (*,*) 'Redshift:'
      read (*,*) z
      write (*,*)  z
      write (*,*) 'Deviation from the mean sigma:'
      read (*,*) sgm
      write (*,*)  sgm
      write (*,*) 'Omega_m:'
      read (*,*) Omega_m
      write (*,*)  Omega_m
      write (*,*) 'Hubble constant H0 (km/s/Mpc):'
      read (*,*) H0
      write (*,*) H0

      i = 0
 1    continue
      i = i + 1
      if (i .gt. len(sprof))  stop 'Empty string!'
      if (sprof(i:i).eq.' ') go to 1
      NFW = (sprof(i:i).eq.'n' .or. sprof(i:i).eq.'N')
      if (.not.NFW .and. sprof(i:i).ne.'b' .and. sprof(i:i).ne.'B') 
     #      stop 'Should be either n or b!'

c Hublle units of s^{-1}:
      H = H0 / pc / 10.d0

c  Valid for flat LCDM:
      Omega_mz = 1.d0/(1.d0 + (1.d0-Omega_m)/Omega_m/(1.d0+z)**3)
      x = Omega_mz - 1.d0
c  Valid for 0.1<Omega_mz<1: [1998ApJ...495...80]
      Delta = (18.d0*pi**2+82.d0*x-39.d0*x**2) / Omega_mz
c  Valid for flat LCDM:
      r_vir = 1.d0/(1.d0+z) * (2.d0*G*m_vir*M_Sun/Delta/H**2/Omega_m)**
     #        (1.d0/3.d0) / kpc
c  Eqs. (14-15) of ApJS 143, 419:
      if (m_vir .ge. 1.0702d11)  then
         c_vir = 10.3d0/(1.d0+z) * (m_vir/1.d13)**(-0.13d0) * 
     #           10.d0**(0.14d0*sgm)
         else
         c_vir = 27d0/(1.d0+z) * (m_vir/1.d9)**(-0.08d0) * 
     #           10.d0**(0.14d0*sgm)
         end if

      r_s = r_vir/c_vir

      V_vir = dsqrt(G*m_vir*M_Sun/r_vir/kpc) / 1.d5
c  Virial temperature for primordial gas with HII and HeII (Barkana & Loeb 2000):
      T_vir = 0.61d0 * mp * (V_vir*1.d5)**2 /2.d0 /boltz

c  Half-mass radius (NFW):
      if (NFW)  then
         A = 0.5d0*(1.d0/(1.d0+c_vir)+dlog(1.d0+c_vir)-1.d0)
         rh = r_vir/c_vir * zriddr(funcn,0.d0,c_vir,1.d-9)
         else
c  Half-mass radius (Burkert):
         A = 0.5d0*(dlog(1.d0+c_vir)+0.5d0*dlog(1.d0+c_vir**2)-
     #        datan(c_vir))
         rh = r_vir/c_vir * zriddr(funcb,0.d0,c_vir,1.d-9)
         end if

c  Crossing time at rh:
      tau = dsqrt(2.d0*(rh*kpc)**3/G/(m_vir*M_sun))/Gyr

c  Scaling mass (Ms = 4*pi * r_s**3 * rho_s):
      if (NFW)  then
         Ms = m_vir/(dlog(1.d0+c_vir)-c_vir/(1.d0+c_vir))
         else
         Ms = 2.d0*m_vir/(dlog(1.d0+c_vir)+0.5d0*dlog(1.d0+c_vir**2)-
     #        datan(c_vir))
         end if

c  Scaling density (the DM density profile is rho=rho_s/x/(1+x)**2 for NFW,
c  and rho=rho_s/(1+x)/(1+x**2) for Burkert; here x=r/r_s):
      if (NFW)  then
         rho_s = m_vir/(4.d0*pi*r_s**3)/
     #         (dlog(1.d0+c_vir)-c_vir/(1.d0+c_vir))
         else
         rho_s = m_vir/(2.d0*pi*r_s**3)/(dlog(1.d0+c_vir)+
     #         0.5d0*dlog(1.d0+c_vir**2)-datan(c_vir))
         end if

c  Scaling velocity (just a scaling parameter; does not correspond to any 
c  physical circular velocity in the halo):
      vs = dsqrt(G*Ms*M_Sun/r_s/kpc) / 1.d5

c  At the scaling radius r_s:
c    Enclosed mass:
      if (NFW)  then
         M0 = Ms * (dlog(2.d0)-0.5d0)
         else
         M0 = Ms * (3.d0*dlog(2.d0)-pi/2.d0)/4.d0
         end if
c    Circular velocity:
      v0 = dsqrt(G*M0*M_Sun/r_s/kpc) / 1.d5

      write (*,*)
      write (*,'(A,1PE13.5,A)') 'r_vir=',r_vir,' kpc'
      write (*,'(A,1PE13.5,A)') 'r_h=',rh,' kpc'
      write (*,'(A,1PE13.5,A)') 'r_s=',r_s,' kpc'
      write (*,'(A,1PE13.5)')   'c_vir=',c_vir
      write (*,'(A,1PE13.5,A)') 'V_vir=',V_vir,' km/s'
      write (*,'(A,1PE13.5,A)') 'T_vir=',T_vir,' K'
      write (*,'(A,1PE12.5,A)') 'tau_h = ',tau,' Gyr'
c      write (*,'(A,1PE13.5,A)') 'M_s=',Ms,' M_Sun'
      write (*,'(A,1PE13.5,A)') 'rho_s=',rho_s,' M_Sun/kpc^3'
c      write (*,'(A,1PE13.5,A)') 'v_s=',vs,' km/s'
      write (*,*)
      write (*,*) '  Values at the scaling radius r_s:'
      write (*,'(A,1PE13.5,A)') 'M(r_s)=',M0,' M_Sun'
      write (*,'(A,1PE13.5,A)') 'V_c(r_s)=',v0,' km/s'
      

      stop
      end
c-------------------------------------------------------------------
      function funcn (r)
      implicit NONE
      real*8 funcn, r, A
      common A

      funcn = 1.d0/(1.d0+r) + dlog(1.d0+r) - 1.d0 - A

      return
      end

c-------------------------------------------------------------------
      function funcb (r)
      implicit NONE
      real*8 funcb, r, A
      common A

      funcb = dlog(1.d0+r)+0.5d0*dlog(1.d0+r**2)-datan(r) - A

      return
      end

c-------------------------------------------------------------------
      FUNCTION zriddr(func,x1,x2,xacc)
      INTEGER MAXIT
      REAL*8 zriddr,x1,x2,xacc,func,UNUSED
      PARAMETER (MAXIT=60,UNUSED=-1.11E30)
      EXTERNAL func
CU    USES func
      INTEGER j
      REAL*8 fh,fl,fm,fnew,s,xh,xl,xm,xnew
      fl=func(x1)
      fh=func(x2)
      if((fl.gt.0..and.fh.lt.0.).or.(fl.lt.0..and.fh.gt.0.))then
        xl=x1
        xh=x2
        zriddr=UNUSED
        do 11 j=1,MAXIT
          xm=0.5*(xl+xh)
          fm=func(xm)
          s=sqrt(fm**2-fl*fh)
          if(s.eq.0.)return
          xnew=xm+(xm-xl)*(sign(1.,fl-fh)*fm/s)
          if (abs(xnew-zriddr).le.xacc) return
          zriddr=xnew
          fnew=func(zriddr)
          if (fnew.eq.0.) return
          if(sign(fm,fnew).ne.fm) then
            xl=xm
            fl=fm
            xh=zriddr
            fh=fnew
          else if(sign(fl,fnew).ne.fl) then
            xh=zriddr
            fh=fnew
          else if(sign(fh,fnew).ne.fh) then
            xl=zriddr
            fl=fnew
          else
            pause 'never get here in zriddr'
          endif
          if(abs(xh-xl).le.xacc) return
11      continue
        pause 'zriddr exceed maximum iterations'
      else if (fl.eq.0.) then
        zriddr=x1
      else if (fh.eq.0.) then
        zriddr=x2
      else
        pause 'root must be bracketed in zriddr'
      endif
      return
      END
c-------------------------------------------------------------------