#!/usr/bin/perl

use strict;

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

use JFile;

use Crystal::AtomType;
use Sci qw($todeg $torad acos);

#my $infile  = "C2H4.dat";
#my $outfile = "C2H4-out.mld";
my $infile  = "C6H5OH.dat";
my $outfile = "C6H5OH-out.mld";

my $in  = new JFile;
$in->Open($infile,  "r") or die "$!: Can not read [$infile]\n";
my $out = new JFile;
$out->Open($outfile, "w") or die "$!: Can not write to [$infile]\n";

my $Condition = $in->ReadLine(1);
$in->ReadLine(1);
my $Title = $in->ReadLine(1);

my (@atom, @d, @idd, @angle1, @idangle1, @angle2, @idangle2, @atom2, @atom3, @atom4);
my (@x, @y, @z);
my $nAtom = 0;
while(1) {
	my $line = $in->ReadLine(1);
	last if($line eq '');
	my $i = $nAtom;
	($atom[$i], $d[$i], $idd[$i], $angle1[$i], $idangle1[$i], $angle2[$i], $idangle2[$i], $atom2[$i], $atom3[$i], $atom4[$i]) 
			= Utils::Split("\\s+", $line);
	$nAtom++;
}

$out->print("\"$Title\"\n");
$out->print("$nAtom\n");
for(my $i = 0 ; $i < $nAtom ; $i++) {
	my $R  = $d[$i];
	my ($x, $y, $z);
	if($i == 0) {
		($x, $y, $z) = (0, 0, 0);
	}
	elsif($i == 1) {
		($x, $y, $z) = ($R, 0, 0);
	}
	elsif($i == 2) {
		$x = $R * cos($angle1[$i] * $torad);
		$y = $R * sin($angle1[$i] * $torad);
		$z = 0.0;
	}
	else {
		my $alpha = $torad * $angle1[$i];
		my $beta  = $torad * $angle2[$i];
		my $sina  = sin($alpha);
		my $cosa  = cos($alpha);
		my $sinb  = sin($beta);
		my $cosb  = cos($beta);
		my $j = $atom2[$i]-1;
		my $k = $atom3[$i]-1;
		my $l = $atom4[$i]-1;

		my @a  = &DifVectors($x[$k], $y[$k], $z[$k], $x[$j], $y[$j], $z[$j]);
		my @b  = &DifVectors($x[$l], $y[$l], $z[$l], $x[$j], $y[$j], $z[$j]);
		my @c  = &OuterProdcut(@b, @a);
#		my @c  = &OuterProdcut(@a, @b);

		my $a  = sqrt(&InnerProdcut(@a, @a));
		my $b  = sqrt(&InnerProdcut(@b, @b));
		my $c  = sqrt(&InnerProdcut(@c, @c));

		my $ab  = &InnerProdcut(@a, @b);

print "$i-$j-$k-$l: $R: $angle1[$i]: $angle2[$i]\n";
printf "a=(%6.3f, %6.3f, %6.3f)\n", @a;
printf "b=(%6.3f, %6.3f, %6.3f)\n", @b;
printf "c=(%6.3f, %6.3f, %6.3f)\n", @c;
#my $ac  = &InnerProdcut(@a, @c);
#my $bc  = &InnerProdcut(@b, @c);
#print "a*c=$ac  b*c=$bc\n";

#my $ab  = &InnerProdcut(@a, @b);
#my $angleab = $todeg * acos($ab / $a / $b);
#print "angle(a-b)=$angleab deg.\n";

		my $B = $R * $a * $sina * $cosb / $c;
		my $A = ($R * $a * $cosa - $B * $ab) / $a / $a;
		my $C  = sqrt($R*$R - $A*$A * $a*$a - $B*$B * $b*$b - 2.0 * $A * $B * $ab) / $c;
		$C = -$C if($angle2[$i] < 0.0);

		($x, $y, $z) = ($x[$j] + ($A*$a[0] + $B*$b[0] + $C*$c[0]), 
						$y[$j] + ($A*$a[1] + $B*$b[1] + $C*$c[1]), 
						$z[$j] + ($A*$a[2] + $B*$b[2] + $C*$c[2]));
	}
	($x[$i], $y[$i], $z[$i]) = ($x, $y, $z);
}

for(my $i = 0 ; $i < $nAtom ; $i++) {
	my ($AtmNum, $name, $charge) = AtomType::GetAtomInformation($atom[$i]);
$x[$i] -= 0.75014;
	$out->printf("%g,%g,%g,%d\n", $x[$i], $y[$i], $z[$i], $AtmNum);
printf("%g,%g,%g,%d\n", $x[$i], $y[$i], $z[$i], $AtmNum);
}
$out->printf("%d\n", $nAtom-1);
for(my $i = 1 ; $i < $nAtom ; $i++) {
	$out->printf("%d,%d,%d\n", $atom2[$i], $i+1, 1);
}

$out->Close();
$in->Close();
exit;

sub OuterProdcut
{
	my (@a) = @_;
	my $dim = 3;
	return ($a[1]*$a[2+$dim]-$a[2]*$a[1+$dim], $a[2]*$a[0+$dim]-$a[0]*$a[2+$dim], $a[0]*$a[1+$dim]-$a[1]*$a[0+$dim]);
}

sub InnerProdcut
{
	my (@a) = @_;
	my $dim = scalar @a / 2;
#print "dim=$dim\n";
	my $c = 0.0;
	for(my $i = 0 ; $i < $dim ; $i++) {
		$c += $a[$i] * $a[$i+$dim];
	}
	return $c;
}

sub DifVectors
{
	my (@a) = @_;
	my $dim = scalar @a / 2;
#print "dim=$dim\n";
	my @c;
	for(my $i = 0 ; $i < $dim ; $i++) {
		$c[$i] = $a[$i] - $a[$i+$dim];
	}
	return @c;
}