#!/usr/bin/perl
#
# A script to calculate the evolution of the scale factor of the
#   universe, given initial parameters Omega, H-nought, and Lambda.
#   Uses equations from Felten and Isaacman, Rev. Mod. Phys. 58, 689 (1986).
#
# MWR 2/27/2003
#
# This version may fix mistakes in the earlier version.  
#   Specifically, I think what I _ought_ to be using as the
#   output time value is the internal time variable, tau,
#   divided by the Hubble constant hnought.
# MWR 3/23/2004
#
# Fixed very minor bug in the forward-stepping portion of the program
#   (and in the backward-stepping portion of the program, too)
#   Was printing current time tau and _next_ time R_new 
#   Now we print current time tau and _current_ time R
# Also, we now expect command-line options
#
#      perl hubble_x.pl   omega_matter  omega_lambda
# 
# and don't bother with the H0 value.
# 
# Changed also so that no output is printed if R(t) < 0.
#
# MWR 3/19/2008
#
# Added new column 11 of output: redshift z 
#     for times before present, actual redshift
#     for times after present, a bogus value of -1
# 
# MWR 3/21/2008
#
# Bug fix: the redshift z wasn't being assigned properly
#     the array was filled correctly, but the age used as index to array
#           wasn't same as age used to look up in the array.
#     new function "find_tau_for_BB()" runs before main loop
#           in order to compute the "tau" value when BB occurs
#           (will be negative number).  We can then use that
#           value of tau_bb, plus main loop's value of tau,
#           to compute the proper age for looking up redshift z in array.
#  
# MWR 5/5/2011




use POSIX;

$debug = 0;

# when to start and end the integration, and the timestep, all in units
#   of the Hubble time (which is very roughly ten billion years)
$start_tau = -5.0;
$end_tau = 5.0;
$dtau = 0.002;

# we keep two global arrays: one to hold the scaled time (tau) at each
#   step of the simulation, and one to hold the scale value (R) at each step.
#   Each array will contain this many elements:
#
#           [ (end_tau - start_tau) / dtau ]  -  1
#
#   (The "minus 1" at the end occurs because we don't include the starting
#   and ending points in the arrays)
#
$tau_array[0] = 0;
$R_array[0] = 0;


#
# first, get input values from command line
#
if ($#ARGV != 1) {
  die("error: usage   perl hubble_x_new.pl  omega_matter omega_lambda \n");
}
$omega_matter = $ARGV[0];
$omega_lambda = $ARGV[1];
$omega_tot = $omega_matter + $omega_lambda;
if ($debug > 0) {
  printf " omega_matter %6.3f  omega_lambda %6.3f  omega_tot  %6.3f \n",
         $omega_matter, $omega_lambda, $omega_tot;
}

# we translate the input "omega_lambda" into the parameter used by
#   Felten and Isaacman, which is exactly three times larger
$lambda = $omega_lambda*3.0;
$omega = $omega_matter;
# this doesn't matter --- will cancel out in end
$hnought = 1.0;

if ($debug > 0) {
  printf " Omega ..$omega..  H-nought ..$hnought..  Lambda ..$lambda.. \n";
}


# we need to know the value of "tau" when the BB occurs in order
#   to look up values for redshift at different tau.
#   Note that if find_tau_for_BB() fails, it will return 0.0
$tau_for_bb = find_tau_for_BB();

# generate an array of redshift-age pairs; that is, for a large range of
#   values of redshift, calculate the age in units of tau.
#   We can later start with a value of tau and use this array to look
#   up the corresponding value of redshift z
#  
# we adopt H0 = 75 as our standard value here 
$H0 = 75;
$hubble_time_in_gyr = compute_age_from_bb_to_z($H0, $omega_matter, $omega_lambda, 0.0);
$z_start = 0.0;
$z_end = 10.0;
$dz = 0.01;
$num_z_array = 0;
for ($z = $z_start; $z <= $z_end; $z += $dz) {
  $age_from_bb = compute_age_from_bb_to_z($H0, $omega_matter, $omega_lambda, $z);

  $z_age_array[$num_z_array] = $age_from_bb / $hubble_time_in_gyr;
  $z_z_array[$num_z_array] = $z;
  if ($debug > 0) {
    printf "  num_z_array %6d  z %8.3f  age %8.3f \n", $num_z_array, $z, 
                  $z_age_array[$num_z_array];
  }
  $num_z_array++;
}





# 
# perform the integration: first, we go backwards in time from
#   the current day (tau = 0) to the starting point.  Then, we
#   go forward in time from the current day to the ending point.
#   We place our results into a pair of arrays: one for the
#   _time_ (tau) at each step, and one for the _scale_ (R) at
#   each step.  We'll print out the values later.
#

# go backwards in time
if (1 == 1) {
$R = 1.0;
$index = int((0.0 - $start_tau)/$dtau);
for ($tau = 0.0; $tau > $start_tau; $tau -= $dtau) {

  # find the redshift for this age
  #    MWR: fix mistake here, wasn't using a proper age to look up
  #         the redshift.  5/5/2011
  # old
  #    $my_age = 1.0 + $tau;
  # new
  $my_age = ($tau - $tau_for_bb) / (0.0 - $tau_for_bb);
  if ($debug > 0) {
    printf " top of main loop, tau is %8.3f   age is %8.3f \n", $tau, $my_age;
  }
  $z = find_z_for_age($num_z_array, $my_age);
  if ($debug > 0) {
    printf "  z is %8.3f \n", $z;
  }
  

  # fill in the current time and size elements
  $tau_array[$index] = $tau;
  $R_array[$index] = $R;

  # calculate the second derivative, d2R/dtau2
  $second_deriv = ((0.0 - ($omega/($R*$R))) + 
                         (2.0/3.0)*($R*$lambda)/($hnought*$hnought));
  $second_deriv /= 2.0;

  # calculate the first derivative, dR/dtau
  $x = big_term($R, $omega, $hnought, $lambda);
  if ($tau == 0.0) {
    $first_deriv = sqrt($x);
  }
  else {
    $first_deriv -= $second_deriv*$dtau;
  }

  # now, we use a second-order method to integrate from current
  #   value of R to the previous one
  $new_R = $R - ($first_deriv*$dtau);

  if ($debug > -1) {
    if ($R >= 0.0) {
      printf "%8.4f  %12.5e  %12.5e = %12.5e  %12.5e = %12.5e  %9.6f  %8.4f  %10.5f \n", 
           $tau, $x, $first_deriv, $first_deriv*$dtau, 
# old version printed new R
#           $second_deriv, 0.5*$second_deriv*$dtau*$dtau, $new_R, $tau/$hnought;
# but we should be printing the current R
           $second_deriv, 0.5*$second_deriv*$dtau*$dtau, $R, $tau/$hnought, $z;
    }
  }

  # prepare for next time through loop
  $R = $new_R;
  $index--;
}
}



# go forwards in time
$R = 1.0;
$index = int(($end_tau - 0)/$dtau);
for ($tau = 0.0; $tau < $end_tau; $tau += $dtau) {
  
  # redshift doesn't make sense now, just set to -1.0
  $z = -1.0;

  # fill in the current time and size elements
  $tau_array[$index] = $tau;
  $R_array[$index] = $R;

  # calculate the second derivative, d2R/dtau2
  $second_deriv = ((0.0 - ($omega/($R*$R))) + 
                         (2.0/3.0)*($R*$lambda)/($hnought*$hnought));
 # $second_deriv /= 2.0*sqrt($x);
  $second_deriv /= 2.0;

  # calculate the first derivative, dR/dtau
  $x = big_term($R, $omega, $hnought, $lambda);
if (0 == 1) {
  if ($x > 0.0) {
    $first_deriv = sqrt($x);
  }
  else {
    $first_deriv += $second_deriv*$dtau;
  }
}

  if ($tau == 0.0) {
    $first_deriv = sqrt($x);
  }
  else {
    $first_deriv += $second_deriv*$dtau;
  }

  # now, we use a second-order method to integrate from current
  #   value of R to the previous one
#  $new_R = $R + ($first_deriv*$dtau) + 0.5*($second_deriv)*($dtau*$dtau);
  $new_R = $R + ($first_deriv*$dtau);
 
  if ($debug > -1) {
    if ($R >= 0.0) {
      printf "%8.4f  %12.5e  %12.5e = %12.5e  %12.5e = %12.5e  %9.6f %8.4f  %10.5f \n", 
           $tau, $x, $first_deriv, $first_deriv*$dtau, 
# this was old version
#           $second_deriv, 0.5*$second_deriv*$dtau*$dtau, $new_R, $tau/$hnought;
# my attempt at fix -- print current tau and current R, not new R
           $second_deriv, 0.5*$second_deriv*$dtau*$dtau, $R, $tau/$hnought, $z;
    }
  }

  # prepare for next time through loop
  $R = $new_R;
  $index++;
}


exit 0;


###########################################################################
# PROCEDURE: big_term
# 
# DESCRIPTION: calculate a mathematical expression which appears in both
#              the first and second derivatives of scale factor with
#              respect to tau.
#
#              Expects four arguments:
#                     R          current value of scale factor
#                     omega      fixed value of overall density factor
#                     hnought    fixed value of Hubble Constant
#                     lambda     fixed value of cosmological constant
#
# RETURNS:
#              a numerical value
#
sub big_term {
   my($R, $omega, $hnought, $lambda);

   $R = $_[0];
   $omega = $_[1];
   $hnought = $_[2];
   $lambda = $_[3];

   $x = ($omega*(1.0 - $R))/($R*$R*$R);
   if ($debug > 0) {
     printf "big_term: x = %12.5e ", $x;
   }
   $x += (1.0/($R*$R));
   $x += ((1.0/3.0)*($lambda/($hnought*$hnought))*(($R*$R - 1.0)/($R*$R)));
   $x *= $R*$R; 
   if ($debug > 0) {
     printf "   final x = %12.5e \n", $x;
   }

   return($x);

}



###############################################################
# Do the computation: find the time from BB to redshift "z"
#
# Input args:
#           H0            =   Hubble const   km/s/Mpc
#           omega_matter  =   matter's contribution to Omega
#           omega_lambda  =   lambda's contribution to Omega
#           z             =   redshift
# 
# We set omega_radiation to a very small value, following Ned
#
# Output:   
#           age            in Gyr
#                 
#
#
sub compute_age_from_bb_to_z {

  my($h, $WR, $WK, $a, $az, $age, $i, $adot);
  my($zage, $n);
  my($Tyr, $age_Gyr);


  $H0 = $_[0];
  $WM = $_[1];
  $WV = $_[2];
  $z  = $_[3];

  $h = $H0/100.0;

  # convert 1/H into Gyr
  $Tyr = 977.8;

  # this is fixed, I don't care about radiation
  # includes 3 massless neutrino species, T0 = 2.72528
  $WR = 4.165E-5/($h*$h);	

  # number of points in the integral
  $n = 1000;

  # prepare to do the calculation
  $h = $H0/100;
  $WK = 1-$WM-$WR-$WV;
  $az = 1.0/(1+1.0*$z);
  $age = 0;
  for ($i = 0; $i < $n; $i++) {
    $a = $az*($i+0.5)/$n;
    $adot = sqrt($WK+($WM/$a)+($WR/($a*$a))+($WV*$a*$a));
    $age = $age + 1/$adot;
  };
  $zage = $az*$age/$n;

  $zage = $zage*pow(10.0, $dzage);
  $age_Gyr = ($Tyr/$H0)*$zage;

  return($age_Gyr);

}

###############################################################
# Do the computation: find the time from redshift "z" to the present 
#
# Input args:
#           H0            =   Hubble const   km/s/Mpc
#           omega_matter  =   matter's contribution to Omega
#           omega_lambda  =   lambda's contribution to Omega
#           z             =   redshift
# 
# We set omega_radiation to a very small value, following Ned
#
# Output:   
#           age            in Gyr
#                 
#
#
sub compute_age_from_z_to_now {

  my($h, $WR, $WK, $a, $az, $age, $i, $adot);
  my($zage, $n);
  my($Tyr, $age_Gyr);
  my($DTT);


  $H0 = $_[0];
  $WM = $_[1];
  $WV = $_[2];
  $z  = $_[3];

  $h = $H0/100.0;

  # convert 1/H into Gyr
  $Tyr = 977.8;

  # this is fixed, I don't care about radiation
  # includes 3 massless neutrino species, T0 = 2.72528
  $WR = 4.165E-5/($h*$h);	

  # number of points in the integral
  $n = 1000;

  # prepare to do the calculation
  $DTT = 0.0;

  #  do integral over a=1/(1+z) from az to 1 in n steps, midpoint rule
  $az = 1.0/(1+1.0*$z);
  for ($i = 0; $i < $n; $i++) {
    $a = $az+(1.0-$az)*($i+0.5)/$n;
    $adot = sqrt($WK+($WM/$a)+($WR/($a*$a))+($WV*$a*$a));
    $DTT = $DTT + 1.0/$adot;
  };
  $DTT = (1.0-$az)*$DTT/$n;
  $age = $DTT+$zage;
  $age_Gyr = $age*($Tyr/$H0);

  return($age_Gyr);

}




###############################################################
# PROCEDURE: find_z_for_age
#
# Given an age, in Gyr, find the redshift z which corresponds
#   to it.  We have a pair of arrays, one with z and one with
#   age; we use a binary search in age to find the right element.
#
# Input args:
#              num_z_array  = number of elements in the two arrays
#              age          = age in Gyr
# 
# Output:   
#              redshift z           
#                 
#
#
sub find_z_for_age {

  my($num, $age, $z);
  my($top, $bot, $mid);
  my($index);

  $num = $_[0];
  $age = $_[1];

  $top = $num - 1;
  $bot = 0;
  
  while ($top - $bot > 1) {

    $mid = floor(($top + $bot)/2);

    if ($debug > 0) {
      printf " want %8.5f   top %4d %8.5f  mid %4d %8.5f  bot %4d %8.5f \n",
                    $age, 
                    $top, $z_age_array[$top],
                    $mid, $z_age_array[$mid],
                    $bot, $z_age_array[$bot],
    }

    $this_age = $z_age_array[$mid];
    if ($this_age < $age) {
      $top = $mid;
    }
    elsif ($this_age >= $age) {
      $bot = $mid;
    }

  }

  # at this point, there are several possibilities
  #
  # the desired age lies above the "bot" element's limit.
  #     Result: return the "bot" element
  if ($age >= $z_age_array[$bot]) {
    $z = $z_z_array[$bot];
    if ($debug > 0) {
      printf "         interp A \n";
    }
  }
  # the desired age lies between "bot" and "mid" elements
  #     Result: linearly interpolate between "bot" and "mid"
  elsif ($age >= $z_age_array[$mid]) {
    $d_age = $z_age_array[$bot] - $z_age_array[$mid];
    $frac = ($age - $z_age_array[$mid]) / $d_age;
    $d_z = ($z_z_array[$bot] - $z_z_array[$mid]);
    $z = $z_z_array[$mid] + $frac*$d_z;
    if ($debug > 0) {
      printf "         interp B \n";
    }
  }
  # the desired age lies between "mid" and "top" elements
  #     Result: linearly interpolate between "bot" and "mid"
  elsif ($age >= $z_age_array[$top]) {
    $d_age = $z_age_array[$mid] - $z_age_array[$top];
    $frac = ($age - $z_age_array[$top]) / $d_age;
    $d_z = $z_z_array[$mid] - $z_z_array[$top];
    $z = $z_z_array[$top] + $frac*$d_z;
    if ($debug > 0) {
      printf "         interp C \n";
    }
  }
  else {
  # the desired age lies below the "top" element's limit
  #    Result: return the "top" element
    $z = $z_z_array[$top];
    if ($debug > 0) {
      printf "         interp D \n";
    }
  }

  return($z);

}





###################################################################
# PROCEDURE: find_tau_for_BB
#
# DESCRIPTION:  This is a copy of the main loop we perform
#               in order to compute scale factor as a function
#               of time before the present.
#
#               We go through the whole loop until we reach the
#               point at which the scale factor becomes zero.
#               We then note the value of tau, and return it.
#
#               If the scale factor reaches SILLY_TAU_VALUE,
#               we halt and return a value of tau = 0.
#
# RETURNS:
#               The value of tau at which the Big Bang occurs
#
sub find_tau_for_BB {

  my($R, $new_R, $index);
  my($tau, $start_tau);
  my($my_age);
  my($z);
  my($second_deriv, $x, $first_deriv);

  # if we go back this far without scale factor becoming zero,
  #    we give up, print error message, return 0.0
  $SILLY_TAU_VALUE = -20.0;
  $start_tau = $SILLY_TAU_VALUE;
  
  $R = 1.0;
  $index = int((0.0 - $start_tau)/$dtau);
  for ($tau = 0.0; $tau > $start_tau; $tau -= $dtau) {
  
    # find the redshift for this age
    $my_age = 1.0 + $tau;
    if ($debug > 0) {
      printf "find_tau: top of main loop, tau is %8.3f   age is %8.3f \n", $tau, $my_age;
    }
    $z = find_z_for_age($num_z_array, $my_age);
    if ($debug > 0) {
      printf "find_tau:  z is %8.3f \n", $z;
    }
    
    # calculate the second derivative, d2R/dtau2
    $second_deriv = ((0.0 - ($omega/($R*$R))) + 
                           (2.0/3.0)*($R*$lambda)/($hnought*$hnought));
    $second_deriv /= 2.0;
  
    # calculate the first derivative, dR/dtau
    $x = big_term($R, $omega, $hnought, $lambda);
    if ($tau == 0.0) {
      $first_deriv = sqrt($x);
    }
    else {
      $first_deriv -= $second_deriv*$dtau;
    }
  
    # now, we use a second-order method to integrate from current
    #   value of R to the previous one
    $new_R = $R - ($first_deriv*$dtau);
  
	 if ($debug > 0) {
      if ($R >= 0.0) {
        printf "find_tau %8.4f  %12.5e  %12.5e = %12.5e  %12.5e = %12.5e  %9.6f  %8.4f  %10.5f \n", 
             $tau, $x, $first_deriv, $first_deriv*$dtau, 
  # old version printed new R
  #           $second_deriv, 0.5*$second_deriv*$dtau*$dtau, $new_R, $tau/$hnought;
  # but we should be printing the current R
             $second_deriv, 0.5*$second_deriv*$dtau*$dtau, $R, $tau/$hnought, $z;
      }
    }
    
    if ($R <= 0.0) {
      if ($debug > 0) {
        printf "find_tau:  R <= 0 at tau = %6.2f \n", $tau;
      }
      return($tau);
    }
  
    # prepare for next time through loop
    $R = $new_R;
    $index--;
  }
  
  # if we get here, we failed to reach scale factor = 0.
  printf STDERR "find_tau_for_BB: reached tau %6.1f without R=0 \n", 
                $SILLY_TAU_VALUE;
  return(0.0);
  
}