#!/usr/bin/perl
#
# Create a distribution of items following some particular power law.
#
# In fact, create a data file with several distributions in different
# columns, for easy comparison.
#
# MWR 6/15/2013

use POSIX;
$debug = 0;

# total number of objects
$total = 1000;
# power law exponent for group A
$a_exp = -1.0;
# power law exponent for group B
$b_exp = 0.0;
# power law exponent for group C
$c_exp = 0.5;


# min mass (solar masses)
$min_mass = 0.05;
# max mass (solar masses)
$max_mass = 10;


# we compute distributions for three power law exponents
$a_exp = -1.0;
$b_exp = 0.0;
$c_exp = 0.5;


# compute the number of stars within each mass bin over the entire range
#   (we'll normalize later)
$d_m = 0.5;
$d_log_m = 0.03;
$m = $min_mass;
$num_bin = 0;
while ($m < $max_mass) {

  $bin_bot = $m;
if (1 == 0) {
  $bin_bot_log = log10($bin_bot);
  $bin_top_log = $bin_bot_log + $d_log_m;
  $bin_top = (10.0**$bin_top_log);
} else {
  $bin_top = $bin_bot + $d_m;
}
  $bin_width = $bin_top - $bin_bot;

  $bin_bot_array[$num_bin] = $bin_bot;
  $bin_top_array[$num_bin] = $bin_top;
  $bin_mid_array[$num_bin] = ($bin_bot + $bin_top)*0.5;

  if ($debug > 1) {
    printf "  m %8.4f   bot %8.4f  top %8.4f  width %8.4f \n", 
         $m, $bin_bot, $bin_top, $bin_width;
  }

  $num_a = ($m**$a_exp)*$bin_width;
  $num_b = ($m**$b_exp)*$bin_width;
  $num_c = ($m**$c_exp)*$bin_width;

  $num_a_array[$num_bin] = $num_a;
  $num_b_array[$num_bin] = $num_b;
  $num_c_array[$num_bin] = $num_c;
  
  if ($debug > 0) {
    printf "#  %4d %8.4f   %8.4f   a %8.4e  b %8.4e  c %8.4e \n", 
         $num_bin, $bin_bot, $bin_top, $num_a, $num_b, $num_c;
  }


  $m = $bin_top;
  $num_bin++;
}


# now normalize the values in the "num_a_array[]" arrays so that
#   the cover the range 0.0 to 1.0
$sum_a = 0;
$sum_b = 0;
$sum_c = 0;
for ($i = 0; $i < $num_bin; $i++) {
  $sum_a += $num_a_array[$i];
  $sum_b += $num_b_array[$i];
  $sum_c += $num_c_array[$i];
}
if ($debug > 0) {
  printf "# sum_a  %10.4f  sum_b %10.4f  sum_c %10.4f \n",
    $sum_a, $sum_b, $sum_c;
}
for ($i = 0; $i < $num_bin; $i++) {
  $num_a_array[$i] *= ($total/$sum_a);
  $num_b_array[$i] *= ($total/$sum_b);
  $num_c_array[$i] *= ($total/$sum_c);
}
  

# print out the results
  for ($i = 0; $i < $num_bin; $i++) {
    printf "  %8.4f  %8.2f %8.2f %8.2f \n",
       $bin_mid_array[$i],
       $num_a_array[$i], $num_b_array[$i], $num_c_array[$i];
  }





exit 0;