#!/usr/bin/perl
#
# Make noisy data with a 1-D gaussian
#    Just one gaussian, at a given location

use POSIX;

$amp = 1.0;
$fwhm = 1.0;
$sigma = $fwhm / 2.35;
$x0 = 12.3333333333;
$noise_stdev = 0.1;



$start_t = 0;
$end_t = 20;
$dt = 0.1;
for ($t = $start_t; $t < $end_t; $t += $dt) {

  $arg = 0.0 - ($t - $x0)*($t - $x0)/(2.0*$sigma*$sigma);
  $y = $amp*exp($arg);
  $noise = $noise_stdev*gaussian_rand();

  printf " %8.3f %8.3f \n", $t, $y + $noise;

}






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

