#!/usr/bin/perl

use lib 'f:/Programs/Perl/lib';
use lib 'd:/Programs/Perl/lib';

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'PerlDir'};
$BaseDir = $ENV{'TkPerlDir'} unless(defined $BaseDir);
print "BaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", @INC);
}

use strict;
#use warnings;

use Math::Matrix;
use Math::MatrixReal;

use Deps;
use Utils;
use JFile;
use CSV;

use MyApplication;
use Sci::Science;
use Sci::Planet;

#===============================================
# デバッグ関係変数
#===============================================
#$PrintLevelが大きいほど、情報が詳しくなる
my $PrintLevel = 0;

#===============================================
# 文字コード関係変数
#===============================================
# sjis, euc, jis, noconv
my $PrintCharCode      = Deps::PrintCharCode();
my $OSCharCode         = Deps::OSCharCode();
my $FileSystemCharCode = Deps::FileSystemCharCode();

my $LF        = Deps::LF();
my $DirSep    = Deps::DirSep();
my $RegDirSep = Deps::RegDirSep();

#===============================================
# Applicationオブジェクト作成
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);
#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);

#===============================================
# スクリプト大域変数
#===============================================
my $InitialDirectory = Deps::GetWorkingDirectory();
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;

#CSV configuration
my $PrintForce = 0;
my $PrintVelocity = 0;
my $PrintPosition = 1;
my $PrintZ = 0;

my $dt = 1.0;
my $nStep = 10000;
my $OutputInterval = 1000;
my $CSVOutputInterval = 10;
my $CriteriaR2ForIdenticalAtom = 0.001;
my $DifferentialEquationSolver = "Verlet";
#my $DifferentialEquationSolver = "Euler";
my $CSVFile = "Planet.csv";

#==========================================
# コマンドラインオプション読み込み
#==========================================
#$App->AddArgument("--Action", "--Action=[Convolute]",       '');
#exit 1 if($App->ReadArgs(0) != 1);
#my $Args = $App->Args();
##my $form = new CGI;
##$Args->SetCGIForm($form);
##$Args->parseInput($WebCharCode);

#==========================================
# メイン関数スタート
#==========================================

&Execute();

exit;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &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 ($pPlanets, $it, $pM, $pTime, $ppfx, $ppfy, $ppfz, $ppvx, $ppvy, $ppvz, $ppx, $ppy, $ppz) = @_;
	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");
	}
}

#Set arrays and initial parameters
sub Initialize
{
	my ($pPlanets, $pM, $pTime, $ppfx, $ppfy, $ppfz, $ppvx, $ppvy, $ppvz, $ppx, $ppy, $ppz) = @_;
	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 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;
#	return (0.0, 0.0) if(abs($r2) == $CriteriaR2ForIdenticalAtom);
	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);
}

sub Execute
{
	my (@M, @Time, @pfx, @pfy, @pfz, @pvx, @pvy, @pvz, @px, @py, @pz);

	&Initialize(\@planets, \@M, \@Time, \@pfx, \@pfy, \@pfz, \@pvx, \@pvy, \@pvz, \@px, \@py, \@pz);
	&Normalize(\@planets, 0, \@M, \@Time, \@pfx, \@pfy, \@pfz, \@pvx, \@pvy, \@pvz, \@px, \@py, \@pz);

#Do simulation
#$nStep=10;
	my (@U, @K, @UTotal);
	$K[0] = 0.0;
	for(my $ip0 = 0 ; $ip0 < $nPlanet ; $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->print("0: $ip0: ($pvx[$ip0][0], $pvy[$ip0][0], $pvz[$ip0][0])\n");
	}
	$UTotal[0] = $U[0] + $K[0];
	$App->print("  U=$U[0] + K=$K[0] = $UTotal[0]\n");
	for(my $it = 1 ; $it < $nStep ; $it++) {
		$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];

			if($it ==1 or $DifferentialEquationSolver eq "Euler") {
				$pvx[$ip0][$it] = $pvx[$ip0][$itm1] + $dt * $pfx[$ip0][$itm1] / $M0;
				$pvy[$ip0][$it] = $pvy[$ip0][$itm1] + $dt * $pfy[$ip0][$itm1] / $M0;
				$pvz[$ip0][$it] = $pvz[$ip0][$itm1] + $dt * $pfz[$ip0][$itm1] / $M0;
				$px[$ip0][$it]  = $px[$ip0][$itm1]  + $dt * $pvx[$ip0][$it];
				$py[$ip0][$it]  = $py[$ip0][$itm1]  + $dt * $pvy[$ip0][$it];
				$pz[$ip0][$it]  = $pz[$ip0][$itm1]  + $dt * $pvz[$ip0][$it];
			}
			elsif($DifferentialEquationSolver eq "Verlet") {
				$px[$ip0][$it]  = 2.0*$px[$ip0][$itm1] - $px[$ip0][$itm2]  + $dt*$dt * $pfx[$ip0][$itm1] / $M0;
				$py[$ip0][$it]  = 2.0*$py[$ip0][$itm1] - $py[$ip0][$itm2]  + $dt*$dt * $pfy[$ip0][$itm1] / $M0;
				$pz[$ip0][$it]  = 2.0*$pz[$ip0][$itm1] - $pz[$ip0][$itm2]  + $dt*$dt * $pfz[$ip0][$itm1] / $M0;
				$pvx[$ip0][$it] = ($px[$ip0][$it]  - $px[$ip0][$itm2]) / 2.0 / $dt;
				$pvy[$ip0][$it] = ($py[$ip0][$it]  - $py[$ip0][$itm2]) / 2.0 / $dt;
				$pvz[$ip0][$it] = ($pz[$ip0][$it]  - $pz[$ip0][$itm2]) / 2.0 / $dt;
			}
			else {
				$App->print("\nError!!: Unrecognized solver [$DifferentialEquationSolver].\n");
				return -2;
			}

			$K[$it] += 0.5 * $M0 *  ( $pvx[$ip0][$it] * $pvx[$ip0][$it]
						+ $pvy[$ip0][$it] * $pvy[$ip0][$it]
						+ $pvz[$ip0][$it] * $pvz[$ip0][$it]);
			if($it % $OutputInterval == 0) {
				$App->print("$it: $ip0: ($pvx[$ip0][$it], $pvy[$ip0][$it], $pvz[$ip0][$it])\n");
			}
		}

		$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");
		}
	}

#Prepare for CSV
	my @LabelList = ("time(t)", "U", "K", "UTotal");
	for(my $ip = 0 ; $ip < $nPlanet ; $ip++) {
		my $n = $planets[$ip]->{Name};
		if($PrintZ) {
			push(@LabelList, "$n fx", "$n fy", "$n fz") if($PrintForce);
			push(@LabelList, "$n vx", "$n vy", "$n vz") if($PrintVelocity);
			push(@LabelList, "$n x",  "$n y",  "$n z")  if($PrintPosition);
		}
		else {
			push(@LabelList, "$n fx", "$n fy") if($PrintForce);
			push(@LabelList, "$n vx", "$n vy") if($PrintVelocity);
			push(@LabelList, "$n x",  "$n y")  if($PrintPosition);
		}
	}
	my @DataArray = (\@Time, \@U, \@K, \@UTotal);
	for(my $ip = 0 ; $ip < $nPlanet ; $ip++) {
		my $n = $planets[$ip]->{Name};
		if($PrintZ) {
			push(@DataArray, $pfx[$ip], $pfy[$ip], $pfz[$ip]) if($PrintForce);
			push(@DataArray, $pvx[$ip], $pvy[$ip], $pvz[$ip]) if($PrintVelocity);
			push(@DataArray, $px[$ip],  $py[$ip],  $pz[$ip])  if($PrintPosition);
		}
		else {
			push(@DataArray, $pfx[$ip], $pfy[$ip]) if($PrintForce);
			push(@DataArray, $pvx[$ip], $pvy[$ip]) if($PrintVelocity);
			push(@DataArray, $px[$ip],  $py[$ip])  if($PrintPosition);
		}
	}

#Write to CSV
	my $CSV = new CSV;
	if($CSV->Save($CSVFile, \@LabelList, \@DataArray, $CSVOutputInterval)) {
		$App->print("Data was saved to [$CSVFile]\n");
	}
	else {
		$App->print("Can not write to [$CSVFile]\n");
		return -1;
	}
	
	return 1;
}