#! /usr/local/bin/perl -w

# This program generates (the frames of) a video of relativistic acceleration
# Written by David A. Madore (<URL: http://www.madore.org/~david/ >)
# around 2014-03-29.

# For more information, see the description on
# <URL: http://www.youtube.com/watch?v=0uunSMipnxA >

# This program is in the Public Domain: you can do whatever you wish
# with it (but it also comes with no guarantee whatsoever).  However,
# I would appreciate being credited if significant parts of it, or
# images produced with it, are being used somewhere.

use strict;
use warnings;
use Math::Trig;
use POSIX qw(floor waitpid);
use File::Temp qw(tempfile);

my $blackbody_octave = 25;
my @blackbody = (
  # (linear) sRGB components of blackbody color ***divided by total power***
  # (up to a constant) for logarithmically evenly spaced temperatures
  # from 1625K to 26000K (so, 25=:$blackbody_octave steps for a factor 2):
  [0.0036, 0.0006, -0.0001],  # T=1625K
  [0.0047, 0.0008, -0.0001],  # T=1671K
  [0.0061, 0.0012, -0.0001],  # T=1718K
  [0.0079, 0.0016, -0.0001],  # T=1766K
  [0.0101, 0.0021, -0.0001],  # T=1816K
  [0.0128, 0.0029, -0.0000],  # T=1867K
  [0.0161, 0.0038, 0.0000],  # T=1919K
  [0.0200, 0.0050, 0.0001],  # T=1973K
  [0.0246, 0.0065, 0.0003],  # T=2029K
  [0.0300, 0.0083, 0.0005],  # T=2086K
  [0.0364, 0.0106, 0.0008],  # T=2144K
  [0.0436, 0.0133, 0.0013],  # T=2204K
  [0.0519, 0.0166, 0.0019],  # T=2266K
  [0.0613, 0.0205, 0.0027],  # T=2330K
  [0.0718, 0.0251, 0.0038],  # T=2396K
  [0.0835, 0.0304, 0.0052],  # T=2463K
  [0.0965, 0.0366, 0.0070],  # T=2532K
  [0.1106, 0.0436, 0.0092],  # T=2603K
  [0.1259, 0.0517, 0.0120],  # T=2677K
  [0.1424, 0.0607, 0.0154],  # T=2752K
  [0.1600, 0.0708, 0.0195],  # T=2829K
  [0.1786, 0.0819, 0.0243],  # T=2909K
  [0.1982, 0.0941, 0.0301],  # T=2991K
  [0.2186, 0.1074, 0.0368],  # T=3075K
  [0.2397, 0.1218, 0.0446],  # T=3161K
  [0.2614, 0.1372, 0.0534],  # T=3250K
  [0.2834, 0.1536, 0.0635],  # T=3341K
  [0.3057, 0.1708, 0.0748],  # T=3435K
  [0.3280, 0.1889, 0.0874],  # T=3532K
  [0.3501, 0.2077, 0.1013],  # T=3631K
  [0.3720, 0.2271, 0.1165],  # T=3733K
  [0.3933, 0.2469, 0.1330],  # T=3838K
  [0.4139, 0.2670, 0.1507],  # T=3946K
  [0.4337, 0.2872, 0.1697],  # T=4057K
  [0.4524, 0.3075, 0.1897],  # T=4171K
  [0.4700, 0.3275, 0.2108],  # T=4288K
  [0.4862, 0.3473, 0.2328],  # T=4409K
  [0.5011, 0.3665, 0.2556],  # T=4533K
  [0.5144, 0.3850, 0.2789],  # T=4660K
  [0.5261, 0.4028, 0.3027],  # T=4791K
  [0.5362, 0.4195, 0.3268],  # T=4926K
  [0.5445, 0.4352, 0.3510],  # T=5065K
  [0.5512, 0.4497, 0.3750],  # T=5207K
  [0.5561, 0.4629, 0.3987],  # T=5353K
  [0.5593, 0.4747, 0.4219],  # T=5504K
  [0.5608, 0.4851, 0.4444],  # T=5659K
  [0.5606, 0.4939, 0.4660],  # T=5818K
  [0.5589, 0.5012, 0.4866],  # T=5981K
  [0.5557, 0.5069, 0.5060],  # T=6149K
  [0.5510, 0.5111, 0.5240],  # T=6322K
  [0.5449, 0.5137, 0.5405],  # T=6500K
  [0.5376, 0.5148, 0.5555],  # T=6683K
  [0.5291, 0.5144, 0.5687],  # T=6871K
  [0.5196, 0.5125, 0.5802],  # T=7064K
  [0.5091, 0.5093, 0.5899],  # T=7262K
  [0.4978, 0.5048, 0.5978],  # T=7467K
  [0.4857, 0.4990, 0.6038],  # T=7676K
  [0.4729, 0.4922, 0.6080],  # T=7892K
  [0.4596, 0.4842, 0.6103],  # T=8114K
  [0.4459, 0.4753, 0.6109],  # T=8342K
  [0.4318, 0.4656, 0.6097],  # T=8577K
  [0.4174, 0.4551, 0.6069],  # T=8818K
  [0.4028, 0.4439, 0.6025],  # T=9066K
  [0.3881, 0.4321, 0.5966],  # T=9321K
  [0.3734, 0.4198, 0.5893],  # T=9583K
  [0.3587, 0.4071, 0.5808],  # T=9852K
  [0.3441, 0.3941, 0.5710],  # T=10129K
  [0.3296, 0.3809, 0.5601],  # T=10414K
  [0.3153, 0.3675, 0.5482],  # T=10707K
  [0.3013, 0.3540, 0.5355],  # T=11008K
  [0.2875, 0.3404, 0.5220],  # T=11317K
  [0.2740, 0.3269, 0.5079],  # T=11635K
  [0.2608, 0.3135, 0.4933],  # T=11962K
  [0.2480, 0.3002, 0.4782],  # T=12299K
  [0.2356, 0.2871, 0.4627],  # T=12645K
  [0.2236, 0.2743, 0.4470],  # T=13000K
  [0.2119, 0.2616, 0.4311],  # T=13365K
  [0.2007, 0.2493, 0.4151],  # T=13741K
  [0.1899, 0.2372, 0.3992],  # T=14128K
  [0.1795, 0.2255, 0.3832],  # T=14525K
  [0.1696, 0.2141, 0.3674],  # T=14933K
  [0.1600, 0.2031, 0.3518],  # T=15353K
  [0.1509, 0.1925, 0.3364],  # T=15785K
  [0.1422, 0.1823, 0.3212],  # T=16228K
  [0.1339, 0.1724, 0.3064],  # T=16685K
  [0.1259, 0.1629, 0.2919],  # T=17154K
  [0.1184, 0.1538, 0.2778],  # T=17636K
  [0.1112, 0.1451, 0.2640],  # T=18132K
  [0.1045, 0.1368, 0.2507],  # T=18641K
  [0.0980, 0.1289, 0.2378],  # T=19165K
  [0.0919, 0.1213, 0.2253],  # T=19704K
  [0.0861, 0.1141, 0.2133],  # T=20258K
  [0.0807, 0.1072, 0.2018],  # T=20828K
  [0.0755, 0.1007, 0.1907],  # T=21413K
  [0.0707, 0.0945, 0.1800],  # T=22015K
  [0.0661, 0.0887, 0.1699],  # T=22634K
  [0.0618, 0.0831, 0.1601],  # T=23271K
  [0.0577, 0.0779, 0.1508],  # T=23925K
  [0.0539, 0.0729, 0.1420],  # T=24597K
  [0.0503, 0.0682, 0.1335],  # T=25289K
  [0.0469, 0.0638, 0.1255],  # T=26000K
);

my $monochrom_octave = 75;
my @monochrom = (
  # (linear) sRGB components of monochromatic color (for fixed power)
  # (up to a constant) for logarithmically evenly spaced frequencies
  # (75=:$monochrom_octave steps for a factor 2):
  [0.0000, 0.0000, 0.0000],  # λ=881nm
  [0.0000, 0.0000, 0.0000],  # λ=873nm
  [0.0000, 0.0000, 0.0000],  # λ=865nm
  [0.0000, 0.0000, 0.0000],  # λ=857nm
  [0.0000, 0.0000, 0.0000],  # λ=849nm
  [0.0000, 0.0000, 0.0000],  # λ=841nm
  [0.0000, 0.0000, 0.0000],  # λ=834nm
  [0.0000, -0.0000, -0.0000],  # λ=826nm
  [0.0000, -0.0000, -0.0000],  # λ=818nm
  [0.0000, -0.0000, -0.0000],  # λ=811nm
  [0.0000, -0.0000, -0.0000],  # λ=803nm
  [0.0000, -0.0000, -0.0000],  # λ=796nm
  [0.0001, -0.0000, -0.0000],  # λ=789nm
  [0.0001, -0.0000, -0.0000],  # λ=782nm
  [0.0002, -0.0000, -0.0000],  # λ=774nm
  [0.0003, -0.0000, -0.0000],  # λ=767nm
  [0.0004, -0.0000, -0.0000],  # λ=760nm
  [0.0007, -0.0001, -0.0000],  # λ=753nm
  [0.0012, -0.0001, -0.0000],  # λ=746nm
  [0.0020, -0.0002, -0.0000],  # λ=739nm
  [0.0033, -0.0004, -0.0000],  # λ=733nm
  [0.0053, -0.0006, -0.0000],  # λ=726nm
  [0.0084, -0.0009, -0.0001],  # λ=719nm
  [0.0133, -0.0014, -0.0001],  # λ=713nm
  [0.0207, -0.0022, -0.0001],  # λ=706nm
  [0.0320, -0.0035, -0.0002],  # λ=699nm
  [0.0495, -0.0054, -0.0003],  # λ=693nm
  [0.0792, -0.0086, -0.0005],  # λ=687nm
  [0.1243, -0.0133, -0.0009],  # λ=680nm
  [0.1828, -0.0194, -0.0013],  # λ=674nm
  [0.2725, -0.0286, -0.0020],  # λ=668nm
  [0.4025, -0.0417, -0.0030],  # λ=662nm
  [0.5707, -0.0578, -0.0044],  # λ=656nm
  [0.7789, -0.0763, -0.0063],  # λ=650nm
  [1.0296, -0.0960, -0.0089],  # λ=644nm
  [1.3074, -0.1134, -0.0123],  # λ=638nm
  [1.5997, -0.1250, -0.0168],  # λ=632nm
  [1.9114, -0.1285, -0.0227],  # λ=626nm
  [2.1994, -0.1160, -0.0300],  # λ=620nm
  [2.4060, -0.0791, -0.0387],  # λ=615nm
  [2.5240, -0.0137, -0.0492],  # λ=609nm
  [2.5438, 0.0833, -0.0615],  # λ=603nm
  [2.4613, 0.2102, -0.0755],  # λ=598nm
  [2.2878, 0.3617, -0.0908],  # λ=592nm
  [2.0372, 0.5315, -0.1067],  # λ=587nm
  [1.7382, 0.7074, -0.1224],  # λ=581nm
  [1.4085, 0.8790, -0.1370],  # λ=576nm
  [1.0686, 1.0389, -0.1498],  # λ=571nm
  [0.7349, 1.1807, -0.1602],  # λ=566nm
  [0.4195, 1.3007, -0.1676],  # λ=560nm
  [0.1301, 1.3953, -0.1716],  # λ=555nm
  [-0.1278, 1.4650, -0.1718],  # λ=550nm
  [-0.3524, 1.5105, -0.1680],  # λ=545nm
  [-0.5430, 1.5289, -0.1590],  # λ=540nm
  [-0.6982, 1.5185, -0.1445],  # λ=535nm
  [-0.8184, 1.4787, -0.1241],  # λ=530nm
  [-0.9016, 1.4064, -0.0978],  # λ=525nm
  [-0.9370, 1.3009, -0.0628],  # λ=520nm
  [-0.9144, 1.1528, -0.0119],  # λ=516nm
  [-0.8480, 0.9845, 0.0532],  # λ=511nm
  [-0.7583, 0.8201, 0.1251],  # λ=506nm
  [-0.6565, 0.6662, 0.1990],  # λ=502nm
  [-0.5657, 0.5375, 0.2846],  # λ=497nm
  [-0.4892, 0.4317, 0.3918],  # λ=492nm
  [-0.4228, 0.3412, 0.5271],  # λ=488nm
  [-0.3610, 0.2613, 0.6977],  # λ=483nm
  [-0.2960, 0.1848, 0.9027],  # λ=479nm
  [-0.2262, 0.1091, 1.1296],  # λ=474nm
  [-0.1522, 0.0367, 1.3662],  # λ=470nm
  [-0.0758, -0.0309, 1.6000],  # λ=466nm
  [-0.0049, -0.0851, 1.7546],  # λ=461nm
  [0.0596, -0.1287, 1.8455],  # λ=457nm
  [0.1155, -0.1633, 1.8928],  # λ=453nm
  [0.1613, -0.1902, 1.9140],  # λ=449nm
  [0.1979, -0.2112, 1.9219],  # λ=445nm
  [0.2225, -0.2236, 1.8946],  # λ=441nm
  [0.2328, -0.2251, 1.8071],  # λ=437nm
  [0.2269, -0.2134, 1.6442],  # λ=433nm
  [0.2038, -0.1875, 1.4018],  # λ=429nm
  [0.1648, -0.1494, 1.0930],  # λ=425nm
  [0.1174, -0.1050, 0.7565],  # λ=421nm
  [0.0786, -0.0696, 0.4965],  # λ=417nm
  [0.0513, -0.0451, 0.3202],  # λ=413nm
  [0.0330, -0.0289, 0.2044],  # λ=409nm
  [0.0207, -0.0180, 0.1265],  # λ=405nm
  [0.0142, -0.0122, 0.0861],  # λ=402nm
  [0.0097, -0.0084, 0.0589],  # λ=398nm
  [0.0060, -0.0052, 0.0363],  # λ=394nm
  [0.0039, -0.0034, 0.0237],  # λ=391nm
  [0.0025, -0.0022, 0.0151],  # λ=387nm
  [0.0016, -0.0014, 0.0098],  # λ=384nm
  [0.0012, -0.0010, 0.0070],  # λ=380nm
  [0.0008, -0.0007, 0.0047],  # λ=377nm
  [0.0005, -0.0004, 0.0030],  # λ=373nm
  [0.0003, -0.0003, 0.0020],  # λ=370nm
  [0.0002, -0.0002, 0.0014],  # λ=366nm
  [0.0002, -0.0001, 0.0009],  # λ=363nm
  [0.0001, -0.0001, 0.0004],  # λ=360nm
  [0.0000, 0.0000, 0.0000],  # λ=356nm
  [0.0000, 0.0000, 0.0000],  # λ=353nm
  [0.0000, 0.0000, 0.0000],  # λ=350nm
);

sub inter_lookup {
    # Lookup an element in a table, performing linear interpolation.
    my $tref = shift;
    my $idx = shift;
    if ( $idx <= 0 ) {
	return $tref->[0];
    } elsif ( $idx >= scalar(@$tref)-1 ) {
	return $tref->[scalar(@$tref)-1];
    } else {
	my $i0 = int(floor($idx));
	my $p = $idx - $i0;  my $q = 1 - $p;
	my @v0 = @{$tref->[$i0]};
	my @v1 = @{$tref->[$i0+1]};
	my @v = ();
	for ( my $k=0 ; $k<scalar(@v0) ; $k++ ) {
	    push @v, $v0[$k] * $q + $v1[$k] * $p;
	}
	return \@v;
    }
}

my $preimgw = 1280;  # Image width
my $preimgh = 720;   # Image height
my $preimgsc = 500;  # Image scale (45° angle mark)
my $aafactor = 2;  # Compute intermediate image this times larger, to anti-alias
#my $preimgw = 640;
#my $preimgh = 360;
#my $preimgsc = 250;
my $imgw = $preimgw * $aafactor;
my $imgh = $preimgh * $aafactor;
my $imgsc = $preimgsc * $aafactor;

my $frcnt = 2500;   # Number of frames
#my $frcnt = 8;

my $accel = 0.2;    # Acceleration value

my $attlength = 1;  # Light attenuation distance

my $nbthreads = 4;  # Number of computing processes to spawn
my @pidtab;
my $thrnum;
for ( $thrnum=0 ; $thrnum<$nbthreads ; $thrnum++ ) {  # Spawn workers
    if ( my $pid = fork ) {
	push @pidtab, $pid;
    } else {
	last;
    }
}

if ( $thrnum == $nbthreads ) {  # Father process just waits for children
    foreach my $pid ( @pidtab ) {
	POSIX::waitpid($pid, 0);
    }
}

my $frnum;  # Frame index
my ($ppmf, $ppmfname);

for ( $frnum = 0 ; $frnum < $frcnt ; $frnum ++ ) {

next unless ($frnum % $nbthreads) == $thrnum;  # Only compute frames we mean to

print "$frnum\n";

$ppmf = File::Temp->new(SUFFIX=>".ppm", UNLINK=>1);
$ppmfname = $ppmf->filename;

printf $ppmf "P3\n%d %d\n255\n", $imgw, $imgh;

my $tim = 20*($frnum / $frcnt);  # Proper time (s), in natural units

my $ctim = cosh($accel*$tim);    # = dz/ds
my $stim = sinh($accel*$tim);    # = dt/ds

my $ot = $stim/$accel;           # t
my $oz = ($ctim - 1)/$accel;     # z

sub min {
    my $ret = shift;
    while ( defined(my $v = shift) ) {
	$ret = $v if $v < $ret;
    }
    return $ret;
}

sub max {
    my $ret = shift;
    while ( defined(my $v = shift) ) {
	$ret = $v if $v > $ret;
    }
    return $ret;
}

sub do_point {
    my $colr = shift;  my $colg = shift;  my $colb = shift;
    my $min = min($colr, $colg, $colb);
    if ( $min < 0 ) {  # Out of color gamut: fix by desaturating
	$colr -= $min;  $colg -= $min;  $colb -= $min;
    }
    my $max = max($colr, $colg, $colb);
    if ( $max > 1 ) {  # Out of luminosity gamut: fix by rescaling
	$colr /= $max;  $colg /= $max;  $colb /= $max;
    }
    # Apply sRGB gamma correction
    $colr = $colr ** (1/2.2);
    $colg = $colg ** (1/2.2);
    $colb = $colb ** (1/2.2);
    my $vr = int(floor($colr*255+0.5));
    my $vg = int(floor($colg*255+0.5));
    my $vb = int(floor($colb*255+0.5));
    printf $ppmf "\t%d %d %d", $vr, $vg, $vb;
}

for ( my $q=0 ; $q<$imgh ; $q++ ) {
  PIXEL:
    for ( my $p=0 ; $p<$imgw ; $p++ ) {
	# Direction of incoming light ray in observer's reference
	my $u = ($p - $imgw/2) / $imgsc;
	my $v = ($q - $imgh/2) / $imgsc;
	my $ttt = sqrt(1 + $u * $u + $v * $v);
	# Direction of incoming light ray in fixed reference
	my $rdz = $ctim - $stim * $ttt;
	my $rdx = $u;  my $rdy = $v;
	my $rdt = - $ttt*$ctim + $stim;
	# Vector of incoming light ray
	my $l;
	my ($rt, $rx, $ry, $rz);
	if ( $rdy > 0 ) {
	    $l = 0;
	    $ry = 0.25;
	} elsif ( $rdy < 0 ) {
	    $l = 1;
	    $ry = -0.25;
	} else {
	    do_point 0,0,0;
	    next PIXEL;
	}
	$rx = $ry * $rdx / $rdy;  $rz = $ry * $rdz / $rdy;
	$rt = $ry * $rdt / $rdy;
	my $dst = - $rt;  # Distance traveled by light (in fixed frame)
	# die unless abs($dst - sqrt($rx*$rx + $ry*$ry + $rz*$rz)) < 1.e-8;
	my $doppler = -$ttt / $rdt;  # Doppler factor
	# Source of light ray
	my ($st, $sx, $sy, $sz);
	$sx = $rx;  $sy = $ry;  $sz = $oz + $rz;  $st = $ot + $rt;
	my ($colr, $colg, $colb);
	if ( $l ) {  # Ceiling tile
	    my $on = (($st/2 - floor($st/2)) > 0.5);  # Blink effect
	    my $fracsx = $sx - floor($sx);
	    my $fracsz = $sz - floor($sz);
	    my $inter = (abs($fracsx-0.5)<0.01)||(abs($fracsz-0.5)<0.01);
	    unless ( $on || $inter ) {
		do_point 0,0,0;
		next PIXEL;
	    }
	    my $magn = ($inter?0.2:0.1);
	    # Doppler effect on intensity is fourth power
	    $magn *= ($doppler*$doppler*$doppler*$doppler);
	    # Attenuation
	    $magn /= exp($dst/$attlength);
	    my $freqidx = 35 + $monochrom_octave*log($doppler)/log(2);
	    my $cref = inter_lookup(\@monochrom, $freqidx);
	    $colr = $magn * $cref->[0];
	    $colg = $magn * $cref->[1];
	    $colb = $magn * $cref->[2];
	} else {
	    my $fracsx = $sx - floor($sx);
	    my $fracsz = $sz - floor($sz);
	    my $bright = (abs($fracsx-0.5)<0.1)||(abs($fracsz-0.5)<0.1);
	    my $magn = ($bright?0.4:0.025);
	    # Doppler effect on intensity is fourth power
	    $magn *= ($doppler*$doppler*$doppler*$doppler);
	    # Attenuation
	    $magn /= exp($dst/$attlength);
	    my $tempidx = ($bright?50:25) + $blackbody_octave*log($doppler)/log(2);
	    my $cref = inter_lookup(\@blackbody, $tempidx);
	    $colr = $magn * $cref->[0];
	    $colg = $magn * $cref->[1];
	    $colb = $magn * $cref->[2];
	}
	do_point $colr, $colg, $colb;
    }
    print $ppmf "\n";
}

close $ppmf;

# And now some ImageMagick!
system ("convert", $ppmfname,
	"-resize", sprintf("%dx%d", $preimgw, $preimgh),
	"-font", "Source-Code-Pro-Regular",
	"-pointsize", "36", "-fill", "white", "-undercolor", "\#00000080",
	"-gravity", "southwest", "-annotate", "+6+6",
	sprintf("v/c=%5.3f", $stim/$ctim),
	sprintf("output%04d.png", $frnum));

}
