From: Cameron McBride <cameron.mcbride@...>
Date: 2006-04-28T07:48:40+09:00
Subject: Cosmological distances with ruby

------=_Part_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

------=_Part_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 = deg.to_f * (Math::PI / 180.0)
  end

  # convert radians to degrees
  def rad_to_deg(rad)
    deg = rad.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 = 1e6 pc; 1pc = 3.26 light years)
#   time:     Myrs (1 mega year = 1 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 = z.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 = 1e20 

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

  # Lookback Time Integrand based on redshift.
  LTI_func = lambda { |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 = h * 100 (km/s/Mpc)
  #   omega_matter = rho_matter / rho_critical -> matter content
  #   omega_lambda = rho_lambda / rho_critical -> dark energy content
  #   omega_k -> curvature
  def initialize( h=0.71, omega_matter=0.27, omega_lambda=0.73, omega_k=0.0 )
    @h       = h
    @omega_m = omega_matter
    @omega_l = omega_lambda
    @omega_k = omega_k

    @z_obs   = 0.0  # redshift we're "observing" from

    @c       = GSL::CONST::MKSA::SPEED_OF_LIGHT / 1000.0 # km/s

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


  # the total energy density of all the respective parts.  
  #   Omega > 1 : closed universe
  #   Omega < 1 : open universe
  #   Omega = 1 : 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) = Ho * E(z) / E(0) 
  def hubble_param(z=0.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 = 1 / Ho (Hogg, EQN 3)
  def time_hubble 
    @t_units / (@h * @h_units) 
  end

  # the Hubble distance (today), in Mpc
  #   Dh = c / 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=@z_obs)
    di = GSL::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=@z_obs)
    dc = dist_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=@z_obs)
    dc = dist_co_los(z,z_obs)
    if @omega_k.zero?
      dm = dc
    else 
      dfrac = dist_hubble / Math.sqrt(@omega_k.abs)
      if @omega_k > 0 
        dm = dfrac * Math.sinh( dc / dfrac )
      elsif @omega_k < 0 
        dm = dfrac * 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 = Dm / (1 + z) from  EQN 18
  def dist_ang(z, z_obs=@z_obs)
    dist_co_trans(z,z_obs) * (1 + z_obs) / (1.0 + z)
  end

  # calculate the luminosity distance, in Mpc
  # (section 7, Hogg)
  #   Dl = Dm * (1 + z) from  EQN 21
  def dist_lum(z, z_obs=@z_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=@z_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=@z_obs)
    ti = GSL::Function.alloc LTI_func
    ti.set_params([@omega_m, @omega_l, @omega_k])

    t  = ti.qag([z_obs.to_f, z.to_f]).first
    t *= time_hubble
  end

  # this is the effective scale factor, really the ratio of a_1 / a_obs
  def scale_factor(z, z_obs=@z_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=@z_obs)
    phys_size = dist_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=@z_obs)
    rad = phys_size.to_f / dist_ang(z,z_obs)
  end

  # physical size of the horizon, in Mpc
  def phys_size_horizon(z,z_obs=@z_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=@z_obs)
    # this really is the ratio of comoving distances here.  
    # in general: 
    #   theta = dist_phys / dist_ang (as seen in ang_size())
    # but the physical distance of the horizon is 
    #   d_phys = r_1 * a_1 -> cordinate size times the scale factor at that point
    # therefore
    #   theta = r_1 * a_1 / r_2 * a_1  = r_1 / r_2 -> ratio of cordinate distance
    # aka comoving distances
    ang = dist_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 = Math.sin(ang)
    end

    ang
  end

end

if $0 == __FILE__
  # Example usage
  
  include AstroUnits
  fmt = "%-.5g"

  cosmo = CosmoDist.new

  z = 1100 # 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 = 10
  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 = 1.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=%g:    %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 = cosmo.dist_co_los(z), pc_to_ly(d)/1e3]
    puts "  angular diameter dist: %7i Mpc (%5.2f billion light years)" % 
            [d = cosmo.dist_ang(z),    pc_to_ly(d)/1e3]
    puts "  luminosity distance:   %7i Mpc (%5.2f billion light years)" % 
            [d = cosmo.dist_lum(z),    pc_to_ly(d)/1e3]
    puts 
  end

  puts 
  z, cosmo.z_obs = 1.0,0.0
  report_stats(cosmo,z)
  z = 2.0
  report_stats(cosmo,z)
  cosmo.z_obs = 1.0
  report_stats(cosmo,z)

end





------=_Part_20169_8026061.1146178118120--