#!/usr/bin/perl
# Compute in a very simple manner the bias incurred by
#    errors in a parallax measurement.
#
# MWR 5/9/2017

use POSIX;

$debug = 0;
my $mypi = 3.14159265;

srand();


# usage: calc_bias.pl      apparent_pi  sigma_pi  num_star
$USAGE = "calc_bias.pl  measured_pi  sigma_pi  num_star ";

if ($#ARGV != 2) {
  printf "usage: %s \n", $USAGE;
  exit(1);
}

$measured_pi = $ARGV[0];
$sigma_pi = $ARGV[1];
$num_star = $ARGV[2];

if ($debug > 0) {
  printf "  measured_pi    %10.2f \n", $measured_pi;
  printf "  sigma_pi       %10.2f \n", $sigma_pi   ;
  printf "  num_star       %10d \n", $num_star;
}


# this first bit is very simple: just compute ratios of volumes
# loop over a range of ratios of uncertainty sigma to measurement mu
my $mu_mas = 10.0;
my $min_sigma_mas = 1.0;
my $max_sigma_mas = 7.0;
my $d_sigma_mas = 1.0;
for (my $sigma_mas = $min_sigma_mas; $sigma_mas <= $max_sigma_mas; 
           $sigma_mas += $d_sigma_mas) {

  my $best_dist = 1000.0 / $mu_mas;
  my $min_dist = 1000.0 / ($mu_mas + $sigma_mas);
  my $max_dist = 1000.0 / ($mu_mas - $sigma_mas);

  my $inner_volume = (4.0*$mypi/3.0)*($best_dist**3.0) - 
                      (4.0*$mypi/3.0)*($min_dist**3.0);

  my $outer_volume = (4.0*$mypi/3.0)*($max_dist**3.0) - 
                      (4.0*$mypi/3.0)*($best_dist**3.0);

  my $ratio = $outer_volume / $inner_volume;

  printf "# mu %7.2f  sig %7.2f   %6.1f %6.1f %6.1f   %8.3e %8.3e   ratio %6.2f \n",
         $mu_mas, $sigma_mas, $min_dist, $best_dist, $max_dist, 
         $inner_volume, $outer_volume, $ratio;


}


# now comes the more complicated part.  Here, we do a detailed
#    statistical modeling
#
#      - given a set of N stars with measured parallax p +/= sigma
#
#      - loop over each star
#           - compute apparent distance d
#           - choose error from gaussian random distribution
#           - add error to measured p to get true angle ptrue
#           - compute true distance dtrue
#           
#      - in the end, we can print/plot distribution of true distances
# 
my $p = $measured_pi;
my $sigma = $sigma_pi;

  for (my $i = 0; $i < $num_star; $i++) {

    my $d = 1000.0 / $p;

    my $error = gaussian_rand()*$sigma;
    my $true_p = $p + $error;
    my $true_d = 1000.0 / $true_p;

    # now accumulate statistics
    $true_p_array[$i] = $true_p;
    $true_d_array[$i] = $true_d;
  

  }


# print out the table 
if ($debug >= 0) {
  for (my $i = 0; $i < $num_star; $i++) {
    printf "  star %5d  true_p %8.2f  true_d %9.3f \n", 
             $i, $true_p_array[$i], $true_d_array[$i];
  }
}


# compute interquartile mean of all true distances
(my $mean_true_d, my $sigma_true_d) = 
       compute_interquartiles($num_star, @true_d_array);
printf "#  true dist: iq mean %8.2f  iq range %8.2f \n",
          $mean_true_d, $sigma_true_d;



exit 0;


######################################################################
# 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: compute_interquartiles
#
# DESCRIPTION: Given an array of numbers, find the interquartile
#              mean and HALF the interquartile range.  
#
# RETURN:      (inter_quartile_mean, half_inter_quartile_range)
#
#
sub compute_interquartiles {

  my $i, $num, $index25, $index50, $index75;
  my $sum, $sum_terms;
  my $iqm, $iqrange;
  my @array, @sorted_array;

  $iqm = 0.0;
  $iqrange = -1.0;

  $num = $_[0];
  for (my $i = 0; $i < $num; $i++) {
    $array[$i] = $_[1 + $i];
  }

  @sorted_array = sort { $a <=> $b } @array;

  if ($num < 4) {
    printf "compute_interquartiles: warning: num %d is small! \n", $num;
  }
  $index25 = floor($num/4.0);
  $index50 = floor($num/2.0);
  $index75 = floor(3.0*$num/4.0);
  if ($debug > 1) {
    printf "   num %4d  index25 %3d  index50 %3d  index75 %3d \n",
          $num, $index25, $index50, $index75;
  }

  # compute interquartile mean
  $sum_terms = 0;
  $sum = 0;
  for (my $i = $index25; $i <= $index75; $i++) {
    if ($debug > 1) {
      printf "  i %5d  sorted_array %7.3f  \n", $i, $sorted_array[$i];
    }
    $sum += $sorted_array[$i];
    $sum_terms++;
  }
  if ($sum_terms > 0) {
    $iqm = $sum / $sum_terms;
    $iqrange = ($sorted_array[$index75] - $sorted_array[$index25]) / 2.0;
  }
  else {
    printf "compute_interquartiles: sum_terms == 0 ?? bogus results \n";
  }

  if ($debug > 1) {
    printf "  iqm %7.3f  iqrange %7.3f \n", $iqm, $iqrange;
  }
  
  return($iqm, $iqrange);
}