#!/usr/bin/perl

$inc_angle = .25;
$two_pi = 6.2831319;

sub vector_product  {
    my($a,$b,$c,$d,$e,$f)=@_;
    # print ($b*$f-$c*$e) . "," . ($c*$d-$a*$f) . "," . ($a*$e-$b*$d) . "\n";
    return($b*$f-$c*$e,$c*$d-$a*$f,$a*$e-$b*$d);
}

sub triangle_normal {
    my($ax,$ay,$az,$bx,$by,$bz,$cx,$cy,$cz)=@_;
    my(@AB)=($bx-$ax,$by-$ay,$bz-$az);
    my(@AC)=($cx-$ax,$cy-$ay,$cz-$az);
    return(vector_product(@AB,@AC));
}

print STDERR "rho = A + B * cos(C * theta) * cos(D * phi)\n";
print STDERR "Enter A:";
$a = <STDIN>;
chomp $a;
print STDERR "Enter B:";
$b = <STDIN>;
chomp $b;
print STDERR "Enter C:";
$c = <STDIN>;
chomp $c;
print STDERR "Enter D:";
$d = <STDIN>;
chomp $d;
print "rho = $a + $b * cos($c * theta) * cos($d * phi)\n";


# determine how many times around we have to go.
$cnr = 1;
$c =~ /(\.[0-9]+)/;
if (defined $1 && $1 != 0) {
	$cnr = 1 / $1;
}

$dnr = 1;
$d =~ /(\.[0-9]+)/;
if (defined $1 && $1 != 0) {
	$dnr = 1 / $1;
}

$theta_inc_angle = $inc_angle;
if ($c >= 1) {
	$theta_inc_angle = $inc_angle / $c;
}

$phi_inc_angle = $inc_angle;
if ($d >= 1) {
	$phi_inc_angle = $inc_angle / $d;
}

open (FACES, ">faces.txt");
open (VERTS, ">verts.txt");
open (NORMS, ">norms.txt");

@verts = ();
@norms = ();
$normals = 1;
$vertices = 1;
for ($theta = 0; $theta < $cnr * $two_pi + $theta_inc_angle; $theta += $theta_inc_angle) {
	
	print "$theta, " . ($cnr * $two_pi + $theta_inc_angle) . "\n";
	$vertcount = 0;
	for ($phi = 0; $phi < $dnr * $two_pi + $phi_inc_angle;  $phi += $phi_inc_angle) {
		$rho = $a + $b * cos($c * $theta) * cos($d * $phi);
		$x = $rho * cos($phi) * cos($theta);
		$y = $rho * cos($phi) * sin($theta);
		$z = $rho * sin($phi);

		# print "theta = $theta ";
		# print "phi = $phi ";
		# print "rho = $rho\n";
		# print "x = $x ";
		# print "y = $y ";
		# print "z = $z\n";
		push @verts, [$x, $y, $z];
		$vertcount++;
	}
	if (defined @oldverts) {
		for ($vc = 0; $vc < $vertcount - 1; $vc++) {
			$verto1 = $oldverts[$vc];
			$verto2 = $oldverts[$vc+1];
			$vertn1 = $verts[$vc];
			$vertn2 = $verts[$vc+1];

			$indo1 = index_vertex($verto1);
			$indo2 = index_vertex($verto2);
			$indn1 = index_vertex($vertn1);
			$indn2 = index_vertex($vertn2);

			$normo2 = index_normal(@$verto1, @$verto2, @$vertn2);
			$normn2 = index_normal(@$verto2, @$vertn2, @$vertn1);
			$normn1 = index_normal(@$vertn2, @$vertn1, @$verto1);
			$normo1 = index_normal(@$vertn1, @$verto1, @$verto2);

			print FACES "f $indo2//$normo2 $indo1//$normo1 $indn1//$normn1\n";
			print FACES "f $indn1//$normn1  $indn2//$normn2  $indo2//$normo2\n";
		}
	}
	if (defined @veryoldverts) {
		foreach $vovert (@veryoldverts) {
			local($pv) = "@$vovert[0],@$vovert[1],@$vovert[2]";
			# $vert{$pv} = undef;
			undef($vert{$pv});
		}
	}
	@veryoldverts = @oldverts;
	@oldverts = @verts;
	@verts = ();
}

@verts = @oldverts;
@oldverts = @veryoldverts;
if (defined @oldverts) {
	for ($vc = 0; $vc < $vertcount - 1; $vc++) {
		$verto1 = $oldverts[$vc];
		$verto2 = $oldverts[$vc+1];
		$vertn1 = $verts[$vc];
		$vertn2 = $verts[$vc+1];

		$indo1 = index_vertex($verto1);
		$indo2 = index_vertex($verto2);
		$indn1 = index_vertex($vertn1);
		$indn2 = index_vertex($vertn2);

		$normo2 = index_normal(@$verto1, @$verto2, @$vertn2);
		$normn2 = index_normal(@$verto2, @$vertn2, @$vertn1);
		$normn1 = index_normal(@$vertn2, @$vertn1, @$verto1);
		$normo1 = index_normal(@$vertn1, @$verto1, @$verto2);

		print FACES "f $indo2//$normo2  $indo1//$normo1  $indn1//$normn1\n";
		print FACES "f $indn1//$normn1  $indn2//$normn2  $indo2//$normo2\n";
	}
}
close FACES;
close VERTS;
close NORMS;

open (VERTS, "verts.txt");
open (NORMS, "norms.txt");
open (FACES, "faces.txt");
open (VENUS, ">venusm.obj");

while (<VERTS>) {
	print VENUS;
}
print VENUS "\n";
while (<NORMS>) {
	print VENUS;
}
print VENUS "\n";
print VENUS "s off\n";
while (<FACES>) {
	print VENUS;
}

close VERTS;
close NORMS;
close FACES;
close VENUS;

sub index_vertex {
	local($vertex) = @_;
	local($pv) = "@$vertex[0],@$vertex[1],@$vertex[2]";
	if (!defined $vert{$pv}) {
		$vert{$pv} = $vertices++;
		print VERTS "v @$vertex[0] " . (@$vertex[1] + 1500) . " @$vertex[2]\n";
	}
	return $vert{$pv};
}

sub index_normal {
	local(@vertices) = @_;
	local(@normal) = triangle_normal(@vertices);
	# print "$normal[0],$normal[1],$normal[2]\n";
	local($pv) = "$normal[0],$normal[1],$normal[2]";
	if (!defined $norm{$pv}) {
		$norm{$pv} = $normals++;
		print NORMS "vn $normal[0] $normal[1] $normal[2]\n";
	}
	return $norm{$pv};
	
}