--- ray/src/util/genklemsamp.pl	2009/06/16 04:54:08	2.2
+++ ray/src/util/genklemsamp.pl	2010/09/26 15:40:53	2.7
@@ -1,12 +1,13 @@
 #!/usr/bin/perl -w
-# RCSid $Id: genklemsamp.pl,v 2.2 2009/06/16 04:54:08 greg Exp $
+# RCSid $Id: genklemsamp.pl,v 2.7 2010/09/26 15:40:53 greg Exp $
 #
 # Sample Klems (full) directions impinging on surface(s)
 #
 #	G. Ward
 #
+use strict;
 if ($#ARGV < 0) {
-	print STDERR, "Usage: genklemsamp [-c N ][-f{a|f|d}] [view opts] [geom.rad ..]\n";
+	print STDERR "Usage: genklemsamp [-c N ][-f{a|f|d}] [view opts] [geom.rad ..]\n";
 	exit 1;
 }
 my $td = `mktemp -d /tmp/genklemsamp.XXXXXX`;
@@ -27,10 +28,10 @@ while ($#ARGV >= 0) {
 		push @vopts, "@ARGV[0..1]";
 		shift @ARGV;
 	} elsif ("$ARGV[0]" =~ /^-v./) {
-		print STDERR, "Unsupported view option: $ARGV[0]\n";
+		print STDERR "Unsupported view option: $ARGV[0]\n";
 		exit 1;
 	} elsif ("$ARGV[0]" =~ /^-./) {
-		print STDERR, "Unknown option: $ARGV[0]\n";
+		print STDERR "Unknown option: $ARGV[0]\n";
 		exit 1;
 	} else {
 		last;
@@ -62,7 +63,7 @@ Kpol = DEGREE * (x1*Kpola(Krow) + (1-x1)*Kpola(Krow-1)
 sin_kpol = sin(Kpol);
 K0 = cos(Kazi)*sin_kpol;
 K1 = sin(Kazi)*sin_kpol;
-K2 = cos(Kpol);
+K2 = sqrt(1 - sin_kpol*sin_kpol);
 ';
 my $ndiv = 145;
 # Do we have any Radiance input files?
@@ -81,13 +82,14 @@ if ($#ARGV >= 0) {
 	push @vopts, ("-vp", $bcube[0]+$bcube[3]/2-$width/2*$vd[0],
 			$bcube[1]+$bcube[3]/2-$width/2*$vd[1],
 			$bcube[2]+$bcube[3]/2-$width/2*$vd[2]);
-	push @vopts, ("-vh", $width, "-vh", $height);
+	push @vopts, ("-vh", $width, "-vv", $height);
 	$vwset = `vwright @vopts V`;
+	$sca = sqrt($nsamp/($width*$height));
 	my $xres;
 	my $yres;
-	my $ntot;
+	my $ntot = 0;
 	# This generally passes through the loop twice to get density right
-	do {
+	while ($ntot < $nsamp) {
 		$xres = int($width*$sca) + 1;
 		$yres = int($height*$sca) + 1;
 		system "vwrays -ff -x $xres -y $yres -pa 0 -pj .7 @vopts " .
@@ -97,15 +99,19 @@ if ($#ARGV >= 0) {
 				q{-e 'nOK=sq(Vdx*$5+Vdy*$6+Vdz*$7)-.999' } .
 				"-if7 -of > $td/origins.flt";
 		$ntot = -s "$td/origins.flt";
-		$ntot /= 3*4;
+		$ntot /= 3*4;		# number of bytes per sample position
 		if ($ntot == 0) {
-			print STDERR, "View direction does not correspond to any surfaces\n";
+			if ($nsamp < 200) {
+				$sca = sqrt(200/($width*$height));
+				redo;
+			}
+			print STDERR "View direction does not correspond to any surfaces\n";
 			exit 1;
 		}
 		$sca *= 1.05 * sqrt($nsamp/$ntot);
-	} while ($ntot < $nsamp);
+	}
 	# All set to produce our samples
-	for ($k = 1; $k <= $ndiv; $k++) {
+	for (my $k = 1; $k <= $ndiv; $k++) {
 		my $rn = rand(10);
 		my $r1 = rand; my $r2 = rand;
 		# Chance of using = (number_still_needed)/(number_left_avail)
@@ -113,9 +119,9 @@ if ($#ARGV >= 0) {
 			"-e 'cond=($nsamp-outno+1)/($ntot-recno+1)-rand($rn+recno)' " .
 			"-e 'Kbin=$k;x1=rand($r1+recno);x2=rand($r2+recno)' " .
 			q{-e '$1=$1+Vo*Vdx; $2=$2+Vo*Vdy; $3=$3+Vo*Vdz' } .
-			q{-e '$4=K0*Vhx+K1*Vvx+K2*Vdx' } .
-			q{-e '$5=K0*Vhy+K1*Vvy+K2*Vdy' } .
-			q{-e '$6=K0*Vhz+K1*Vvz+K2*Vdz' } .
+			q{-e '$4=K0*Vhx-K1*Vvx+K2*Vdx' } .
+			q{-e '$5=K0*Vhy-K1*Vvy+K2*Vdy' } .
+			q{-e '$6=K0*Vhz-K1*Vvz+K2*Vdz' } .
 			"-if3 $td/origins.flt";
 	}
 } else {
@@ -123,7 +129,7 @@ if ($#ARGV >= 0) {
 	my $xres = int($width*$sca) + 1;
 	my $yres = int($height*$sca) + 1;
 	my $ntot = $xres * $yres;
-	for ($k = 1; $k <= $ndiv; $k++) {
+	for (my $k = 1; $k <= $ndiv; $k++) {
 		my $rn = rand(10);
 		my $r1 = rand; my $r2 = rand;
 		my $r3 = rand; my $r4 = rand;
@@ -137,9 +143,9 @@ if ($#ARGV >= 0) {
 			q{-e '$1=Vpx+Vo*Vdx+hpos*Vhx+vpos*Vvx' } .
 			q{-e '$2=Vpy+Vo*Vdy+hpos*Vhy+vpos*Vvy' } .
 			q{-e '$3=Vpz+Vo*Vdz+hpos*Vhz+vpos*Vvz' } .
-			q{-e '$4=K0*Vhx+K1*Vvx+K2*Vdx' } .
-			q{-e '$5=K0*Vhy+K1*Vvy+K2*Vdy' } .
-			q{-e '$6=K0*Vhz+K1*Vvz+K2*Vdz' } ;
+			q{-e '$4=K0*Vhx-K1*Vvx+K2*Vdx' } .
+			q{-e '$5=K0*Vhy-K1*Vvy+K2*Vdy' } .
+			q{-e '$6=K0*Vhz-K1*Vvz+K2*Vdz' } ;
 	}
 }
 system "rm -rf $td";