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

# Generate POVray files of an animation representing 4-dimensional
# regular polytopes seen as a tesselation of the 3-sphere (through
# which the observer travels).

# Written by David A. Madore <URL: http://www.madore.org/~david/ >
# on 2010-11-12.

# This program is in the Public Domain: do what you will with it, but
# I won't take any responsability.

# How to use this script: choose a regular solid by setting the
# @vertices and $edgelensq variables below, and possibly adjust
# $nbframes.  Then run the script to produce that many .pov files
# named frame0000.pov to frame1199.pov (assuming 1200 frames, say).
# It is probably desirable to rerun the script if the chosen travel
# path crosses a beam (a warning will then be printed).  Then use
# POVray to render each of these .pov files (e.g.,
# for i in `seq 0 1199` ; do povray +I `printf frame%04d.pov $i` render.ini ; done
# where render.ini is your favorite POVray INI config file), and
# create an animation from the resulting images.

# Quick mathematical notes: the way this works is by using gnomonic
# projection, centered on the observer, from the 3-sphere to Euclidean
# 3-space.  The reason it works is that gnomonic projection preserves
# alignment, and preserves angles at the center of projection.  In
# essence, this shows that what an observer in the 3-sphere sees of a
# configuration of lines is the same thing as what an observer in
# Euclidean space sees of the gnomonic projection of this
# configuration.  There is a problem, however, in that the 3-sphere
# has two hemispheres: the "near" or "cis" hemisphere being that
# centered around the observer whereas the "far" or "trans" hemisphere
# is centered around its antipodal point; points in the "cis"
# hemisphere pose no particular difficulty, but points in the "trans"
# hemisphere obey a kind of reverse perspective which it is hard to
# persuade an ordinary raytracer to render.  A hack around this
# consists of drawing untextured 2D lines for the "trans" lines (by
# drawing them on a distant plane parallel to the plane of
# projection).  Set $use_trans_hack to 0 to deactivate this (and
# everything in the the "trans" hemisphere will appear to disappear in
# a kind of black fog).

# Notes for myself:

# My render.ini contains:

# Width = 640
# Height = 480
# Display = off
# Pause_when_Done = off
# Bounding_Threshold = 3
# Test_Abort=On
# Test_Abort_Count=100
# Quality=11
# Antialias_Depth=3
# Antialias=On
# Antialias_Threshold=0.1
# Jitter_Amount=0.5
# Jitter=On

# Encoding to x264 using ffmpeg:

# ffmpeg -metadata title="120-cell navigation" \
# -metadata author="David A. Madore" -metadata copyright="Public Domain" \
# -r 25 -s 640x480 -aspect 4:3 -i frame%04d.png -vcodec libx264 \
# -vpre slower -vpre main -refs 8 -bf 2 -flags2 dct8x8+wpred \
# -me_method umh -subq 7 -trellis 2 -directpred 3 -crf 22.5 \
# -aspect 4:3 120cell.avi

use strict;
use warnings;

use constant SQRT5 => sqrt(5);
use constant PHI => (1+SQRT5)/2;
use constant PI => 4 * atan2(1, 1);

# For each regular polytope (inscribed in a sphere centered at the
# origin and with some unspecified radius), the list of vertices, and
# the square of the edge length.

my @vertices_16cell = (
    [ 0, 0, 0, 2 ], [ 2, 0, 0, 0 ], [ 0, 0, 0, -2 ], [ 0, 2, 0, 0 ],
    [ -2, 0, 0, 0 ], [ 0, 0, 2, 0 ], [ 0, -2, 0, 0 ], [ 0, 0, -2, 0 ],
);
my $edgelensq_16cell = 8;

my @vertices_tesseract = (
    [ 1, 1, 1, 1 ], [ -1, 1, 1, 1 ], [ 1, -1, 1, 1 ], [ 1, 1, -1, 1 ],
    [ 1, 1, 1, -1 ], [ -1, -1, 1, 1 ], [ -1, 1, -1, 1 ], [ -1, 1, 1, -1 ],
    [ 1, -1, -1, 1 ], [ 1, -1, 1, -1 ], [ 1, 1, -1, -1 ], [ -1, -1, -1, 1 ],
    [ -1, -1, 1, -1 ], [ -1, 1, -1, -1 ], [ 1, -1, -1, -1 ], [ -1, -1, -1, -1 ],
);
my $edgelensq_tesseract = 4;

my @vertices_24cell = (
    [ 1, 1, 1, 1 ], [ -1, 1, 1, 1 ], [ 1, -1, 1, 1 ], [ 1, 1, -1, 1 ],
    [ 1, 1, 1, -1 ], [ -1, -1, 1, 1 ], [ -1, 1, -1, 1 ], [ -1, 1, 1, -1 ],
    [ 1, -1, -1, 1 ], [ 1, -1, 1, -1 ], [ 1, 1, -1, -1 ], [ -1, -1, -1, 1 ],
    [ -1, -1, 1, -1 ], [ -1, 1, -1, -1 ], [ 1, -1, -1, -1 ], [ -1, -1, -1, -1 ],
    [ 0, 0, 0, 2 ], [ 2, 0, 0, 0 ], [ 0, 0, 0, -2 ], [ 0, 2, 0, 0 ],
    [ -2, 0, 0, 0 ], [ 0, 0, 2, 0 ], [ 0, -2, 0, 0 ], [ 0, 0, -2, 0 ],
);
my $edgelensq_24cell = 4;

my @vertices_600cell = (
    [ 1, 1, 1, 1 ], [ -1, 1, 1, 1 ], [ 1, -1, 1, 1 ], [ 1, 1, -1, 1 ],
    [ 1, 1, 1, -1 ], [ -1, -1, 1, 1 ], [ -1, 1, -1, 1 ], [ -1, 1, 1, -1 ],
    [ 1, -1, -1, 1 ], [ 1, -1, 1, -1 ], [ 1, 1, -1, -1 ], [ -1, -1, -1, 1 ],
    [ -1, -1, 1, -1 ], [ -1, 1, -1, -1 ], [ 1, -1, -1, -1 ], [ -1, -1, -1, -1 ],
    [ 0, 0, 0, 2 ], [ 2, 0, 0, 0 ], [ 0, 0, 0, -2 ], [ 0, 2, 0, 0 ],
    [ -2, 0, 0, 0 ], [ 0, 0, 2, 0 ], [ 0, -2, 0, 0 ], [ 0, 0, -2, 0 ],
    [ 1, PHI, 1/PHI, 0 ], [ 1/PHI, 1, PHI, 0 ], [ 1, 0, PHI, 1/PHI ], [ -1, PHI, 1/PHI, 0 ],
    [ 1, -(PHI), 1/PHI, 0 ], [ 1, PHI, -1/PHI, 0 ], [ PHI, 1/PHI, 1, 0 ], [ 1/PHI, 0, 1, PHI ],
    [ -1/PHI, 1, PHI, 0 ], [ 1/PHI, -1, PHI, 0 ], [ 1/PHI, 1, -(PHI), 0 ], [ PHI, 1, 0, 1/PHI ],
    [ 1, 1/PHI, 0, PHI ], [ -1, 0, PHI, 1/PHI ], [ 1, 0, -(PHI), 1/PHI ], [ 1, 0, PHI, -1/PHI ],
    [ -1, -(PHI), 1/PHI, 0 ], [ -1, PHI, -1/PHI, 0 ], [ 1, -(PHI), -1/PHI, 0 ], [ PHI, 0, 1/PHI, 1 ],
    [ -(PHI), 1/PHI, 1, 0 ], [ PHI, -1/PHI, 1, 0 ], [ PHI, 1/PHI, -1, 0 ], [ 1/PHI, PHI, 0, 1 ],
    [ -1/PHI, 0, 1, PHI ], [ 1/PHI, 0, -1, PHI ], [ 1/PHI, 0, 1, -(PHI) ], [ -1/PHI, -1, PHI, 0 ],
    [ -1/PHI, 1, -(PHI), 0 ], [ 1/PHI, -1, -(PHI), 0 ], [ 0, PHI, 1, 1/PHI ], [ -(PHI), 1, 0, 1/PHI ],
    [ PHI, -1, 0, 1/PHI ], [ PHI, 1, 0, -1/PHI ], [ 0, 1, 1/PHI, PHI ], [ -1, 1/PHI, 0, PHI ],
    [ 1, -1/PHI, 0, PHI ], [ 1, 1/PHI, 0, -(PHI) ], [ -1, 0, -(PHI), 1/PHI ], [ -1, 0, PHI, -1/PHI ],
    [ 1, 0, -(PHI), -1/PHI ], [ -1, -(PHI), -1/PHI, 0 ], [ -(PHI), 0, 1/PHI, 1 ], [ PHI, 0, -1/PHI, 1 ],
    [ PHI, 0, 1/PHI, -1 ], [ -(PHI), -1/PHI, 1, 0 ], [ -(PHI), 1/PHI, -1, 0 ], [ PHI, -1/PHI, -1, 0 ],
    [ 0, 1/PHI, PHI, 1 ], [ -1/PHI, PHI, 0, 1 ], [ 1/PHI, -(PHI), 0, 1 ], [ 1/PHI, PHI, 0, -1 ],
    [ -1/PHI, 0, -1, PHI ], [ -1/PHI, 0, 1, -(PHI) ], [ 1/PHI, 0, -1, -(PHI) ], [ -1/PHI, -1, -(PHI), 0 ],
    [ 0, -(PHI), 1, 1/PHI ], [ 0, PHI, -1, 1/PHI ], [ 0, PHI, 1, -1/PHI ], [ -(PHI), -1, 0, 1/PHI ],
    [ -(PHI), 1, 0, -1/PHI ], [ PHI, -1, 0, -1/PHI ], [ 0, -1, 1/PHI, PHI ], [ 0, 1, -1/PHI, PHI ],
    [ 0, 1, 1/PHI, -(PHI) ], [ -1, -1/PHI, 0, PHI ], [ -1, 1/PHI, 0, -(PHI) ], [ 1, -1/PHI, 0, -(PHI) ],
    [ -1, 0, -(PHI), -1/PHI ], [ -(PHI), 0, -1/PHI, 1 ], [ -(PHI), 0, 1/PHI, -1 ], [ PHI, 0, -1/PHI, -1 ],
    [ -(PHI), -1/PHI, -1, 0 ], [ 0, -1/PHI, PHI, 1 ], [ 0, 1/PHI, -(PHI), 1 ], [ 0, 1/PHI, PHI, -1 ],
    [ -1/PHI, -(PHI), 0, 1 ], [ -1/PHI, PHI, 0, -1 ], [ 1/PHI, -(PHI), 0, -1 ], [ -1/PHI, 0, -1, -(PHI) ],
    [ 0, -(PHI), -1, 1/PHI ], [ 0, -(PHI), 1, -1/PHI ], [ 0, PHI, -1, -1/PHI ], [ -(PHI), -1, 0, -1/PHI ],
    [ 0, -1, -1/PHI, PHI ], [ 0, -1, 1/PHI, -(PHI) ], [ 0, 1, -1/PHI, -(PHI) ], [ -1, -1/PHI, 0, -(PHI) ],
    [ -(PHI), 0, -1/PHI, -1 ], [ 0, -1/PHI, -(PHI), 1 ], [ 0, -1/PHI, PHI, -1 ], [ 0, 1/PHI, -(PHI), -1 ],
    [ -1/PHI, -(PHI), 0, -1 ], [ 0, -(PHI), -1, -1/PHI ], [ 0, -1, -1/PHI, -(PHI) ], [ 0, -1/PHI, -(PHI), -1 ],
);
my $edgelensq_600cell = 1.52786404500042;

my @vertices_120cell = (
    [ 0, 0, 2, 2 ], [ 2, 0, 0, 2 ], [ 0, 0, -2, 2 ], [ 0, 0, 2, -2 ],
    [ -2, 0, 0, 2 ], [ 2, 2, 0, 0 ], [ 0, 2, 0, 2 ], [ 2, 0, 0, -2 ],
    [ 0, 0, -2, -2 ], [ 2, -2, 0, 0 ], [ 0, -2, 0, 2 ], [ -2, 0, 0, -2 ],
    [ -2, 2, 0, 0 ], [ 0, 2, 2, 0 ], [ 2, 0, 2, 0 ], [ 0, 2, 0, -2 ],
    [ -2, -2, 0, 0 ], [ 0, 2, -2, 0 ], [ 2, 0, -2, 0 ], [ 0, -2, 0, -2 ],
    [ 0, -2, 2, 0 ], [ -2, 0, 2, 0 ], [ 0, -2, -2, 0 ], [ -2, 0, -2, 0 ],
    [ 1, 1, 1, SQRT5 ], [ -1, 1, 1, SQRT5 ], [ SQRT5, 1, 1, 1 ], [ 1, -1, 1, SQRT5 ],
    [ 1, 1, -1, SQRT5 ], [ 1, 1, 1, -(SQRT5) ], [ SQRT5, -1, 1, 1 ], [ -1, -1, 1, SQRT5 ],
    [ -1, 1, -1, SQRT5 ], [ -1, 1, 1, -(SQRT5) ], [ -(SQRT5), 1, 1, 1 ], [ 1, SQRT5, 1, 1 ],
    [ SQRT5, 1, -1, 1 ], [ SQRT5, 1, 1, -1 ], [ 1, -1, -1, SQRT5 ], [ 1, -1, 1, -(SQRT5) ],
    [ 1, 1, -1, -(SQRT5) ], [ -(SQRT5), -1, 1, 1 ], [ 1, SQRT5, -1, 1 ], [ -1, SQRT5, 1, 1 ],
    [ SQRT5, -1, -1, 1 ], [ SQRT5, -1, 1, -1 ], [ -1, -1, -1, SQRT5 ], [ -1, -1, 1, -(SQRT5) ],
    [ -1, 1, -1, -(SQRT5) ], [ 1, -(SQRT5), 1, 1 ], [ -(SQRT5), 1, -1, 1 ], [ -(SQRT5), 1, 1, -1 ],
    [ 1, 1, SQRT5, 1 ], [ 1, SQRT5, 1, -1 ], [ SQRT5, 1, -1, -1 ], [ 1, -1, -1, -(SQRT5) ],
    [ 1, -(SQRT5), -1, 1 ], [ -1, -(SQRT5), 1, 1 ], [ -(SQRT5), -1, -1, 1 ], [ -(SQRT5), -1, 1, -1 ],
    [ -1, SQRT5, -1, 1 ], [ 1, 1, SQRT5, -1 ], [ 1, SQRT5, -1, -1 ], [ 1, -1, SQRT5, 1 ],
    [ -1, SQRT5, 1, -1 ], [ SQRT5, -1, -1, -1 ], [ -1, -1, -1, -(SQRT5) ], [ 1, 1, -(SQRT5), 1 ],
    [ 1, -(SQRT5), 1, -1 ], [ -(SQRT5), 1, -1, -1 ], [ -1, 1, SQRT5, 1 ], [ -1, -(SQRT5), -1, 1 ],
    [ 1, 1, -(SQRT5), -1 ], [ 1, -(SQRT5), -1, -1 ], [ 1, -1, -(SQRT5), 1 ], [ -1, -(SQRT5), 1, -1 ],
    [ -(SQRT5), -1, -1, -1 ], [ 1, -1, SQRT5, -1 ], [ -1, SQRT5, -1, -1 ], [ -1, 1, SQRT5, -1 ],
    [ -1, -1, SQRT5, 1 ], [ -1, 1, -(SQRT5), 1 ], [ 1, -1, -(SQRT5), -1 ], [ -1, -(SQRT5), -1, -1 ],
    [ -1, 1, -(SQRT5), -1 ], [ -1, -1, -(SQRT5), 1 ], [ -1, -1, SQRT5, -1 ], [ -1, -1, -(SQRT5), -1 ],
    [ 1/PHI**2, PHI, PHI, PHI ], [ -1/PHI**2, PHI, PHI, PHI ], [ PHI, 1/PHI**2, PHI, PHI ], [ 1/PHI**2, -(PHI), PHI, PHI ],
    [ 1/PHI**2, PHI, -(PHI), PHI ], [ 1/PHI**2, PHI, PHI, -(PHI) ], [ PHI, -1/PHI**2, PHI, PHI ], [ -1/PHI**2, -(PHI), PHI, PHI ],
    [ -1/PHI**2, PHI, -(PHI), PHI ], [ -1/PHI**2, PHI, PHI, -(PHI) ], [ -(PHI), 1/PHI**2, PHI, PHI ], [ PHI, PHI, 1/PHI**2, PHI ],
    [ PHI, 1/PHI**2, -(PHI), PHI ], [ PHI, 1/PHI**2, PHI, -(PHI) ], [ 1/PHI**2, -(PHI), -(PHI), PHI ], [ 1/PHI**2, -(PHI), PHI, -(PHI) ],
    [ 1/PHI**2, PHI, -(PHI), -(PHI) ], [ -(PHI), -1/PHI**2, PHI, PHI ], [ PHI, PHI, -1/PHI**2, PHI ], [ PHI, -1/PHI**2, -(PHI), PHI ],
    [ PHI, -1/PHI**2, PHI, -(PHI) ], [ -1/PHI**2, -(PHI), -(PHI), PHI ], [ -1/PHI**2, -(PHI), PHI, -(PHI) ], [ -1/PHI**2, PHI, -(PHI), -(PHI) ],
    [ PHI, -(PHI), 1/PHI**2, PHI ], [ -(PHI), 1/PHI**2, -(PHI), PHI ], [ -(PHI), 1/PHI**2, PHI, -(PHI) ], [ -(PHI), PHI, 1/PHI**2, PHI ],
    [ PHI, PHI, PHI, 1/PHI**2 ], [ PHI, PHI, 1/PHI**2, -(PHI) ], [ PHI, 1/PHI**2, -(PHI), -(PHI) ], [ 1/PHI**2, -(PHI), -(PHI), -(PHI) ],
    [ PHI, -(PHI), -1/PHI**2, PHI ], [ -(PHI), -1/PHI**2, -(PHI), PHI ], [ -(PHI), -1/PHI**2, PHI, -(PHI) ], [ -(PHI), PHI, -1/PHI**2, PHI ],
    [ PHI, PHI, PHI, -1/PHI**2 ], [ PHI, PHI, -1/PHI**2, -(PHI) ], [ PHI, -1/PHI**2, -(PHI), -(PHI) ], [ -1/PHI**2, -(PHI), -(PHI), -(PHI) ],
    [ -(PHI), -(PHI), 1/PHI**2, PHI ], [ PHI, PHI, -(PHI), 1/PHI**2 ], [ PHI, -(PHI), 1/PHI**2, -(PHI) ], [ -(PHI), 1/PHI**2, -(PHI), -(PHI) ],
    [ PHI, -(PHI), PHI, 1/PHI**2 ], [ -(PHI), PHI, 1/PHI**2, -(PHI) ], [ -(PHI), PHI, PHI, 1/PHI**2 ], [ -(PHI), -(PHI), -1/PHI**2, PHI ],
    [ PHI, PHI, -(PHI), -1/PHI**2 ], [ PHI, -(PHI), -1/PHI**2, -(PHI) ], [ -(PHI), -1/PHI**2, -(PHI), -(PHI) ], [ PHI, -(PHI), PHI, -1/PHI**2 ],
    [ -(PHI), PHI, -1/PHI**2, -(PHI) ], [ -(PHI), PHI, PHI, -1/PHI**2 ], [ PHI, -(PHI), -(PHI), 1/PHI**2 ], [ -(PHI), -(PHI), 1/PHI**2, -(PHI) ],
    [ -(PHI), PHI, -(PHI), 1/PHI**2 ], [ -(PHI), -(PHI), PHI, 1/PHI**2 ], [ PHI, -(PHI), -(PHI), -1/PHI**2 ], [ -(PHI), -(PHI), -1/PHI**2, -(PHI) ],
    [ -(PHI), PHI, -(PHI), -1/PHI**2 ], [ -(PHI), -(PHI), PHI, -1/PHI**2 ], [ -(PHI), -(PHI), -(PHI), 1/PHI**2 ], [ -(PHI), -(PHI), -(PHI), -1/PHI**2 ],
    [ 1/PHI, 1/PHI, 1/PHI, PHI**2 ], [ -1/PHI, 1/PHI, 1/PHI, PHI**2 ], [ PHI**2, 1/PHI, 1/PHI, 1/PHI ], [ 1/PHI, -1/PHI, 1/PHI, PHI**2 ],
    [ 1/PHI, 1/PHI, -1/PHI, PHI**2 ], [ 1/PHI, 1/PHI, 1/PHI, -(PHI)**2 ], [ PHI**2, -1/PHI, 1/PHI, 1/PHI ], [ -1/PHI, -1/PHI, 1/PHI, PHI**2 ],
    [ -1/PHI, 1/PHI, -1/PHI, PHI**2 ], [ -1/PHI, 1/PHI, 1/PHI, -(PHI)**2 ], [ -(PHI)**2, 1/PHI, 1/PHI, 1/PHI ], [ 1/PHI, PHI**2, 1/PHI, 1/PHI ],
    [ PHI**2, 1/PHI, -1/PHI, 1/PHI ], [ PHI**2, 1/PHI, 1/PHI, -1/PHI ], [ 1/PHI, -1/PHI, -1/PHI, PHI**2 ], [ 1/PHI, -1/PHI, 1/PHI, -(PHI)**2 ],
    [ 1/PHI, 1/PHI, -1/PHI, -(PHI)**2 ], [ -(PHI)**2, -1/PHI, 1/PHI, 1/PHI ], [ 1/PHI, PHI**2, -1/PHI, 1/PHI ], [ -1/PHI, PHI**2, 1/PHI, 1/PHI ],
    [ PHI**2, -1/PHI, -1/PHI, 1/PHI ], [ PHI**2, -1/PHI, 1/PHI, -1/PHI ], [ -1/PHI, -1/PHI, -1/PHI, PHI**2 ], [ -1/PHI, -1/PHI, 1/PHI, -(PHI)**2 ],
    [ -1/PHI, 1/PHI, -1/PHI, -(PHI)**2 ], [ 1/PHI, -(PHI)**2, 1/PHI, 1/PHI ], [ -(PHI)**2, 1/PHI, -1/PHI, 1/PHI ], [ -(PHI)**2, 1/PHI, 1/PHI, -1/PHI ],
    [ 1/PHI, 1/PHI, PHI**2, 1/PHI ], [ 1/PHI, PHI**2, 1/PHI, -1/PHI ], [ PHI**2, 1/PHI, -1/PHI, -1/PHI ], [ 1/PHI, -1/PHI, -1/PHI, -(PHI)**2 ],
    [ 1/PHI, -(PHI)**2, -1/PHI, 1/PHI ], [ -1/PHI, -(PHI)**2, 1/PHI, 1/PHI ], [ -(PHI)**2, -1/PHI, -1/PHI, 1/PHI ], [ -(PHI)**2, -1/PHI, 1/PHI, -1/PHI ],
    [ -1/PHI, PHI**2, -1/PHI, 1/PHI ], [ 1/PHI, 1/PHI, PHI**2, -1/PHI ], [ 1/PHI, PHI**2, -1/PHI, -1/PHI ], [ 1/PHI, -1/PHI, PHI**2, 1/PHI ],
    [ -1/PHI, PHI**2, 1/PHI, -1/PHI ], [ PHI**2, -1/PHI, -1/PHI, -1/PHI ], [ -1/PHI, -1/PHI, -1/PHI, -(PHI)**2 ], [ 1/PHI, 1/PHI, -(PHI)**2, 1/PHI ],
    [ 1/PHI, -(PHI)**2, 1/PHI, -1/PHI ], [ -(PHI)**2, 1/PHI, -1/PHI, -1/PHI ], [ -1/PHI, 1/PHI, PHI**2, 1/PHI ], [ -1/PHI, -(PHI)**2, -1/PHI, 1/PHI ],
    [ 1/PHI, 1/PHI, -(PHI)**2, -1/PHI ], [ 1/PHI, -(PHI)**2, -1/PHI, -1/PHI ], [ 1/PHI, -1/PHI, -(PHI)**2, 1/PHI ], [ -1/PHI, -(PHI)**2, 1/PHI, -1/PHI ],
    [ -(PHI)**2, -1/PHI, -1/PHI, -1/PHI ], [ 1/PHI, -1/PHI, PHI**2, -1/PHI ], [ -1/PHI, PHI**2, -1/PHI, -1/PHI ], [ -1/PHI, 1/PHI, PHI**2, -1/PHI ],
    [ -1/PHI, -1/PHI, PHI**2, 1/PHI ], [ -1/PHI, 1/PHI, -(PHI)**2, 1/PHI ], [ 1/PHI, -1/PHI, -(PHI)**2, -1/PHI ], [ -1/PHI, -(PHI)**2, -1/PHI, -1/PHI ],
    [ -1/PHI, 1/PHI, -(PHI)**2, -1/PHI ], [ -1/PHI, -1/PHI, -(PHI)**2, 1/PHI ], [ -1/PHI, -1/PHI, PHI**2, -1/PHI ], [ -1/PHI, -1/PHI, -(PHI)**2, -1/PHI ],
    [ 0, 1/PHI**2, 1, PHI**2 ], [ 1, 0, 1/PHI**2, PHI**2 ], [ 0, -1/PHI**2, 1, PHI**2 ], [ 0, PHI**2, 1/PHI**2, 1 ],
    [ 0, 1/PHI**2, -1, PHI**2 ], [ 0, 1/PHI**2, 1, -(PHI)**2 ], [ -1, 0, 1/PHI**2, PHI**2 ], [ 1/PHI**2, 1, 0, PHI**2 ],
    [ 1, PHI**2, 0, 1/PHI**2 ], [ 1, 0, -1/PHI**2, PHI**2 ], [ 1, 0, 1/PHI**2, -(PHI)**2 ], [ 0, PHI**2, -1/PHI**2, 1 ],
    [ 0, -1/PHI**2, -1, PHI**2 ], [ 0, -1/PHI**2, 1, -(PHI)**2 ], [ 1/PHI**2, 0, PHI**2, 1 ], [ 0, -(PHI)**2, 1/PHI**2, 1 ],
    [ 0, 1, PHI**2, 1/PHI**2 ], [ 0, PHI**2, 1/PHI**2, -1 ], [ 0, 1/PHI**2, -1, -(PHI)**2 ], [ 1/PHI**2, -1, 0, PHI**2 ],
    [ -1, PHI**2, 0, 1/PHI**2 ], [ -1, 0, -1/PHI**2, PHI**2 ], [ -1, 0, 1/PHI**2, -(PHI)**2 ], [ -1/PHI**2, 1, 0, PHI**2 ],
    [ 1/PHI**2, PHI**2, 1, 0 ], [ 1/PHI**2, 1, 0, -(PHI)**2 ], [ 1, -(PHI)**2, 0, 1/PHI**2 ], [ 1, 1/PHI**2, PHI**2, 0 ],
    [ 1, PHI**2, 0, -1/PHI**2 ], [ 1, 0, -1/PHI**2, -(PHI)**2 ], [ -1/PHI**2, 0, PHI**2, 1 ], [ 0, -(PHI)**2, -1/PHI**2, 1 ],
    [ 0, 1, PHI**2, -1/PHI**2 ], [ 0, PHI**2, -1/PHI**2, -1 ], [ 0, -1/PHI**2, -1, -(PHI)**2 ], [ PHI**2, 1/PHI**2, 0, 1 ],
    [ 1/PHI**2, 0, -(PHI)**2, 1 ], [ 1/PHI**2, 0, PHI**2, -1 ], [ 0, 1, -(PHI)**2, 1/PHI**2 ], [ 0, -(PHI)**2, 1/PHI**2, -1 ],
    [ PHI**2, 0, 1, 1/PHI**2 ], [ 0, -1, PHI**2, 1/PHI**2 ], [ -1/PHI**2, -1, 0, PHI**2 ], [ 1/PHI**2, PHI**2, -1, 0 ],
    [ 1/PHI**2, -1, 0, -(PHI)**2 ], [ -1, -(PHI)**2, 0, 1/PHI**2 ], [ -1, 1/PHI**2, PHI**2, 0 ], [ -1, PHI**2, 0, -1/PHI**2 ],
    [ -1, 0, -1/PHI**2, -(PHI)**2 ], [ -1/PHI**2, PHI**2, 1, 0 ], [ -1/PHI**2, 1, 0, -(PHI)**2 ], [ 1/PHI**2, -(PHI)**2, 1, 0 ],
    [ 1, 1/PHI**2, -(PHI)**2, 0 ], [ 1, -(PHI)**2, 0, -1/PHI**2 ], [ PHI**2, 1, 1/PHI**2, 0 ], [ 1, -1/PHI**2, PHI**2, 0 ],
    [ PHI**2, -1/PHI**2, 0, 1 ], [ -1/PHI**2, 0, -(PHI)**2, 1 ], [ -1/PHI**2, 0, PHI**2, -1 ], [ 0, 1, -(PHI)**2, -1/PHI**2 ],
    [ 0, -(PHI)**2, -1/PHI**2, -1 ], [ PHI**2, 0, 1, -1/PHI**2 ], [ 0, -1, PHI**2, -1/PHI**2 ], [ -(PHI)**2, 1/PHI**2, 0, 1 ],
    [ PHI**2, 1/PHI**2, 0, -1 ], [ 1/PHI**2, 0, -(PHI)**2, -1 ], [ -(PHI)**2, 0, 1, 1/PHI**2 ], [ 0, -1, -(PHI)**2, 1/PHI**2 ],
    [ PHI**2, 0, -1, 1/PHI**2 ], [ -1/PHI**2, PHI**2, -1, 0 ], [ -1/PHI**2, -1, 0, -(PHI)**2 ], [ 1/PHI**2, -(PHI)**2, -1, 0 ],
    [ -1, 1/PHI**2, -(PHI)**2, 0 ], [ -1, -(PHI)**2, 0, -1/PHI**2 ], [ PHI**2, -1, 1/PHI**2, 0 ], [ -1, -1/PHI**2, PHI**2, 0 ],
    [ -1/PHI**2, -(PHI)**2, 1, 0 ], [ -(PHI)**2, 1, 1/PHI**2, 0 ], [ 1, -1/PHI**2, -(PHI)**2, 0 ], [ PHI**2, 1, -1/PHI**2, 0 ],
    [ -(PHI)**2, -1/PHI**2, 0, 1 ], [ PHI**2, -1/PHI**2, 0, -1 ], [ -1/PHI**2, 0, -(PHI)**2, -1 ], [ -(PHI)**2, 0, 1, -1/PHI**2 ],
    [ 0, -1, -(PHI)**2, -1/PHI**2 ], [ PHI**2, 0, -1, -1/PHI**2 ], [ -(PHI)**2, 1/PHI**2, 0, -1 ], [ -(PHI)**2, 0, -1, 1/PHI**2 ],
    [ -1/PHI**2, -(PHI)**2, -1, 0 ], [ -(PHI)**2, -1, 1/PHI**2, 0 ], [ -1, -1/PHI**2, -(PHI)**2, 0 ], [ PHI**2, -1, -1/PHI**2, 0 ],
    [ -(PHI)**2, 1, -1/PHI**2, 0 ], [ -(PHI)**2, -1/PHI**2, 0, -1 ], [ -(PHI)**2, 0, -1, -1/PHI**2 ], [ -(PHI)**2, -1, -1/PHI**2, 0 ],
    [ 0, 1/PHI, PHI, SQRT5 ], [ PHI, 0, 1/PHI, SQRT5 ], [ 0, -1/PHI, PHI, SQRT5 ], [ 0, SQRT5, 1/PHI, PHI ],
    [ 0, 1/PHI, -(PHI), SQRT5 ], [ 0, 1/PHI, PHI, -(SQRT5) ], [ -(PHI), 0, 1/PHI, SQRT5 ], [ 1/PHI, PHI, 0, SQRT5 ],
    [ PHI, SQRT5, 0, 1/PHI ], [ PHI, 0, -1/PHI, SQRT5 ], [ PHI, 0, 1/PHI, -(SQRT5) ], [ 0, SQRT5, -1/PHI, PHI ],
    [ 0, -1/PHI, -(PHI), SQRT5 ], [ 0, -1/PHI, PHI, -(SQRT5) ], [ 1/PHI, 0, SQRT5, PHI ], [ 0, -(SQRT5), 1/PHI, PHI ],
    [ 0, PHI, SQRT5, 1/PHI ], [ 0, SQRT5, 1/PHI, -(PHI) ], [ 0, 1/PHI, -(PHI), -(SQRT5) ], [ 1/PHI, -(PHI), 0, SQRT5 ],
    [ -(PHI), SQRT5, 0, 1/PHI ], [ -(PHI), 0, -1/PHI, SQRT5 ], [ -(PHI), 0, 1/PHI, -(SQRT5) ], [ -1/PHI, PHI, 0, SQRT5 ],
    [ 1/PHI, SQRT5, PHI, 0 ], [ 1/PHI, PHI, 0, -(SQRT5) ], [ PHI, -(SQRT5), 0, 1/PHI ], [ PHI, 1/PHI, SQRT5, 0 ],
    [ PHI, SQRT5, 0, -1/PHI ], [ PHI, 0, -1/PHI, -(SQRT5) ], [ -1/PHI, 0, SQRT5, PHI ], [ 0, -(SQRT5), -1/PHI, PHI ],
    [ 0, PHI, SQRT5, -1/PHI ], [ 0, SQRT5, -1/PHI, -(PHI) ], [ 0, -1/PHI, -(PHI), -(SQRT5) ], [ SQRT5, 1/PHI, 0, PHI ],
    [ 1/PHI, 0, -(SQRT5), PHI ], [ 1/PHI, 0, SQRT5, -(PHI) ], [ 0, PHI, -(SQRT5), 1/PHI ], [ 0, -(SQRT5), 1/PHI, -(PHI) ],
    [ SQRT5, 0, PHI, 1/PHI ], [ 0, -(PHI), SQRT5, 1/PHI ], [ -1/PHI, -(PHI), 0, SQRT5 ], [ 1/PHI, SQRT5, -(PHI), 0 ],
    [ 1/PHI, -(PHI), 0, -(SQRT5) ], [ -(PHI), -(SQRT5), 0, 1/PHI ], [ -(PHI), 1/PHI, SQRT5, 0 ], [ -(PHI), SQRT5, 0, -1/PHI ],
    [ -(PHI), 0, -1/PHI, -(SQRT5) ], [ -1/PHI, SQRT5, PHI, 0 ], [ -1/PHI, PHI, 0, -(SQRT5) ], [ 1/PHI, -(SQRT5), PHI, 0 ],
    [ PHI, 1/PHI, -(SQRT5), 0 ], [ PHI, -(SQRT5), 0, -1/PHI ], [ SQRT5, PHI, 1/PHI, 0 ], [ PHI, -1/PHI, SQRT5, 0 ],
    [ SQRT5, -1/PHI, 0, PHI ], [ -1/PHI, 0, -(SQRT5), PHI ], [ -1/PHI, 0, SQRT5, -(PHI) ], [ 0, PHI, -(SQRT5), -1/PHI ],
    [ 0, -(SQRT5), -1/PHI, -(PHI) ], [ SQRT5, 0, PHI, -1/PHI ], [ 0, -(PHI), SQRT5, -1/PHI ], [ -(SQRT5), 1/PHI, 0, PHI ],
    [ SQRT5, 1/PHI, 0, -(PHI) ], [ 1/PHI, 0, -(SQRT5), -(PHI) ], [ -(SQRT5), 0, PHI, 1/PHI ], [ 0, -(PHI), -(SQRT5), 1/PHI ],
    [ SQRT5, 0, -(PHI), 1/PHI ], [ -1/PHI, SQRT5, -(PHI), 0 ], [ -1/PHI, -(PHI), 0, -(SQRT5) ], [ 1/PHI, -(SQRT5), -(PHI), 0 ],
    [ -(PHI), 1/PHI, -(SQRT5), 0 ], [ -(PHI), -(SQRT5), 0, -1/PHI ], [ SQRT5, -(PHI), 1/PHI, 0 ], [ -(PHI), -1/PHI, SQRT5, 0 ],
    [ -1/PHI, -(SQRT5), PHI, 0 ], [ -(SQRT5), PHI, 1/PHI, 0 ], [ PHI, -1/PHI, -(SQRT5), 0 ], [ SQRT5, PHI, -1/PHI, 0 ],
    [ -(SQRT5), -1/PHI, 0, PHI ], [ SQRT5, -1/PHI, 0, -(PHI) ], [ -1/PHI, 0, -(SQRT5), -(PHI) ], [ -(SQRT5), 0, PHI, -1/PHI ],
    [ 0, -(PHI), -(SQRT5), -1/PHI ], [ SQRT5, 0, -(PHI), -1/PHI ], [ -(SQRT5), 1/PHI, 0, -(PHI) ], [ -(SQRT5), 0, -(PHI), 1/PHI ],
    [ -1/PHI, -(SQRT5), -(PHI), 0 ], [ -(SQRT5), -(PHI), 1/PHI, 0 ], [ -(PHI), -1/PHI, -(SQRT5), 0 ], [ SQRT5, -(PHI), -1/PHI, 0 ],
    [ -(SQRT5), PHI, -1/PHI, 0 ], [ -(SQRT5), -1/PHI, 0, -(PHI) ], [ -(SQRT5), 0, -(PHI), -1/PHI ], [ -(SQRT5), -(PHI), -1/PHI, 0 ],
    [ 1/PHI, 1, PHI, 2 ], [ -1/PHI, 1, PHI, 2 ], [ PHI, 1/PHI, 1, 2 ], [ 1/PHI, -1, PHI, 2 ],
    [ 1/PHI, 2, 1, PHI ], [ 1/PHI, 1, -(PHI), 2 ], [ 1/PHI, 1, PHI, -2 ], [ PHI, -1/PHI, 1, 2 ],
    [ -1/PHI, -1, PHI, 2 ], [ -1/PHI, 2, 1, PHI ], [ -1/PHI, 1, -(PHI), 2 ], [ -1/PHI, 1, PHI, -2 ],
    [ -(PHI), 1/PHI, 1, 2 ], [ 1, PHI, 1/PHI, 2 ], [ PHI, 2, 1/PHI, 1 ], [ PHI, 1/PHI, -1, 2 ],
    [ PHI, 1/PHI, 1, -2 ], [ 1/PHI, 2, -1, PHI ], [ 1/PHI, -1, -(PHI), 2 ], [ 1/PHI, -1, PHI, -2 ],
    [ 1, 1/PHI, 2, PHI ], [ 1/PHI, -2, 1, PHI ], [ 1/PHI, PHI, 2, 1 ], [ 1/PHI, 2, 1, -(PHI) ],
    [ 1/PHI, 1, -(PHI), -2 ], [ -(PHI), -1/PHI, 1, 2 ], [ 1, PHI, -1/PHI, 2 ], [ PHI, 2, -1/PHI, 1 ],
    [ PHI, -1/PHI, -1, 2 ], [ PHI, -1/PHI, 1, -2 ], [ -1/PHI, 2, -1, PHI ], [ -1/PHI, -1, -(PHI), 2 ],
    [ -1/PHI, -1, PHI, -2 ], [ 1, -1/PHI, 2, PHI ], [ -1/PHI, -2, 1, PHI ], [ -1/PHI, PHI, 2, 1 ],
    [ -1/PHI, 2, 1, -(PHI) ], [ -1/PHI, 1, -(PHI), -2 ], [ 1, -(PHI), 1/PHI, 2 ], [ -(PHI), 2, 1/PHI, 1 ],
    [ -(PHI), 1/PHI, -1, 2 ], [ -(PHI), 1/PHI, 1, -2 ], [ -1, PHI, 1/PHI, 2 ], [ 1, 2, PHI, 1/PHI ],
    [ 1, PHI, 1/PHI, -2 ], [ PHI, -2, 1/PHI, 1 ], [ PHI, 1, 2, 1/PHI ], [ PHI, 2, 1/PHI, -1 ],
    [ PHI, 1/PHI, -1, -2 ], [ -1, 1/PHI, 2, PHI ], [ 1/PHI, -2, -1, PHI ], [ 1/PHI, PHI, 2, -1 ],
    [ 1/PHI, 2, -1, -(PHI) ], [ 1/PHI, -1, -(PHI), -2 ], [ 2, 1, 1/PHI, PHI ], [ 1, 1/PHI, -2, PHI ],
    [ 1, 1/PHI, 2, -(PHI) ], [ 1/PHI, PHI, -2, 1 ], [ 1/PHI, -2, 1, -(PHI) ], [ 2, 1/PHI, PHI, 1 ],
    [ 1/PHI, -(PHI), 2, 1 ], [ 1, -(PHI), -1/PHI, 2 ], [ -(PHI), 2, -1/PHI, 1 ], [ -(PHI), -1/PHI, -1, 2 ],
    [ -(PHI), -1/PHI, 1, -2 ], [ -1, PHI, -1/PHI, 2 ], [ 1, 2, PHI, -1/PHI ], [ 1, PHI, -1/PHI, -2 ],
    [ PHI, -2, -1/PHI, 1 ], [ PHI, 1, 2, -1/PHI ], [ PHI, 2, -1/PHI, -1 ], [ PHI, -1/PHI, -1, -2 ],
    [ -1, -1/PHI, 2, PHI ], [ -1/PHI, -2, -1, PHI ], [ -1/PHI, PHI, 2, -1 ], [ -1/PHI, 2, -1, -(PHI) ],
    [ -1/PHI, -1, -(PHI), -2 ], [ 2, 1, -1/PHI, PHI ], [ 1, -1/PHI, -2, PHI ], [ 1, -1/PHI, 2, -(PHI) ],
    [ -1/PHI, PHI, -2, 1 ], [ -1/PHI, -2, 1, -(PHI) ], [ 2, -1/PHI, PHI, 1 ], [ -1/PHI, -(PHI), 2, 1 ],
    [ -1, -(PHI), 1/PHI, 2 ], [ 1, 2, -(PHI), 1/PHI ], [ 1, -(PHI), 1/PHI, -2 ], [ -(PHI), -2, 1/PHI, 1 ],
    [ -(PHI), 1, 2, 1/PHI ], [ -(PHI), 2, 1/PHI, -1 ], [ -(PHI), 1/PHI, -1, -2 ], [ -1, 2, PHI, 1/PHI ],
    [ -1, PHI, 1/PHI, -2 ], [ 1, -2, PHI, 1/PHI ], [ PHI, 1, -2, 1/PHI ], [ PHI, -2, 1/PHI, -1 ],
    [ 2, PHI, 1, 1/PHI ], [ PHI, -1, 2, 1/PHI ], [ 2, -1, 1/PHI, PHI ], [ -1, 1/PHI, -2, PHI ],
    [ -1, 1/PHI, 2, -(PHI) ], [ 1/PHI, PHI, -2, -1 ], [ 1/PHI, -2, -1, -(PHI) ], [ 2, 1/PHI, PHI, -1 ],
    [ 1/PHI, -(PHI), 2, -1 ], [ -2, 1, 1/PHI, PHI ], [ 2, 1, 1/PHI, -(PHI) ], [ 1, 1/PHI, -2, -(PHI) ],
    [ -2, 1/PHI, PHI, 1 ], [ 1/PHI, -(PHI), -2, 1 ], [ 2, 1/PHI, -(PHI), 1 ], [ -1, -(PHI), -1/PHI, 2 ],
    [ 1, 2, -(PHI), -1/PHI ], [ 1, -(PHI), -1/PHI, -2 ], [ -(PHI), -2, -1/PHI, 1 ], [ -(PHI), 1, 2, -1/PHI ],
    [ -(PHI), 2, -1/PHI, -1 ], [ -(PHI), -1/PHI, -1, -2 ], [ -1, 2, PHI, -1/PHI ], [ -1, PHI, -1/PHI, -2 ],
    [ 1, -2, PHI, -1/PHI ], [ PHI, 1, -2, -1/PHI ], [ PHI, -2, -1/PHI, -1 ], [ 2, PHI, 1, -1/PHI ],
    [ PHI, -1, 2, -1/PHI ], [ 2, -1, -1/PHI, PHI ], [ -1, -1/PHI, -2, PHI ], [ -1, -1/PHI, 2, -(PHI) ],
    [ -1/PHI, PHI, -2, -1 ], [ -1/PHI, -2, -1, -(PHI) ], [ 2, -1/PHI, PHI, -1 ], [ -1/PHI, -(PHI), 2, -1 ],
    [ -2, 1, -1/PHI, PHI ], [ 2, 1, -1/PHI, -(PHI) ], [ 1, -1/PHI, -2, -(PHI) ], [ -2, -1/PHI, PHI, 1 ],
    [ -1/PHI, -(PHI), -2, 1 ], [ 2, -1/PHI, -(PHI), 1 ], [ -1, 2, -(PHI), 1/PHI ], [ -1, -(PHI), 1/PHI, -2 ],
    [ 1, -2, -(PHI), 1/PHI ], [ -(PHI), 1, -2, 1/PHI ], [ -(PHI), -2, 1/PHI, -1 ], [ 2, -(PHI), 1, 1/PHI ],
    [ -(PHI), -1, 2, 1/PHI ], [ -1, -2, PHI, 1/PHI ], [ -2, PHI, 1, 1/PHI ], [ PHI, -1, -2, 1/PHI ],
    [ 2, PHI, -1, 1/PHI ], [ -2, -1, 1/PHI, PHI ], [ 2, -1, 1/PHI, -(PHI) ], [ -1, 1/PHI, -2, -(PHI) ],
    [ -2, 1/PHI, PHI, -1 ], [ 1/PHI, -(PHI), -2, -1 ], [ 2, 1/PHI, -(PHI), -1 ], [ -2, 1, 1/PHI, -(PHI) ],
    [ -2, 1/PHI, -(PHI), 1 ], [ -1, 2, -(PHI), -1/PHI ], [ -1, -(PHI), -1/PHI, -2 ], [ 1, -2, -(PHI), -1/PHI ],
    [ -(PHI), 1, -2, -1/PHI ], [ -(PHI), -2, -1/PHI, -1 ], [ 2, -(PHI), 1, -1/PHI ], [ -(PHI), -1, 2, -1/PHI ],
    [ -1, -2, PHI, -1/PHI ], [ -2, PHI, 1, -1/PHI ], [ PHI, -1, -2, -1/PHI ], [ 2, PHI, -1, -1/PHI ],
    [ -2, -1, -1/PHI, PHI ], [ 2, -1, -1/PHI, -(PHI) ], [ -1, -1/PHI, -2, -(PHI) ], [ -2, -1/PHI, PHI, -1 ],
    [ -1/PHI, -(PHI), -2, -1 ], [ 2, -1/PHI, -(PHI), -1 ], [ -2, 1, -1/PHI, -(PHI) ], [ -2, -1/PHI, -(PHI), 1 ],
    [ -1, -2, -(PHI), 1/PHI ], [ -2, -(PHI), 1, 1/PHI ], [ -(PHI), -1, -2, 1/PHI ], [ 2, -(PHI), -1, 1/PHI ],
    [ -2, PHI, -1, 1/PHI ], [ -2, -1, 1/PHI, -(PHI) ], [ -2, 1/PHI, -(PHI), -1 ], [ -1, -2, -(PHI), -1/PHI ],
    [ -2, -(PHI), 1, -1/PHI ], [ -(PHI), -1, -2, -1/PHI ], [ 2, -(PHI), -1, -1/PHI ], [ -2, PHI, -1, -1/PHI ],
    [ -2, -1, -1/PHI, -(PHI) ], [ -2, -1/PHI, -(PHI), -1 ], [ -2, -(PHI), -1, 1/PHI ], [ -2, -(PHI), -1, -1/PHI ],
);
my $edgelensq_120cell = 0.583592135001262;

# Which regular polytope to draw (replace @vertices_foobar and
# $edgelensq_foobar for some "foobar" in "16cell", "tesseract",
# "24cell", "600cell" or "120cell").
my @vertices = @vertices_120cell;
my $edgelensq = $edgelensq_120cell;

# Number of frames to compute
my $nbframes = 1200;

# Whether to use the hack to display lines of the "trans" hemisphere.
my $use_trans_hack = 1;

# (No user-serviceable parts below.)

# Compute the list of edges.
my @edges;
for ( my $i=0 ; $i<scalar(@vertices) ; $i++ ) {
    for ( my $j=$i+1 ; $j<scalar(@vertices) ; $j++ ) {
	my $vi = $vertices[$i];
	my $vj = $vertices[$j];
	my $dsq = (($$vi[0]-$$vj[0])*($$vi[0]-$$vj[0])
		   + ($$vi[1]-$$vj[1])*($$vi[1]-$$vj[1])
		   + ($$vi[2]-$$vj[2])*($$vi[2]-$$vj[2])
		   + ($$vi[3]-$$vj[3])*($$vi[3]-$$vj[3]));
	push @edges, [$i, $j] if ( abs($dsq - $edgelensq) < 1.e-9 );
    }
}
printf STDERR "%d vertices, %d edges\n", scalar(@vertices), scalar(@edges);

# A generator of gaussian numbers
my $gaussian_store;
sub gaussian {
    my $val;
    if ( defined($val = $gaussian_store) ) {
	$gaussian_store = undef;
	return $val;
    }
    my $u0 = rand;
    my $u1 = rand;
    $val = sqrt(-2*log($u0)) * cos(2*PI*$u1);
    $gaussian_store = sqrt(-2*log($u0)) * sin(2*PI*$u1);
    return $val;
}

# Perform a Gram-Schmidt orthonormalization.
sub gram_schmidt {
    my $b = shift;
    for ( my $i=0 ; $i<4 ; $i++ ) {
	for ( my $j=0 ; $j<$i ; $j++ ) {
	    my $dot = 0;
	    for ( my $k=0 ; $k<4 ; $k++ ) {
		$dot += $$b[$i][$k] * $$b[$j][$k];
	    }
	    for ( my $k=0 ; $k<4 ; $k++ ) {
		$$b[$i][$k] -= $dot * $$b[$j][$k];
	    }
	}
	my $norm = 0;
	for ( my $k=0 ; $k<4 ; $k++ ) {
	    $norm += $$b[$i][$k] * $$b[$i][$k];
	}
	$norm = sqrt($norm);
	for ( my $k=0 ; $k<4 ; $k++ ) {
	    $$b[$i][$k] /= $norm;
	}
    }
}

# The initial orthonormal matrix
my @orthog0;
@orthog0 = ([gaussian, gaussian, gaussian, gaussian],
	    [gaussian, gaussian, gaussian, gaussian],
	    [gaussian, gaussian, gaussian, gaussian],
	    [gaussian, gaussian, gaussian, gaussian]);
gram_schmidt \@orthog0;

# The main loop, over all frames:
for ( my $fnum=0 ; $fnum<$nbframes ; $fnum++ ) {

    my @orthog;
    my $ctheta = cos(2*$fnum*PI/$nbframes);
    my $stheta = sin(2*$fnum*PI/$nbframes);
    $orthog[0] = [ $ctheta*$orthog0[0][0] - $stheta*$orthog0[3][0],
		   $ctheta*$orthog0[0][1] - $stheta*$orthog0[3][1],
		   $ctheta*$orthog0[0][2] - $stheta*$orthog0[3][2],
		   $ctheta*$orthog0[0][3] - $stheta*$orthog0[3][3] ];
    $orthog[1] = $orthog0[1];
    $orthog[2] = $orthog0[2];
    $orthog[3] = [ $ctheta*$orthog0[3][0] + $stheta*$orthog0[0][0],
		   $ctheta*$orthog0[3][1] + $stheta*$orthog0[0][1],
		   $ctheta*$orthog0[3][2] + $stheta*$orthog0[0][2],
		   $ctheta*$orthog0[3][3] + $stheta*$orthog0[0][3] ];

    my @rotated = ();
    for ( my $n=0 ; $n<scalar(@vertices) ; $n++ ) {
	my @l = ();
	for ( my $i=0 ; $i<4 ; $i++ ) {
	    my $v = 0;
	    for ( my $j=0 ; $j<4 ; $j++ ) {
		$v += $orthog[$i][$j] * $vertices[$n][$j];
	    }
	    $l[$i] = $v;
	}
	push @rotated, \@l;
    }

    open my $f, ">", sprintf("frame%04d.pov", $fnum)
	or die "can't open output file";

    # Print the orthogonal projection matrix as a comment, so it can
    # be reused (to reuse it, set @orthog0 above to the matrix found
    # in frame0000.pov).
    print $f "// Orthogonal projection matrix used for this frame:\n";
    for ( my $i=0 ; $i<4 ; $i++ ) {
	print $f "// orth[$i] = (";
	for ( my $j=0 ; $j<4 ; $j++ ) {
	    printf $f "%+11.9f", $orthog[$i][$j];
	    print $f "," if $j<3;
	}
	print $f ")";
	print $f "  // <-- current point" if $i==3;
	print $f "  // <-- direction of view" if $i==0;
	print $f "  // <-- right (?)" if $i==1;
	print $f "  // <-- up (?)" if $i==2;
	print $f "\n";
    }

    print $f <<'__EOF__';
camera { location <0,0,0> look_at <1,0,0> up <0,0,1> right <0,1.5,0> }
light_source { <0,0,0>, 1 }
#declare Std = texture { pigment { color <1,1,1> } finish { ambient .25 diffuse .4 phong .35 reflection .1 } }
#declare Far = texture { pigment { color <1,0.5,0.5> } finish { ambient .4 diffuse .0 } }
__EOF__

    foreach my $e ( @edges ) {
	my $n0 = $$e[0];
	my $n1 = $$e[1];
	my $v0 = $rotated[$n0];
	my $v1 = $rotated[$n1];
	my @p0 = ($$v0[0]/$$v0[3], $$v0[1]/$$v0[3], $$v0[2]/$$v0[3]);
	my @p1 = ($$v1[0]/$$v1[3], $$v1[1]/$$v1[3], $$v1[2]/$$v1[3]);
	my @dir = ($p1[0]-$p0[0], $p1[1]-$p0[1], $p1[2]-$p0[2]);
	my $dirnorm = sqrt($dir[0]*$dir[0] + $dir[1]*$dir[1] + $dir[2]*$dir[2]);
	@dir = ($dir[0]/$dirnorm, $dir[1]/$dirnorm, $dir[2]/$dirnorm);
	if ( $use_trans_hack ) {
	    # This ugly hack will display the lines of the "far"
	    # ("trans") hemisphere in the 3-sphere, i.e., the one away
	    # from the observer, as untextured red lines in the
	    # background.
	    my @q0 = (1, $$v0[1]/$$v0[0], $$v0[2]/$$v0[0]);
	    my @q1 = (1, $$v1[1]/$$v1[0], $$v1[2]/$$v1[0]);
	    if ( ($$v0[0]<0)^($$v1[0]<0) ) {
		my @qdir = (1, $q1[1]-$q0[1], $q1[2]-$q0[2]);
		my $qdirnorm = sqrt($qdir[1]*$qdir[1] + $qdir[2]*$qdir[2]);
		@qdir = ($qdir[0]/$qdirnorm, $qdir[1]/$qdirnorm, $qdir[2]/$qdirnorm);
		printf $f "cylinder { <10000,%.5f,%.5f>, <10000,%.5f,%.5f>, 20 texture { Far } }\n", 1.e4*$q0[1], 1.e4*$q0[2], 1.e4*($q0[1]-10*$qdir[1]), 1.e4*($q0[2]-10*$qdir[2]);
		printf $f "cylinder { <10000,%.5f,%.5f>, <10000,%.5f,%.5f>, 20 texture { Far } }\n", 1.e4*$q1[1], 1.e4*$q1[2], 1.e4*($q1[1]+10*$qdir[1]), 1.e4*($q1[2]+10*$qdir[2]);
	    } else {
		printf $f "cylinder { <10000,%.5f,%.5f>, <10000,%.5f,%.5f>, 20 texture { Far } }\n", 1.e4*$q0[1], 1.e4*$q0[2], 1.e4*$q1[1], 1.e4*$q1[2];
	    }
	}
	next if $$v0[3] < 0 && $$v1[3] < 0;  # Line is entirely "trans".
	if ( $$v0[3] < 0 ) {  # The v0 edge is "trans" whereas v1 is "cis".
	    @p1 = ($p0[0]-1.e4*$dir[0], $p0[1]-1.e4*$dir[1], $p0[2]-1.e4*$dir[2]);
	} elsif ( $$v1[3] < 0 ) {  # The v1 edge is "trans" whereas v0 is "cis".
	    @p0 = ($p1[0]+1.e4*$dir[0], $p1[1]+1.e4*$dir[1], $p1[2]+1.e4*$dir[2]);
	}
	# Draw the actual "cis" beam.
	printf $f "cylinder { <%.5f,%.5f,%.5f>, <%.5f,%.5f,%.5f>, 0.01 texture { Std } }\n", $p0[0], $p0[1], $p0[2], $p1[0], $p1[1], $p1[2];
	my @cross = ( $p1[1]*$dir[2]-$p1[2]*$dir[1],
		      $p1[2]*$dir[0]-$p1[0]*$dir[2],
		      $p1[0]*$dir[1]-$p1[1]*$dir[0] );
	printf STDERR "warning: camera crosses beam at frame %d\n", $fnum
	    if ( $cross[0]*$cross[0] + $cross[1]*$cross[1] + $cross[2]*$cross[2] < 1.2e-4
		 && ( ( $p0[0]*$dir[0]+$p0[1]*$dir[1]+$p0[2]*$dir[2]<0 )
		      ^ ( $p1[0]*$dir[0]+$p1[1]*$dir[1]+$p1[2]*$dir[2]<0 ) ) );
    }

    close $f;
}