;################ @COPYRIGHT  C.brook -  A. di cintio  ############################################
;################ cbabrook@gmail.com - dicintio.arianna@gmail.com #################################
;this IDL code computes the density profiles and velocity curves of a
;distribution of DM haloes. it uses 2 density profiles for DM
;haloes,the universal NFW one (http://adsabs.harvard.edu/abs/1997ApJ...490..493N)
; and the mass-dependent DC14 profile (http://adsabs.harvard.edu/abs/2014MNRAS.441.2986D).
; The DC14 profile takes into account the effect of baryonic feedback on the DM distribution within
;galaxies (see http://adsabs.harvard.edu/abs/2014MNRAS.437..415D for reference). 
; Feel free to use this code  but please cite the DC14 profile  and acknowledge C.Brook A.Di Cintio
;####################################################################################################

FUNCTION r2density, X   ; gives r^2*density of DC14 profile
common share,a,b,g,cn,c_dcb,Rv,rscale,r_core,rs
RETURN, x^2/((c_dcb*x/Rv)^g*(1+(c_dcb*x/Rv)^a)^((b-g)/a))
end

FUNCTION r2nfw, X   ; gives r^2*density  of NFW profile            
common share,a,b,g,cn,c_dcb,Rv,rscale,r_core,rs         
RETURN, x^2/((cn*x/Rv)*(1+(cn*x/Rv))^2)                     
end              



pro DC14_NFW_curves
common share,a,b,g,cn,c_dcb,Rv,rscale,r_core,rs 
                ;parameters in Planck cosmology
                H = 0.068 ; km/s kpc-1
                Grav = 4.302e-6  ; kpc (km/s)2 Msun-1
                rho_c = (3.* H^2)/(8.*!pi*Grav)
                omega_m=0.3
                delta_c=96.; this could be a fixed overdensity - like 100 or 200- or a cosmology dependent one, a la  Bryan & Norman 1998

loadct,39

;create distribution of haloes with halo mass 10^8 to 10^12 -> dwarfs to Milky Way scales
logM=8+(indgen(41)*.1)
M_cdm=10^logM
;use abundance matching relation - here Guo+11 - to derive the corresponding M* for each Mhalo
MsMhG=0.129*((M_cdm/(10^11.4))^(-.926)+(M_cdm/(10^11.4))^.261)^(-2.44)
Ms=MsMhG*M_cdm


;define efficiency - in the DC14 formalism, this parameter determines if  a galaxy forms a
;central DM core or it retains the primordial cuspy NFW density  
;**********IMPORTANT : if you wish to plot nfw curves for every halo,rather than dc14 ones, simply add  -1000. in efficiency definition,
;so that every halo has a nfw profile with nfw concentration***************************

eff=alog10(Ms/M_cdm); -1000. 
print,eff


; concentration-mass relation from planck cosmology, Dutton & Maccio 14
logC=1.025-0.097*alog10(M_cdm*0.68/1e12)
C=10^logC
;define virial radius, virial velocity  and scale radius
Rv_cdm= (2*Grav*M_cdm/(delta_c*H^2))^(1/3.)  ;
Vvir=sqrt(Grav*M_cdm/Rv_cdm)
Rs_cdm=Rv_cdm/C


bins=500L
dir='./'
set_plot,'ps'
device,filename=dir+'DC14curves.eps',/encapsulated,$
/color,/inches,xsize=3.4,ysize=2.6
plot,[.1,.1],[.1,.1],/ylog,/xlog,$
xrange=[0.01,300],/xstyle,yrange=[1,300],/ystyle,xtit='Radius [kpc]',ytit='V!dcirc!n(r)[km/s]',$
position=[0.16,0.18,.968,.97],charthick=3




for i = 0L,n_elements(M_cdm)-1L do begin
Rv=Rv_cdm[i]
radius=(indgen(long(bins*Rv),/long)+1L)/(1.*bins)
;define concentration in DC14 formalism.this follows eq 6 in DC14
;paper and includes contraction of haloes at large masses, i.e. at the
;Milky Way mass                     
cn=(1.0+0.0004*exp(2.4*(eff[i]+4.5)))*C[i]
rs=Rv/cn 

;if the efficiency parameter is within the range explored in DC14,then
;compute the DC14 profile and corresponding velocity...
if eff[i] ge -4.11 and eff[i] le -1.3 then begin            
A= 2.94182 - alog10( (10^(eff[i]+2.3268)^(-1.08))+(10^(eff[i]+2.3268)^(2.29)))
B= 4.23 + 1.34*eff[i] + 0.26*eff[i]^2                                    
G =-0.06 + alog10( (10^(eff[i]+2.56496))^(-0.6776) + (10^(eff[i]+2.56496)))
c_dcb=cn*((g-2)/(2-b))^(1./a)                 

r_dcb=Rv/c_dcb                                               
nrad=n_elements(radius)                                                    
rho_s=M_cdm[i]/(4*!pi*qsimp('r2density',.2,Rv,eps=1e-5))                    
den_dcb=rho_s/((radius/r_dcb)^g*(1+(radius/r_dcb)^a)^((b-g)/a))               
mass_dcb=4/3.*!pi*((radius+1./bins/2.)^3-(radius-1./bins/2.)^3)*$   
den_dcb                                          
cummass=cumulate(mass_dcb)       
;take into account the baryon fraction before plotting V(r) curve of DM halo                                        
Vc_dcb=sqrt(.856*grav*cummass/radius)                               

endif    

;...else, compute a NFW profile and corresponding velocity
if (eff[i] lt -4.11 or eff[i] gt -1.3) then begin
cn=C[i]
rho_0_NFW= M_cdm[i]/(4.*!pi* qsimp('r2nfw',.01,Rv,eps=1e-5))
den_nfw=rho_0_NFW/((cn*radius/Rv) * (1+(cn*radius/Rv)^2)) 
mass_nfw=4/3.*!pi*((radius+1/bins/2.)^3-(radius-1./bins/2.)^3)*den_nfw
cummass=cumulate(mass_nfw) 
Vc_DCB=sqrt(.856*grav*cummass/radius)   
endif


oplot,radius,Vc_DCB,color=240-6*i,thick=3
endfor

device,/close
set_plot,'x'
stop
end