yorick banner

Home

Manual

Packages

Global Index

Keywords

Quick Reference


H0_cst=67;


func cosmo_param (&Omega_m,&Omega_lambda,&Omega_k,&q0)
/* DOCUMENT cosmo_param(Omega_m,Omega_lambda,Omega_k,q0)

   Same as IDL routine:
   Given any two of the four input parameters
      (1) the normalized matter density Omega_m
      (2) the normalized cosmolgical constant, Omega_lambda
      (3) the normalized curvature term, Omega_k
      (4) the deceleration parameter q0
   this program will derive the remaining two.
   Here "normalized" means divided by the closure density
   so that Omega_m + Omega_lambda + Omega_k = 1. For a more
   precise definition see Caroll, Press, & Turner (1992, ArAA, 30, 499).
   
   If less than two parameters are defined, this procedure sets default 
   values of Omega_k=0 (flat space), Omega_lambda = 0.7, Omega_m = 0.3
   and q0 = -0.5 in these order (by avoiding parameter defined) until
   two parameters are defined.

   If more than two parameters are defined upon input (overspecification), 
   then the first two defined parameters in the ordered list Omega_m, 
   Omega_lambda, Omega_k, q0 are used to define the cosmology.

   INPUTS-OUTPUTS:
      Omega_m     : Normalized matter energy density
                    non-negative numeric scalar
      Omega_lambda: Normalized cosmological constant
                    numeric scalar
      Omega_k     : Normalized curvature parmeter
                    numeric scalar. This is  zero
                    for a flat universe
      q0          : Deceleration parameter
                    numeric scalar
                       = 0.5*Omega_m - Omega_lambda

*/
{
  local Nk,Nl,Nm,Nq;
  Nk=numberof(Omega_k)!=0;
  Nl=numberof(Omega_lambda)!=0;
  Nm=numberof(Omega_m)!=0;
  Nq=numberof(q0)!=0;
  NSum=Nk+Nl+Nm+Nq;
  if(NSum==0)
    {
      Omega_k=0.;
      Omega_lambda=0.7;
    }

  if(NSum==1)
    {
      if(Nk) Omega_lambda=min(0.7,1-Omega_k);
      if(Nl) Omega_k=0;
      if(Nm)
        {
          Omega_k=0;
          Omega_lambda=1-Omega_m;
        }
      if(Nq)
        {
          Omega_k=0;
          Omega_lambda=(1-2*q0)/3.;
        }
    }

  if(NSum==2)
    {
      if(Nk&&Nl) {} // do nothing!
      if(Nk&&Nq) Omega_lambda=(1-Omega_k-2*q0)/3.;
      if(Nk&&Nm) Omega_lambda=1-Omega_m-Omega_k;

      if(Nl&&Nq) Omega_k=1-2*q0-3*Omega_lambda;
      if(Nl&&Nm) Omega_k=1-Omega_m-Omega_lambda;

      if(Nq&&Nm)
        {
          Omega_lambda=0.5*Omega_m-q0;
          Omega_k=1-Omega_m-Omega_lambda;
        }
    }
  if(NSum==3)
    {
      if(Nk) Omega_k=1-Omega_m-Omega_lambda;
      if(Nl) Omega_lambda=1-Omega_m-Omega_k;
    }

  Omega_k=double(Omega_k);
  Omega_lambda=double(Omega_lambda);
  Omega_m=double(1-Omega_k-Omega_lambda);
  q0=double(0.5*Omega_m-Omega_lambda);
  return;
}

func computeH (z,Om,Ok,Ol,H0=)
{
  if(is_void(H0)) H0=H0_cst;

  fct=sqrt(Om*(1+z)^3+Ok*(1+z)^2+Ol);
  return H0*fct;
}
func ldist (z,q0,lambda0)
/* DOCUMENT  ldist(z,q0=,lambda0=)

   Compute the term to integrate to compute the luminosity
   distance. z can be a VECTOR of values.
   
   KEYWORDS: q0      : Deceleration parameter,numeric scalar
                       -a*(a'')/(a')^2 (Omega_m/2-lambda0)
                        
             lambda0 : Cosmological constant, normalized to
                       th closure density.
   
   SEE ALSO: lumdist
*/
{
  local good,bad,denom;
  local term1,term2,term3;
  local out;
  
  term1=(1+z)^2;
  term2=1+2*(q0+lambda0)*z; 
  term3=z*(2+z)*lambda0;
  denom=term1*term2-term3;
  out=array(structof(z(1)),dimsof(z));
  if(numberof((good=where(denom>0)))) out(good)=1/sqrt(denom(good));
  return out;
}


func uage (z,q0,lambda0)
/* DOCUMENT  univAge(z,q0=,lambda0=)

   Compute the term to integrate to compute the age of the univers.
   z can be a VECTOR of values.
   
   KEYWORDS: q0      : Deceleration parameter,numeric scalar
                       -a*(a'')/(a')^2 (Omega_m/2-lambda0)
                        
             lambda0 : Cosmological constant, normalized to
                       th closure density.
   
   SEE ALSO: lumdist
*/
{
  local good,bad,denom;
  local term1,term2,term3;
  local out;
  
  term1=(1+z)^2;
  term2=1+2*(q0+lambda0)*z; 
  term3=z*(2+z)*lambda0;
  denom=term1*term2-term3;
  out=array(structof(z(1)),dimsof(z));
  if(numberof((good=where(denom>0)))) out(good)=1/((1+z(good))*sqrt(denom(good)));
  return out;
}


func lumdist (z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)
/* DOCUMENT lumdist(z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)

   See IDL routine...
   PURPOSE: 
     Calculate luminosity distance (in Mpc) of an object given its redshift 
   EXPLANATION:
     The luminosity distance in the Friedmann-Robertson-Walker model is 
     taken from  Caroll, Press, and Turner (1992, ARAA, 30, 499), p. 511
     Uses a closed form (Mattig equation) to compute the distance when the 
     cosmological constant is zero.   Otherwise integrates the function using
     simpson_cosmo.	
   EXAMPLE:
     Plot the distance of a galaxy in Mpc as a function of redshift out
     to z = 5.0, assuming the default cosmology (Omega_m=0.3, Lambda = 0.7,
     H0 = 70 km/s/Mpc)
     
     z = span(0,5,50);
     plg,lumdist(z),z;
     xytitles,"z","Distance (Mpc)"
  
   SEE ALSO: ldist
 */
{
  local SPEEDOFLIGHT;
  SPEEDOFLIGHT=2.9979e5;
  cosmo_param,Omega_m,lambda0,k,q0;
  if(is_void(h0)) h0=H0_cst;
  if(is_void(silent)||!silent)
    write,format="H0:%3i Omega_m:%5.2f Lambda0:%5.2f q0:%5.2f k:%5.2f\n",h0,
      Omega_m,lambda0,q0,k;

  if(lambda0==0)
    {
      denom=sqrt(1+2*q0*z)+1+q0*z;
      return SPEEDOFLIGHT*z/h0*(1+z*(1-q0)/denom);
    }

  dlum=array(0.,dimsof(z));
  for(i=1;i<=numberof(z);i++)
    if(z(i)<=0)
      dlum(i)=0;
    else
      dlum(i)=simpson_cosmo(ldist,0,z(i),q0,lambda0);

  if(k>0)
    dlum = sinh(sqrt(k)*dlum)/sqrt(k);
  else if(k<0)
    dlum = sin(sqrt(-k)*dlum)/sqrt(-k);

  return dlum*SPEEDOFLIGHT*(1+z)/h0;
}

func univAge (z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)
/* DOCUMENT univAge(z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)

   See IDL routine...
   PURPOSE: 
     Calculate age of the universe in s
   EXPLANATION:
     The luminosity distance in the Friedmann-Robertson-Walker model is 
     taken from  Caroll, Press, and Turner (1992, ARAA, 30, 499), p. 511
     Uses a closed form (Mattig equation) to compute the distance when the 
     cosmological constant is zero.   Otherwise integrates the function using
     simpson_cosmo.	
   EXAMPLE:
     Plot the distance of a galaxy in Mpc as a function of redshift out
     to z = 5.0, assuming the default cosmology (Omega_m=0.3, Lambda = 0.7,
     H0 = 70 km/s/Mpc)
     
     z = span(0,5,50);
     plg,lumdist(z),z;
     xytitles,"z","Distance (Mpc)"
  
   SEE ALSO: ldist
 */
{
  local SPEEDOFLIGHT;
  SPEEDOFLIGHT=2.9979e5;
  cosmo_param,Omega_m,lambda0,k,q0;
  if(is_void(h0)) h0=H0_cst;
  if(is_void(silent)||!silent)
    write,format="H0:%3i Omega_m:%5.2f Lambda0:%5.2f q0:%5.2f k:%5.2f\n",h0,
      Omega_m,lambda0,q0,k;

  auni=array(0.,dimsof(z));
  for(i=1;i<=numberof(z);i++)
    if(z(i)<=0)
      auni(i)=0;
    else
      auni(i)=simpson_cosmo(uage,0,z(i),q0,lambda0);

//   if(k>0)
//     dlum = sinh(sqrt(k)*dlum)/sqrt(k);
//   else if(k<0)
//     dlum = sin(sqrt(-k)*dlum)/sqrt(-k);

  return auni/h0*3.0856775876793e+19; //(Mpc/km)
}


func angdist (z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)
/* DOCUMENT angdist(z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)

   Compute the Angular diameter distance (lumdist/(1+z)^2)
     
   SEE ALSO: lumdist,ldist,cosmo_param
 */
{
  return lumdist(z,h0=h0,k=k,lambda0=lambda0,Omega_m=Omega_m,q0=q0,silent=silent)/(1+z)^2;
}


func angsize (dl,z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)
/* DOCUMENT  angsize(dl,z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)

   Compute the apparent angular size of an object as a function of
   redshift.
   INPUTS:   dl - proper size of the object (in kpc)
              z - redhsift values 

   EXAMPLE: plot the angular size of a 50 kpc diameter galaxy as
            a function of redshift for the default cosmology
            (Lambda = 0.7, Omega_m=0.3) up to z = 5
            
        z = span(0.1,5,100)    //Angular size undefined at z = 0
        pg,angsize(50,z),z;
        xytitles,"z","Angular Size (\")"
   
              
   SEE ALSO: lumdist,ldist,cosmo_param,angdist,propersize
 */
{
  local factor;
  local RADEG3600;
  RADEG3600=206264.806247;
  factor=RADEG3600*(dl/1000.); //lumdist return Mpc...
  return factor/angdist(z,h0=h0,k=k,lambda0=lambda0,Omega_m=Omega_m,q0=q0,silent=silent);
}

func propersize (dtheta,z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)
/* DOCUMENT  propersize(dl,z,h0=,k=,lambda0=,Omega_m=,q0=,silent=)

   Compute the seperation (in kpc) of two line of sight sparated
   by dtheta arcseconds.
   INPUTS:  dtheta - angular separation (in ")
                 z - redhsift values 

   EXAMPLE: plot the separation in kpc beetween two line of sight
            separated with 2" as a function of redshift for the
            default cosmology (Lambda = 0.7, Omega_m=0.3) up to
            z = 5
            
        z = span(0,5,100)
        plg,propersize(2,z),z;
        xytitles,"z","Proper Size (kpc)"
   
              
   SEE ALSO: lumdist,ldist,cosmo_param,angdist,angsize
 */
{
  local factor;
  local RADEG3600;
  RADEG3600=206264.806247;
  factor=1000*dtheta/RADEG3600; //1000 because angdist return Mpc...
  return factor*angdist(z,h0=h0,k=k,lambda0=lambda0,Omega_m=Omega_m,q0=q0,silent=silent);
}


func simpson_cosmo (function, a, b, q0,lambda0,epsilon,notvector=)
/* DOCUMENT integral= simpson(function, a, b)
         or integral= simpson(function, a, b, epsilon)
     returns the integral of FUNCTION(x) from A to B.  If EPSILON is
     given, Simpson's rule is refined until that fractional accuracy
     is obtained.  EPSILON defaults to 1.e-6.

     If the notvector= keyword is supplied and non-zero, then FUNCTION
     may not be called with a list of x values to return a list of
     results.  By default, FUNCTION is assumed to be a vector function.

     If the function is very smooth, romberg may work better.

   SEE ALSO: romberg, max_doublings
 */
{
  if (!epsilon || epsilon<0.) epsilon= 1.e-6;
  a= double(a);
  b= double(b);
  ost= os= -1.e-30;
  for (i=1 ; i<=max_doublings ; ++i) {
    st= trapezoid_cosmo(function, a, b, i,q0,lambda0, notvector);
    s= (4.*st - ost)/3.;
    if (abs(s-os) < epsilon*abs(os)) return s;
    os= s;
    ost= st;
  }
  error, "no convergence after "+pr1(2^i)+" function evaluations";
}

local max_doublings ;
/* DOCUMENT max_doublings= 20
     is the maximum number of times romberg or simpson will split the
     integration interval in half.  Default is 20.
 */
max_doublings= 20;

func trapezoid_cosmo (function, a, b, n, q0,lambda0,notvector)
{
  extern _trapezoid_i, _trapezoid_s;
  if (n==1) {
    _trapezoid_i= 1;
    return _trapezoid_s= 0.5*(b-a)*(function(a,q0,lambda0)+function(b,q0,lambda0));
  }
  dx= (b-a)/_trapezoid_i;
  if (!notvector) {
    /* conserve memory-- dont try more than 8192 points at a time */
    if (_trapezoid_i <= 8192) {
      s= sum(function(span(a,b,_trapezoid_i+1)(zcen),q0,lambda0));
    } else {
      x= a+(indgen(8192)-0.5)*dx;
      s= sum(function(x,q0,lambda0));
      dx2= 8192*dx;
      for (i=8193 ; i<=_trapezoid_i ; i+=8192) s+= sum(function((x+=dx2),q0,lambda0));
    }
  } else {
    x= a+0.5*dx;
    s= function(x,q0,lambda0);
    for (i=2 ; i<=_trapezoid_i ; ++i) s+= function((x+=dx),q0,lambda0);
  }
  _trapezoid_i*= 2;
  return _trapezoid_s= 0.5*(_trapezoid_s + dx*s);
}