# -*-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 POSIX qw();
# 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;
print "Module ".$name_str.": ",$module_val,"\n";
print "Local ".$name_str.": ",$local_val,"\n";
my $adiff = abs($local_val-$module_val);
print "Absolute value of difference: ",
sprintf("%.3g",$adiff),"\n";
if ($adiff < 1e-5) {
print "Absolute value of difference within tolerance for $name_str\n";
} else {
print "!!**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});
print $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");
print "\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);