view release on metacpan or  search on metacpan

t/cosmology2.t  view on Meta::CPAN

  $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

t/cosmology2.t  view on Meta::CPAN

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