#!/usr/bin/env wxPerl
use Wx;

# ============================================================
#  A class representing the main canvas to plot the points on
# ============================================================
package MyCanvas;
use base qw(Wx::Window);
use Wx::Event qw(EVT_PAINT EVT_LEFT_DOWN EVT_MOTION EVT_LEFT_UP);
sub new {
    my $ref = shift;
    my $self = $ref->SUPER::new(@_);
    #  Declare the events
    EVT_PAINT($self, \&OnPaint);
    EVT_LEFT_DOWN($self, \&OnMouseDown);
    EVT_MOTION($self, \&OnMouseMove);
    EVT_LEFT_UP($self, \&OnMouseUp);
    return $self;
}

sub GetParentFrame {
    my $self = shift;
    my ($win, $win1);
    $win = $self;
    while (defined($win1 = $win->GetParent())) {
        $win = $win1;
    }
    return $win;
}

sub OnPaint {
    my $self = shift;
    #  Creating a PaintDC is required in every OnPaint handler
    #  ($dc will be automatically DESTROY'ed when it becomes out of scope)
    my $dc = Wx::PaintDC->new($self);
    $dc->BeginDrawing();
    $self->GetParentFrame()->Draw($dc, $self->GetClientSize());
    $dc->EndDrawing();
}

sub OnMouseDown {
    my ($self, $mevent) = @_;
    my $pt = $mevent->GetPosition();
#    print "pt = (", $pt->x(), ", ", $pt->y(), ")\n";
    $self->GetParentFrame()->SetDragStartPosition($pt);
}

sub OnMouseMove {
    my ($self, $mevent) = @_;
    if ($mevent->Dragging()) {
        my $pt = $mevent->GetPosition();
        $self->GetParentFrame()->SetDragPosition($pt);
    }
}


sub OnMouseUp {
    my ($self, $mevent) = @_;
    my $pt = $mevent->GetPosition();
    $self->GetParentFrame()->EndDrag($pt);
}

# ============================================================
#  A class representing the main window
# ============================================================
package MyFrame;
use base qw(Wx::Frame);
use Wx qw(:sizer);             #  Import symbols related to wxSizer
use Wx qw(wxRA_SPECIFY_ROWS);  #  Import symbols related to wxRadioBox
use Wx::Event qw(EVT_BUTTON);
use Wx qw(:pen :brush);
use Wx qw(:filedialog wxID_OK wxID_CANCEL);

sub new {
    my ($ref) = @_;

    #  Call the constructor of the superclass (wxFrame)
    my $self = $ref->SUPER::new(undef, -1, 'NPointsOnSphere with wxPerl', [-1, -1], [560, 432]);

    #  Create a horizontal sizer to adjust the sizes of the subwindows
    my $sizer = Wx::BoxSizer->new(wxHORIZONTAL);
    $self->SetSizer($sizer);
    
    #  Create a first panel (to contain a MyCanvas window and a StaticText)
    my $panel = Wx::Panel->new($self, -1, [0, 0], [400, 432]);
    $sizer->Add($panel, 1, wxEXPAND);  # Horizontal weight = 1, vertically expands

    #  Create a vertical sizer to adjust sizes of MyCanvas window and a StaticText
    my $vsizer = Wx::BoxSizer->new(wxVERTICAL);
    $panel->SetSizer($vsizer);
    
    #  Create a MyCanvas window and add to the sizer
    $self->{_canvas} = MyCanvas->new($panel, -1, [0, 0], [400, 400]);
    $vsizer->Add($self->{_canvas}, 1, wxEXPAND); # Vertical weight = 1, horizontally expands
    
    #  Create a StaticText control
    $self->{_caption} = Wx::StaticText->new($panel, -1, "", [0, 400], [400, 32]);
    $vsizer->Add($self->{_caption}, 0, wxEXPAND);

    #  Create a second panel (to contain controls)
    $panel = Wx::Panel->new($self, -1, [400, 400], [160, 400]);
    $sizer->Add($panel, 0, wxEXPAND);  # Horizontal weight = 0, vertically expands

    #  Create another vertical sizer to adjust the sizes and positions of the controls
    $vsizer = Wx::BoxSizer->new(wxVERTICAL);
    $panel->SetSizer($vsizer);

    #  Create controls in the panel
    $vsizer->Add(
        Wx::StaticText->new($panel, -1, "Number of points"),
        0,      #  No vertical stretching
        wxALL,  #  Make border all around
        4);    #  Border width
    $vsizer->Add(
        $self->{_npointsSpin} = Wx::SpinCtrl->new($panel, -1, "20"),
        0, wxALL, 4);
    $self->{_npointsSpin}->SetRange(2, 100);
    $vsizer->Add(
        Wx::StaticText->new($panel, -1, "Start positions"),
        0, wxALL, 4);
    $vsizer->Add(
        $self->{_initposRadio} = Wx::RadioBox->new($panel, -1, "", [-1, -1], [-1, -1], ["Spiral", "Random", "C3 symmetric" ], 0, wxRA_SPECIFY_ROWS),
        0, wxALL, 4);
    $vsizer->Add(
        Wx::StaticText->new($panel, -1, "Max cycle"),
        0, wxALL, 4);
    $vsizer->Add(
        $self->{_maxcycleText} = Wx::TextCtrl->new($panel, -1, "1000"),
        0, wxALL, 4);
    $vsizer->AddSpacer(8);
    $vsizer->Add(
        $self->{_runButton} = Wx::Button->new($panel, -1, "Run"),
        0, wxALL | wxALIGN_CENTER_HORIZONTAL, 4);
    $vsizer->Add(
        $self->{_stopButton} = Wx::Button->new($panel, -1, "Stop"),
        0, wxALL | wxALIGN_CENTER_HORIZONTAL, 4);
    $vsizer->Add(
        $self->{_exportButton} = Wx::Button->new($panel, -1, "Export"),
        0, wxALL | wxALIGN_CENTER_HORIZONTAL, 4);
    $vsizer->AddStretchSpacer();
    $vsizer->Add(
        $self->{_quitButton} = Wx::Button->new($panel, -1, "Quit"),
        0, wxALL | wxALIGN_BOTTOM | wxALIGN_CENTER_HORIZONTAL, 4);
    $vsizer->AddSpacer(14);
    $self->{_stopButton}->Enable(0);
    $self->{_exportButton}->Enable(0);
    
    #  Bind event handling functions
    EVT_BUTTON($self, $self->{_runButton}->GetId(), \&OnRun);
    EVT_BUTTON($self, $self->{_stopButton}->GetId(), \&OnStop);
    EVT_BUTTON($self, $self->{_exportButton}->GetId(), \&OnExport);
    EVT_BUTTON($self, $self->{_quitButton}->GetId(), \&OnQuit);
    
    #  Initialize instance variables
    $self->{_quat} = [0, 0, 0, 1];
    $self->{_tempQuat} = [0, 0, 0, 1];
    $self->{_depth} = 8.0;
    $self->{_height} = 2.4;

#    $sizer->SetSizeHints($self);
#    my $panel = Wx::Panel->new($self, -1);
#    my $button = Wx::Button->new($panel, -1, 'Click me!', [30, 20], [-1, -1]);
    return $self;
}

sub OnRun {
    my $self = shift;
    $self->{_npoints} = $self->{_npointsSpin}->GetValue();
    $self->{_maxcycle} = $self->{_maxcycleText}->GetValue();
    $self->{_initpos} = $self->{_initposRadio}->GetSelection();
    if ($self->{_initpos} == 2) {
        $self->{_points} = c3SphereSet($self->{_npoints});
    } elsif ($self->{_initpos} == 1) {
        $self->{_points} = randomSphereSet($self->{_npoints});
    } else {
        $self->{_points} = generalizedSpiralSet($self->{_npoints});
    }
    $self->UpdateDisplayList(1);
    print "$self->{_npoints} $self->{_maxcycle} $self->{_initpos}\n";
    $self->{_runButton}->Enable(0);
    $self->{_stopButton}->Enable(1);
    $self->{_exportButton}->Enable(0);
    $self->{_quitButton}->Enable(0);
    $self->UniformSphere();
}

sub OnStop {
    my ($self) = @_;
    $self->{_stopFlag} = 1;
}

sub OnExport {
    my ($self) = @_;
    my $dialog = Wx::FileDialog->new($self, "Save coordinates to:", "", "export.txt", "*.txt", wxSAVE|wxOVERWRITE_PROMPT);
    if ($dialog->ShowModal() == wxID_OK) {
        my $file = $dialog->GetPath();
        if ($file ne "" && open($file, ">$file")) {
            my $i;
            for ($i = 0; $i < $self->{_npoints}; $i++) {
                printf $file "%-5d %11.8f %11.8f %11.8f\n", $i+1, @{$self->{_points}->[$i]}[0..2];
            }
            close($file);
        }
    }
    $dialog->Destroy();
}

sub OnQuit {
    my ($self) = @_;
    $self->{_stopFlag} = 1;
    $self->Close();
}

sub Draw {
    my ($self, $dc, $size) = @_;
    my $width = $size->GetWidth();
    my $height = $size->GetHeight();
    my $screenSize = ($width > $height ? $height : $width);
    $dc->SetBackground(wxWHITE_BRUSH);
    $dc->Clear();
    $self->{_screenWidth} = $width;
    $self->{_screenHeight} = $height;
    $self->{_screenSize} = $screenSize;
    my ($i, $n, $xc, $yc, $r, $k, $tw, $th, $label);
#    print "$dc $width $height\n";
    my $list = $self->{_objlist};
    my $depth = $self->{_depth};
    my $realheight = $self->{_height};
    $dc->SetPen(wxBLACK_PEN);
    $dc->SetBrush(wxCYAN_BRUSH);
    for ($i = 0; $i <= $#{$list}; $i++) {
        $k = ($depth + $list->[$i]->[3]) / ($depth * $realheight) * $screenSize;
        $r = $k * 0.08;
        $xc = $list->[$i]->[1] * $k + $width / 2;
        $yc = $list->[$i]->[2] * $k + $height / 2;
#        print "entry = @{$list->[$i]} DrawCircle($xc, $yc, $r)\n";
        $dc->DrawCircle($xc, $yc, $r);
        $label = $list->[$i]->[0] + 1;
        ($tw, $th) = $dc->GetTextExtent($label);
        $dc->DrawText($label, $xc - $tw / 2, $yc - $th / 2);
    }
}

sub SetDragStartPosition {
    my ($self, $pt) = @_;
    my ($x, $y);
    $x = ($pt->x() - $self->{_screenWidth} * 0.5) / $self->{_screenSize};
    $y = ($pt->y() - $self->{_screenHeight} * 0.5) / $self->{_screenSize};
    $self->{_xstart} = $x;
    $self->{_ystart} = $y;
}

sub SetDragPosition {
    my ($self, $pt) = @_;
    my ($x, $y);
    return if (!defined($self->{_xstart}) || !defined($self->{_ystart}));
    $x = ($pt->x() - $self->{_screenWidth} * 0.5) / $self->{_screenSize};
    $y = ($pt->y() - $self->{_screenHeight} * 0.5) / $self->{_screenSize};
    calcTrackballRotation($self->{_xstart}, $self->{_ystart}, $x, $y, $self->{_tempQuat});
    $self->UpdateDisplayList(0);
    $self->{_canvas}->Update();
}

sub EndDrag {
    my ($self, $pt) = @_;
    return if (!defined($self->{_xstart}) || !defined($self->{_ystart}));
    quatMul($self->{_quat}, $self->{_tempQuat}, $self->{_quat});
    @{$self->{_tempQuat}} = (0, 0, 0, 1);
    $self->UpdateDisplayList(0);
    $self->{_canvas}->Update();
    $self->{_xstart} = $self->{_ystart} = undef;
}

#  Update $self->{_objlist} and redraw the canvas
#  $self->{_objlist} = [[i1, x1, y1, z1], [i2, x2, y2, z2],...];
#  i: the index of the point, (x, y, z): the screen coordinates
#  The list is sorted by the descending order of z coordinates
sub UpdateDisplayList {
    my ($self, $create) = @_;
    my ($i, $k, $xx, $yy, $rr, $xs, $ys, $xe, $ye, $points, $list);
    my $conv = [];
    my $myQuat = quatMul($self->{_quat}, $self->{_tempQuat});
    my $myQuatConj = quatConjugate($myQuat);
    my @z;
    $points = $self->{_points};
    if ($create) {
        #  Allocate the memory
        $self->{_objlist} = [];
        for ($i = 0; $i <= $#{$points}; $i++) {
            $self->{_objlist}->[$i] = [0, 0, 0, 0];
        }
    }
    my $list = $self->{_objlist};
    for ($i = 0; $i <= $#{$points}; $i++) {
        quatMul($myQuatConj, quatMul($points->[$i], $myQuat, $conv), $conv);
#        printf "conv($i) = [%f %f %f %f]\n", @$conv;
        @{$list->[$i]}[0..3] = ($i, @{$conv}[0..2]);
    }
    @{$list} = sort { $a->[3] <=> $b->[3] } @{$list};
    $self->{_canvas}->Refresh();
}

sub UniformSphere {
    my ($self) = @_;
    my $points = $self->{_points};
    my $npoints = $#{$points} + 1;
    my @force;
    my $r = 1.0;
    my $iter = 0;
    my $maxIter = $self->{_maxcycle};
    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];
    }
    $self->{_stopFlag} = 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;
        }
        $self->{_caption}->SetLabel(sprintf("Cycle %d, repulsion energy %f, max force %f,\ndamp %f, max move %f%s\n", $iter, $eng, $maxf, $damp, $maxw, ($converged ? " --- converged." : "")));
        $self->UpdateDisplayList(0);
        $self->{_canvas}->Update();
    #    while ($main::app->Pending()) {
    #        $main::app->Dispatch();
    #    }
        $main::app->Yield();
        last if $self->{_stopFlag};
        last if $converged || $iter >= $maxIter;
    }
    $self->{_runButton}->Enable(1);
    $self->{_stopButton}->Enable(0);
    $self->{_exportButton}->Enable(1);
    $self->{_quitButton}->Enable(1);
    return $eng;
}

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 c3SphereSet {
    my ($npoints) = @_;
    my $results = [];
    my ($i, $j, $n, $m, $x, $y, $z);
    $n = int($npoints / 3);
    $m = $npoints % 3;
    for ($i = 0; $i < $n; $i++) {
        my $v = randomUnitVector();
        $results->[$i * 3] = $v;
        $x = $v->[0] * (-0.5) + $v->[1] * 0.866025;
        $y = $v->[0] * (-0.866025) + $v->[1] * (-0.5);
        $results->[$i * 3 + 1] = [$x, $y, $v->[2]];
        $x = $v->[0] * (-0.5) + $v->[1] * (-0.866025);
        $y = $v->[0] * (0.866025) + $v->[1] * (-0.5);
        $results->[$i * 3 + 2] = [$x, $y, $v->[2]];
    }
    if ($m > 0) {
        unshift @{$results}, [0, 0, 1];
        if ($m > 1) {
            push @{$results}, [0, 0, -1];
        }
    }
    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 ($MyFrame::haveNextNextGaussian) {
        $MyFrame::haveNextNextGaussian = 0;
        return $MyFrame::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);
        $MyFrame::nextNextGaussian = $v2 * $mul;
        $MyFrame::haveNextNextGaussian = 1;
        return $v1 * $mul;
    }
}

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

# ============================================================
#  The application class
# ============================================================

package MyApp;
use base qw(Wx::App);
sub OnInit {
    my $frame = MyFrame->new;
    $frame->Show(1);
}

package main;

$main::app = MyApp->new;
$main::app->MainLoop;