#!/usr/bin/perl
#
# Compute the scale factor as a function of time,
#    using starting point at current time and
#    working forward or backward via differential equation
#
# MWR 5/4/2011

use POSIX;

$debug = 0;


$pi = 3.14159;
$G = 1.0;
$rho_0 = 1.0;
$c = 1.0;
$k = 1.0;


# starting values
$a = 1.0;


$start_t = 1.0;
$end_t = 100.0;
$d_t = 0.01;

for ($t = $start_t; $t <= $end_t; $t += $d_t) {

  # first, we compute the value of da_dt
  $x1 = ((8.0*$pi*$G*$rho_0)/3.0)*(1.0/$a);
  $x2 = $k*$c*$c;
  if ($debug > 0) {
    printf " x1 %12.5f  x2 %12.5f \n", $x1, $x2;
  } 
  if ($x1 > $x2) {
    $x3 = $x1 - $x2;
    $da_dt = sqrt( $x3 );
  }
  else {
    $x3 = $x2 - $x1;
    $da_dt = -sqrt( $x3 );
  }

  # next, we compute the next value of a
  $a_new = $a + $da_dt*$d_t;

  # print the current values
  printf " t  %12.5e   a %12.5e   da_dt %12.5e \n",
          $t, $a, $da_dt;

  $a = $a_new;
}








exit 0;