;;;;; powspec.pro - Computes starting and stopping redshift for
;;;;;               cosmological simulations
;;;;;
;;;;; Manodeep Sinha v1.0 23rd March, 2011. 
;;;;; email: manodeep@gmail.com
;;;;;
;;;;;
;;;;; To run:        Change the ngrid, boxsize, rms_ratio_start, delta_sqr_stop
;;;;;                and z_obs to values of your liking. The code will
;;;;;                compute the rms displacement etc and hopefully give 
;;;;;                (correct) useful answers. 
;;;;;
;;;;; Bugs:           Please email me - manodeep@gmail.com
;;;;;
;;;;; Known Issues:   Code is only tested for reasonable values of
;;;;;                 ngrid/boxsize - you will have to check that
;;;;;                 nothing funky is happening, e.g., if you are trying 
;;;;;                 to run a 1 Mpc box with 4096^3 particles..
;;;;;
;;;;; Revisions:      v1.1 19th June, 2012: Made the main routine a
;;;;;                 procedure and created a main level code that
;;;;;                 runs the example code; i.e., made it slightly
;;;;;                 more user-friendly. - MS


function solve_for_growth,z
compile_opt idl2,strictarrsubs
common cosmology_tf, omega_m,omega_lambda,omega_b,ns,ns_bar,sigma8,h,delta_h,norm,t_cmb
common root_solve,req_growth_fac

return, growth_d(z)-req_growth_fac

end

function tophatsigma2
compile_opt idl2,strictarrsubs
common sigma8_norm,r_tophat
r_tophat = 8.0 ;; Mpc/h

return, qromb('sigma2_int',0.0,500.0d/r_tophat,/double)
end

function sigma2_int,k
compile_opt idl2,strictarrsubs
common sigma8_norm

kr = r_tophat * k
kr2 = kr * kr
kr3 = kr2 * kr

if kr lt 1d-8 then $
  return, 0.0

w = 3 * (sin(kr) / kr3 - cos(kr) / kr2)
x = 4 * !pi * k * k * w * w * powspec_eh(k)

return, x


end

function growth_d,z
compile_opt idl2,strictarrsubs
common cosmology_tf
if n_elements(z) eq 0 then stop
theta = t_cmb/2.7d
zeq = 2.5d4*omega_m*h^2.0*theta^(-4.0)
gz_sqr  = omega_m*(1+z)^3.0 + (1-omega_m-omega_lambda)*(1+z)^2 + omega_lambda
omega_m_z = omega_m*(1+z)^3/gz_sqr
omega_lambda_z = omega_lambda/gz_sqr

return, (1+zeq)/(1+z)*2.5*omega_m_z/(omega_m_z^(4./7.) - omega_lambda_z + (1+omega_m_z/2.0)*(1+omega_lambda_z/70.0))

end




function powspec_eh, k
compile_opt idl2,strictarrsubs
common cosmology_tf
tk = tk_eh_fit(k)
return, norm*k*tk^2.0
end

function tk_eh_fit, k
compile_opt idl2,strictarrsubs
common cosmology_tf

hubble = h
omegam = omega_m
ombh2 = omega_b *hubble^2.0

theta = t_cmb / 2.7d
ommh2 = omegam * hubble * hubble
s = 44.5 * alog(9.83 / ommh2) / sqrt(1. + 10. * exp(0.75 * alog(ombh2))) * hubble
a = 1. - 0.328 * alog(431. * ommh2) * ombh2 / ommh2 + 0.380 * alog(22.3 * ommh2) * (ombh2 / ommh2) * (ombh2 / ommh2)
gamma = a + (1. - a) / (1. + exp(4 * alog(0.43 * k * s)))
gamma *= (omegam * hubble)
q = k * theta * theta / gamma
L0 = alog(2.0 * exp(1.) + 1.8 * q)
C0 = 14.2 + 731.d / (1. + 62.5 * q)
tmp = L0 / (L0 + C0 * q * q)

return, tmp

end


;;; this will compute some 'useful' things using the powerspectrum.
pro powspec,ngrid,boxsize,delta_sqr_stop,rms_ratio_start

compile_opt idl2,strictarrsubs
common cosmology_tf
common root_solve

if n_elements(ngrid) eq 0 or n_elements(boxsize) eq 0 then begin
    print,'ERROR: Provide both N (for a N^3 simulation) and the boxsize [in Mpc/h] '
    return
endif

boxsize = double(boxsize) ;;; ensure that boxsize is double -> it does imply that the calling routine boxsize variable will become double as well. 


;;; I am not actually using the CMB normalisation -- using sigma8 instead
delta_h  = 1.94d-5*(omega_m)^(-0.785 -0.05d*alog(omega_m))*exp(-1.02d*ns_bar - 1.70*ns_bar^2.0d)       ;;; eqn. 32 in EH 97

res = tophatsigma2()
norm = sigma8^2.0/res


;;;; how much power are you willing to have in delta_sqr for deciding a stopping z. 
;;;; Delta_sqr = 1 means definitely non-linear (however, mode coupling actually causes 
;;;; non-linearity to set in at a lower value of delta_sqr)
if n_elements(delta_sqr_stop) eq 0 then $
  delta_sqr_stop = 0.1d 

;;;; this is the initial displacement. If it is > 1.0, then
;;;; shell-crossing occurs.
if n_elements(rms_ratio_start) eq 0 then $
  rms_ratio_start = 0.25d ;;; ratio of initial rms/delta_p desired. Max. disp ~ 3 \times rms. rms_ratio_start = 0.25 => shell-crossing is a 4-sigma event

z_obs = 20.0 ;; where the first real objects should be somewhat reliable
z_start = [65.0,80.0,99.0,199.0,249.0,299.0,399.0,499.0,799.0] ;;; do the rms calculations at least at these z's

;;;;; that's it. and ensure that kmin and kmax are relatively large
;;;;; and encompass your box_kmin and box_kmax (they should, unless
;;;;; you have a really weird box/N combo)

totn  = ngrid^3.0
delta_p = boxsize/ngrid
box_kmin = 2.0*!pi/boxsize ;;; fundamental mode
box_kmax = !pi/delta_p ;;; nyquist frequency

;;;; now generate pk
kmin = 1d-5
kmax = 1d4
nkbins = 1d6 

dlogk = (alog10(kmax)-alog10(kmin))/(nkbins-1.0)
all_log_k = dindgen(nkbins)*dlogk + alog10(kmin)
all_pk    = dblarr(nkbins)
all_k     = dblarr(nkbins)
counter=0
for log10k=alog10(kmin),alog10(kmax),dlogk do begin
    k = 10.d^log10k
    all_pk[counter] = powspec_eh(k)*(k^(ns-1.0))
    all_k[counter]  = k
    counter=counter+1
endfor

a_obs = 1.0/(1+z_obs)
a_min_start = a_obs/10.0
z_min_start = 1/a_min_start-1.0

zmin = 50.0d > z_min_start
dz = 50.0
xx = fix(zmin/dz)
if xx gt 0 then zmin = xx*dz
zmax = 800.0d

growth_0 = growth_d(0.0)
valid_inds = where(all_k ge box_kmin and all_k le box_kmax,cnt)
if cnt eq 0 then begin
    print,'Something is wrong: I can not find any k that fits inside the box'
    stop
endif else begin
    ;;; integrate between kmin and kmax (if the integral is to computed
    ;;; between box_kmin and box_kmax just comment the following line)
;    valid_inds = lindgen(n_elements(all_k))
    integral = int_tabulated(all_k[valid_inds],all_pk[valid_inds],/double)
endelse

print,z_obs,z_min_start,format='("Assuming at least 10 scale factors between z and z_obs = ",F6.1," gives a min. starting z = ",F6.1)'
mpc_to_kpc = 1d3
for iz = zmin,zmax,dz do begin
    growth_factor = growth_d(iz)
    rms = sqrt(4.0*!pi*integral/3.0d)*(growth_factor/growth_0) ;;; p(k) contains D1(z)^2/D1(0)^2 -> they come out of the integral as a linear ratio
    print,iz,mpc_to_kpc*rms,rms/delta_p*100.0,$
      format='("At z = ",F7.1," rms = ",F6.2," kpc/h   rms/delta_p = ",F6.1," %")'
    if rms gt boxsize then stop
endfor

;;;; Now figure out what z_init should really be used
print,''
req_growth_fac = sqrt(rms_ratio_start^2.0*3.0*growth_0^2.0*delta_p^2.0/(4.0*!pi*integral))
print,rms_ratio_start,format='("Now calculating starting z for rms/delta_p = ",F10.3,"  ...", $)'
z_guess = [10.0]
z_init = newton(z_guess,'solve_for_growth',itmax=10000,tolx=0.1,/double)
print,"..done"
print,z_init,format='("The simulation should be started at at least z = ",F8.1)'

;;; Also report for whatever redshift is in z_start
print,''
for i=0l,n_elements(z_start)-1 do begin
    iz = z_start[i]
    growth_factor = growth_d(iz)
    rms = sqrt(4.0*!pi/3.0*integral)*(growth_factor/growth_0)
    print,iz,1d3*rms,rms/delta_p*100.0,$
      format='("At z = ",F7.1," rms = ",F6.2," kpc/h   rms/delta_p = ",F6.1," %")'

endfor


;;; Now find the stopping redshift -- let's find the power
;;; in the fundamental mode.

delta_sqr = 4*!pi*all_k^3.0*all_pk
xx = min(abs(all_k - box_kmin),index)
k_stop = all_k[index]

if delta_sqr[index] le delta_sqr_stop then begin
    print,''
    print,delta_sqr_stop,format='("This simulation can proceed to z=0 with delta_stop =  ",F10.2)'
    z_stop = 0.0
endif else begin
    print,''
    print,format='("This box should not be evolved to z=0 ")'
    print,''
    print,delta_sqr_stop,format='("Now calculating stopping z for delta_stop = ",F10.2," ...", $)'

    req_growth_fac = sqrt(delta_sqr_stop*growth_0^2.0/delta_sqr[index])
    z_guess = [0.0]
    z_stop = newton(z_guess,'solve_for_growth',itmax=1000,tolx=0.1,/double) ;;; don't really care about z being more precise than 0.1

    print,"..done"
Endelse

print,''
print,'Simulation params:'
if boxsize gt 1 then $
  print,boxsize,format='("Boxsize = ",I6," Mpc/h")' $
else $
  print,long(1d3*boxsize),format='("Boxsize = ",I6," kpc/h")' 

print,ngrid,  format='("Ngrid   = ",I6)'
print,rms_ratio_start,format='("Initial (max) displacement from grid (rms/delta_p) = ",F10.3)'
print,delta_sqr_stop,format='("Max. power in \Delta^2 for the fundamental mode at final z = ",F10.3)'

print,''
print,'Suggested params:'
print,z_init,format='("IC redshift = ",F10.1)'
print,1d3*delta_p/40.0,1d3*delta_p/20.0,format='("Softening length = [",F0.1," - ",F0.1,"] kpc/h ")'
print,z_stop,format='("Final redshift = ",F6.1 )'

end


;;; end of all the actual routines to compute the starting and
;;; stopping redshifts. Below is the example code to run this idl
;;; routine. 








;;;; Change ngrid and boxsize in the section below. Additionally you
;;;; can set the cosmological parameters to your favourite model. I
;;;; have noted down ngrid, boxsize and starting redshifts for some of
;;;; the more famous simulations. The astute reader will find out that
;;;; some of these simulations were not started at a high enough
;;;; redshift or proceeded to z=0 when the fundamental mode was
;;;; already non-linear. 
;;;;
;;;;
;;;; Main-level routine to run an example code. Here I have chosen the
;;;; parameters for the Full Universe Run. - MS 06/19/2012.

compile_opt idl2,strictarrsubs
common cosmology_tf

;;;; Set ngrid and boxsize according to your simulation



;;; The DEUS Full Universe run with z_init = 100
ngrid = 8192   ;; 8192^3 simulation.
boxsize = 21d3 ;; in mpc/h -> 21 Gpc/h boxsize full universe run. 

;;; Bolshoi simulation with z_init = 80. 
; ngrid   = 2048 
; boxsize = 250.0 

;;; MultiDark Simulation with z_init = 65
; ngrid = 2048
; boxsize = 1d3 


;;; Millenium simulation with z_init = 127
; ngrid = 2160
; boxsize = 500.0 

;;; Millenium-II simulation with z_init = 127
; ngrid = 2160
; boxsize = 100.0 

;;; MXXL simulation with z_init = 63 
; ngrid = 6720
; boxsize = 4.1d3 

;;; Horizon Run 3 with z_init = 27
; ngrid = 7210
; boxsize = 10815.0

;;; Horizon Run 2 with z_init = 32
; ngrid = 6000
; boxsize = 7200.0

;;; Horizon Run 1 with z_init = 23
; ngrid = 4120
; boxsize = 6592.0

;;;; how much power are you willing to have in delta_sqr for deciding a stopping z. 
;;;; Delta_sqr = 1 means definitely non-linear (however, mode coupling actually causes 
;;;; non-linearity to set in at a lower value of delta_sqr)
delta_sqr_stop = 0.1


;;; ratio of initial rms displacement/delta_p(=boxsize/ngrid)
;;; desired. rms_ratio_start = 0.25 => shell-crossing is a 4-sigma
;;; event. You can be even more conservative and choose a smaller
;;; value of rms_ratio_start but then the starting redshift will be
;;; very high and you have to worry about integration errors in the
;;; simulation code itself. (My recommendation is 0.2-0.25)
rms_ratio_start = 0.25

;;; Set your favourite cosmological model.

;;;; WMAP-5 Cosmology
;     print,"Setting `WMAP-5' cosmology"
;     omega_m = 0.258d
;     omega_b = 0.044d
;     omega_lambda = 0.742d
;     ns = 0.96d
;     ns_bar = ns - 1.0d
;     sigma8 = 0.796d
;     norm = 1.0d
;     h = 0.719d
;     t_cmb = 2.725d
    

;;;; LAS DAMAS cosmology
print,"Setting `LAS DAMAS' cosmology"
omega_m = 0.25d
omega_b = 0.04d
omega_lambda = 0.75d
ns = 0.96d
ns_bar = ns - 1.0d
sigma8 = 0.8d
norm = 1.0d
h = 0.7d
t_cmb = 2.725d

;;; call the powspec routine
powspec,ngrid,boxsize,delta_sqr_stop,rms_ratio_start


end