/*

Code for calculating the mean initial and final density profiles
around a virialized dark matter halo. The initial profile is derived
from the statistics of the initial Gaussian random field, accounting
for the problem of peaks within peaks using the extended
Press-Schechter model. Spherical collapse then yields the typical
density and velocity profiles of the gas and dark matter that
surrounds the final, virialized halo. In additional to the mean
profile, +/- 1-sigma profiles are calculated and can be used as an
estimate of the scatter. The code was developed by Rennan Barkana,
2002/3.

The corresponding paper [MNRAS 347, 59 (2004)], entitled " A model for
infall around virialized halos", is available on the astro-ph archive,
at: http://arxiv.org/abs/astro-ph/0212458 .

*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "Infall.h"

#define IntN 400
#define FMT "%-11.6e "
#define FMT4 "%-11.6e %-11.6e %-11.6e %-11.6e "
#define FMT8 "%-11.6e %-11.6e %-11.6e %-11.6e %-11.6e %-11.6e %-11.6e %-11.6e "

/* 
   Compile:
   gcc -O Infall.c -o Infall.out -lm

   Files:
   Infall.c InfallSub.c Infall.h 
*/

long idum,iWhite,ifilter;

double OMEGA0,OMEGAB,OMEGALAMBDA,OMEGAR,hubble,sigma8,
  THETA27,RUNNINGnS,RUNNINGk0,RUNNINGdn,betac,alpha_gamma,
  sound_horizon_fit,omhh,omega_matter,omega_lambda,tilt,tiltk0,
  tiltdn,omega_baryon,snorm,omega_curv,scale;  

double ageFac1,ageFac2,ageFac3,zhalo,omega,omegal,nk0,ndndlk,
  omegab,hub,sig8,npower;
  
int IntProf,iInfall,ido,InfSig;

double Mhalo,SphCz1,SphCRi,MProfR,dpp,SphCtf,SphCHi,SphCk,
  SphCnorm,SphCz2,Pass1,Pass2,Pass3,MProfx[IntN3],MProfy[IntN3],
  MProfy2[IntN3],PassXi1,PassXi2,PassXi4;

main()
{
  double t1,t2,z2,x1,x2,M,t3,t4,t5,t6,t7,t9,t10,y1,t11,t12,t13,t15,
    t20,maxRm1,Mtop,t21;
  int i,menu;

  iInfall=1;
  IntProf=IntN3;

  printf("Enter Cosmological Parameters:\n");
  printf("Note!! Flat Lambda-CDM universe assumed in various places.\n");
  printf("Parameter Menu:\n"); 
  printf("0: Enter all parameters.\n");
  printf("1: Use parameters from the 5-yr WMAP paper (WMAP5+SN+BAO).\n");
  printf("\n");
  scanf("%i",&menu);
  if (menu == 0)
    {
      printf("Enter Omega_0, Om_L, Om_b, h, sigma_8, tilt\n");
      scanf("%lg %lg %lg %lg %lg %lg",&omega,&omegal,&omegab,&hub,
	    &sig8,&npower);
    }
  else 
    {
      /* WMAP5 (+SN+BAO) power-law LCDM from Komatsu et al. 2008: */
      omega=.28;
      omegal=.72;
      omegab=.046;
      hub=.7;
      sig8=.82;
      npower=0.96;
      printf("WMAP5: LCDM(Om=0.28), Om_b=.046, h=.7, sig8=.82, n=.96\n");
    } 
  nk0=0.05; 
  ndndlk=0.; /* Assume no running of the spectral index. */
  Initialize();
  ageFac1=omega_lambda / omega_matter;
  ageFac2=(9.7776e9/hubble)*2./(3.*sqrt(omega_lambda));
  ageFac3=9.7776e9/hubble;

  printf("Choose which initial profile to use: (0:mean, -1:-1 sigma, 1:+1 sigma)\n");
  scanf("%i",&ido);
  printf("Enter desired virial redshift, and halo virial mass M [solar units]\n");
  scanf("%lg %lg",&z2,&Mhalo);
  maxRm1=300.;
  printf("Now finding initial M_tophat that yields final M_vir...\n");
  InfSig=0;
  if (ido == -1)
    InfSig=-1;
  else if (ido == 1)
    InfSig=1;
  M=Mhalo;
  Mtop=M;
  (z2 > 1.) ? GenMeanProf(RofM(M),.00099,&maxRm1,z2) :
    GenMeanProf(RofM(M),.000099,&maxRm1,z2);
  InitSphCol(z2,M,&t2,&t3,&t4);
  t5=FindMvir(z2);
  Mtop=M/t5;
  maxRm1=300.;
  (z2 > 1.) ? GenMeanProf(RofM(Mtop),.00099,&maxRm1,z2) :
    GenMeanProf(RofM(Mtop),.000099,&maxRm1,z2);
  InitSphCol(z2,Mtop,&t2,&t3,&t4);
  t6=FindMvir(z2);
  /* Tried two top-hat masses, got corresponding virial masses. Now linearly interpolate 
     (in the log) between the two guesses to get an accurate final top-hat value: */
  t3=log(t5)/log(SQR(t5)/t6);
  t4=log(M)-t3*log(M*t5);
  Mtop=exp(t4+t3*log(M));
  printf("M_tophat = "FMT"\n",Mtop);  

  printf("\nMain Menu: \n");
  printf("Find initial infall profile                          [1]\n");
  printf("Find final profile after spherical collapse          [2]\n");
  scanf("%i",&menu);
  printf("\n");
  if (menu==1)
    {
      M=Mtop;
      zhalo=z2;
      iInfall=0;
      t1=RofM(M);
      t7=sigmaOfM(M);
      t9=wofz(zhalo);
      maxRm1=300.;
      GenMeanProf(t1,.00099,&maxRm1,zhalo);
      printf("Initial rM [Mpc], z=0 sigma(rM), critical threshold nu = \n");
      printf(FMT FMT FMT"\n",t1,t7,t9);
      printf("r/rM, delta(r)/nu (asymptotic), delta(r)/nu =\n");
      x1=.001;
      x2=DMIN(maxRm1-.01,100.);
      for (i=0; i<136; i++)
	{
	  t2=t1*(1.+exp(log(x1)+log(x2/x1)*(i/(IntN3-1.))));
	  M=MofR(t2);
	  t6=XiOfR1R2(t1,t2); 
	  t3=sigmaOfM(M);
	  t10=SQR(t9/t7)*(t6/SQR(t3))*(1.-t6/SQR(t7))/(1.-SQR(t6/(t7*t3)));
	  if (InfSig == -1) t11=NewInfallm1sig(t6/SQR(t7),t10);
	  else if (InfSig == 1) t11=NewInfallp1sig(t6/SQR(t7),t10);
	  else t11=NewInfall(t6/SQR(t7),t10);
	  printf(FMT FMT FMT"\n",t2/t1,MeanProf(t2),t11);
	  /* Initial overdensity is 1.69/D(z).           */
	}
      }
  else
    {
      printf("Running spherical collapse...\n");
      maxRm1=300.;
      (z2 > 1.) ? GenMeanProf(RofM(Mtop),.00099,&maxRm1,z2) :
	GenMeanProf(RofM(Mtop),.000099,&maxRm1,z2);
      InitSphCol(z2,Mtop,&t2,&t3,&t4);
      t15=t3*cbrt(Mhalo/Mtop);
      t21=6.558e-5*sqrt(Mhalo/t15);
      printf("Initial rM [com Mpc], final R_vir [phys Mpc], z=0 sigma(rM), virial circular velocity V_c [km/s] =\n");
      printf(FMT4"\n",SphCRi,t15,sigmaOfM(Mtop),t21);
      printf("Mass shells:\nt=0 r [com Mpc], M [solar], del_bar_i, delta_i, R/");
      printf("R_vir, del_bar, delta, V_pec/V_c, V_tot/V_c =\n");
      x1=(z2 > 1.) ? .001 : 0.0001;
      x2=maxRm1-.01;
      for (i=0; i<IntN3; i++)
	{
	  y1=1.+exp(log(x1)+log(x2/x1)*(i/(IntN3-1.)));
	  /* Now integrate the differential equations, in order to include the cosmological constant */
	  SphColODE(y1,Mtop,&t1,&t7,&t2,&t5,&t11,&t12,&t3,&t6,&t13,&t20,1.);
	  printf(FMT8 FMT"\n",
		 t1,t2,t7,t12,t5/t15,t11,t3,t6/t21,t13/t21);
	}  
    }
  /* End of main program */
}

#include "InfallSub.c"