#!/usr/bin/perl
#
# Compute the photon-counting luminosity distance 
#   for a range of values of redshift, given some
#   initial cosmological parameters.
#
# We call the "phot_dist" program in the 
#   ~/z1/snap/zero_point directory to do the work.
#
#
# MWR 5/10/2011

use POSIX;

$debug = 0;

$phot_dist_exe = "~/z1/snap/zero_point/progs/phot_dist";
$lambda = 0.0;
$H0 = 70.0;
$omega_m = 1.00;

$start_z = 0.001;
$end_z = 10.0;
$mult_z = 1.05;


#@H0_list = (30, 50, 70, 90, 110);
#foreach $H0 (@H0_list) {

@omega_lambda_list = (0.40, 0.50, 0.60);
@omega_m_list = (0.20, 0.30, 0.4);
foreach $omega_m (@omega_m_list) {

 foreach $lambda (@omega_lambda_list) {
    
    $outfile = sprintf "lumdist_%05.2f_%05.2f_%06.1f.dat",
       $omega_m, $lambda, $H0;
    open(OUTFILE, ">$outfile") || die("can't open file $outfile");
    
    printf OUTFILE "# %8s  %8s %8s %6s   %6s \n",
         "z", "omega_m", "lambda", "H0", "(m-M)";
    
    for ($z = $start_z; $z < $end_z; $z *= $mult_z) {
    
      $cmd = "$phot_dist_exe $omega_m $lambda $z $H0 ";
      $ret = exec_cmd($cmd);
    
      chomp($ret);
      @words = split(/\s+/, $ret);
      $dist_mod = $words[1];

      # compute the luminosity distance in Mpc
      $dist_pc = 10**(0.2*($dist_mod + 5.0));
      $dist_Mpc = $dist_pc/1.0E6;

      # correct for my error -- calc photon-counting dist properly
      $dist_pc *= sqrt(1.0 + $z);
      $dist_Mpc *= sqrt(1.0 + $z);
      $dist_mod = 5.0*log10($dist_pc) - 5.0;

    
      printf OUTFILE " %10.4e  %8.4f %8.4f %6.1f   %6.3f  %10.4e\n",
         $z, $omega_m, $lambda, $H0, $dist_mod, $dist_Mpc;
    
    }
    close(OUTFILE);

  }
    
}





exit 0;



##############################################################################
# PROCEDURE: exec_cmd
#
# DESCRIPTION: Execute the given shell command line.  If the $debug flag
#              is set, we print to stdout the command line, and
#              also print out the result string.
#
# RETURNS:
#              the result of the command
#
#
sub exec_cmd {

  my($cmd, $ret);

  $cmd = $_[0];

  if ($debug > 0) {
    printf "cmd is ..$cmd.. \n";
  }
  $ret = `$cmd`;
  if ($debug > 0) {
    printf "ret is ..$ret.. \n";
  }

  return($ret);
}


