#!/usr/bin/perl
#
# Calculate the scale factor, R, versus time, t, for simple
# models of the universe.  Print the results in the form
#
#         time      closed   critical  open  empty
#
# MWR 3/16/2002

use POSIX;

$omega_empty = 0.001;
$omega_open = 0.5;
$omega_closed = 2.0;
$h0 = 1.0;

$t_start = 0;
$t_end = 10.0;
$t_h = 1.0;
$dt = ($t_end - $t_start)/1000.0;


for ($t = $t_start; $t <= $t_end; $t += $dt) {

  # this is the critical density universe
  $crit_r = pow(1.5, (2.0/3.0)) * pow(($t/$t_h), (2.0/3.0));

  $closed_r = calc_closed_r($omega_closed, $h0, $t);

  $open_r = calc_open_r($omega_open, $h0, $t);

  $empty_r = calc_open_r($omega_empty, $h0, $t);

  printf "%8.3f  %8.3f  %8.3f  %8.3f  %8.3f \n", 
           $t, $crit_r, $closed_r, $open_r, $empty_r;
}


exit 0;


##########################################################################
# Given a particular time 't', calculate the scale factor of the
# CLOSED universe at that time.
#
# We use parametric equations for both R and t as a function of
# another variable, 'x', so this is a little bit indirect.
# It's also inefficient, but who cares?
#
# Usage:   calc_closed_r    omega h0 t
#
# where
#              omega          is the value of density
#              h0             is the hubble constant at current time
#              desired_t      is the desired time
#
# Returns
#       scale factor at the given time

sub calc_closed_r {
  my($omega, $h0, $desired_t);
  my($r, $t);
  my($x, $dx, $i, $j);
  my($top, $bot);

  $omega = $_[0];
  $h0 = $_[1];
  $desired_t = $_[2];
  
  $dx = 0.01;

  for ($x = 0; $x < 10; $x += $dx) {
 
    # determine the time for this value of 'x'
    $top = $omega*($x - sin($x));
    $bot = 2*$h0*pow(($omega - 1.0), 1.5);
    $t = $top/$bot;

    # determine the scale factor for this value of 'x'
    $top = $omega*(1.0 - cos($x));
    $bot = 2*($omega - 1.0);
    $r = $top/$bot;

    if ($t > $desired_t) {
      return($r);
    }
    # printf "%8.3f %8.3f \n", $t, $r;
  }

  # we should never get here
  printf "ERROR ran out of range of calc_closed_r\n";
  return(0);

}


##########################################################################
# Given a particular time 't', calculate the scale factor of the
# OPEN universe at that time.
#
# We use parametric equations for both R and t as a function of
# another variable, 'x', so this is a little bit indirect.
# It's also inefficient, but who cares?
#
# Usage:   calc_open_r    omega h0 t
#
# where
#              omega          is the value of density
#              h0             is the hubble constant at current time
#              desired_t      is the desired time
#
# Returns
#       scale factor at the given time

sub calc_open_r {
  my($omega, $h0, $desired_t);
  my($r, $t);
  my($x, $dx, $i, $j);
  my($top, $bot);

  $omega = $_[0];
  $h0 = $_[1];
  $desired_t = $_[2];
  
  $dx = 0.001;

  for ($x = 0; $x < 50; $x += $dx) {
 
    # determine the time for this value of 'x'
    $top = $omega*(sinh($x) - $x);
    $bot = 2*$h0*pow((1.0 - $omega), 1.5);
    $t = $top/$bot;

    # determine the scale factor for this value of 'x'
    $top = $omega*(cosh($x) - 1.0);
    $bot = 2*(1.0 - $omega);
    $r = $top/$bot;

    if ($t > $desired_t) {
      return($r);
    }
    # printf "%8.3f %8.3f \n", $t, $r;
  }

  # we should never get here
  printf "ERROR ran out of range of calc_open_r\n";
  return(0);

}
