#!/usr/bin/perl
#
# Generate synthetic populations of planetary nebulae,
#    as they would be observed in some distant galaxy.
#
# Print apparent magnitudes to stdout.
#
# MWR 4/23/2011

use POSIX;

$debug = 0;

if (1 == 1) {
  $t = time;
  srand($t);
}

# Set A: no incompleteness, dist_mod = 26.7
#$dist_mod = 26.7;
# Set B: incompleteness, dist_mod = 34.2
$dist_mod = 34.2;
# how many objects in the sample?
$num_obj = 1000;
# should we add incompleteness?
$do_incomplete = 1;
$lim_mag = 27.20;
$alpha = 1.8;

# most probable absolute mag
$m_peak = -7.1;
# width of distribution in magnitude
$mag_sigma = 1.1;
# how far away from peak do we go in each direction?
$mag_span = 4.0;
# at what intervals of magnitude do we sample the population?
$mag_interval = 0.05;
# arbitrary scaling constant
$K = 100.0;

# first, we build a population of objects which have the proper
#    distribution; we'll draw randomly from them later
$num_pop = 0.0;
for ($delta_mag = 0.0 - $mag_span; $delta_mag < $mag_span; $delta_mag += $mag_interval) {

  # the absolute magnitude of objects in this little bin
  $mag = $m_peak + $delta_mag;
  
  # the number of objects in this little bin (which we will turn into
  #     a _relative_ number later)
  $num = $K*exp(- (($mag - $m_peak)**2.0) / (2.0*$mag_sigma*$mag_sigma) );
  $num = $num + 1.0;
  if ($debug > 0) {
    printf "   delta_mag %6.2f  mag %6.2f  num %12.3f \n",
        $delta_mag, $mag, $num;
  }
  
  # fill the cumulative array with "num" entries, all of which have
  #    the current magnitude
  for ($j = 0; $j < $num; $j++) {
    $mag_array[$num_pop + $j] = $mag;
  }

  $num_pop += $num;

}
if ($debug > 0) {
  printf "# population has %d objects \n", $num_pop;
  for ($i = 0; $i < $num_pop; $i++) {
    printf "# %6d  %6.3f \n", $i, $mag_array[$i];
  }
}



for ($i = 0; $i < $num_obj; $i++) {

  # select a random number uniformly distributed between 0 and 1
  $r = rand();
  if ($debug > 1) {
    printf "r is %12.5f \n", $r;
  }

  # translate it into a number between 0 and num_pop
  $r_int = floor($r*$num_pop);

  # use it to pick out one of the objects
  $obj_mag = $mag_array[$r_int];

  $obj_mag += $dist_mod;

  # now, determine if we should include this or not,
  #   due to incompleteness
  if ($do_incomplete == 1) {
    if (include_item($obj_mag, $lim_mag, $alpha) == 0) {
      $i--;
      next;
    }
  }

  printf " r  %12.5f   r_int %12d   mag %6.3f \n", 
          $r, $r_int, $obj_mag;
  
}

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


##########################################################################
# PROCEDURE: include_item
#
# DESCRIPTION: Given the magnitude of an object, and a "limiting" magnitude,
#              decide if the object should be detected or not.
#
#              See Fleming et al., AJ 109, 1044 (1995)
#
#              Arguments:   
#                     mag          magnitude of item
#                     lim_mag      mag at which 50% objects detected
#                     alpha        param used for sharpness of decline
#
#
# RETURNS:     
#                   1 if the object will be detected
#                   0 if the object will NOT be detected
#
sub include_item {

  my($mag, $lim_mag, $alpha);
  my($f, $r);
  my($top, $bot);

  $mag = $_[0];
  $lim_mag = $_[1];
  $alpha = $_[2];

  # generate the likelihood that we detect object (range 1.0 to 0.0)
  $top = $alpha*($mag - $lim_mag);
  $bot = sqrt(   1.0 + $alpha*$alpha*( ($mag - $lim_mag)**2.0 )  );
  $f = 0.5 * (  1.0 - ($top/$bot) );

  # now generate a random number between 0.0 and 1.0
  $r = rand();

  # we detect the object if f > r
  if ($f > $r) {
    return(1);
  }
  else {
    return(0);
  }
  
}

