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);
}
|