————# -*-perl-*-
#
# thanks to Brad Holden for providing an independent check
# of the Astro::Cosmology code. Since the algorithm is
# essentially taken from the same source it's not completely
# independent, but was written by someone else and
# uses a different integration scheme for the calculations.
#
use
strict;
$|++;
use
Test;
plan
tests
=> 16;
use
PDL;
use
PDL::Math;
use
Astro::Cosmology;
# time for some actual comments
# This "module" is written so Doug has something to test against
# for his Astro::Cosmology package
my
$nstep
= 1e6;
my
$_tol
= 1/
$nstep
;
# this routine determines whether or not you have flat universe or
# if you have to use the correction when going between the comoving
# distance and the interesting distances (lum, ang, etc.)
sub
_not_flat {
my
$matter
=
shift
;
my
$lambda
=
shift
;
my
$tot
=
$matter
+
$lambda
;
if
(
abs
(1-
$tot
) <
$_tol
) {
return
(0);
}
elsif
(
$tot
> 1) {
return
(-1);
# negative curvature
}
else
{
return
(1);
# positive curvature
}
}
=pod
=head1 B<d_c>
Inputs
=over 4
=item $z
Input redshift
=item $matter
Omega matter - defaults to 0.3
=item $lambda
Omega lambda - defaults to 0.7
=item $nsteps
Number of integration steps - defaults to 10000
=back
Outputs
=over 4
=item comoving distance in units of the horizon distance
=back
Computation is done via direct integration using PDL's intover.
=cut
sub
d_c {
my
$z
=
shift
;
my
$matter
= 0.3;
my
$lambda
= 0.7;
my
$nsteps
= 10000;
$matter
=
shift
if
(
@_
);
$lambda
=
shift
if
(
@_
);
$nsteps
=
shift
if
(
@_
);
my
$curve
= 1.0 - (
$lambda
+
$matter
);
# here I do the actual integration, slowly, in PDL
my
$zs
= xvals
$nsteps
;
$zs
*=
$z
/
$nsteps
;
$zs
+= 0.5
*$z
/
$nsteps
;
# below is the E(z) term from David Hogg's writeup
my
$ezs
=
$matter
*(1+
$zs
)**3;
$ezs
+=
$curve
*(1+
$zs
)**2;
$ezs
+=
$lambda
;
# below is integrand for the comoving distance
$ezs
=
sqrt
(
$ezs
);
$ezs
= 1.0/
$ezs
;
$ezs
*=
$z
/
$nsteps
;
return
(intover(
$ezs
));
}
=pod
=head1 B<d_m>
Inputs
=over 4
=item $z
Input redshift
=item $matter
Omega matter - defaults to 0.3
=item $lambda
Omega lambda - defaults to 0.7
=item $nsteps
Number of integration steps - defaults to 10000
=back
Outputs
=over 4
=item transverse comoving distance in units of the horizon distance
=back
Calls B<d_c()> to get the comoving distance, then computes curvature.
Uses formula from David Hogg's astro-ph/9905116 paper.
=cut
sub
d_m {
my
$z
=
shift
;
my
$matter
= 0.3;
my
$lambda
= 0.7;
my
$nsteps
= 10000;
$matter
=
shift
if
(
@_
);
$lambda
=
shift
if
(
@_
);
$nsteps
=
shift
if
(
@_
);
# here I have to "correct" the comoving distance depending on the
# the curvature of the universe. So, I calculate the curvature,
# calculate the comoving distance along the line of sight and
# only then do compute whether the curvature is positive, negatuve
# or "zero", which means less than the tolerance variable $_tol
#
# _not_flat is the routine that checks the curvature, is uses the
# input lambda and matter, not the computed curvature
my
$curve
= 1.0- (
$lambda
+
$matter
);
my
$dc
= d_c(
$z
,
$matter
,
$lambda
,
$nsteps
);
if
(!_not_flat(
$matter
,
$lambda
)) {
return
(
$dc
);
}
elsif
(_not_flat(
$matter
,
$lambda
) < 0) {
# negative curvature
$curve
=
abs
(
$curve
);
return
(
sin
(
$dc
*sqrt
(
$curve
))/
sqrt
(
$curve
));
}
else
{
return
(sinh(
$dc
*sqrt
(
$curve
))/
sqrt
(
$curve
));
}
}
=pod
=head1 B<look_back_z>
Inputs
=over 4
=item $z
Input redshift
=item $matter
Omega matter - defaults to 0.3
=item $lambda
Omega lambda - defaults to 0.7
=item $nsteps
Number of integration steps - defaults to 10000
=back
Outputs
=over 4
=item lookback time in units of the Hubble time
=back
Computes by direct integration using PDL's intover routine, like
B<d_c()>
=cut
sub
look_back_z {
my
$z
=
shift
;
my
$matter
= 0.3;
my
$lambda
= 0.7;
my
$nsteps
= 10000;
$matter
=
shift
if
(
@_
);
$lambda
=
shift
if
(
@_
);
$nsteps
=
shift
if
(
@_
);
$nsteps
= POSIX::floor(
$nsteps
);
my
$curve
= 1.0 - (
$lambda
+
$matter
);
# here I do the actual integration, slowly, in PDL
my
$zs
= xvals
$nsteps
;
$zs
*=
$z
/
$nsteps
;
$zs
+= 0.5
*$z
/
$nsteps
;
# below is the E(z) term from David Hogg's writeup
my
$ezs
=
$matter
*(1+
$zs
)**3;
$ezs
+=
$curve
*(1+
$zs
)**2;
$ezs
+=
$lambda
;
$ezs
=
sqrt
(
$ezs
);
# below is integrand for the lookback time
$ezs
*= (1+
$zs
);
$ezs
= 1.0/
$ezs
;
$ezs
*=
$z
/
$nsteps
;
return
(intover(
$ezs
));
}
=pod
=head1 B<comoving_volume_z>
Inputs
=over 4
=item $z
Input redshift
=item $matter
Omega matter - defaults to 0.3
=item $lambda
Omega lambda - defaults to 0.7
=item $nsteps
Number of integration steps - defaults to 10000
=back
Outputs
=over 4
=item comoving volume in units of the Hubble volume
=back
Computes using an analytical formula from Carrol, Press and Turner, ARA&Ap, 1992, 30, 499
=cut
sub
comoving_volume_z {
my
$z
=
shift
;
my
$matter
= 0.3;
my
$lambda
= 0.7;
my
$nsteps
= 10000;
$matter
=
shift
if
(
@_
);
$lambda
=
shift
if
(
@_
);
$nsteps
=
shift
if
(
@_
);
$nsteps
= POSIX::floor(
$nsteps
);
my
$curve
= 1.0 - (
$lambda
+
$matter
);
my
$sacurve
=
sqrt
(
abs
(
$curve
));
my
$dm
= d_m(
$z
,
$matter
,
$lambda
,
$nsteps
);
my
$vc
;
# this formula is from Carrol, Press and Turner, ARA&Ap, 1992, 30, 499
# however, it seems odd, as the volume will be negative for negative
# curvature.
if
(!_not_flat(
$matter
,
$lambda
)) {
$vc
=
$dm
**3;
$vc
/= 3.0;
}
elsif
(_not_flat(
$matter
,
$lambda
) < 0) {
# negative curvature
$vc
=
$dm
*sqrt
(1+
$curve
*$dm
**2) - asin(
$dm
*$sacurve
)/
$sacurve
;
$vc
/= 2
*$curve
;
}
else
{
$vc
=
$dm
*sqrt
(1+
$curve
*$dm
**2) - asinh(
$dm
*$sacurve
)/
$sacurve
;
$vc
/= 2
*$curve
;
}
return
(
$vc
);
}
sub
print_output {
my
$module_val
=
shift
;
my
$local_val
=
shift
;
my
$name_str
=
shift
;
"Module "
.
$name_str
.
": "
,
$module_val
,
"\n"
;
"Local "
.
$name_str
.
": "
,
$local_val
,
"\n"
;
my
$adiff
=
abs
(
$local_val
-
$module_val
);
"Absolute value of difference: "
,
sprintf
(
"%.3g"
,
$adiff
),
"\n"
;
if
(
$adiff
< 1e-5) {
"Absolute value of difference within tolerance for $name_str\n"
;
}
else
{
"!!**Absolute value of difference exceeds tolerance for $name_str **!!\n"
;
}
ok(
$adiff
< 1.0e-5 );
}
# Tests Doug's code by comparing the local subroutines to values returned
# by dougs Module.
sub
run_test {
my
$z
=
shift
;
my
$om
=
shift
;
my
$ol
=
shift
;
my
$nstep
=
shift
;
my
$ac_obj
= Astro::Cosmology->new({
Matter
=>
$om
,
Lambda
=>
$ol
,
H0
=>0});
$ac_obj
,
"\n"
;
my
$m_ld
=
$ac_obj
->lum_dist(
$z
);
my
$m_ad
=
$ac_obj
->adiam_dist(
$z
);
my
$m_lbt
=
$ac_obj
->lookback_time(
$z
);
my
$m_vc
=
$ac_obj
->comov_vol(
$z
);
my
$l_dm
= d_m(
$z
,
$om
,
$ol
,
$nstep
);
my
$l_ld
=
$l_dm
*(1+
$z
);
my
$l_ad
=
$l_dm
/(1+
$z
);
my
$l_lbt
= look_back_z(
$z
,
$om
,
$ol
,
$nstep
);
my
$l_vc
= comoving_volume_z(
$z
,
$om
,
$ol
,
$nstep
);
print_output(
$m_ld
,
$l_ld
,
"luminosity distance"
);
print_output(
$m_ad
,
$l_ad
,
"angular diameter distance"
);
print_output(
$m_lbt
,
$l_lbt
,
"lookback time"
);
print_output(
$m_vc
,
$l_vc
,
"comoving volume"
);
"\n"
;
}
# Canonical flat universe
run_test(1.235,0.3,0.7,1e6);
# Open universe using best fit values from my thesis
run_test(1.235,0.32,0.0,1e6);
# weird closed universe
run_test(1.235,1.32,0.0,1e6);
# wacky open universe
run_test(1.235,1.32,-1.0,1e6);