#!/usr/bin/perl
#
# Compute numerically the integral of a function 
#
#              dr 
#     ---------------------
#         1         1   
#       ----   -  ----- 
#         r       r_max
#
# from r=0 to r=r_max
#
# MWR 4/7/2011

$debug = 0;
$tiny = 1e-6;


for ($r_max = 1.0; $r_max < 100.0; $r_max += 1.0) {
#for ($r_max = 3.0; $r_max <= 3.0; $r_max += 1.0) {

  
  $d_r = 0.001;
  $start_r = $d_r;
  $end_r = $r_max;
  
  $sum = 0;
  for ($r = $start_r; $r < $end_r; $r += $d_r) {

    if (($end_r - $r) < $tiny) {
      last;
    }
  
    $func = eval_func($r, $r_max);
    $this_piece = $func*$d_r;
    $sum += $this_piece;
  
    if ($debug > 0) {
      printf " %9.4e  %9.4e  %9.4e  %9.4e  \n", 
         $r, $func, $this_piece, $sum;
    }
  
  
  }

  printf " %9.4e  %9.4e \n", $r_max, $sum;

}
  

##########################################
# PROCEDURE: evaluate the function we are
#            going to integrate, given 
#            arguments "r" and "r_max".
#
#
sub eval_func {

  my($r, $r_max);
  my($retval);

  $r = $_[0];
  $r_max = $_[1];

  $retval = (  (1.0/$r) - (1.0/$r_max) );
  $retval = 1.0 / $retval;

  return($retval);
}
