#!/usr/bin/perl
#
# generate a list of pseudo-observed magnitudes,
#    using an input list of real star magnitudes
#    and some model for the color terms
#    and some random scatter
#
# MWR 9/7/2011

$debug = 1;

# read in the real star data
$landolt_file = "./pg1633_landolt.dat";
$num_stars = 0;
open(LANDOLT_FILE, "$landolt_file") || die("can't open $landolt_file");
while (<LANDOLT_FILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  @words = split(/\s+/, $line);
  if ($#words != 5) {
    printf " bad input line in %s \n", $landolt_file;
    printf " %s \n", $line;
    exit(1);
  }
  $landolt_name_array[$num_stars] = $words[0];
  $landolt_umag_array[$num_stars] = $words[1];
  $landolt_bmag_array[$num_stars] = $words[2];
  $landolt_vmag_array[$num_stars] = $words[3];
  $landolt_rmag_array[$num_stars] = $words[4];
  $landolt_imag_array[$num_stars] = $words[5];

  $num_stars++;
}
close(LANDOLT_FILE);

if ($debug > 1) {
  for ($i = 0; $i < $num_stars; $i++) {
    printf " star %3d has name %10s and V %7.4f \n", 
          $i, $landolt_name_array[$i], $landolt_vmag_array[$i];
  }
}


# now, compute the "observed" magnitude of each star in "b" and "v"
$b_zero = 2.58;
$v_zero = 1.35;
$b_term = 0.24;
$v_term = 0.0;
$b_scatter = 0.05;
$v_scatter = 0.03;

for ($i = 0; $i < $num_stars; $i++) {
  $color = $landolt_bmag_array[$i] - $landolt_vmag_array[$i];

  $b = $landolt_bmag_array[$i] 
          + $b_zero 
          + $b_term*($color)
          + $b_scatter*gaussian_rand();

  $v = $landolt_vmag_array[$i] 
          + $v_zero 
          + $v_term*($color)
          + $v_scatter*gaussian_rand();

  if ($debug > 0) {
     printf " %9s  %7.4f %7.4f %6.4f   %7.4f %7.4f %6.4f  \n", 
         $landolt_name_array[$i], 
         $landolt_bmag_array[$i],
         $b, $b_scatter,
         $landolt_vmag_array[$i],
         $v, $v_scatter;
  }
}


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);
}




######################################################################
# PROCEDURE: gaussian_rand
#
# DESCRIPTION: Returns a random number drawn from a gaussian (normal)
#              distribution with mean 0.0 and standard deviation 1.0.
#
#              Stolen from 
#
#                 http://www.unix.com.ua/orelly/perl/cookbook/ch02_11.htm
#
#
sub gaussian_rand {
    my ($u1, $u2);  # uniformly distributed random numbers
    my $w;          # variance, then a weight
    my ($g1, $g2);  # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w = $u1*$u1 + $u2*$u2;
    } while ( $w >= 1 );

    $w = sqrt( (-2 * log($w))  / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;
    # could return both if wanted
    #   return wantarray ? ($g1, $g2) : $g1;
    # but I need just one
    return($g1);
}

