```------art_20169_8026061.1146178118120
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
Content-Disposition: inline

I'm doubtful of any interest, but in the odd chance that someone else
might want to play, I'm posting this code.  If someone does look, all
comments are welcome (be it ruby style or corrections)!

Attached is a small library that calculates comological distances
(like light travel time, comoving distances, age of the universe,
angular size of the horizon) given a specific cosmology.  This is
something I threw together for personal use in solving a homework
problem or two and making some quick estimates.  References are
included in the comments for further information.

It requires the ruby GSL bindings (http://rb-gsl.rubyforge.org) for
numerical integration.

The library also includes some sample usage at the end (which is
invoked if the script is just run from the command line).

Cameron

p.s. to interested parties, the science / calculations are pretty well
described by Ned Wright's page:
http://www.astro.ucla.edu/~wright/CosmoCalc.html

------art_20169_8026061.1146178118120
Content-Type: application/x-ruby; name=cosmo_dist.rb
Content-Transfer-Encoding: 7bit
X-Attachment-Id: f_emjomqbj
Content-Disposition: attachment; filename="cosmo_dist.rb"

#!/usr/bin/env ruby

# These are some basic routines for cosmology distance and time estimates.
#   Cameron McBride (cameron.mcbride / gmail.com)
#   April 2006

# use the GSL integration routines and some physical constants
require 'rbgsl'

# Simple module to handle angle conversions for some standard astronomical units
module AstroUnits

rad  eg.to_f * (Math::PI / 180.0)
end

end

end

end

# convert parsecs to light-years
def pc_to_ly(pc)
pc.to_f * 3.26163
end

# convert light-years to parsecs
def ly_to_pc(ly)
ly.to_f / 3.26163
end

end

# Implementation of some cosmological distance and time estimators, based on the
# content of the universe SINCE the last scattering surface (CMB, z ~ 1090).
# Small contributions are neglected, and radiation content of the universe is ignored
# (which is an acceptable assumption since recombination).
#
# Although this is standard cosmology, these routines are based on arguments
# presented in:
# Hogg, D. 2000:: "Distance Measures in cosmology", http://arxiv.org/abs/astro-ph/9905116
# Weinberg, S., 1972:: _Gravitation_ _and_ _Cosmology_, John Wiley & Sons, New York.
#
# Results were checked / compared to Ned Wright's cosmology calculator
# * Ned Wright's http://www.astro.ucla.edu/~wright/CosmoCalc.html
# (if you look at the javascript, he is much more careful about small details)
#
# All input "times" are in numerical values of redshift, which is related to the scale factor.
# See the above references for further definitions.
#
# Although the units can be changed (see @h_units, @t_units, @pc_to_l), the default output
# is in:
#   distance: Mpc  (1 mega parsec  e6 pc; 1pc  .26 light years)
#   time:     Myrs (1 mega year   million years)
#
# Given the simplicity of the calculations, none of these estimates should be taken as more
# accurate than around 2-3 significant digits (~1% estimates).
class CosmoDist

# this is E(z) as described by Hogg EQN 14
def self.chi(z, om, ol, ok)
zf  .to_f # ensure it's a float
Math.sqrt( om*(1.0+zf)**3 + ol + ok*(1.0+zf)**2 )
end

# effective infinite redshift
Z_INF  e20

# Comoving distance integrand.
CDI_func  ambda { |z, omegas| 1.0 / chi(z, *omegas) }

# Lookback Time Integrand based on redshift.
LTI_func  ambda { |z, omegas| 1.0 / chi(z, *omegas) / (1.0 + z) }

attr_accessor :h, :omega_m, :omega_l, :omega_k, :z_obs
attr_accessor :c, :h_units, :t_units, :pc_to_l

# initialize content of the universe.
#   Ho   * 100 (km/s/Mpc)
#   omega_matter  ho_matter / rho_critical -> matter content
#   omega_lambda  ho_lambda / rho_critical -> dark energy content
#   omega_k -> curvature
def initialize( h71, omega_matter27, omega_lambda73, omega_k0 )
@h
@omega_m  mega_matter
@omega_l  mega_lambda
@omega_k  mega_k

@z_obs    .0  # redshift we're "observing" from

@c        SL::CONST::MKSA::SPEED_OF_LIGHT / 1000.0 # km/s

@h_units  00.0                    # Ho   * h_units (km/s/Mpc)
@t_units  3.09e19 / 3.15e7) / 1e6 # convert to Myrs from s * (Mpc / km)
@pc_to_l  e6                      # parsecs to Mpc
end

# the total energy density of all the respective parts.
#   Omega > 1 : closed universe
#   Omega < 1 : open universe
#   Omega   : critical density (what we currently believe is true)
def omega
@omega_m + @omega_l + @omega_k
end

# the Hubble parameter as a function of redshift.
#   H(z)  o * E(z) / E(0)
def hubble_param(z0)
@h * @h_units *
self.class.chi(z  , @omega_m, @omega_l, @omega_k) /
self.class.chi(0.0, @omega_m, @omega_l, @omega_k)
end

# the Hubble time (today), in Myrs
#   th   / Ho (Hogg, EQN 3)
def time_hubble
@t_units / (@h * @h_units)
end

# the Hubble distance (today), in Mpc
#   Dh   / Ho (Hogg, EQN 4)
def dist_hubble
@c / (@h * @h_units)
end

# integrate 1/E(z) without multiplying by the Hubble distance
def dist_co_frac(z, z_obs_obs)
di  SL::Function.alloc CDI_func
di.set_params([@omega_m, @omega_l, @omega_k])

di.qag([z_obs.to_f,z.to_f]).first
end

# the comoving line-of-sight distance, in Mpc
#   Dc in Hogg, EQN 15
def dist_co_los(z, z_obs_obs)
dc  ist_hubble * dist_co_frac(z,z_obs)
end

# the comoving transverse distance, in Mpc
# This is also called the proper motion distance.
# (Hogg, section 6: Dm in EQN 16).
# This value is identical to the LOS comoving distance without any curvature.
def dist_co_trans(z, z_obs_obs)
dc  ist_co_los(z,z_obs)
if @omega_k.zero?
dm  c
else
dfrac  ist_hubble / Math.sqrt(@omega_k.abs)
if @omega_k > 0
dm  frac * Math.sinh( dc / dfrac )
elsif @omega_k < 0
dm  frac * Math.sin(  dc / dfrac )
else
raise "Don't know what to do with Omega_k  {omega_k}"
end
end

dm
end

# calculate the angular diameter distance, in Mpc
# (section 6, Hogg)
#   Da  m / (1 + z) from  EQN 18
def dist_ang(z, z_obs_obs)
dist_co_trans(z,z_obs) * (1 + z_obs) / (1.0 + z)
end

# calculate the luminosity distance, in Mpc
# (section 7, Hogg)
#   Dl  m * (1 + z) from  EQN 21
def dist_lum(z, z_obs_obs)
dist_co_trans(z,z_obs) * (1.0 + z) / (1.0 + z_obs)
end

# calculate the distance modulus, output is difference of magnitudes
# (aka just a number based on the log of luminosity, normalize so it's related
#  to distance)
# (EQN 25, Hogg)
def dist_mod(z,z_obs_obs)
5.0 * Math.log10( dist_lum(z,z_obs) / (10*@pc_to_l)  )
end

# calculate the lookback time, in Myrs
# (section 10, Hogg)
#   t_l in EQN 30
def time_lookback(z, z_obs_obs)
ti  SL::Function.alloc LTI_func
ti.set_params([@omega_m, @omega_l, @omega_k])

t   i.qag([z_obs.to_f, z.to_f]).first
t * ime_hubble
end

# this is the effective scale factor, really the ratio of a_1 / a_obs
def scale_factor(z, z_obs_obs)
(1.0 + z_obs) / (1.0 + z)
end

# calculate the age of the universe at a given redshift, in Myrs
def age(z)
time_lookback(Z_INF,0.0) - time_lookback(z,0)
end

# calculate the physical size of an object based on it's angular size (in radians)
# at a given redshift.  assumes small angle approximation.
end

# angular size of an object (in radians) given a physical size (in Mpc)
# defined from the angular diameter distance, and discussed in
# Weinberg, pg 422 (EQN 14.4.16)
def ang_size(phys_size, z, z_obs_obs)
end

# physical size of the horizon, in Mpc
def phys_size_horizon(z,z_obs_obs)
dist_co_trans(Z_INF,z) * scale_factor(z,z_obs)
end

# the angular size of the horizon at the given redshift (z) when observed from
# a specific redshift (z_obs).
def ang_horizon(z, z_obs_obs)
# this really is the ratio of comoving distances here.
# in general:
#   theta  ist_phys / dist_ang (as seen in ang_size())
# but the physical distance of the horizon is
#   d_phys  _1 * a_1 -> cordinate size times the scale factor at that point
# therefore
#   theta  _1 * a_1 / r_2 * a_1   _1 / r_2 -> ratio of cordinate distance
# aka comoving distances
ang  ist_co_frac(Z_INF,z) / dist_co_frac(z,z_obs)

# some sanity checks
if(ang > Math::PI)
# we know this formulation is divergent, throw an error if we go there
raise "Horizon angle divergent! (for zz} at z_obsz_obs})"
elsif(ang > 0.01)
ang  ath.sin(ang)
end

ang
end

end

if \$0 __FILE__
# Example usage

include AstroUnits
fmt  %-.5g"

cosmo  osmoDist.new

z  100 # CMB: surface of last scattering
puts "Angular size of the horizon at z  {fmt} is #{fmt} degrees" % [z, rad_to_deg(cosmo.ang_horizon(z))]
puts "  physical size of horizon:  #{fmt} Mpc" % cosmo.phys_size_horizon(z)
puts "  angular diameter distance: #{fmt} Mpc" % cosmo.dist_ang(z)

puts
puts "Age of universe at z  {fmt} #{fmt} Myrs" % [z, cosmo.age(z)]
puts "Age of universe now:        #{fmt} Myrs" % [cosmo.age(0.0)]

puts
z  0
puts "Angular size of the horizon at z  {fmt} (from zfmt}) is #{fmt} degrees" %
cosmo.z_obs  .0
puts "Angular size of the horizon at z  {fmt} (from zfmt}) is #{fmt} degrees" %

def report_stats(cosmo,z)
puts "Stats for z  g -> %g:" % [cosmo.z_obs,z]
puts "  universe age at z:    %6i Myrs" % [z,cosmo.age(z)]
puts "  lookback time:          %6i Myrs" % cosmo.time_lookback(z)
puts "  comoving LOS distance: %7i Mpc (%5.2f billion light years)" %
[d  osmo.dist_co_los(z), pc_to_ly(d)/1e3]
puts "  angular diameter dist: %7i Mpc (%5.2f billion light years)" %
[d  osmo.dist_ang(z),    pc_to_ly(d)/1e3]
puts "  luminosity distance:   %7i Mpc (%5.2f billion light years)" %
[d  osmo.dist_lum(z),    pc_to_ly(d)/1e3]
puts
end

puts
z, cosmo.z_obs  .0,0.0
report_stats(cosmo,z)
z  .0
report_stats(cosmo,z)
cosmo.z_obs  .0
report_stats(cosmo,z)

end

------art_20169_8026061.1146178118120--

```