#!/usr/bin/perl

use Tcl::Tk;
$interp = new Tcl::Tk;

#  Set up the main window with TCL
$interp->Eval(qq{
    wm title . "N Points on Sphere"
    wm geometry . "640x480"
    frame .frame -width 160 -height 480
    pack .frame -expand 1 -side right -fill y
    label .caption -relief sunken -text " \n " -width 2048 -justify left -anchor w -font {Helvetica 12}
    pack .caption -expand 1 -side bottom -fill x
    canvas .c -width 2048 -height 2048 -background AntiqueWhite -relief sunken
    pack .c -expand 1 -fill both -side left
    label .frame.npointsLabel -text "Number of points"
    spinbox .frame.npoints -from 2 -to 100 -width 5 -validate key -vcmd {string is integer %P}
    .frame.npoints set 20
    pack .frame.npointsLabel .frame.npoints -side top -anchor w
    label .frame.initLabel -text "Start positions"
    set initpos 0
    radiobutton .frame.initpos1 -text "Spiral" -variable initpos -value 0
    radiobutton .frame.initpos2 -text "Random" -variable initpos -value 1
    pack .frame.initLabel -side top -pady {16 0} -anchor w
    pack .frame.initpos1 .frame.initpos2 -side top -anchor w
    label .frame.cycleLabel -text "Max cycle"
    set maxCycle 1000
    entry .frame.cycle -textvariable maxCycle -width 6 -validate key -vcmd {string is integer %P}
    pack .frame.cycleLabel -side top -pady {16 0} -anchor w
    pack .frame.cycle -side top -anchor w
    button .frame.run -text "Run" -width 10
    button .frame.stop -text "Stop" -width 10 -state disabled 
    button .frame.quit -text "Quit" -width 10
    pack .frame.run -side top -pady {40 4}
    pack .frame.stop -side top -pady 4
    pack .frame.quit -side bottom -pady 20
});
$canvas = $interp->widget('.c');
$caption = $interp->widget('.caption');

#  Bind the perl callback procedure to Tk events
#  Note: Tcl::Ev() must appear first in the argument list
$canvas->bind('<1>', [\&mouseDown, Tcl::Ev('%x', '%y'), $canvas]);
$canvas->bind('<B1-Motion>', [\&mouseMove, Tcl::Ev('%x', '%y'), $canvas]);
$canvas->bind('<ButtonRelease-1>', [\&mouseUp, Tcl::Ev('%x', '%y'), $canvas]);
$canvas->bind('<Configure>', [\&canvasConfigure, Tcl::Ev('%w', '%h'), $canvas]);
$interp->widget('.')->bind('<Destroy>', [\&windowDestroyed]);
$interp->widget('.frame.run')->configure(-command=>[\&runSimulation]);
$interp->widget('.frame.stop')->configure(-command=>[\&stopSimulation]);
$interp->widget('.frame.quit')->configure(-command=>[\&quit]);

$npoints = 0;

$depth = 8.0;
$height = 2.4;
$screenWidth = $screenHeight = $screenSize = 480;
$maxIter = 1000;

$xstart = 0;
$ystart = 0;
$quat = [0, 0, 0, 1];
$tempQuat = [0, 0, 0, 1];

$interp->MainLoop;

exit(0);

sub generalizedSpiralSet {
    #  Generate $npoints points uniformly on a unit sphere
    my ($npoints) = @_;
    my @col;
    my @lon;
    my ($sn, $i, $xi, $yi, $sc);
    my $results = [];
    $col[0] = 3.14159265358979;
    $lon[0] = 0.0;
    $col[$npoints - 1] = $lon[$npoints - 1] = 0.0;
    $sn = sqrt($npoints);
    for ($i = 1; $i < $npoints - 1; $i++) {
        $xi = -1.0 + 2.0 * $i / ($npoints - 1);
        $yi = sqrt(1.0 - $xi * $xi);
        $col[$i] = atan2($yi, $xi);
        $lon[$i] = $lon[$i - 1] + 3.6 / ($sn * $yi);
    }
    for ($i = 0; $i < $npoints; $i++) {
        $sc = sin($col[$i]);
        $results->[$i] = [$sc * cos($lon[$i]), $sc * sin($lon[$i]), -cos($col[$i])];
    }
    return $results;
}

sub randomSphereSet {
    my ($npoints) = @_;
    my $results = [];
    my $i;
    for ($i = 0; $i < $npoints; $i++) {
	$results->[$i] = randomUnitVector();
    }
    return $results;
}

sub randomUnitVector {
    my ($x, $y, $z, $w);
    $x = gaussianRand();
    $y = gaussianRand();
    $z = gaussianRand();
    $w = 1.0/sqrt($x*$x + $y*$y + $z*$z);
    return [$x*$w, $y*$w, $z*$w];
}

sub gaussianRand {
    if ($haveNextNextGaussian) {
	$haveNextNextGaussian = 0;
	return $nextNextGaussian;
    } else {
	my ($v1, $v2, $s, $mul);
	while ($s >= 1.0 || $s == 0.0) {
	    $v1 = 2 * rand() - 1;
	    $v2 = 2 * rand() - 1;
	    $s = $v1*$v1 + $v2*$v2;
	}
	$mul = sqrt(-2*log($s)/$s);
	$nextNextGaussian = $v2 * $mul;
	$haveNextNextGaussian = 1;
	return $v1 * $mul;
    }
}

sub uniformSphere {
    my ($points) = @_;
    my $npoints = $#{$points} + 1;
    my @force;
    my $r = 1.0;
    my $iter = 0;
    my $size = sqrt($npoints) * 5;
    my ($i, $j, $pi, $fi, $pj, $fx, $fy, $fz, $x, $y, $z, $dx, $dy, $dz, $w, $ww);
    my ($eng, $maxw, $maxf, $damp);
    my $converged = 0;
    for ($i = 0; $i < $npoints; $i++) {
	$force[$i] = [0, 0, 0];
    }
    while (1) {
	#  Clear the force
	for ($i = 0; $i < $npoints; $i++) {
	    $force[$i]->[0] = $force[$i]->[1] = $force[$i]->[2] = 0.0;
	}
	$eng = 0.0;
	$maxf = 0;
	for ($i = 0; $i < $npoints; $i++) {
	    $pi = $points->[$i];
	    for ($j = $i + 1; $j < $npoints; $j++) {
		$pj = $points->[$j];
		$fx = $pi->[0] - $pj->[0];
		$fy = $pi->[1] - $pj->[1];
		$fz = $pi->[2] - $pj->[2];
		$w = sqrt($fx*$fx + $fy*$fy + $fz*$fz);
		if ($w < 1e-3) {
		    $w = 1e-3;
		}
		$ww = $r / ($w ** 3);
		$eng = $r / $w;
		$fx *= $ww;
		$fy *= $ww;
		$fz *= $ww;
		$force[$i]->[0] += $fx;
		$force[$i]->[1] += $fy;
		$force[$i]->[2] += $fz;
		$force[$j]->[0] -= $fx;
		$force[$j]->[1] -= $fy;
		$force[$j]->[2] -= $fz;
	    }
	    $pi = $force[$i];
	    $w = $pi->[0]*$pi->[0] + $pi->[1]*$pi->[1] + $pi->[2]*$pi->[2];
	    if ($maxf < $w) {
		$maxf = $w;
	    }
	}
#    print "maxf = $maxf\n";
	if ($maxf > 10) {
	    $damp = 10 / $maxf;
	} else {
	    $damp = 1.0;
	}
	$maxw = 0.0;
	for ($i = 0; $i < $npoints; $i++) {
	    $pi = $points->[$i];
	    $fi = $force[$i];
	    #    printf "%d [%f,%f,%f][%f,%f,%f]\n", $i+1, @$pi, @$fi;
	    $x = $pi->[0] + $damp * $fi->[0];
	    $y = $pi->[1] + $damp * $fi->[1];
	    $z = $pi->[2] + $damp * $fi->[2];
	    $w = 1.0/sqrt($x*$x + $y*$y + $z*$z);
	    $x *= $w;
	    $y *= $w;
	    $z *= $w;
	    $dx = $pi->[0] - $x;
	    $dy = $pi->[1] - $y;
	    $dz = $pi->[2] - $z;
	    $w = sqrt($dx*$dx + $dy*$dy + $dz*$dz);
	    if ($w > $maxw) {
		$maxw = $w;
	    }
	    $pi->[0] = $x;
	    $pi->[1] = $y;
	    $pi->[2] = $z;
	}
	$iter++;
	if ($maxw < 1e-6) {
	    $converged = 1;
	}
	$caption->configure(-text=>sprintf("Cycle %d, repulsion energy %f, max force %f,\ndamp %f, max move %f%s", $iter, $eng, $maxf, $damp, $maxw, $converged ? " --- converged." : ""));
	updateDisplay();
	$interp->update();
	last if $stopFlag;
	last if $converged || $iter >= $maxIter;
    }
    return $eng;
}

sub updateDisplay {
    my ($create) = @_;
    my ($i, $k, $xx, $yy, $rr, $xs, $ys, $xe, $ye);
    my $conv = [];
    my $myQuat = quatMul($quat, $tempQuat);
    my $myQuatConj = quatConjugate($myQuat);
    my @z;
    if ($create) {
	$canvas->delete("all");
	undef @id, @label_id;
    }
    for ($i = 0; $i <= $#points; $i++) {
	quatMul($myQuatConj, quatMul($points[$i], $myQuat, $conv), $conv);
#    printf "conv($i) = [%f %f %f %f]\n", @$conv;
	$z[$i] = $conv->[2];
	$k = ($depth + $conv->[2]) / ($depth * $height) * $screenSize;
	$xx = $conv->[0] * $k + $screenSize / 2;
	$yy = $conv->[1] * $k + $screenSize / 2;
	$rr = 0.05 * $k;
	$xs = $xx - $rr;
	$xe = $xx + $rr;
	$ys = $yy - $rr;
	$ye = $yy + $rr;
	if ($create) {
	    $id[$i] = $canvas->createOval($xs, $ys, $xe, $ye, -fill=>'cyan');
	    $label_id[$i] = $canvas->createText($xx, $yy, -text=>sprintf("%d",$i+1));
	} else {
	    $canvas->coords($id[$i], $xs, $ys, $xe, $ye);
	    $canvas->coords($label_id[$i], $xx, $yy);
	}
    }
    my @order = sort { $z[$a] <=> $z[$b] } (0..$#points);
    foreach $i (@order) {
	$interp->Eval(".c raise $id[$i]");
	$interp->Eval(".c raise $label_id[$i]");
    }
}

sub vecLength {
    my ($a) = @_;
    return sqrt($a->[0]*$a->[0] + $a->[1]*$a->[1] + $a->[2]*$a->[2]);
}

sub vecDot {
    my ($a, $b) = @_;
    return $a->[0]*$b->[0] + $a->[1]*$b->[1] + $a->[2]*$b->[2];
}

sub vecCross {
    my ($a, $b, $c) = @_;
    $c = [] if !defined($c);
    @{$c}[0..2] = ($a->[1]*$b->[2]-$a->[2]*$b->[1], $a->[2]*$b->[0]-$a->[0]*$b->[2], $a->[0]*$b->[1]-$a->[1]*$b->[0]);
    return $c;
}

sub vecNormalize {
    my ($a, $b) = @_;
    my $w = 1.0 / vecLength($a);
    $b = [] if !defined($b);
    @{$b}[0..2] = ($a->[0]*$w, $a->[1]*$w, $a->[2]*$w);
    return $b;
}

#  Quaternion operations
#  Note: The expression of quaternion used here is unconventional; the imaginary 
#  elements are placed first, and the real element is placed last.
#  Hence the rotation of ANG along the axis (X, Y, Z) becomes
#  [sin(ANG/2)*X, sin(ANG/2)*Y, sin(ANG/2)*Z, cos(ANG/2)].
sub quatMul {
    my ($q1, $q2, $q3) = @_;
    $q3 = [] if !defined($q3);
    @{$q3}[0..3] = (
		    $q1->[3]*$q2->[0] + $q1->[0]*$q2->[3] + $q1->[1]*$q2->[2] - $q1->[2]*$q2->[1],
		    $q1->[3]*$q2->[1] + $q1->[1]*$q2->[3] + $q1->[2]*$q2->[0] - $q1->[0]*$q2->[2],
		    $q1->[3]*$q2->[2] + $q1->[2]*$q2->[3] + $q1->[0]*$q2->[1] - $q1->[1]*$q2->[0],
		    $q1->[3]*$q2->[3] - $q1->[0]*$q2->[0] - $q1->[1]*$q2->[1] - $q1->[2]*$q2->[2]
		    );
    return $q3;
}

sub quatConjugate {
    my ($q1, $q2) = @_;
    $q2 = [] if !defined($q2);
    @{$q2}[0..3] = (-$q1->[0], -$q1->[1], -$q1->[2], $q1->[3]);
    return $q2;
}

sub calcTrackballRotation {
    #  Handling trackball: calculate the rotation between the two
    # points (xs, ys) to (xe, ye) within the unit circle in the XY plane.
    my ($xs, $ys, $xe, $ye, $quat) = @_;
    my ($r, $r2);
    $r = $xs * $xs + $ys * $ys;
    if ($r > 1.0) {
	$zs = 0.0;
	$r = 1.0 / sqrt($r);
	$xs *= $r;
	$ys *= $r;
    } else {
	$zs = sqrt(1.0 - $r);
    }
    $r = $xe * $xe + $ye * $ye;
    if ($r > 1.0) {
	$ze = 0.0;
	$r = 1.0 / sqrt($r);
	$xe *= $r;
	$ye *= $r;
    } else {
	$ze = sqrt(1.0 - $r);
    }
    $vs = [$xs, $ys, $zs];
    $ve = [$xe, $ye, $ze];
    $cs = vecDot($vs, $ve);
    if ($cs < -0.9999) {
	#  180-degree rotation around (-y, x, 0)
	$r = 1.0/sqrt($xs*$xs + $ys*$ys);
	$quat->[0] = -$ys * $r;
	$quat->[1] = $xs * $r;
	$quat->[2] = 0.0;
	$quat->[3] = 0.0;
    } elsif ($cs > 0.9999) {
	#  No rotation
	$quat->[0] = $quat->[1] = $quat->[2] = 0.0;
	$quat->[3] = 1.0;
    } else {
	$vx = [0, 0, 0];
	vecNormalize(vecCross($ve, $vs, $vx), $vx);
	$sn = sqrt(0.5 - 0.5 * $cs);
	$cs = sqrt(0.5 + 0.5 * $cs);
	$quat->[0] = $vx->[0] * $sn;
	$quat->[1] = $vx->[1] * $sn;
	$quat->[2] = $vx->[2] * $sn;
	$quat->[3] = $cs;
    }
    return $quat;
}

sub runSimulation {
    $interp->update();
    $interp->Eval('.frame.run configure -state disable; .frame.stop configure -state normal');
    $stopFlag = 0;
    $npoints = $interp->widget('.frame.npoints')->get();
    if ($npoints < 2 || $npoints > 99999) {
	$npoints = ($npoints < 2 ? 2 : 99999);
	$interp->widget('.frame.npoints')->set($npoints);
    }
    $maxIter = $interp->widget('.frame.cycle')->get();
    if ($maxIter < 2 || $maxIter > 99999) {
	$maxIter = ($maxIter < 2 ? 2 : 99999);
	$interp->widget('.frame.cycle')->configure(-text=>$maxIter);
    }
    $initpos = $interp->Eval('expr {$initpos}');
    if ($initpos) {
	@points = @{randomSphereSet($npoints)};
    } else {
	@points = @{generalizedSpiralSet($npoints)};
    }
    updateDisplay(1);
    uniformSphere(\@points);
    if (!$windowDestroyed) {
	$interp->Eval('.frame.run configure -state normal; .frame.stop configure -state disable');
    }
}

sub stopSimulation {
    $stopFlag = 1;
}

sub quit {
    $interp->widget('.')->destroy();
}

sub windowDestroyed {
    $windowDestroyed = 1;
    $stopFlag = 1;
}

sub mouseDown {
    my ($x, $y, $widget) = @_;
    $xstart = $x;
    $ystart = $y;
}

sub mouseMove {
    my ($x, $y, $widget) = @_;
    my ($xs, $ys, $xe, $ye);
    $xs = ($xstart - $screenWidth * 0.5) / $screenSize;
    $xe = ($x - $screenWidth * 0.5) / $screenSize;
    $ys = ($ystart - $screenHeight * 0.5) / $screenSize;
    $ye = ($y - $screenHeight * 0.5) / $screenSize;
    calcTrackballRotation($xs, $ys, $xe, $ye, $tempQuat);
    updateDisplay();
}

sub mouseUp {
    quatMul($quat, $tempQuat, $quat);
    $tempQuat->[0] = $tempQuat->[1] = $tempQuat->[2] = 0.0;
    $tempQuat->[3] = 1.0;
    updateDisplay();
}

sub canvasConfigure {
    my ($width, $height, $widget) = @_;
    $screenWidth = $width;
    $screenHeight = $height;
    if ($width < $height) {
	$screenSize = $width;
    } else {
	$screenSize = $height;
    }
    updateDisplay();
}