#===============================================
# TkPlanet
#===============================================
package TkPlanet;
use clib::TkFittingCommon;
@ISA = qw(TkFittingCommon);

#公開したいサブルーチン
#@EXPORT = qw(DelSpace Reduce01 MakePath RegExpQuote);

use strict;
use Math::Spline;

use Utils;
use IniFile;
use CSV;

use MultiColumnData;
use MyTk::GraphFrameArray;

use Sci qw($pi $kB $R $NA $c $e $e0 $u0 $me $mp $mn $h $hbar $G $DayToSecond $SecondToDay $AU);
use Sci::Planet;

#my $G = Sci::G();
#my $DayToSecond = Sci::DayToSecond();
#my $SecondToDay = Sci::SecondToDay();
#my $AU = Sci::AstronomicalUnit();
my $G1 = $G * $DayToSecond * $DayToSecond / $AU / $AU / $AU;

my $Sun     = new Planet("太陽");
my $Mercury = new Planet("水星");
my $Venus   = new Planet("金星");
my $Earth   = new Planet("地球");
my $Mars    = new Planet("火星");
my $Jupiter = new Planet("木星");
#$Jupiter->{M} = $Sun->{M} / 30.0;
my $Saturn  = new Planet("土星");
my $Uranus  = new Planet("天王星");
my $Neptune = new Planet("海王星");
my $Pluto   = new Planet("冥王星");
my @planets = ($Sun, $Mercury, $Venus, $Earth, $Mars, $Jupiter,
		$Saturn, $Uranus, $Neptune, $Pluto);
my $nPlanet = @planets;

my $PrintForce = 0;
my $PrintVelocity = 0;
my $PrintPosition = 1;
my $PrintZ = 0;

my $dt                         = 1.0;
my $nStep                      = 1000000;
my $OutputInterval             = 1000;
my $GraphOutputInterval        = 10;
my $GraphPlotInterval          = 100;

#============================================================
#　変数等取得関数
#============================================================

#============================================================
#　コンストラクタ、デストラクタ
#============================================================
#呼び出されることはない
sub new
{
	my ($module, $app) = @_;
	my $this = {};
	bless $this;

	return $this;
}

sub DESTROY
{
	my $this = shift;
	$this->SUPER::DESTROY(@_);
}

#============================================================
#　継承クラスで定義しなおす関数
#============================================================
sub InsertToolBar
{
	my ($this, $frame) = @_;
}

sub CreateLeftFrame
{
	my ($this, $ConfigSide, @args) = @_;
#	return $this->SUPER::CreateLeftFrame();
	return undef;
}

sub CreateSelectFilePane
{
	my ($this) = @_;
#	return $this->SUPER::CreateSelectFilePane();
	return undef;
}

sub CreateFileContentPane
{
	return undef;
}

sub CreateWidgets
{
	my ($this, $DoSuperOnly) = @_;
	$DoSuperOnly = 0 if(!defined $DoSuperOnly);

	my $App    = $this->App();
	my $args   = $App->Args();
	my $mw     = $this->mw();
	my $Program = "Planet2007";

	$this->{ini}       = new IniFile;
	$this->{DataFile}  = new MultiColumnData;
	$this->{Executing} = 0;

	$this->SUPER::CreateWidgets();
	return if($DoSuperOnly);

	$mw->SetTitle("$Program / TkPlot");
	$this->App()->SetProgram("$Program / TkPlot");

#===================================================
# ペインを作製
#===================================================
	my $EntryWidth =  12;

	my $Frame = $mw->MyFrame()->pack(-anchor => 'nw');
	$this->MakeLabelEntry($Frame, "dt",                  "dt:",         \$dt,                               "%g", $EntryWidth, "", 0, 0);
	$this->MakeLabelEntry($Frame, "nStep",               "nStep:",      \$nStep,                            "%d", $EntryWidth, "", 0, 1);
	$this->MakeLabelEntry($Frame, "OutputInterval",      "OutputItvl:", \$OutputInterval,                   "%d", $EntryWidth, "", 0, 2);
	$this->MakeLabelEntry($Frame, "GraphOutputInterval", "GraphIntvl:", \$GraphOutputInterval,              "%d", $EntryWidth, "", 0, 3);
	$this->MakeLabelEntry($Frame, "GraphPlotInterval",   "PlotIntvl:",  \$GraphPlotInterval,                "%d", $EntryWidth, "", 0, 4);
	$this->MakeLabelEntry($Frame, "GraphRange",          "GraphRange:", $this->{ini}->pVariable("GraphRange", 3.0), "%g", $EntryWidth, "", 0, 5);

	$Frame = $mw->MyFrame()->pack(-anchor => 'nw');
	$this->{ChooseSolverListBox} = $Frame->MyBrowseEntry(
			-label      => "Differential Equation Solver:",
			-state      => "readonly",
			-takefocus  => 1,
			-width      => 15,
			-Selections => ["Euler", "Verlet"],
			-SelIndex   => 1,
		)->pack(-side => 'left', -expand => 'yes');

	for(my $i = 0 ; $i < $nPlanet ; $i++) {
		my $p = $planets[$i];
		$Frame = $mw->MyFrame()->pack(-anchor => 'nw');
		$this->MakeLabelEntry($Frame, "Name$i", "Name:", $p->pVariable("Name", ""),  "%s", $EntryWidth, "",      0, 0);
		$this->MakeLabelEntry($Frame, "M$i",    "M:",    $p->pVariable("M", 1.0),    "%g", $EntryWidth, "kg",    0, 1); # Mass
		$this->MakeLabelEntry($Frame, "RevR$i", "RevR:", $p->pVariable("RevR", 1.0), "%g", $EntryWidth, "m",     0, 2); # Revolution Radius
		$this->MakeLabelEntry($Frame, "RevV$i", "RevV:", $p->pVariable("RevV", 1.0), "%g", $EntryWidth, "m/s",   0, 3); # Revolution Velocity
		$this->MakeLabelEntry($Frame, "R$i",    "R:",    $p->pVariable("R", 1.0),    "%g", $EntryWidth, "m",     0, 4); # Rdius
		$this->MakeLabelEntry($Frame, "V$i",    "V:",    $p->pVariable("V", 1.0),    "%g", $EntryWidth, "m3",    0, 5); # Volume
		$this->MakeLabelEntry($Frame, "d$i",    "d:",    $p->pVariable("d", 1.0),    "%g", $EntryWidth, "g/cm3", 0, 6); # Density
#		$this->MakeLabelEntry($Frame, "RotV$i", "RotV:", $p->pVariable("RotV", 1.0), "%g", $EntryWidth, "m/s",   0, 7); # Rotation Velocity
	}

	$this->UpdateParameters();
	$Frame = $mw->MyFrame()->pack(-anchor => 'nw');
	$this->{RecalcButton} = $Frame->MyButton(
			-text => '&Recalc',
			-takefocus => 1,
			-command => sub { $this->ButtonPressed('RButtonDown', 'Recalc'); },
		)->pack(-side => 'left');

#ツールバーのOpenボタンをファイル読み込みにバインドする
	$mw->{OpenButton}->configure(
			-command  => sub { $this->ChooseFile($this->{SampleCSVEntry}); },
		) if($mw->{OpenButton});
}

#Set arrays and initial parameters
sub Initialize
{
#	my ($this, $InitializeParameterFile, $InitializeData) = @_;
	my ($this, $pPlanets, $pM, $pTime, $ppfx, $ppfy, $ppfz, $ppvx, $ppvy, $ppvz, $ppx, $ppy, $ppz) = @_;

	my $App = $this->App();

	my $nPlanet = @$pPlanets;
	$pTime->[0] = 0.0;
	for(my $ip = 0 ; $ip < $nPlanet ; $ip++) {
		$pM->[$ip] = $pPlanets->[$ip]->{M};
		my $r0  = $pPlanets->[$ip]->{RevR} / $AU;
		my $v0  = $pPlanets->[$ip]->{RevV} * $DayToSecond / $AU;

$App->printf("%8s: M=%8.3ekg R=%8.3ekm RevR=%8.3e RevV=%8.3e\n",
	$pPlanets->[$ip]->{Name}, $pM->[$ip], $pPlanets->[$ip]->{R}*1e-3, $r0, $v0);

		my (@fx, @fy, @fz, @vx, @vy, @vz, @x, @y, @z);
		$x[0]  = $r0;
		$y[0]  = 0.0;
		$z[0]  = 0.0;
		$vx[0] = 0.0;
		$vy[0] = $v0;
		$vz[0] = 0.0;
		$ppfx->[$ip] = \@fx;
		$ppfy->[$ip] = \@fy;
		$ppfz->[$ip] = \@fz;
		$ppvx->[$ip] = \@vx;
		$ppvy->[$ip] = \@vy;
		$ppvz->[$ip] = \@vz;
		$ppx->[$ip]  = \@x;
		$ppy->[$ip]  = \@y;
		$ppz->[$ip]  = \@z;
	}
}

sub SelChangeListBox
{
	my ($this, $type, $lb) = @_;
#print "this:$this  lb=$lb\n";
	my $s = $lb->GetText();
	if($type eq 'X' or $type eq 'Y1') {
	}
	$this->AssignGraphData();
	$this->Draw();
}

sub ButtonPressed
{
	my ($this, $event, $type) = @_;
	my $App = $this->App();
#print "this=$this\n";
	if($type eq 'Recalc') {
		$this->Recalc();
	}
	elsif($type eq 'SaveParam') {
		$this->SaveParameterFile($this->{ini}->{ParameterFile});
		return;
	}
	elsif($type eq 'EditParam') {
		my $Editor  = $this->App()->{EditorPath};
		my $File = $this->{ini}->{ParameterFile};
		my $command = "$Editor $File";
		system($command);
		return;
	}
	elsif($type eq 'EditSampleCSV') {
		my $Editor  = $this->App()->{EditorPath};
		my $File = $this->{ini}->{SampleCSVFile};
		my $command = "$Editor $File";
		system($command);
		return;
	}
	elsif($type eq 'EditSubstrateCSV') {
		my $Editor  = $this->App()->{EditorPath};
		my $File = $this->{ini}->{SubstrateCSVFile};
		my $command = "$Editor $File";
		system($command);
		return;
	}
}

#===================================================
# データ処理
#===================================================

sub Recalc
{
	my ($this, $ResetViewRange) = @_;
	$ResetViewRange = 1 if(!defined $ResetViewRange);

	my $App = $this->App();
	if(!$this->{Executing}) {
		$this->{Executing} = 1;
		$this->{RecalcButton}->SetText("&Stop");
		$this->mw()->update();
	}
	else {
		$this->{Executing} = 0;
		$this->{RecalcButton}->SetText("&Recalc");
		return;
	}

	my $DifferentialEquationSolver = $this->{ChooseSolverListBox}->GetText();

	my (@M, @Time, @pfx, @pfy, @pfz, @pvx, @pvy, @pvz, @px, @py, @pz);
	$this->Initialize(\@planets, \@M, \@Time, \@pfx, \@pfy, \@pfz, \@pvx, \@pvy, \@pvz, \@px, \@py, \@pz);
	$this->Normalize(\@planets, 0, \@M, \@Time, \@pfx, \@pfy, \@pfz, \@pvx, \@pvy, \@pvz, \@px, \@py, \@pz);

	my (@x, @y, @z);

	my (@U, @K, @UTotal);
	$K[0] = 0.0;
	my $c = 0;
	for(my $ip0 = 0 ; $ip0 < $nPlanet ; $ip0++) {
		$x[$ip0] = [];
		$y[$ip0] = [];
		$z[$ip0] = [];

		my $M0 = $M[$ip0];
		($U[0], $pfx[$ip0][0], $pfy[$ip0][0], $pfz[$ip0][0]) = CalUFi($ip0, 0, $nPlanet, \@M, \@px, \@py, \@pz);
		$K[0] += 0.5 * $M0 * ( $pvx[$ip0][0] * $pvx[$ip0][0]
							 + $pvy[$ip0][0] * $pvy[$ip0][0]
							 + $pvz[$ip0][0] * $pvz[$ip0][0]);
		$App->printf("%05d: %2d %10s: (%12.5g, %12.5g, %12.5g)\n", 
					0, $ip0, $planets[$ip0]->{Name}, $pvx[$ip0][0], $pvy[$ip0][0], $pvz[$ip0][0]);

		$x[$ip0][$c] = $px[$ip0][0];
		$y[$ip0][$c] = $py[$ip0][0];
		$z[$ip0][$c] = $pz[$ip0][0];
#print "c=$c: $px[$ip0][0],$py[$ip0][0],$pz[$ip0][0]\n";
	}
	$c++;

	$this->{pXArray} = \@x;
	$this->{pYArray} = \@y;
	$this->{pZArray} = \@z;
	$this->CreateGraphFrame() if(!defined $this->{GraphFrameArray});
	$this->AssignGraphData(1);

	my $Invdt2   = 1.0 / ($dt + $dt);
	my $dtSqr = $dt * $dt;
	$UTotal[0] = $U[0] + $K[0];
	$App->print("  U=$U[0] + K=$K[0] = $UTotal[0]\n");
	for(my $it = 1 ; $it < $nStep ; $it++) {
		last if(!$this->{Executing});

		$Time[$it] = $dt * $it;
		my $itm1 = $it-1;
		my $itm2 = $it-2;
		$K[$it] = 0.0;
		for(my $ip0 = 0 ; $ip0 < $nPlanet ; $ip0++) {
			my $M0 = $M[$ip0];
			my $pxi  = $px[$ip0];
			my $pyi  = $py[$ip0];
			my $pzi  = $pz[$ip0];
			my $pvxi = $pvx[$ip0];
			my $pvyi = $pvy[$ip0];
			my $pvzi = $pvz[$ip0];

			if($it ==1 or $DifferentialEquationSolver eq "Euler") {
				$pvxi->[$it] = $pvxi->[$itm1] + $dt * $pfx[$ip0][$itm1] / $M0;
				$pvyi->[$it] = $pvyi->[$itm1] + $dt * $pfy[$ip0][$itm1] / $M0;
				$pvzi->[$it] = $pvzi->[$itm1] + $dt * $pfz[$ip0][$itm1] / $M0;
				$pxi->[$it]  = $pxi->[$itm1]  + $dt * $pvxi->[$it];
				$pyi->[$it]  = $pyi->[$itm1]  + $dt * $pvyi->[$it];
				$pzi->[$it]  = $pzi->[$itm1]  + $dt * $pvzi->[$it];
			}
			elsif($DifferentialEquationSolver eq "Verlet") {
				$pxi->[$it]  = 2.0*$pxi->[$itm1] - $pxi->[$itm2]  + $dtSqr * $pfx[$ip0][$itm1] / $M0;
				$pyi->[$it]  = 2.0*$pyi->[$itm1] - $pyi->[$itm2]  + $dtSqr * $pfy[$ip0][$itm1] / $M0;
				$pzi->[$it]  = 2.0*$pzi->[$itm1] - $pzi->[$itm2]  + $dtSqr * $pfz[$ip0][$itm1] / $M0;
				$pvxi->[$it] = ($pxi->[$it]  - $pxi->[$itm2]) * $Invdt2;
				$pvyi->[$it] = ($pyi->[$it]  - $pyi->[$itm2]) * $Invdt2;
				$pvzi->[$it] = ($pzi->[$it]  - $pzi->[$itm2]) * $Invdt2;
			}
			else {
				$App->print("\nError!!: Unrecognized solver [$DifferentialEquationSolver].\n");
				return -2;
			}

			$K[$it] += 0.5 * $M0 *  ( $pvxi->[$it] * $pvxi->[$it] + $pvyi->[$it] * $pvyi->[$it] + $pvzi->[$it] * $pvzi->[$it]);
			$px[$ip0][$it]  = $pxi->[$it];
			$py[$ip0][$it]  = $pyi->[$it];
			$pz[$ip0][$it]  = $pzi->[$it];
			$pvx[$ip0][$it] = $pvxi->[$it];
			$pvy[$ip0][$it] = $pvyi->[$it];
			$pvz[$ip0][$it] = $pvzi->[$it];

			if($it % $OutputInterval == 0) {
				$App->printf("%05d: %2d %10s: (%12.5g, %12.5g, %12.5g)\n", 
							$it, $ip0, $planets[$ip0]->{Name}, $pvx[$ip0][$it], $pvy[$ip0][$it], $pvz[$ip0][$it]);
			}
			if($it % $GraphOutputInterval == 0) {
				$x[$ip0][$c] = $px[$ip0][$it];
				$y[$ip0][$c] = $py[$ip0][$it];
				$z[$ip0][$c] = $pz[$ip0][$it];
			}
		}

		$U[$it] = 0.0;
		for(my $ip0 = 0 ; $ip0 < $nPlanet ; $ip0++) {
			my $M0 = $M[$ip0];
			($U[$it], $pfx[$ip0][$it], $pfy[$ip0][$it], $pfz[$ip0][$it]) = CalUFi($ip0, $it, $nPlanet, \@M, \@px, \@py, \@pz);
		}
		$UTotal[$it] = $U[$it] + $K[$it];

		if($it % $OutputInterval == 0) {
			$App->print("  U=$U[$it] + K=$K[$it] = $UTotal[$it]\n");
		}
		if($it % $GraphOutputInterval == 0) {
			$c++;
		}
		if($it % $GraphPlotInterval == 0) {
			$this->Draw();
			$this->mw()->update();
		}
	}

	$this->Draw();
}

sub CreateGraphFrame
{
	my ($this, $canvas, $TargetData) = @_;
	$canvas = $this->Canvas();
	my $App = $this->App();
	my $font = $App->{'GraphFrameFont'};
	my @font = split(/,/, $font) if($font);
	my $w = $canvas->width();
	my $h = $canvas->height();

	my $GraphFrameArray = $this->{GraphFrameArray} = new GraphFrameArray($this->mw());

	$GraphFrameArray->SetCanvasSize($w, $h);

	$GraphFrameArray->AddGraphFrame();
	my $GraphFrame0  = $GraphFrameArray->GetGraphFrame(0);
	my $FramePosStr0 = $App->{"GraphFrame0Position"};
	$GraphFrame0->SetPositionByStr($FramePosStr0);

	my $XScale0      = $GraphFrame0->GetXScale(0);
	my $YScale0      = $GraphFrame0->GetYScale(0);
	$XScale0->SetScaleStringVisible(1);
	$XScale0->SetCaptionVisible(1);

	$GraphFrame0->SetViewRange(0, 0, 1, 1);
}

sub AssignGraphData
{
	my ($this, $ResetViewRange) = @_;
	$ResetViewRange = 1 if(!defined $ResetViewRange);

	my $GraphFrameArray = $this->GetGraphFrameArray();
	my @GraphFrame;
	$GraphFrame[0]  = $GraphFrameArray->GetGraphFrame(0);
	$GraphFrame[0]->ClearAllData();

	my $pxArray = $this->{pXArray};
	my $pyArray = $this->{pYArray};
	my $pzArray = $this->{pZArray};
	my $nArray = @$pxArray;

	for(my $i = 0 ; $i < $nArray ; $i++) {
		$GraphFrame[0]->AddGraphData($pxArray->[$i], $pyArray->[$i],
				1, "black", "", 6, "red", 0, "red", "XAutoSkip", "Energy", $planets[$i]->Name()) if($pyArray->[$i]);
	}

	$GraphFrame[0]->SetYScalePlotType('x');
	if($ResetViewRange) {
		$GraphFrame[0]->CalMinMax();
		my $r = $this->{ini}->{GraphRange};
		if($r <= 0.0) {
			$GraphFrame[0]->AdjustViewRange(0.05, 0.05, 0.05, 0.05);
		}
		else {
			$GraphFrame[0]->SetViewRange(-$r, -$r, $r, $r);
		}
	}
}

#==========================================
# &Subroutines
#==========================================
sub CalTotalMass
{
	my ($nPlanet, $pM) = @_;

	my $MSum = 0.0;
	for(my $ip = 0 ; $ip < $nPlanet ; $ip++) {
		my $M = $pM->[$ip];
		$MSum += $M;
	}
	return $MSum;
}

sub CalTotalMomentum
{
	my ($nPlanet, $it, $pM, $ppvx, $ppvy, $ppvz) = @_;

	my ($Px, $Py, $Pz) = (0.0, 0.0, 0.0);
	for(my $ip = 0 ; $ip < $nPlanet ; $ip++) {
		my $M = $pM->[$ip];
		$Px += $M * $ppvx->[$ip][$it];
		$Py += $M * $ppvy->[$ip][$it];
		$Pz += $M * $ppvz->[$ip][$it];
	}
	return ($Px, $Py, $Pz);
}

#Normalize total momentum to zero
sub Normalize
{
	my ($this, $pPlanets, $it, $pM, $pTime, $ppfx, $ppfy, $ppfz, $ppvx, $ppvy, $ppvz, $ppx, $ppy, $ppz) = @_;

	my $App = $this->App();
	my $nPlanet = @$pPlanets;

	my $MSum = &CalTotalMass($nPlanet, $pM);
	my ($Px, $Py, $Pz) = &CalTotalMomentum($nPlanet, $it, $pM, $ppvx, $ppvy, $ppvz);
	$App->print("Initial momentum: ($Px, $Py, $Pz)\n");

	$App->print("Normalized velocities\n");
	for(my $ip = 0 ; $ip < $nPlanet ; $ip++) {
		my $p = $pPlanets->[$ip];
		my $M = $p->{M};
		$ppvx->[$ip][$it] -= $Px / $MSum;
		$ppvy->[$ip][$it] -= $Py / $MSum;
		$ppvz->[$ip][$it] -= $Pz / $MSum;
		$App->print("  $p->{Name}: (", 
				$ppvx->[$ip][$it], ", ", 
				$ppvy->[$ip][$it], ", ", 
				$ppvz->[$ip][0], ")\n");
	}
}

sub CalUFij
{
	my ($m0, $x0, $y0, $z0, $m1, $x1, $y1, $z1) = @_;

	my $dx = $x1 - $x0;
	my $dy = $y1 - $y0;
	my $dz = $z1 - $z0;
	my $r2 = $dx*$dx + $dy*$dy + $dz*$dz;
	my $r = sqrt($r2);

	my $g = $G1 * $m0 * $m1;
	my $u  = -$g / $r;
	my $f  =  $g / $r2;
	my $fx = $f * $dx / $r;
	my $fy = $f * $dy / $r;
	my $fz = $f * $dz / $r;
	return ($r2, $r, $u, $f, $fx, $fy, $fz);
}

sub CalUFi
{
	my ($ip0, $it, $nPlanet, $pM, $ppx, $ppy, $ppz) = @_;

	my $M0 = $pM->[$ip0];
	my ($fx, $fy, $fz) = (0.0, 0.0, 0.0);
	my $U = 0.0;
	for(my $ip1 = 0 ; $ip1 < $nPlanet ; $ip1++) {
		next if($ip1 == $ip0);

		my $M1 = $pM->[$ip1];
		my ($r2, $r, $Uij, $Fij, $Fxij, $Fyij, $Fzij)
			= &CalUFij($M0, $ppx->[$ip0][$it], $ppy->[$ip0][$it], $ppz->[$ip0][$it],
				   $M1, $ppx->[$ip1][$it], $ppy->[$ip1][$it], $ppz->[$ip1][$it]);

		$U += $Uij;
		$fx += $Fxij;
		$fy += $Fyij;
		$fz += $Fzij;
#print "$it: $ip0 - $ip1: dx=($dx, $dy, $dz)\n";
#print "$it: x=($px[$ip1][$itm1], $py[$ip1][$itm1], $pz[$ip1][$itm1]) r=$r\n";
#print "$it: f/m=($pfx[$ip0][$itm1], $pfy[$ip0][$itm1], $pfz[$ip0][$itm1])\n";
	}
	return ($U, $fx, $fy, $fz);
}

1;