------ 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
# convert degrees to radians
def deg_to_rad(deg)
rad eg.to_f * (Math::PI / 180.0)
end
# convert radians to degrees
def rad_to_deg(rad)
deg ad.to_f * (180.0/Math::PI)
end
# convert radians to arcminuts
def rad_to_arcmin(rad)
rad_to_deg(rad) * 60.0
end
# convert radians to arcseconds
def rad_to_arcsec(rad)
rad_to_deg(rad) * 3600.0
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( h 71, omega_matter 27, omega_lambda 73, omega_k 0 )
@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(z 0)
@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.
def phys_size(rad, z, z_obs _obs)
phys_size ist_ang(z,z_obs) * rad.to_f
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)
rad hys_size.to_f / dist_ang(z,z_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 z z} at z_obs z_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 z fmt}) is #{fmt} degrees" %
[z, cosmo.z_obs, rad_to_deg(cosmo.ang_horizon(z))]
cosmo.z_obs .0
puts "Angular size of the horizon at z {fmt} (from z fmt}) is #{fmt} degrees" %
[z, cosmo.z_obs, rad_to_deg(cosmo.ang_horizon(z))]
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--