#!/usr/bin/perl
#
# Compute the number of galaxies predicted by a real luminosity function
#
# MWR 4/6/2011

use POSIX;

$debug = 0;

# absolute magnitude of the Sun in V-ish
$magsun = 4.9;
# these are parameters from 
#    http://adsabs.harvard.edu/abs/2009MNRAS.399.1106M
$magstar = -19.59;
$phi = 0.0140;
$alpha = -1.10;

$lumstar = 10**(-0.4*($magstar - $magsun));
if ($debug > 0) {
  printf " magstar %6.3f  magsun %6.3f  diff %6.3f lumstar %8.3e \n", 
       $magstar, $magsun, ($magstar - $magsun), $lumstar;
}

$start_mag = -5;
$end_mag = -25;
$d_mag = -0.1;

for ($mag = $start_mag; $mag >= $end_mag; $mag += $d_mag) {

  if ($debug > 0) {
     printf " mag is %6.3f \n", $mag;
  }

  # compute luminosity in solar units
  $lum = 10**(-0.4*($mag - $magsun));

  # compute interval in luminosity
  $next_mag = $mag + $d_mag;
  $next_lum = 10**(-0.4*($next_mag - $magsun));
  $d_lum = $next_lum - $lum;
  if ($d_lum < 0) {
    $d_lum = 0.0 - $d_lum;
  }
  
  # compute number of galaxies within this mag interval
  $x1 = ($phi / $lumstar);
  $x2 = (($lum/$lumstar)**$alpha);
  $x3 = exp(-$lum/$lumstar);
  $x4 = $d_lum;
  if ($debug > 0) {
    printf " %9.4e  %9.4e  ( %8.3e %8.3e %8.3e)   %9.4e  \n",
             $x1, $x2, $lum, $lumstar, $x3, $x4;
  }

  $num = ($phi/$lumstar) * (($lum/$lumstar)**$alpha) * 
                  (exp(-$lum/$lumstar)) * $d_lum;

  printf "  %6.3f  %9.4e    %9.4e \n", $mag, $lum, $num;

}


exit 0;
