#!/usr/bin/perl -w
use strict;
use Config;

use constant bits => $Config{randbits};
use constant b2 => int(bits/2);
use constant b1 => bits - b2;

#sub rand2 {return rand 1, rand 1}
sub rand2 {
  my $r = int rand(2**bits);
  my $r2 = $r >> b1;
  return((($r & ((1<<b1)-1))+0.5)/(1<<b1), ($r2 + 0.5)/(1<<b2));
}

#my @r;
#@r = rand2, print "@r\n" for 1..40;

# INT dr r exp(-r^2/2) = 1 - exp(-r^2/2) = 1 - u
# r = sqrt(-2 log u)

my $prev_n;
sub normal_dist () {
  if (defined $prev_n) {
    my $v = $prev_n;
    undef $prev_n;
    return $v;
  }
  my ($u, $u1) = rand2;		# Random from 0 to 1
  my $r = sqrt(-2 * log $u);
  $u1 *= 2*atan2(0,-1);
  $prev_n = $r * sin $u1;
  return $r * cos $u1;
}

#my $s = 0;
#my $l = 9999;
#$s += normal_dist**2 for 1..$l;
#print sqrt($s) - sqrt($l);		# Approx normal with sigma=1

# Old formula Y' = 0.299 R' + 0.587 G' + 0.114 B'.
# Modern HDTV systems use a different definition of RGB
# phosphor characteristics and the ITU-R BT.709 formula
# Y' =  0.2126  R'  + 0.7152 G' + 0.0722 B'.

sub vector_product ($$$$$$) {
  my ($x1, $y1, $z1, $x2, $y2, $z2) = @_;
  ($y1 * $z2 - $y2 * $z1, $z1 * $x2 - $z2 * $x1, $x1 * $y2 - $x2 * $y1)
}

sub vector_norm ($$$) {
  my ($x, $y, $z) = @_;
  sqrt($x**2 + $y**2 + $z**2);
}

sub normalize_vector {
  my ($x, $y, $z) = @_;
  my $n = vector_norm($x, $y, $z);
  $x/$n, $y/$n, $z/$n;
}

# Given a vector, returns two orthonormal vector perpendicular to it;
# Assume that the initial vector is not proportional to 1,0,0
sub normal_vectors ($$$) {
  my ($x, $y, $z) = normalize_vector @_;
  my ($x1, $y1, $z1) = normalize_vector vector_product 1,0,0, $x, $y, $z;
  $x1, $y1, $z1, vector_product $x, $y, $z, $x1, $y1, $z1;
}

my ($x1, $y1, $z1, $x2, $y2, $z2);
BEGIN {
 ($x1, $y1, $z1, $x2, $y2, $z2) = normal_vectors 0.2126, 0.7152, 0.0722;
 #my @v = vector_product($x1, $y1, $z1, $x2, $y2, $z2);
 #print "@v\n";
}
use constant {R1 => $x1, G1 => $y1, B1 => $z1, R2 => $x2, G2 => $y2, B2 => $z2};

use constant color_noise => $ARGV[0];
use constant color_base => 128;

sub noise_color () {
  my ($r1, $r2) = (normal_dist, normal_dist);
  int(0.5 + color_base + color_noise * R1 * $r1 + color_noise * R2 * $r2),
  int(0.5 + color_base + color_noise * G1 * $r1 + color_noise * G2 * $r2),
  int(0.5 + color_base + color_noise * B1 * $r1 + color_noise * B2 * $r2)
}

use constant {w => 200, h => 400};

use GD;

die "Usage: $0 COLORNOISE\n" unless @ARGV == 1;

my $im = GD::Image->newTrueColor(2*w, h);
my $black = $im->colorAllocate(0,0,0);
for my $x (0..w-1) {
  for my $y (0..h-1) {
    my (@c) = noise_color;
	#print "@c\n";
    my $i = $im->colorExact(@c);
    $i = $im->colorAllocate(@c) unless $i >= 0;
    $im->setPixel($x, $y, $i);
    my $g = int(0.5 + color_base + color_noise/(h-1)*$y*normal_dist);
    $i = $im->colorExact($g, $g, $g);
    $i = $im->colorAllocate($g, $g, $g) unless $i >= 0;
    $im->setPixel($x+w, $y, $i);
  }
}
for my $frac (0..10) {
  my $h = int(0.5 + (h - 1) * $frac/10);
  $im->line(int(w*1.9), $h, 2*w-1, $h, $black);
}
binmode STDOUT;
print STDOUT $im->png;